Solved – Three-level partially nested model

lme4-nlmemultilevel-analysisnested datarrepeated measures

I am modeling change over time in group psychotherapy subjects using R and lme4.
My data have the following structure:

subject (id)
time (code 1-10 for equally spaced repeated measures)
outcome (for every repeated measure)
treatment (0/1 for psychotherapy/waiting list control)

My first two-level model with random slopes and intercepts works well and is simple:

lmer(outcome ~ time * treatment + (time | subject), data=data, REML=FALSE)

Now I was wondering if I should use a three-level partially nested model because the group psychotherapy subjects are nested within therapist (there where several therapist providing the treatment), but controls were non-nested. I am guessing I should account at least for main effects of therapists as argued by De Jong, Moerbeek & Van der Leeden (2010).

De Jong, K., Moerbeek, M., & Van der Leeden, R. (2010). A priori power analysis in longitudinal three-level multilevel models: an example with therapist effects. Psychotherapy Research, 20(3), 273-284.

I have found a very useful resource in the following link under "Partially nested models":
http://rpsychologist.com/r-guide-longitudinal-lme-lmer

The author gives the following code for a model virtually identical to the one I would like to test:

lmer(outcome ~ time * treatment + (1 | group:subject) + (0 + time | therapist:subject) + (0 + time:treatment | group) + (0 + treatment | group), data=data)

The data he presents is virtually identical but he adds in the model the "group" variable. I don't understand why this is because treatment/control grouping is the same as treatment/nontreatment grouping. If a subject received treatment than he is in the experimental group, if not he in the control group. How would you write this three-level partially nested model? I feel confused. Thank you.

Best Answer

Question was eloquently answered by Thierry on Stackoverflow.

I post his response, maybe it will be helpful for some one.

Keep it simple. Just add a dummy therapist none to the subjects without therapist. Then fit the model below.
lmer(outcome ~ time * treatment + (time | therapist/subject), data=data)
therapist none is confounded with treatment waiting list. Therapist is a random effect and thus penalised. Treatment is a fixed effect and not penalised. Hence all information will go to the treatment effect for waiting list and the therapist effect for none will be zero.

Related Solutions

r – Correct Specification of Longitudinal Model in lme4 for R

I'll imagine a concrete example, with more context, to make things easy. Assume you measure the score on test of 3k students of 200 schools and you measured each student at 4 time points (say, at each quarter). You have a covariate at student level that doesn't vary by time (like sex), that you called pred1.obs and a covariate by school that vary by time (say the number of meetings between teachers and parents until that moment in time). If this example resembles your study, than I think you have to set up a three level model (individual level, group level and time level for the groups): i = 1 ... 3000 individuals t = 1... 4 periods g = 1... 200 groups

The model would be:

y_i ~ N(a + b_[groups_g] + b.ind*pred.obs1_i, sigma^2) # 1st level
b_g = N(gamma + gamma_[time] + gamma.g[time_t]*pred2.grp, sigma.b^2) # 2nd level
gamma.g_t = N(0, sigma.gamma^2) # 3rd level

Note that you would have the slope at the second level (group level) varying by time, which makes sense, since you expect that the effect of schools on the perfomance of students may vary by time, depending of the valu of the covariate at the level of schools. I'm not that sure how to estimate this with lmer (I know how to estimate a Bayesian model using WinBugs or Jags, calling them by R). In any case, here is my suggestion.

In lme4, I'd try: First, expand pred2.grp (the covariate at group level that vary by time) to the individual level, then you would have repetead measures by individuals at the group and time level. Then:

lmer(outcome ~ pred1.obs + pred2.grp + (1|group))

Mixed Effects Models – Analyzing Longitudinal Data with Nested Random Effects

Assuming that each subject is only in one group, you have a nested design. Conceptually, it makes more sense to treat group as a fixed effect. As you have only three groups, it wouldn't make much sense statistically speaking anyways because you'd be asking the software to estimate a variance for group assuming a normal distribution from only three observations.

Your sample size is quite low which might be a problem. Nevertheless, I suggest starting with a simple random intercept model:

lmer(outcome~day + group + (1|subject), data = dat)

This model fits a global intercept which is simply the intercept for the reference group, deviations from that intercept for the remaining groups, a single slope for the effect of day and a random intercept for subject. Hence, this model assumes that each group has the same trajectory over time, the same slope but different intercepts. To allow each group their own slope, you could fit the following model with an interaction between day and group:

lmer(outcome~day*group + (1|subject), data = dat)

Lastly, you could also allow for random slopes:

lmer(outcome~day*group + (day|subject), data = dat)

All these models assume a linear relationship between the outcome and day which might be unrealistic. If you suspect nonlinear relationships, you could easily accomodate those by using polynomials or (my recommendation) restricted cubic splines (aka natural splines) with, say, 3 knots. Finally, have a look at this post which goes into more details about such longitudinal models.

Best Answer

Related Solutions

r – Correct Specification of Longitudinal Model in lme4 for R

Mixed Effects Models – Analyzing Longitudinal Data with Nested Random Effects

Related Question