2
$\begingroup$

I wrote a model to look at the effect of species + treatment on swimming speed in copepods. I had a bunch of videos and pulled the average swimming speed from each video. Here's a barplot of the data.

> model <- lm(log(Avg_Speed) ~ Species + Treatment + Species:Treatment, 
+               data = total)
> summary(model)

Call:
lm(formula = log(Avg_Speed) ~ Species + Treatment + Species:Treatment, 
    data = total)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.70807 -0.29267 -0.06822  0.26581  1.77693 

Coefficients:
                               Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      1.1753     0.1050  11.192  < 2e-16 ***
SpeciesAcartia                  -0.5187     0.1560  -3.324  0.00107 ** 
SpeciesParvo                    -1.6376     0.1560 -10.494  < 2e-16 ***
SpeciesOithona                  -2.1682     0.1560 -13.894  < 2e-16 ***
TreatmentFicoll                 -0.3131     0.1512  -2.071  0.03978 *  
SpeciesAcartia:TreatmentFicoll  -0.2102     0.2430  -0.865  0.38803    
SpeciesParvo:TreatmentFicoll     0.3127     0.2213   1.413  0.15943    
SpeciesOithona:TreatmentFicoll   0.5356     0.2213   2.420  0.01650 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5655 on 184 degrees of freedom
Multiple R-squared:  0.6696,    Adjusted R-squared:  0.6571 
F-statistic: 53.28 on 7 and 184 DF,  p-value: < 2.2e-16

Calanus Control is the intercept in the model. Now, this looks ok except for when you actually look at what R is saying about the groups - there's a significant effect of treatment for Oithona but not Acartia? Huh? Take a look at the barplot. There's no way that's accurate. So I calculated the means by hand (I can share the data over Google Drive if you want):

> mean(log(subset(total, subset = Species == "Parvo" & Treatment == "Control")[,3]))
[1] -0.4623169
> mean(log(subset(total, subset = Species == "Parvo" & Treatment == "Ficoll")[,3]))
[1] -0.4627922
> mean(log(subset(total, subset = Species == "Oithona" & Treatment == "Control")[,3]))
[1] -0.9929218
> mean(log(subset(total, subset = Species == "Oithona" & Treatment == "Ficoll")[,3]))
[1] -0.7704915
> mean(log(subset(total, subset = Species == "Acartia" & Treatment == "Control")[,3]))
[1] 0.6565436
> mean(log(subset(total, subset = Species == "Acartia" & Treatment == "Ficoll")[,3]))
[1] 0.1331612
> mean(log(subset(total, subset = Species == "Calanus" & Treatment == "Control")[,3]))
[1] 1.175278
> mean(log(subset(total, subset = Species == "Calanus" & Treatment == "Ficoll")[,3]))
[1] 0.862128

Note that these means haven't been backtransformed so they're actually log(swimming speed). This is so they can be directly compared to the model output. By comparing the group means calculated by hand, you can see that the interaction coefficients are way off. I checked the difference between various group means and it seems like R is comparing random/incorrect groups to get these coefficients. For example, the SpeciesParvo:TreatmentFicoll coefficient is actually the difference between the Oithona Ficoll group and the Parvocalanus Control group (or Parvocalanus Ficoll, they have the same mean). So... what gives? Any help is much appreciated!

$\endgroup$
1
  • 1
    $\begingroup$ The interaction coefficient estimates (and associated standard errors) are still conditional on the reference levels used, even if it's a reference level for the interaction. So while the model output appears to say there's no evidence for a treatment effect for Acartia, what it's actually saying is there's no evidence of a difference between the Acartia treatment effect and the Calanus treatment effect. $\endgroup$ – Matt Tyers Feb 24 '17 at 0:28
0
$\begingroup$

Here's an example of how I would use the coefficient estimates to arrive at the group means:

Picking on Parvo/Ficoll, the mean should be [intercept] + [Parvo] + [Ficoll] + [Parvo:Ficoll]. And, in fact,

> 1.1753-1.6376-0.3131+0.3127
[1] -0.4627
$\endgroup$
1
  • $\begingroup$ THANK YOU it turns out my reference for interaction coefficients was worded in a very misleading way and implied that Parvo/Ficoll should be calculated as intercept+parvo+interaction (ignoring ficoll) $\endgroup$ – housetyrell Feb 24 '17 at 15:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.