day, and their interaction. It makes sense to model time as a continuous variable. The response and fixed effects parts of the model formula should be
weight ~ trt + day + trt:day.
planton the group side.
Notice we use
cs(day | plant) to specify that it is a
compound-symmetry error structure, and the repeated measures are by day
within each individual plant.
fit_CS <- glmmTMB(weight ~ trt + day + trt:day + cs(day | plant), data = pederson.lettuce.repeated)
residuals_CS <- simulateResiduals(fit_CS) plot(residuals_CS)
The residual diagnostics do not look good. The residual versus predicted plot has a distinct U-shaped trend. For the assumption of homogeneous variance to be met, there shouldn’t be any trend in this plot.
This is identical to the syntax above except for the addition of the
family = lognormal argument.
fit_CS_lognormal <- glmmTMB(weight ~ trt + day + trt:day + cs(day | plant), data = pederson.lettuce.repeated, family = lognormal) residuals_CS_lognormal <- simulateResiduals(fit_CS_lognormal) plot(residuals_CS_lognormal)
The residual diagnostics look much better. There is no trend in the residual versus predicted plot.
Notice we use
ar1(0 + factor(day) | plant). The
0 is because you must explicitly remove the intercept when
fitting AR(1) error structure, and
AR(1) is intended for equally sized discrete time steps.
fit_AR1_lognormal <- glmmTMB(weight ~ trt + day + trt:day + ar1(0 + factor(day) | plant), data = pederson.lettuce.repeated, family = lognormal) residuals_AR1_lognormal <- simulateResiduals(fit_AR1_lognormal) plot(residuals_AR1_lognormal) AIC(fit_CS_lognormal, fit_AR1_lognormal)
The residual diagnostics look good, and the information criterion comparison indicates that the AR(1) fit is far superior. Despite having roughly similar number of parameters to the compound symmetry fit, it is a much better fit to the data. We will use it in the means comparison in the following exercise.
trt_means_by_day <- emmeans(fit_AR1_lognormal, ~ trt | day, at = list(day = c(1, 21, 40)), type = 'response') cld(trt_means_by_day, adjust = 'tukey', Letters = letters)
Even though it looks like treatment 1 weight gain is slightly higher, this difference is not statistically significant for any of the three days tested. There is a lot of variation around the mean trend within each of the three treatments, as can be seen from the plot we made earlier.