How to choose random effect in a repeated measures mixel model

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!

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.

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.

• This is a great post, Isabella (+1)! Can you clarify why you suggest using the gamlss package? – Erik Ruzek Oct 5 '20 at 21:59
• Your post is great too, Erik! One of the reasons I suggested the gamlss package is for pedagogical reasons - unlike the lmer() function in the lme4 package, the gamlss() function in the gamlss package can fit a series of models and allow the user to compare them without fretting about the warning boundary (singular) fit posted by lmer. – Isabella Ghement Oct 5 '20 at 22:10
• Also, if the study design used Subject as a random grouping factor, I think it is best that the model formulated for the data reflects that design by including (at a minimum) a random Subject effect, as well as fixed effects for the predictor variables of interest in the study. – Isabella Ghement Oct 5 '20 at 22:12
• @Isabella I didn't know the gamless& package, it's so interesting! So useful the possibility to compare different models! With lme4 package to obtain p-value for a variable I used the anova function comparing different models with/without the variable of interest. With gamless& how can I do this? Thank you again!!! – ArTu Oct 6 '20 at 11:16
• @IsabellaGhement You're an inexhaustible source of invaluable information!!!! Can I ask you the difference in p-value obtained with anova(getSmo(cv6)) and with summary(getSmo(cv6))? I obtain overlapping p-value except for Day . Thank you so so much! You've made my day! – ArTu Oct 7 '20 at 12:35

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.