Multi-level / mixed-effect / hierarchical linear regression. Model and interpretation advice? 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. 
 A: 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.

A: 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. 



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)


