#-------------------------------------------------------------------------# # Illustration of different approaches to repeated measures analysis # #-------------------------------------------------------------------------# #-------------------------------------------------------------------------# # Overview of approaches: # 1) a multivariate analysis (MANOVA) # 2) use split-plot analysis with potentially downweighted split plot degrees of freedom # 3) define contrasts among time levels, reducing complexity # 4) model the covariance structure # 5) use random coefficient models # # data from microbial CO2 growth over time under different moisture conditions from # "Design of Experiments: Statistical Principles of Research Design and Analysis" # by Robert O. Kuehl, second edition, page 504. #-------------------------------------------------------------------------# data <- scan(what = list(moisture="",container=0,day=0,CO2=0), multi.line=F) Control 1 2 .22 Control 1 4 .56 Control 1 6 .66 Control 1 8 .89 Control 2 2 .68 Control 2 4 .91 Control 2 6 1.06 Control 2 8 .80 Control 3 2 .68 Control 3 4 .45 Control 3 6 .72 Control 3 8 .89 P24 4 2 2.53 P24 4 4 2.7 P24 4 6 2.1 P24 4 8 1.5 P24 5 2 2.59 P24 5 4 1.43 P24 5 6 1.35 P24 5 8 0.74 P24 6 2 0.56 P24 6 4 1.37 P24 6 6 1.87 P24 6 8 1.21 P26 7 2 0.22 P26 7 4 0.22 P26 7 6 0.20 P26 7 8 0.11 P26 8 2 0.45 P26 8 4 0.28 P26 8 6 1.24 P26 8 8 0.86 P26 9 2 0.22 P26 9 4 0.33 P26 9 6 0.34 P26 9 8 0.20 P28 10 2 0.22 P28 10 4 0.80 P28 10 6 0.80 P28 10 8 0.37 P28 11 2 0.22 P28 11 4 0.62 P28 11 6 0.89 P28 11 8 0.95 P28 12 2 0.22 P28 12 4 0.56 P28 12 6 0.69 P28 12 8 0.63 data table15.6 <- data.frame(data) # using the within function to convert variables to factor table15.6 <- within(table15.6, { moisture <- factor(moisture) container <- factor(container) day <- factor(day) } ) # get an interaction plot with(table15.6, interaction.plot(day,moisture,CO2) ) #-------------------------------------------------------------------------# # afex can be used for the first, second, and fifth approaches # #-------------------------------------------------------------------------# # use the afex package to get several types of analyses: # standard split-plot analysis, assuming the HF condition is valid library(afex) library(emmeans) table15.6.afex0 <- mixed(CO2 ~ moisture + (1|moisture:container) +day +moisture:day, data = table15.6) nice(table15.6.afex0) # using aov_ez to get the Mauchley test, adjusted degrees of freedom approach, # and the MANOVA approach table15.6.afex1 <- aov_ez(id="container",dv="CO2",data=table15.6, between=c("moisture"),within=c("day")) # the summary function produces Type III anova results assuming sphericity, # and includes the Mauchly tests for sphericity summary(table15.6.afex1) #-------------------------------------------------------------------------# # Second approach to repeated measures # #-------------------------------------------------------------------------# # the nice function produces Type III anova results correcting for # possible lack of sphericity by downweighting split plot degrees of freedom nice(table15.6.afex1) #-------------------------------------------------------------------------# # First approach to repeated measures # #-------------------------------------------------------------------------# # the Anova component of the aov_ez object produces the MANOVA results # for the split plot tests (not discussed in the text) table15.6.afex1$Anova # an interaction plot lsmip(table15.6.afex1, moisture ~ day) #-------------------------------------------------------------------------# # Fifth approach to repeated measures # #-------------------------------------------------------------------------# # must convert day back to being numeric to fit a polynomial pattern table15.6$day <- as.numeric(table15.6$day) table15.6.afexrc1 <- mixed(CO2 ~ moisture +day +moisture:day +(1 + day|container), data = table15.6) summary(table15.6.afexrc1) nice(table15.6.afexrc1) qqnorm(summary(table15.6.afexrc1)$residuals) qqline(summary(table15.6.afexrc1)$residuals) # We can use the nlme package to compare covariance structures with repeated # measures data. # Comparing the results here to those from Proc MIXED, the loglikelihood values # and AIC values are consistent with each other. The BIC values are not quite # consistent, but here lead to the same best model. # In the ANOVA results, the F values for the split plot tests match those # from Proc MIXED, but R does not seem to use the correct degrees of freedom library(nlme) # first we use the groupedData function to create a data structure approriate for # use with the repeated measures analysis t15.6.longplot <- groupedData(CO2 ~ as.numeric(day) | container, outer = ~ moisture, data = table15.6) plot(t15.6.longplot,outer = ~ moisture) # now a structure for analysis # t15.6.long2 <- groupedData(CO2 ~ as.factor(day) | container, outer = ~ moisture, # data = table15.6) t15.6.long <- groupedData(CO2 ~ as.numeric(moisture)*as.numeric(day) | container, data = table15.6) #-------------------------------------------------------------------------# # Fourth approach to repeated measures # #-------------------------------------------------------------------------# # For some of these structures, the F values and P values do not agree # with SAS Proc Mixed results using the DDFM = KR option. # Instead they agree with Proc Mixed results using other options such # as DDFM = BW, CON, or RES # compound symmetric covariance structure table15.6.cs <- gls(CO2 ~ moisture * day, data = t15.6.long, corr = corCompSymm(, form= ~ 1 | container) ) summary(table15.6.cs) anova(table15.6.cs) # order 1 autoregressive covariance structure table15.6.ar1 <- gls(CO2 ~ moisture * day, data = t15.6.long, corr = corAR1(, form= ~ 1 | container)) summary(table15.6.ar1) anova(table15.6.ar1) # order 1 heterogeneous autoregressive covariance structure table15.6.arh1 <- gls(CO2 ~ moisture * day, data = t15.6.long, corr = corAR1(, form = ~ 1 | container), weight = varIdent(form = ~ 1 | day)) summary(table15.6.arh1) anova(table15.6.arh1) # an unstructured covariance matrix table15.6.un <- gls(CO2 ~ moisture * day, data = t15.6.long, corr=corSymm(form = ~ 1 | container), weights = varIdent(form = ~ 1 | day)) summary(table15.6.un) anova(table15.6.un)