0
$\begingroup$

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!

$\endgroup$
7
  • $\begingroup$ In your code, you have the same model 3 times and you don't have all teh models that are in your output, so something is wrong there. $\endgroup$
    – Peter Flom
    Sep 4, 2018 at 0:12
  • $\begingroup$ Building your model in such a stepwise fashion will invalidate p-values, confidence intervals etc. from the final model. Why not prespecify the standard MMRM for this type of situation? I.e. Change or value at each post baseline visit = treatment + visit + baseline + visit*treatment + visit *baseline with an unstructured (simplifies to just a random effect + allowing for different variance at each visit when you just have two part baseline visits) covariance matrix describing the correlation of visits within a patient? $\endgroup$
    – Björn
    Sep 4, 2018 at 6:47
  • $\begingroup$ @Bjorn I see. The model you are describing is that a ANCOVA where baseline measures are used as a fixed effect? Can this be carried out within lme? $\endgroup$
    – DiscoR
    Sep 4, 2018 at 16:14
  • $\begingroup$ @Björn I am trying to run the analysis you suggested. However I am not sure how to use only post baseline values as outcome (outcome at time 6 weeks and 12 weeks). What would the mixed effect code look like for this? $\endgroup$
    – DiscoR
    Sep 5, 2018 at 16:39
  • 1
    $\begingroup$ In my experience there usually is. $\endgroup$
    – Björn
    Sep 5, 2018 at 17:04

1 Answer 1

3
$\begingroup$

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().

$\endgroup$
3
  • $\begingroup$ @DimitrisRizopoulas Thanks for this response. I understand this better now. Thanks for this clarification. What do you mean when you say to use time as a categorical variable in the fixed effects, wouldn't time * diet automatically include in the individual effects of time and diet separately? Please let me know if I am missing something here. Also I would appreciate any pointers towards how to carry out the second part of your suggestion 'fitting a completely unstructured covariance matrix using gas()' $\endgroup$
    – DiscoR
    Sep 4, 2018 at 16:10
  • $\begingroup$ If the time variable in dat2 is numeric (you can see that by doing str(dat2$time)), then by including in your model you assume a linear relationship between time and your outcome weight. If you transform it to a factor (e.g., using the factor() function), and include it in the model, then it will be treated as a categorical variable, and you will obtain a different coefficient per time point. The completely unstructured covariance matrix can be fitted by gm <- gls(weight ~ time * diet, data = dat2, correlation = corSymm(form = ~ 1 | code), weights = varIdent(form = ~ 1 | time)) $\endgroup$ Sep 4, 2018 at 19:26
  • $\begingroup$ thank that worked! my final question is if the output of the gls model can be compared with lme4 by anova? To assess which provides better fit? $\endgroup$
    – DiscoR
    Sep 4, 2018 at 19:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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