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!
