chickdata<- read.csv("chickdata.csv")Last assignment (optional)
Repeated measures
As discussed in class, repeated measures can be on of the most complex experimental designs to analyze. There are many reasons for this, but the main one is that time is hard to model, and that there is temporal autocorrelation.
Chicken data
We will do the same example as in class. Make sure to look at the slides! Download the chickdata file and read it
As a reminder, this is what the data looks like:
| Chick | Diet | Day 0 | Day 2 | Day 4 | Day 6 | Day 8 | Day 10 | Day 12 | Day 14 | Day 16 | Day 18 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1 | 41 | 48 | 56 | 68 | 75 | 86 | 102 | 107 | 135 | 155 |
| 2 | 1 | 42 | 43 | 48 | 60 | 79 | 107 | 142 | 164 | 197 | 199 |
| 10 | 2 | 40 | 49 | 63 | 85 | 124 | 163 | 217 | 240 | 277 | 308 |
| 11 | 2 | 43 | 55 | 64 | 77 | 88 | 96 | 108 | 110 | 132 | 148 |
| 19 | 3 | 40 | 51 | 63 | 74 | 85 | 103 | 123 | 138 | 170 | 203 |
| 20 | 3 | 42 | 48 | 64 | 81 | 105 | 132 | 159 | 178 | 221 | 263 |
| 28 | 4 | 40 | 49 | 62 | 84 | 104 | 125 | 160 | 172 | 203 | 234 |
| 29 | 4 | 43 | 56 | 69 | 95 | 132 | 158 | 184 | 187 | 197 | 197 |
The data that you downloaded is in the long format. What you see in the table is the same data in wide format. We can use the tidyr function pivot_longer and pivot_wider to go from wide to long and long to wide. This can be super useful.
Most analyses in R require our data to be in the long format.
This is some of our data in graphical form:
Additive model (split plot)
In the additive model, we can use the following equation:
\[ y_{ijk} = \mu + \alpha_i + \beta_j + \alpha\beta_{ij} + \delta_{ik} + \epsilon_{ijk} \]
where,
\(\mu\) is the grand mean, \(\alpha_i\) is the effect of diet, \(\beta_j\) is the effect of time, \(\alpha \beta_{ij}\) is the interaction of time and diet, \(\delta_ik\) is the effect of the kth subject receiving the ith treatment, and \(\epsilon_{ijk}\) is the residual error.
To fit this model, we use the followinf:
modelaov<-aov(Weight ~ Diet * Time + Error(Chick),data=chickdata)
summary(modelaov)
##
## Error: Chick
## Df Sum Sq Mean Sq
## Diet 1 21915 21915
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 1 9771 9771 12.33 0.000504 ***
## Time 1 918896 918896 1159.18 < 2e-16 ***
## Diet:Time 1 13628 13628 17.19 4.23e-05 ***
## Residuals 355 281412 793
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1While this might be ok, there is an issue. Our errors are not independent, so we might be underestimating our p-value. This is because of an an issue called temporal autocorrelation. How much a chick grew in time 2 will affect its growth in time 4, and that growth will affect time 6 and so on. We need to take this into account by adjusting our p-values up. You should ONLY do the following, if the original ANOVA was significant. If it wasn’t, there is no need, as the p-value will only go up.
In order to adjust our p-values, we estimate the sphericity of the data. There are two methods to do this, but the good news is that R runs both.
First off, we need “wide” data:
chickdata_wide <- tidyr::pivot_wider(chickdata, names_from = Time,
names_prefix = "Day", values_from = Weight)Check the new dataset.
To run a spherical test, we run a manova :
man<-manova(cbind(Day0, Day2, Day4, Day6, Day8,Day10,Day12,Day14,Day16,Day18) ~ Diet, data = chickdata_wide)This model is already taking into account both the effects of diet, time and the interactions. However, it is running it a bit differently. The results from this are:
anova(man, X = ~1, test = "Spherical")
## Analysis of Variance Table
##
##
## Contrasts orthogonal to
## ~1
##
## Greenhouse-Geisser epsilon: 0.1362
## Huynh-Feldt epsilon: 0.1387
##
## Df F num Df den Df Pr(>F) G-G Pr H-F Pr
## (Intercept) 1 220.2712 9 306 0.00000000 0.000000 0.000000
## Diet 1 3.2681 9 306 0.00081711 0.069891 0.068993
## Residuals 34Where Intercept is the effect of time, and Diet is the interaction of time and diet.
In this case, our p-values for the interaction went from p<0.0001 to p>0.05, therefore they are not significant anymore. This is a result of the p-value adjustment.
Profile analysis
Probably the better way to do it is with a profile analysis. In this case, rather than looking at the wight of the chickens, we look at the weight difference between measurements.
To do this, we do the following:
man2<-manova(cbind(Day2-Day0,Day4- Day2, Day6-Day4, Day8-Day6,Day10- Day8,Day12-Day10,Day14-Day12,Day16-Day14,Day18-Day16) ~ Diet, data = chickdata_wide)See how we have only 9 intervals, despite having 10 times.
summary.aov(man2)
## Response 1 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 1 86.806 86.806 12.199 0.001348 **
## Residuals 34 241.944 7.116
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response 2 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 1 81.339 81.339 11.837 0.001555 **
## Residuals 34 233.633 6.872
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response 3 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 1 256.81 256.806 14.913 0.0004805 ***
## Residuals 34 585.50 17.221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response 4 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 1 325.36 325.36 4.2696 0.04648 *
## Residuals 34 2590.87 76.20
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response 5 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 1 53.36 53.356 0.7516 0.392
## Residuals 34 2413.53 70.986
##
## Response 6 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 1 224.5 224.45 1.6374 0.2093
## Residuals 34 4660.5 137.07
##
## Response 7 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 1 39.2 39.200 0.5153 0.4778
## Residuals 34 2586.7 76.079
##
## Response 8 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 1 180 180.00 0.7873 0.3811
## Residuals 34 7773 228.62
##
## Response 9 :
## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 1 680.6 680.56 2.8453 0.1008
## Residuals 34 8132.4 239.19To plot these results, we do the following:
Weightm<-chickdata_wide[,3:12]
weightdif<-Weightm[2:10]-Weightm[1:9]
colnames(weightdif) <- paste("interval", 1:9, sep=".")
kable(head(weightdif))| interval.1 | interval.2 | interval.3 | interval.4 | interval.5 | interval.6 | interval.7 | interval.8 | interval.9 |
|---|---|---|---|---|---|---|---|---|
| 7 | 8 | 12 | 7 | 11 | 16 | 5 | 28 | 20 |
| 1 | 5 | 12 | 19 | 28 | 35 | 22 | 33 | 2 |
| 7 | 11 | 14 | 23 | 26 | 18 | 6 | 8 | 4 |
| 6 | 9 | 17 | 16 | 24 | 31 | 30 | 45 | 30 |
| 11 | 7 | 11 | 15 | 14 | -9 | 3 | -2 | 9 |
| 4 | 10 | 7 | 12 | 7 | 7 | 9 | 3 | 10 |
This is the data that we have now:
This is a complex analysis, so, it’s OK if you get a bit lost.
We need to estimate the means for each week, for each diet:
Diet1<-colMeans(weightdif[1:9,])
Diet2<-colMeans(weightdif[10:18,])
Diet3<-colMeans(weightdif[19:27,])
Diet4<-colMeans(weightdif[28:36,])Now, we calculate the Standard Errors:
SE <- sqrt(diag(stats:::vcov.mlm(man2)))
SE <- SE[names(SE)==":(Intercept)"] # Only use "intercept" SEs
unname(SE) ## Ignore the names
## [1] 1.089037 1.070169 1.694136 3.563752 3.439629 4.779717 3.560877 6.172758
## [9] 6.313868Finally, we can put everything together, and plot:
growthDF <- data.frame(interval = rep(1:9, 4),
Diet = rep(c("1", "2", "3", "4"), each = 9),
growth = c(Diet1, Diet2, Diet3, Diet4),
SE = rep(SE, 2))
ggplot(growthDF, aes(x = interval, y = growth, color = Diet)) +
geom_path() +
geom_point() +
geom_errorbar(aes(ymin = growth - SE, ymax = growth + SE), width = 0.1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
scale_x_continuous("Time interval") +
scale_y_continuous("Growth rate", limits = c(-1, 40))+
theme_classic()
Mixed model with correlation
Finally, you can run it as a mixed model with correlated data:
library(nlme)
lmm2 <- lme(fixed = Weight ~ Diet *Time,
random = ~ (Time-1)|Chick,
data = chickdata,
correlation = corAR1())
summary(lmm2)
## Linear mixed-effects model fit by REML
## Data: chickdata
## AIC BIC logLik
## 2716.979 2744.104 -1351.49
##
## Random effects:
## Formula: ~(Time - 1) | Chick
## Time Residual
## StdDev: 2.771936 15.63854
##
## Correlation Structure: AR(1)
## Formula: ~1 | Chick
## Parameter estimate(s):
## Phi
## 0.8287057
## Fixed effects: Weight ~ Diet * Time
## Value Std.Error DF t-value p-value
## (Intercept) 37.27489 6.217745 322 5.994922 0.0000
## Diet -0.58195 2.270399 34 -0.256320 0.7992
## Time 6.17670 1.218546 322 5.068908 0.0000
## Diet:Time 0.99140 0.444950 322 2.228125 0.0266
## Correlation:
## (Intr) Diet Time
## Diet -0.913
## Time -0.243 0.221
## Diet:Time 0.221 -0.243 -0.913
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.7100799 -0.8975167 -0.2215318 0.2858934 2.3128864
##
## Number of Observations: 360
## Number of Groups: 36
anova(lmm2)
## numDF denDF F-value p-value
## (Intercept) 1 322 357.1048 <.0001
## Diet 1 34 0.0859 0.7713
## Time 1 322 302.7071 <.0001
## Diet:Time 1 322 4.9645 0.0266This is what I would recommend. And it actually shows an effect!
Power analysis
This is a very straightforward analysis.
There are four parameters in a model that are dependent on each other:
- Degrees of freedom
- Effect (how big of a difference?)
- \(\alpha\) generally 0.05
- Power
If you have one, you can estimate the rest. For effect, you can choose an effect size metric or you can also give exact differences (and standard deviation)
For example, if we are comparing the effects of fertilizer on plant growth. We believe a growth of 8 cm would be biologically important, we want an \(\alpha\) of 0.05 and a power of 0.5 (generally pretty good power). We have observed these populations have a standard deviation of ~3 we can do the following:
power.t.test(delta=5,power=.80,sd=3)
##
## Two-sample t test power calculation
##
## n = 6.76095
## delta = 5
## sd = 3
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* groupEssentially, we only need 7 plants in each group!
If we want to find:
1) smaller differences we need more plants:
power.t.test(delta=2,power=.80,sd=3)
##
## Two-sample t test power calculation
##
## n = 36.3058
## delta = 2
## sd = 3
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* group37 plants per group to be able to find these super small differences!
Power analysis ANOVA
In Anova, we are able to use the package pwr. Because in ANOVA, oftentimes we have complex models with multiple factors, it is hard to completely provide an effect size in absolute terms. Instead of this, we can use Cohen’s conventional effect size.
For our example, we will do a 2 factorial effect. The Cohen’s value would be:
library(pwr)
eff<-cohen.ES("f2", size = "medium")$effect.size
eff
## [1] 0.15For a small value it would be:
eff<-cohen.ES("f2", size = "small")$effect.size
eff
## [1] 0.02In this case, I will use 0.9, which is a value between both of them.
Let’s assume in this case we are doing a full factorial 3 by 3 ANOVA
pwr.f2.test(u = 1, # I'm assuming you plan to test a 1 df effect
f2 = 0.09, # We're using Cohen's effect size guidelines, shown above
sig.level = .05, # Our alpha
power = .80) # Our desired power
##
## Multiple regression power calculation
##
## u = 1
## v = 87.17156
## f2 = 0.09
## sig.level = 0.05
## power = 0.8In this case the degrees of freedom of the “error term” are 87.17. If we add the 8 df’s from the 3*3 design plus the intercept, we get that we need 97 observations TOTAL. In this case 11 observations per group would be enough.