Assessing the need for multi level model

Question

I am trying to asses the need for multi level model. I know have to use a multi level model but I am doing this any way to include it as background information.

My study was a 12 week long study involving two diets (diet). Outcome measurements such as weight, waist circumference etc. were taken at time (time) = 0,6,and 12 week intervals. (Code) stands for subjects.

I am using Discovering Statistics Using R by Andy Field, Jeremy Miles and Zoe Fields. This book suggests to build a base model that is "intercept only" and then a model that is "random intercept only" and to compare the outputs in anova to check whether the random intercept model improves the model.

These are the models I created to asses the need for multi level model.

1. Predict weight from intercept only
2. Predict weight from intercept only but let intercepts vary across code (subjects)
3. Predict weight from intercept only but let intercepts vary across diet
4. Predict weight from intercept only but let intercepts vary across time
5. Predict weight from intercept only but let intercepts vary across time and code (subject)

6. Predict weight from intercept only but let intercepts vary across diet and code (subject)

   interceptOnly <-gls (weight ~ 1, data = dat2, method = "ML")
   randomInterceptOnly <-lme(weight ~ 1, data = dat2, random = ~1|code, method = "ML")
   randomInterceptOnly <-lme(weight ~ 1, data = dat2, random = ~1|diet, method = "ML")
   randomInterceptOnly <-lme(weight ~ 1, data = dat2, random = ~1|time, method = "ML")
   randomInterceptOnlytimecode <- lme(weight ~ 1, data = dat2, random = ~time|code, method = "ML")
   randomInterceptOnlydietcode <- lme(weight ~ 1, data = dat2, random = ~diet|code, method = "ML")
   

I then used anova() to determine which model provides an improvement.

anova(interceptOnly,randomInterceptOnly,randomInterceptOnlycode,randomInterceptOnlydiet,randomInterceptOnlytime,randomInterceptOnlytimecode,randomInterceptOnlycodediet)

                           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
interceptOnly                   1  2 905.1530 910.5900 -450.5765                        
randomInterceptOnly             2  3 691.2088 699.3643 -342.6044 1 vs 2 215.9442  <.0001
randomInterceptOnlycode         3  3 691.2088 699.3643 -342.6044                        
randomInterceptOnlydiet         4  3 890.7966 898.9521 -442.3983                        
randomInterceptOnlytime         5  3 907.1530 915.3085 -450.5765                        
randomInterceptOnlytimecode     6  5 648.5658 662.1583 -319.2829 5 vs 6 262.5873  <.0001
randomInterceptOnlycodediet     7  5 695.1268 708.7193 -342.5634                      

(Models 2 and 3 are redundant)

Based on this it seems that 'random intercept only code' and 'random intercept only time | code' provide significant improvement to the model. Therefore, a mixed effect model is warranted.

Is this a correct interpretation?

Thanks for taking the time to read this!

Answer 1

A couple of points:

1. Because the choice of the random-effects can be affected by the chosen fixed effects, what it is typically done is that you select your random effect with a general/flexible model for your fixed effects. In your case, I'd say that this entails including in the fixed effects the time and diet effect, and their interaction.

2. Random effects are typically included to account for the correlation of the measurements within a group/cluster. In your case, the groups/clusters are the subjects for whom you have repeated measurements over time. Hence, it would make more sense to include random effects only for the code variable. A potential series of models to consider is:

   fm0 <- gls(weight ~ time * diet, data = dat2) fm1 <- lme(weight ~ time * diet, data = dat2, random = ~ 1 | code) fm2 <- lme(weight ~ time * diet, data = dat2, random = ~ time | code)

However, because in your case you only have three time points, you may also consider including time as a categorical variable in your fixed effects and/or fitting a completely unstructured covariance matrix using gls().