How to choose random effect in a repeated measures mixel model

Question

I'm trying to analyze these data with a mixed effect model.

I would like to evaluate the association between Score (an ecographic score) and val 1 and val2 (clinical parameters - continuous variables) in different subjects and times. The three variables (Score, val1, val2) have been measured at day 0, 1 and 2.

I'm trying this model:

• Day as fixed effect
• Subject as random effect
• Val1 and Val2 as random effect considered variable between subjects

I've considered correlated random intercept and slope.

Models using lme4: First model:

m1<-lmer(Score ~ Day + (1 + Day |Subject) + (1 + val1 |Subject)+(1 + val2 |Subject), mydata, REML = FALSE)
summary(m1)

First output:

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Score ~ Day + (1 + Day | Subject) + (1 + val1 | Subject) + (1 +      val2 | Subject)
   Data: mydata

     AIC      BIC   logLik deviance df.resid 
    93.3    104.0    -34.6     69.3        6 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1801 -0.4187  0.1480  0.6878  1.6237 

Random effects:
 Groups    Name        Variance  Std.Dev.  Corr 
 Subject   (Intercept) 0.000e+00 0.0000000      
           Day         5.897e-06 0.0024285  NaN 
 Subject.1 (Intercept) 4.397e-07 0.0006631      
           val1        2.713e-08 0.0001647 -1.00
 Subject.2 (Intercept) 6.871e+00 2.6211653      
           val2        1.305e-01 0.3612815 -1.00
 Residual              2.528e+00 1.5899007      
Number of obs: 18, groups:  Subject, 6

Fixed effects:
            Estimate Std. Error t value
(Intercept)  5.42316    0.61129   8.872
Day          0.04776    0.46063   0.104

Correlation of Fixed Effects:
    (Intr)
Day -0.739
convergence code: 0
boundary (singular) fit: see ?isSingular

Second model:

m2<-lmer(Score ~ Day + (1 + Day |Subject) + (1 |val1)+(1 |val2), mydata, REML = FALSE)
summary(m2)

Second output:

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: Score ~ Day + (1 + Day | Subject) + (1 | val1) + (1 | val2)
   Data: mydata

     AIC      BIC   logLik deviance df.resid 
    85.6     92.7    -34.8     69.6       10 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0400 -0.6711  0.3046  0.6761  1.6340 

Random effects:
 Groups   Name        Variance  Std.Dev.  Corr
 Subject  (Intercept) 8.794e-03 9.378e-02     
          Day         4.501e-02 2.122e-01 1.00
 val1     (Intercept) 0.000e+00 0.000e+00     
 val2     (Intercept) 2.714e-10 1.647e-05     
 Residual             2.674e+00 1.635e+00     
Number of obs: 18, groups:  Subject, 6; val1, 4; val2, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept)  5.36111    0.61064   8.779
Day          0.08333    0.47995   0.174

Correlation of Fixed Effects:
    (Intr)
Day -0.749
convergence code: 0
boundary (singular) fit: see ?isSingular

My data:

mydata<- data.frame(Subject  = c(1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6), 
                    Day=c(0,0,0,0,0,0,1,1,1,1,1,1,2,2,2,2,2,2), 
                    val1=c(5,5,5,1, 6, 5, 5, 5, 5, 1, 5, 6, 5, 5, 5, 1, 4, 6), 
                    val2=c(7,7, 7, 4, 6, 7, 7, 7, 7, 4, 6, 6, 7, 7, 7, 4, 5, 6), 
                    Score=c(8, 6, 6, 4, 6, 6, 2, 3, 3, 4, 6, 7, 5, 5, 7, 5, 8, 7))

I've considere a linear model. Here are the plots:

What do you think? How I can modify my code for a better definition of the model and for a better answer to my question? Thank you!

Answer 1

Fascinating thread. Further to Erik's excellent response, I would like to suggest that you also look into fitting your linear mixed effects models using the gamlss function in the gamlss package in R.

Here are some examples of models you can consider fitting:

install.packages("gamlss")

library(gamlss)


## linear mixed effects model with  
## a fixed intercept and 
## a random intercept for Subject 
cv0 <- gamlss(Score ~ 1 + 
          re(random = ~1|Subject), 
          data = mydata) 

cv0
summary(cv0)
getSmo(cv0)

## linear mixed effects model with  
## a fixed effect for Day and 
## a random intercept for Subject 
cv1 <- gamlss(Score ~ 1 + Day + 
              re(random=~1|Subject), 
              data = mydata) 
cv1
summary(cv1)
getSmo(cv1)

## linear mixed effects model with  
## fixed effects for Day, val1 and val2 and 
## a random intercept for Subject 
cv2 <- gamlss(Score ~ 1 + Day + val1 + val2 + 
              re(random=~1|Subject), 
              data = mydata) 

cv2 
summary(cv2)
getSmo(cv2)

