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
> 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
>