R for SAS Users Lesson 3 Exercise Answers

Exercise 1

• 1. The response variable is weight.
• 2. The fixed effects are treatment, trt, time, 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.
• 3. The repeated measures are grouped by individual plant. Any random effect should have plant on the group side.

Exercise 2

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)

Exercise 3

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.

Exercise 4

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.

Exercise 5

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 factor(day) because 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.

Exercise 6

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.