## linear mixed effects model with  
## fixed effects for Day, val1 and val2,  
## a random intercept for Subject and
## a random slope for Day
cv3 <- gamlss(Score ~ 1 + Day + val1 + val2 + 
              re(random=~ 1 + Day|Subject, 
                 opt="optim", numIter=100), 
              data = mydata) 

cv3 
summary(cv3)
getSmo(cv3)

## linear mixed effects model with  
## fixed effects for Day, val1 and val2,  
## a random intercept for Subject,
## a random slope for Day and a random slope for val1
cv4 <- gamlss(Score ~ 1 + Day + val1 + val2 + 
              re(random=~Day + val1|Subject, 
                 opt="optim", numIter=100), 
              data = mydata) 

cv4 
summary(cv4)
getSmo(cv4)
 


## linear mixed effects model with  
## fixed effects for Day, val1 and val2,  
## a random intercept for Subject,
## a random slope for Day and a random slope for val2
cv5 <- gamlss(Score ~ 1 + Day + val1 + val2 + 
              re(random=~Day + val2|Subject, 
                 opt="optim", numIter=100), 
              data = mydata) 

cv5 
summary(cv5)
getSmo(cv5)

Note that you will NOT be able to fit the model:

## linear mixed effects model with  
## fixed effects for Day, val1 and val2,  
## a random intercept for Subject,
## a random slope for Day,
## a random slope for val1 and 
## a random slope for val2
cv6 <- gamlss(Score ~ 1 + Day + val1 + val2 + 
              re(random=~Day + val1 + val2|Subject, 
                 opt="optim", numIter=100), 
              data = mydata) 

cv6 
summary(cv6)

as R will throw an error: fewer observations than random effects in all level 1 groups.

On the plus side, you will be able to fit the sequence of models cv0 to cv5 and compare these models in terms of their (generalized) AIC value to determine which model is preferred for your data:

 GAIC(cv0, cv1, cv2, cv3, cv4, cv5)

 >  GAIC(cv0, cv1, cv2, cv3, cv4, cv5)
      df      AIC
cv0 1.000000 71.63080   <- lowest AIC value
cv5 4.000062 72.46398      so model cv0 is preferred
cv4 4.000027 72.46398
cv3 4.000031 72.46398
cv2 4.000000 72.46398
cv1 2.000000 73.60104

You can also use:

GAIC.table(cv0, cv1, cv2, cv3, cv4, cv5)

to compute the (generalized) AIC with different penalties k.

You can interrogate the preferred model cv0 using commands such as:

## summary for preferred model
summary(cv0)

## random effects for preferred model
getSmo(cv0)

## model diagnostic plots for preferred model
plot(cv0)

 ## worm plot for preferred model
 wp(cv0)

There is not much difference in (generalized) AIC models between cv0 and the next best models, cv2 to cv5 - less than 1 point on the (generalized) AIC scale. So you may decide to still report one of these models, while being transparent about the fact that it has the next lowest (generalized) AIC value and very small standard deviations for the random effect(s) it includes. (You would probably be most interested in reporting model cv3?)

PS: Your val1 and val2 values look a bit odd to me in the data.

library(dplyr)
library(magrittr)

mydata <- mydata %>% 
          arrange(Subject)

mydata

mydata$Subject <- factor(mydata$Subject)

> mydata
Subject Day val1 val2 Score
1        1   0    5    7     8
2        1   1    5    7     2
3        1   2    5    7     5

4        2   0    5    7     6
5        2   1    5    7     3
6        2   2    5    7     5

7        3   0    5    7     6
8        3   1    5    7     3
9        3   2    5    7     7

10       4   0    1    4     4
11       4   1    1    4     4
12       4   2    1    4     5

13       5   0    6    6     6
14       5   1    5    6     6
15       5   2    4    5     8

16       6   0    5    7     6
17       6   1    6    6     7
18       6   2    6    6     7

The values of val1 and val2 are constant across Days for Subject 1, Subject 2, Subject 3 and Subject 4 (suggesting that val1 and val2 are Subject-level predictor variables, whose values remain constant across Days for the same subject).

However, the values of val1 and val2 for Subject 5 and Subject 6 do change across Days within subject. Not sure this is by accident or by design? A quick quality check of the data might be warranted before proceeding with your modelling.

For more information on the gamlss& package of R, see: https://www.gamlss.com/the-r-packages/

Addendum

