I am using linear mixed effect models to study the effect of a physical therapy treatment on two different populations of patients (Adolescents and Adults). The treatment is the same for both populations, and for each subject the (same) treatment is applied to both legs (Left/Right). The question is whether the treatment has a different effect on the two populations and the two legs (I am also interested in their potential interaction).

I therefore have two factors as fixed effects: age (factors with levels: Adolescents and Adults), and leg (factors with levels: Left and Right), with possibile interaction (age*leg). The random effect is subject-id. However, I am uncertain on how to design the model: The two populations are independent to each other, but the two legs are not (i.e. two legs within the same subject).

So far I used a model where I nested leg within subject-id:

options(contrasts = c("contr.sum","contr.poly"))
lCtr  <- lmeControl(maxIter = 1000, niterEM = 500, msVerbose = FALSE, opt = 'optim')

linM <- lme(dat ~ leg*age, random = ~1|sbjID/leg, data=dat_trf, na.action=na.omit, method = "ML", control=lCtr )

Is this design correct?

The summary provides the following results:

> summary(linM)

Linear mixed-effects model fit by maximum likelihood
  Data: dat_trf 
       AIC      BIC    logLik
  118.2452 142.7825 -52.12259

Random effects:
 Formula: ~1 | sbjID
StdDev:   0.3601248

 Formula: ~1 | leg %in% sbjID
        (Intercept)   Residual
StdDev:    0.144943 0.08273855

Fixed effects:  dat ~ leg * age 
                  Value  Std.Error  DF  t-value p-value
(Intercept)   1.6142702 0.03496504 121 46.16813  0.0000
leg1          0.0232176 0.01088832 121  2.13234  0.0350
age1         -0.0514979 0.03496504 121 -1.47284  0.1434
leg:age1     -0.0089302 0.01088832 121 -0.82017  0.4137
             (Intr)   leg1  age1  
leg1          0.000              
age1         -0.171  0.000       
leg1:age1     0.000 -0.171  0.000

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-0.942245288 -0.249664440 -0.002136358  0.199769397  1.004654898 

Number of Observations: 246
Number of Groups: 
            sbjID leg %in% sbjID 
              123            246 

The number of groups seems to be correct: There are indeed 123 subjects (51 adolescents and 72 adults), for a total of 246 measurements (two leg per subjects). However, I suspect that something is wrong, because I obtain a relatively small p-value for leg1 when the data look like in the following boxplots:

enter image description here

This weird result is confirmed by anova:

> anova.lme(linM,type="marginal")
            numDF denDF   F-value p-value
(Intercept)     1   121 2131.4958  <.0001
leg             1   121    4.5469  0.0350
age             1   121    2.1693  0.1434
leg:age         1   121    0.6727  0.4137

Is there anything wrong?

The assumptions of normality and constant variance of residuals are met. Thanks in advance!


2 Answers 2


This all seems sensible to me. I would explain the surprising result of a significant (at $\alpha = 0.05$) $p$-value by noting that your plot shows the whole population. It's not unreasonable that average effect of going from left to right leg (0.02 * 2 = 0.04 in whatever units your response has), while small, is still significant across the population since your data set is reasonably large. Squinting at the picture, I can convince myself that more of the lines are decreasing from left to right than increasing (since you are using sum-to-zero contrasts, a parameter estimate of 0.02 means the response for left legs is 0.02 units higher than the population mean on average).

Also note that these effects, while significant, are small relative to the standard deviations at all levels (ranging from 0.08 to 0.36), consistent with the noisy-looking picture.

To help convince yourself, try drawing boxplots of (left-right) for each subject, subdivided by age, and including the boxplot notches that indicate approximate 95% CIs for the median. Computing (left-right) will eliminate the variation in the intercept; showing the notches will contrast the inference on the average difference between groups (analogous to a standard error of the mean) with the overall variation shown by the boxes ($\approx$ IQR, analogous to the standard deviation).

  • $\begingroup$ Thanks! Would running multiple-comparison post-hoc tests be wrong in this case (because I should just conclude that there is a significant effect of leg, independently of age)? $\endgroup$
    – Cristiano
    May 9, 2023 at 14:18
  • $\begingroup$ I wouldn't bother with post hoc tests, but that's a matter of taste. Follow the conventions of your field. $\endgroup$
    – Ben Bolker
    May 9, 2023 at 14:29

It seems like you have only one measurement for each leg of every patient. Hence the leg %in% sbjID matches with the residual. Use only (1|sbjID) as random effect.

The fixed effect structure assumes that there might be a systematic difference between the right and left legs. Is there any biological reason why you would expect that?

  • $\begingroup$ Thanks! I am not sure I follow why I should use (1|sbjID) as random effect. Could you elaborate? To answer your question: Yes, it is because of dominance! $\endgroup$
    – Cristiano
    May 9, 2023 at 15:33
  • $\begingroup$ @Thierry, good catch about confounding. (Fixed-effect estimates are probably fine though.) (lmer would have complained about this; in lme it will only become apparent if you try to run intervals() on the fit ... $\endgroup$
    – Ben Bolker
    May 10, 2023 at 16:08
  • $\begingroup$ @Cristiano: you need repeated measures for each leg of the subject in order to fit ~1|sbjID/leg, hence you must simplify the random effects to ~1|sbjID. If dominance is important, then code the legs as dominant/non dominant instead of left/right. That makes more sense IMHO. $\endgroup$
    – Thierry
    May 11, 2023 at 7:27
  • $\begingroup$ @Thierry: My understanding was that, theoretically, I should account for the fact that the measurements for left and right legs come from the same subject. Are you saying that I should do that, but I cannot do it because I do not have enough data (only one measurements per leg per subject)? Or you are saying that I should not do it? Also, I am still missing why, in practice, I cannot fit ~1|sbjID/leg if I have only one measurement per leg (@Ben Bolker: intervals works fine). Apologies: I realize I am may be missing something basic here. $\endgroup$
    – Cristiano
    May 11, 2023 at 15:10
  • $\begingroup$ I would expect intervals() to either give you NA values for the random-effects variances, or unrealistically/ridiculously wide confidence intervals, if you only have one measurement per subject per leg. Why this doesn't/shouldn't work: there is a unique random effect for every observation (subject × leg) specified by the random effects component, and a unique random effect for every observation specified by the residual variance term of a Gaussian, so these should be confounded/jointly unidentifiable. $\endgroup$
    – Ben Bolker
    May 11, 2023 at 16:29

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.