Jonas Lang

News Short Vita Research Team Publications Contact

18.06.2021 Is it possible to fit group-specific random slopes in nlme and lme4?

A question that recently surfaced in discussions is whether it is possible to fit group-specific random slopes in nlme and lme4. The researcher in the example was interested in testing the idea that two groups would show different amounts of slope change over time. The model can be fit in both R packages using the code below.


library(lme4)
sleepstudy$Dummy<-ifelse(as.numeric(sleepstudy$Subject)>9,1,0)
testmod<-lmer(Reaction~0
+dummy(Dummy,"0")+dummy(Dummy,"1")+
+Days:dummy(Dummy,"0")+Days:dummy(Dummy,"1")
+(0+dummy(Dummy,"0")+Days:dummy(Dummy,"0")|Subject)
+(0+dummy(Dummy,"1")+Days:dummy(Dummy,"1")|Subject),
data=sleepstudy)

testmod

library(nlme)
testmodb<-lme(Reaction~0
+dummy(Dummy,"0")+dummy(Dummy,"1")+
+Days:dummy(Dummy,"0")+Days:dummy(Dummy,"1"),
random=list(Subject=pdBlocked(list(
pdSymm(~ 0+dummy(Dummy,"0")+Days:dummy(Dummy,"0")),
pdSymm(~ 0+dummy(Dummy,"1")+Days:dummy(Dummy,"1"))
))),data=sleepstudy)

testmodb


When one runs this code in R, one gets the following results (identical in both packages):


> library(lme4)
> sleepstudy$Dummy<-ifelse(as.numeric(sleepstudy$Subject)>9,1,0)
> testmod<-lmer(Reaction~0
+ +dummy(Dummy,"0")+dummy(Dummy,"1")+
+ +Days:dummy(Dummy,"0")+Days:dummy(Dummy,"1")
+ +(0+dummy(Dummy,"0")+Days:dummy(Dummy,"0")|Subject)
+ +(0+dummy(Dummy,"1")+Days:dummy(Dummy,"1")|Subject),
+ data=sleepstudy)
>
> testmod
Linear mixed model fit by REML ['lmerMod']
Formula:
Reaction ~ 0 + dummy(Dummy, "0") + dummy(Dummy, "1") + +Days:dummy(Dummy,
"0") + Days:dummy(Dummy, "1") + (0 + dummy(Dummy, "0") +
Days:dummy(Dummy, "0") | Subject) + (0 + dummy(Dummy, "1") +
Days:dummy(Dummy, "1") | Subject)
Data: sleepstudy
REML criterion at convergence: 1726.288
Random effects:
Groups Name Std.Dev. Corr
Subject dummy(Dummy, "0") 28.080
dummy(Dummy, "0"):Days 6.439 0.09
Subject.1 dummy(Dummy, "1") 23.203
dummy(Dummy, "1"):Days 3.567 0.06
Residual 25.592
Number of obs: 180, groups: Subject, 18
Fixed Effects:
dummy(Dummy, "0") dummy(Dummy, "1")
252.292 250.519
dummy(Dummy, "0"):Days dummy(Dummy, "1"):Days
7.388 13.546
>
> library(nlme)
> testmodb<-lme(Reaction~0
+ +dummy(Dummy,"0")+dummy(Dummy,"1")+
+ +Days:dummy(Dummy,"0")+Days:dummy(Dummy,"1"),
+ random=list(Subject=pdBlocked(list(
+ pdSymm(~ 0+dummy(Dummy,"0")+Days:dummy(Dummy,"0")),
+ pdSymm(~ 0+dummy(Dummy,"1")+Days:dummy(Dummy,"1"))
+ ))),data=sleepstudy)
>
> testmodb
Linear mixed-effects model fit by REML
Data: sleepstudy
Log-restricted-likelihood: -863.1442
Fixed: Reaction ~ 0 + dummy(Dummy, "0") + dummy(Dummy, "1") + +Days:dummy(Dummy, "0") + Days:dummy(Dummy, "1")
dummy(Dummy, "0") dummy(Dummy, "1")
252.291581 250.518629
dummy(Dummy, "0"):Days dummy(Dummy, "1"):Days
7.388489 13.546083

Random effects:
Composite Structure: Blocked

Block 1: dummy(Dummy, "0"), dummy(Dummy, "0"):Days
Formula: ~0 + dummy(Dummy, "0") + Days:dummy(Dummy, "0") | Subject
Structure: General positive-definite
StdDev Corr
dummy(Dummy, "0") 28.080130 d(D,"0
dummy(Dummy, "0"):Days 6.439404 0.094

Block 2: dummy(Dummy, "1"), dummy(Dummy, "1"):Days
Formula: ~0 + dummy(Dummy, "1") + Days:dummy(Dummy, "1") | Subject
Structure: General positive-definite
StdDev Corr
dummy(Dummy, "1") 23.209858 d(D,"1
dummy(Dummy, "1"):Days 3.567597 0.062
Residual 25.591791

Number of Observations: 180
Number of Groups: 18
>