The article Analysis of longitudinal multilevel experiments using GAMLSSs by Thomas et al. (https://arxiv.org/abs/1810.03085) mentions that it is possible to get the anova() function to work with a gamlss model provided the model is specified in a particular way using the re() function. For example, if the linear mixed effects model is specified like so:

library(gamlss)

fm <- gamlss(Score ~ re(fixed = ~ 1 + Day + val1 + val2, 
                        random = ~ 1|Subject, 
                        method = "ML"), 
       data = mydata)
fm

then applying the function getSmo() to the model m will convert m into an lme object so that the anova() function from the nlme package will work on it:

library(lme)

anova(getSmo(fm))

The anova() function will produce output like the one shown below:

> anova(getSmo(fm))
            numDF denDF   F-value p-value
(Intercept)     1     9 0.0000000  1.0000
Day             1     9 0.0308174  0.8645
val1            1     9 1.6835327  0.2267
val2            1     9 2.9404394  0.1205

It is also possible to use the Anova() function from the car package on the model m:

 library(car)

 Anova(getSmo(fm))

I tried to see if anova() works with two nested linear mixed effects models fitted with the ML method (i.e., Maximum Likelihood), but that unfortunately doesn't seem to work:

fm1 <- gamlss(Score ~ re(fixed = ~ 1 + Day, 
                        random = ~ 1|Subject, 
                        method = "ML"), 
       data = mydata)
fm1

fm2 <- gamlss(Score ~ re(fixed = ~ 1 + Day + val1, 
                        random = ~ 1|Subject, 
                        method = "ML"), 
       data = mydata)
fm2

fm1.lme <- getSmo(fm1)
fm2.lme <- getSmo(fm2)

anova(fm1.lme, fm2.lme)

Indeed, the last anova() command produces the following error: Error in eval(x$call$fixed) : object 'fix.formula' not found.

The gamlss package has an LR.test() function for comparing nested models, but that doesn't seem to work for models with random effects specified as above. Trying something like:

LR.test(fm1, fm2)

produces the error message: Error in LR.test(fm1, fm2) : The difference in df's is negative. Are the models nested?.

Perhaps it's worth e-mailing the gamlss package authors to ask for further clarifications on whether it is possible to compare linear mixed effects models via the anova() function in the nlme package or the LR.test() function in the gamlss package?

If you just want to compare a model with and without a variable of interest, perhaps use the drop1() function in the gamlss package? Something like this (note that the model is now specified with the fixed effects outside of the re() function for the drop1() function to work):

library(gamlss)

fm <- gamlss(Score ~  Day + val1 + val2 + re(random = ~ 1|Subject), 
       data = mydata)
fm

drop1(fm, what="mu")

In case this helps, see also the handout available https://docs.ufpr.br/~taconeli/CE06218/Ex5.pdf. Other commands mentioned in this handout that will work on the getSmo() version of a gamlss model with random effects are as follows:

## summary
summary(getSmo(fm))

## residual level=1 plot
plot(getSmo(fm))

## fixed effect estimates
fixef(getSmo(fm))

## random effect estimates
ranef(getSmo(fm))

## plot of random effect estimates
plot(ranef(getSmo(fm)))

## fitted coefficients
coef(getSmo(fm))

## confidence intervals
intervals(getSmo(fm))

Maybe check this reference out too:
https://teses.usp.br/teses/disponiveis/11/11134/tde-06042018-150012/publico/Gustavo_Thomas.pdf.

Answer 2

Your models are over-parameterized presently. Start simple and then build up. The first question is whether you have variation between subjects on the outcome. Said differently, are Scores correlated within Subjects? To answer this, look at the simplest longitudinal model that only has a random intercept for Subject and no predictors:

cv0 <- lmer(Score ~ 1 + (1|Subject), mydata) 

REML criterion at convergence: 69.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0007 -0.6819  0.3054  0.7500  1.4960 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 0.02593  0.161   
 Residual             2.94444  1.716   
Number of obs: 18, groups:  Subject, 6

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)   5.4444     0.4098 5.0000   13.29 4.32e-05 ***
---

This is a big clue to what is going on. You have almost no between Subject variation in your outcome. The proportion of variance in Score that is between subjects is less than 1%:

.02593/(.02593+2.9444) # Formula for calculating an ICC (intraclass correlation coefficient)
[1] 0.00872967

Once you add Day as a fixed effect and val1 and val2, you remove all between-subject variance in Score. This means that you do not need a multilevel model. You can use a simpler OLS model:

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Score ~ 1 + Day + val1 + val2 + (1 | Subject)
   Data: mydata

REML criterion at convergence: 64.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.86667 -0.35792 -0.03183  0.55532  1.77156 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 0.000    0.000   
 Residual             2.704    1.644   
Number of obs: 18, groups:  Subject, 6

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)   
(Intercept)  7.92843    2.64103 14.00000   3.002  0.00951 **
Day         -0.01724    0.48084 14.00000  -0.036  0.97190   
val1         0.92173    0.42918 14.00000   2.148  0.04973 * 
val2        -1.06432    0.62068 14.00000  -1.715  0.10843   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr) Day    val1  
Day  -0.330              
val1  0.495 -0.103       
val2 -0.880  0.154 -0.827
convergence code: 0
boundary (singular) fit: see ?isSingular

One caveat - if you actually have more subjects than you shared here, then you may have more between-Subject variance to play with. But at least with what you've provided here, there is an almost 0 correlation of Score within Subject, thus obviating the need for a mixed effects model.