Patrick Forscher (https://forscher.uark.edu/) recently noted in two tweets (https://twitter.com/psforscher/status/1087458806928420866 and https://twitter.com/psforscher/status/1087458806928420866) that social psychologists are regularly interested in comparing the size of variance components in multilevel models in their research (e.g., Hehman, Sutherland, Flake, & Slepian, 2017; Snijders & Kenny, 1999).
One common scenario is that a psychologists seeks to test whether the variance of two random effects in a multilevel model differ from each other. This scenario especially occurs in the well-known social relations model (Snijders & Kenny, 1999).
A second scenario described in the tweets is that a psychologist seeks to test whether either the random effect, the residual variance, or both differ between two conditions or subgroups in a dataset.
Earlier research has typically used software like SAS or MLwin to
model these types of scenarios. However, the R environement is
becoming increasingly popular among psychologists.
Patrick Forscher was therefore interested in doing these analyses in R. I responsed with R code but David Kenny (http://davidakenny.net/) noted an issue with my initial solution because of the default setting of the as.factor function in R and sent SAS output to check (thank you!)
Below I therefore provide thoroughly checked code in the lme4 (Bates, Mächler, Bolker, & Walker, 2015) and nlme (Pinheiro & Bates, 2000) packages in R for both scenarios that gives similar results as SAS. I use the classic sleepstudy data
in the lme4 package. This dataset was not designed for this type of
research question but is sufficient to illustrate the methodology.
In addition, I also use a very small typical social relations scenario.
library(nlme)
library(lme4)
# cross-classified model in lmer
lmer(Reaction~1+(1|Days)+(1|Subject),data=sleepstudy)
# cross-classified model with constrained variances in lme
sleepstudy$one<-1
mdays<-model.matrix(~as.factor(Days)-1, data = sleepstudy)
msubject<-model.matrix(~Subject-1, data = sleepstudy)
m1 <- lme(Reaction~1,
random=list(one=pdIdent(~mdays+msubject-1))
,data=sleepstudy)
# lmer model above in nlme
m2 <- lme(Reaction~1,
random=list(one=pdBlocked(list(pdIdent(~ as.factor(Days)-1),
pdIdent(~ as.factor(Subject)-1))))
,data=sleepstudy)
summary(m1)
summary(m2)
anova(m1,m2)
library(nlme)
newt<-data.frame(matrix(c(
1, 1, 1,
1, 2, 3,
1, 3, 3,
2, 1, 3,
2, 2, 4,
2, 3, 5,
3, 1, 4,
3, 2, 5,
3, 3, 7),ncol=3,byrow=T))
names(newt)<-c("actor","partner","dv")
newt$group<-1
mactor<-model.matrix(~as.factor(actor)-1, data = newt)
mpartner<-model.matrix(~as.factor(partner)-1, data = newt)
m1 <- lme(dv~1,
random=list(group=pdIdent(~mactor+mpartner-1))
,data=newt)
VarCorr(m1)
#unequal variances
m2 <- lme(dv~1,
random=list(group=pdBlocked(list(pdIdent(~ as.factor(actor)-1),
pdIdent(~ as.factor(partner)-1))))
,data=newt)
VarCorr(m2)
anova(m1,m2)
# compare random effects variance across two parts of a dataset
stacked<-data.frame(DV=c(sleepstudy$Reaction,sqrt(sleepstudy$Reaction)*100),
Subject= c(sleepstudy$Subject,paste(sleepstudy$Subject,"2",sep="_")),
dummy12=rep(c(0,1),each=nrow(sleepstudy)))
(m1<-lmer(DV~0+as.factor(dummy12)+(1|Subject),stacked))
(m2<-lmer(DV~0+as.factor(dummy12)+(0+dummy(dummy12,"0")|Subject)+
(0+dummy(dummy12,"1")|Subject),stacked))
anova(m1,m2,refit=FALSE)
# same in lme
(m1<-lme(DV~as.factor(dummy12)-1,
random=~1|Subject,data=stacked))
(m2<-lme(DV~as.factor(dummy12)-1,
random=list(Subject=pdDiag(~as.factor(dummy12)-1)),
data=stacked))
anova(m1,m2)
# compare random effects variances across subdimensions from the same
# subject
stacked<-data.frame(DV=c(sleepstudy$Reaction,sqrt(sleepstudy$Reaction)*100),
Subject= rep(sleepstudy$Subject,2),
dummy12=rep(c(0,1),each=nrow(sleepstudy)))
(m1<-lme(DV~as.factor(dummy12)-1,
random=list(Subject=pdCompSymm(~as.factor(dummy12)-1)),data=stacked))
(m2<-lme(DV~as.factor(dummy12)-1,
random=list(Subject=pdSymm(~as.factor(dummy12)-1)),
data=stacked))
anova(m1,m2)
# compare error variance
m1<-lme(DV~0+as.factor(dummy12),random=~1|Subject,stacked)
m2<-lme(DV~0+as.factor(dummy12),random=~1|Subject,
data=stacked,weights=varIdent(form=~1|as.factor(dummy12)))
summary(m1)
summary(m2)
anova(m1,m2)
# compare error variance with heterogenous random effects variance
m1<-lme(DV~as.factor(dummy12)-1,
random=list(Subject=pdSymm(~as.factor(dummy12)-1)),stacked)
m2<-lme(DV~as.factor(dummy12)-1,
random=list(Subject=pdSymm(~as.factor(dummy12)-1)),
data=stacked,weights=varIdent(form=~1|as.factor(dummy12)))
summary(m1)
summary(m2)
anova(m1,m2)
Literature
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting Linear mixed-effects models using lme4. Journal of Statistical Software, 67 (1), https://dx.doi.org/10.18637/jss.v067.i01
Hehman, E., Sutherland, C. A. M., Flake, J. K., & Slepian, M. L. (2017, May 8). The Unique Contributions of Perceiver and Target Characteristics in Person Perception. Journal of Personality and Social Psychology. Advance online publication. https://dx.doi.org/10.1037/pspa0000090
Pinheiro, J. C., & Bates, D. M. (2000). Mixed-effects models in S and S-Plus. New York: Springer.
Snijders, T. A. B., & Kenny, D. A. (1999). The social relations model for family data: A multilevel approach. Personal Relationships, 6(4), 471-486. https://dx.doi.org/10.1111/j.1475-6811.1999.tb00204.x