3
$\begingroup$

I want to run and then report a multi-level model on the ChickWeight data in R.

I have done several analyses of which I have a few questions. I've been told the AIC is a good way to test a model's predictive out-of-sample performance.

I have considered the following models, and the bottom one has the lowest AIC, which if i'm not wrong, is a good thing?

library(tidyverse)
library(lme4)
data("ChickWeight", package = "datasets")

chickm <- lmer(weight ~ Time + (1 | Chick), data = ChickWeight)
chickm2 <- lmer(weight ~ Time + (Time | Chick), data = ChickWeight)
chickm3 <- lmer(weight ~ Time + Diet + (Time | Chick), data = ChickWeight)
chickm4 <- lmer(weight ~ Time * Diet + (Time | Chick), data = ChickWeight)

I want to report chickm4 since the interaction effect seems to improve the model. However chickm3 would also be acceptable since I would like Diet in there.

Next question is about the intercept that the model(s) produce. The actual intercept (weight of chick when days = 0 aka it is born, is ~ 41g). Yet chickm3 gives the intercept as 26g and chickm4 as 33g. You can see in the graph below that the linear lines-of-best-fit cause this. The intercepts of the lines are all below the true points.

library(tidyverse)

ggplot(ChickWeight, aes(x = Time, y = weight, colour = Diet)) + geom_point() +
  stat_smooth(method = 'lm', se = F) + theme_minimal()

Is this a problem or an example of a bad model, or just a downside to the data and/or linear models compared to non-linear models? I do have to report a multi-level linear regression though, so I can't change this.

Final question. chickm3 shows 'Model failed to converge with max|grad| = 0.00225537 (tol = 0.002, component 1)' when I summary() it. What does this mean and is it bad? My chickm4 model doesn't do this.

$\endgroup$

2 Answers 2

4
$\begingroup$

A couple of points:

  • These models are nested and therefore you can also use formal hypothesis testing using F-test to compare them. These F-test are provided by the lmerTest package.
  • In general it is not good to only look at measures like AIC or p-values to decide which terms to include in the model, especially the fixed effects in this case. You should also see how substantive is the difference between the additive model chickm3 and the interaction model chickm4. You could more easily assess this using an effects plots - see, e.g., the effects or ggeffects packages.
  • Indeed the effect that the intercept of the model does not correspond to the average value you have for Time equal to zero and Diet equal to the reference level is because you fit a line, and the estimate of the slope also affects the estimate of the intercept. You could mitigate this affect but fitting a nonlinear model for the Time variable using, e.g., splines.
  • Finally, regarding the convergence warning message, you could ignore it in this case because it is very close to the tolerance value.
$\endgroup$
3
$\begingroup$

I've made a few example models which should be nested if you just test them in order. The model I believe best is chickm6 which:

  • removes the baseline Time 0 Diet effect, which is intuitive since at Time 0 (I would imagine at birth), the Diet should not have any effect on chick weight
  • added a linear two-piece spline, with the knot placed at Time 6 (somewhat determined by eyeballing the loess curve). This allows the modelled baseline Time 0 response (39.8) to be much closer to observed mean weight (41), while keeping the model simple, explainable and still linear. Of course you may try with more knots.

enter image description here

library(stats)
library(tidyverse)
library(lme4)
data("ChickWeight", package = "datasets")

ChickWeight <- as_tibble(ChickWeight) %>%
  mutate(Time_6p = pmax(0, Time - 6))

chickm4 <- lmer(weight ~ Time * Diet + (Time | Chick), data = ChickWeight)
chickm5 <- lmer(weight ~ Time + Time_6p + Time:Diet + (Time | Chick), data = ChickWeight)
chickm6 <- lmer(weight ~ Time + Time_6p + Time:Diet + Time_6p:Diet + (Time | Chick), data = ChickWeight)
chickm7 <- lmer(weight ~ Time*Diet + Time_6p*Diet + (Time | Chick), data = ChickWeight)

anova(chickm4, chickm5, chickm6, chickm7) # not sure if it actually does a refit with ML

pred_df <- ChickWeight %>%
  distinct(Time, Time_6p, Diet) %>%
  mutate(pred4 = predict(chickm4, newdata = ., re.form = ~0)) %>%
  mutate(pred6 = predict(chickm6, newdata = ., re.form = ~0))

ChickWeight %>%
  mutate(Chick = forcats::fct_reorder(as.character(Chick), sample(row_number()))) %>% 
  ggplot(aes(x = Time)) + 
  geom_line(aes(y = weight, color = Chick)) + 
  geom_line(data = pred_df, aes(y = pred4, linetype = 'pred4'), size = 1) +
  geom_line(data = pred_df, aes(y = pred6, linetype = 'pred6'), size = 1) + 
  geom_smooth(aes(y = weight), method = 'loess') +
  scale_color_discrete(guide = FALSE) +
  facet_wrap(~Diet)

$\endgroup$
1
  • $\begingroup$ thanks @user3408956 i really appreciate the time you put into this. The thing is this is related to my course work and whilst i understand why the weight of the chick between 0 to ~6 days isnt much affected by Diet, and your 'knotted' line is much better, I don't have the skills to do these sorts of plots (nor is it currently expected of me, hence why I unfortunately can't really use the great graph you've produced.) I suppose the best thing for me to do is point out the limitations of my chickm4 model, which is presumably the poor intercept and poor predictions for Time between 0 and 6, ish? $\endgroup$ Commented Jan 2, 2020 at 16:02

Not the answer you're looking for? Browse other questions tagged or ask your own question.