2

I have data on survival (time in minutes) for animals (salt-marsh hoppers Orchestia gammarellus) from saline and brackish plots, and my model shows a statistically clear interaction between salinity and dry weight (see image below). (actually my model is a GLM with gamma distribution).

I would like to calculate how large this effect is, but I know I cannot use something like Cohen's d, which would partition out the covariate. I was thinking about partial correlation, $\eta^2_{partial}$, which is calculated as $\eta^2_{partial} = \frac{SS_{Treat}}{SS_{Treat} + SS_{Res}}$ (according to this post).

If we look at the image below and focus on the Brackish treatment (red line), my gut feeling is that the effect size is the ratio between the sum of squares of only the red points to the red line versus the sum of squares of all points to the red line (not the overall average). Is this indeed what partial correlation does? Does this not create a smaller number when the difference is larger? Or is $SS_{Res}$ for all points against an overall trendline?

How do I calculate effect size for ANCOVA in R?

I have tried running pcor.test from the ppcor library, but that accepts only continuous explanatory variables.

Sample data:

survival_time dry_weight salinity
552           6.7542     B
554           4.9427     A

Raoul Van Oosten
  • 45
  • 6
  • 1
    What do you think of a ratio of the SSE of the intercept-and-covariate-only model and the SSE of the full model? A large value would indicate a major reduction in variability once you account for the groups, and this ratio does not consider sample size. – Dave Dec 11 '19 at 09:43
  • Sorry for not responding earlier, but I think this is actually exactly what I was looking for. I think this is what eta squared in ANOVA actually is: the effect of a particular variable. Since I have a GLM, it took me a while to figure out I can use the residual deviance as a measure of SSE. – Raoul Van Oosten Dec 16 '19 at 16:11

1 Answers1

2

As you state, partial eta-squared is the SS for that treatment divided by the sum of SS for that treatment and the residual. You just take these values out of the anova table and calculate. You have to be a little careful if you are using type I or type II sums of squares, so that the anova and the partial eta-squared match.

The following data is taken from the Handbook of Biological Statistics.

Input = ("
Species   Temp   Pulse
 ex       20.8   67.9
 ex       20.8   65.1
 ex       24     77.3
 ex       24     78.7
 ex       24     79.4
 ex       24     80.4
 ex       26.2   85.8
 ex       26.2   86.6
 ex       26.2   87.5
 ex       26.2   89.1
 ex       28.4   98.6
 ex       29    100.8
 ex       30.4   99.3
 ex       30.4  101.7
 niv      17.2   44.3
 niv      18.3   47.2
 niv      18.3   47.6
 niv      18.3   49.6
 niv      18.9   50.3
 niv      18.9   51.8
 niv      20.4   60
 niv      21     58.5
 niv      21     58.9
 niv      22.1   60.7
 niv      23.5   69.8
 niv      24.2   70.9
 niv      25.9   76.2
 niv      26.5   76.1
 niv      26.5   77
 niv      26.5   77.7
 niv      28.6   84.7
")

Data = read.table(textConnection(Input),header=TRUE)

my_aov = aov (Pulse ~ Temp + Species + Temp:Species, data = Data)

summary(my_aov)

   ###              Df Sum Sq Mean Sq  F value   Pr(>F)    
   ###     Temp          1   7895    7895 2505.583  < 2e-16 ***
   ###     Species       1    598     598  189.789 9.91e-14 ***
   ###     Temp:Species  1      4       4    1.357    0.254    
   ###     Residuals    27     85       3

7895 / (7895 + 85)

   ### 0.9893484

if(!require(lsr)){install.packages("lsr")}

library(lsr)

etaSquared(my_aov)

   ###                    eta.sq eta.sq.part
   ### Temp         0.5099036842  0.98093004
   ### Species      0.0696797682  0.87545480
   ### Temp:Species 0.0004982162  0.04785434
Sal Mangiafico
  • 7,128
  • 2
  • 10
  • 24
  • Thanks! Now, I have a GLM with gamma distribution because I have time-to-death data, so there is no SS given by summary. Is it still possible to calculate this? [Here](https://stats.stackexchange.com/questions/46345/how-to-calculate-goodness-of-fit-in-glm-r) I read `1 - (Residual Deviance/Null Deviance)` is a good statistic for model fit. Is $SS_{partial}$ equivalent to residual deviance if I run the model with only that treatment included? – Raoul Van Oosten Dec 13 '19 at 08:22
  • My main model has a residual deviance of 1.6209, and the model with only brackish points of 0.96414. So can I say $\eta^2_{partial} = 0.96414 / (0.96414 + 1.6209) = 0.37$? – Raoul Van Oosten Dec 13 '19 at 10:46
  • I also wonder: when I take the saline points I obviously get a different RSS (i.e. 0.65680), and therefore a different eta squared. How should I interpret and report this? I would like to get an effect of the variable in general, not specific levels. – Raoul Van Oosten Dec 13 '19 at 16:34
  • I'm really not sure a standard method to do this with a glm. It might make sense to get a measure of fit for increasing complex models, and then present the difference as an increase in, say, r-squared. This is probably not ideal, but should be readily understandable for your audience. Using the data from my answer, the following will output Efron's R-square for models. – Sal Mangiafico Dec 13 '19 at 16:57
  • `if(!require(rcompanion)){install.packages("rcompanion")}; model0 = glm(Pulse ~ 1, data=Data); model1 = glm(Pulse ~ Temp, data=Data); model2 = glm(Pulse ~ Temp + Species, data=Data); library(rcompanion); accuracy(list(model0, model1, model2))` – Sal Mangiafico Dec 13 '19 at 16:57
  • Thanks! I agree that a difference in something like r squared is ideal. As far as I know eta squared is more or less the same. So I could compare residual deviance of the simplified model (1.6209) to the fuller model (1.3993) as treatment and residuals, respectively, to get $\eta^2_{partial} = \frac{1.62}{(1.62 + 1.40)} = 0.54$. This can then be added to the p-value of 0.02. That makes a lot of sense to me. – Raoul Van Oosten Dec 16 '19 at 10:38