NOTE: This program is also available as a graphic user interface and its use does not require knowledge of R: https://aguinis.shinyapps.io/jomr/
Below is the R code used for analyses described in the following article:
Aguinis, H., Gottfredson, R. K., & Culpepper, S. A. (2013). Best-practice recommendations for estimating cross-level interaction effects using multilevel modeling. Journal of Management, 39, 1490-1528. [copyright notice: You may download this article for one-time personal use only; please obtain publisher permission for any further distribution, publication, or commercial use.] [ pdf ]
#Setting Working Directory and Reading Data File
setwd('C:/Documents/JOM')
exdata <- read.csv("JOM example data.csv",
header = TRUE,
sep = ",")
library(lme4)
library(RLRsim)
#STEP 1: Null Model
lmm.fit1=lmer(Y ~(1|l2id),data=exdata,REML=F)
summary(lmm.fit1)
# Compute ICC
iccy=VarCorr(lmm.fit1)$l2id[1,1]/(VarCorr(lmm.fit1)$l2id[1,1]+attr(VarCorr(lmm.fit1),'sc')^2)
iccy
#STEP 2: Random Intercept and Fixed Slope Model
lmm.fit2=lmer(Y ~(1|l2id)+Xc+I(Wj-mean(Wj) ),data=exdata,REML=F)
summary(lmm.fit2)
# Computing pseudo R-squared
yhat2=model.matrix(lmm.fit2)%*%fixef(lmm.fit2)
cor(yhat2,exdata$Y)^2
#STEP 3: Random Intercept and Random Slope model
lmm.fit3=lmer(Y ~Xc+(Xc|l2id)+I(Wj-mean(Wj) ),data=exdata,REML=F)
summary(lmm.fit3)
# Print VC Estimates
VarCorr(lmm.fit3)
# Computing pseudo R-squared
yhat3=model.matrix(lmm.fit3)%*%fixef(lmm.fit3)
cor(yhat3,exdata$Y)^2
# Crainceanu & Ruppert (2004) Test of Slope Variance Component
obs.LRT <- 2*(logLik(lmm.fit3)-logLik(lmm.fit2))[1]
X <- getME(lmm.fit3, "X")
Z <- t(as.matrix(getME(lmm.fit3,"Zt")))
sim.LRT <- LRTSim(X, Z, 0, diag(ncol(Z)))
(pval <- mean(sim.LRT > obs.LRT))
# Nonparametric Bootstrap Function
REMLVC=VarCorr(lmer(Y ~Xc+(Xc|l2id)+I(Wj-mean(Wj) ),data=exdata,REML=T))$l2id[1:2,1:2]
U.R=chol(REMLVC)
REbootstrap=function(Us,es,X,gs){
nj=nrow(Us)
idk=sample(1:nj,size=nj,replace=T)
Usk=as.matrix(Us[idk,])
esk=sample(es,size=length(es),replace=T)
S=t(Usk)%*%Usk/nj
U.S = chol(S)
A=solve(U.S)%*%U.R
Usk = Usk%*%A
datk=expand.grid(l1id = 1:6,l2id = 1:nj)
colnames(X)=c('one','Xc','Wjc')
datk=cbind(datk,X)
datk$yk = X%*%gs + Usk[datk$l2id,1]+Usk[datk$l2id,2]*X[,2]+esk
lmm.fitk=lmer(yk ~Xc+(Xc|l2id)+Wjc,data=datk,REML=F)
tau11k = VarCorr(lmm.fitk)$l2id[2,2]
tau11k
}
# Implementing Bootstrap
confint(lmm.fit3,method="boot")
#STEP 4: Cross-Level Interaction Model
lmm.fit4=lmer(Y ~(Xc|l2id)+Xc*I(Wj-mean(Wj) ),data=exdata,REML=F)
summary(lmm.fit4)
# Print VC Estimates
VarCorr(lmm.fit4)
# Computing pseudo R-squared
yhat4=model.matrix(lmm.fit4)%*%fixef(lmm.fit4)
cor(yhat4,exdata$Y)^2
#Interaction Plots
#Code creates graphs in pdf format in the same directory as the data file
gammas=fixef(lmm.fit4)
pdf('intplot.xw.pdf',width=10,height=8)
par(mar=c(3.25,3.25,.5,.5),cex=2,bty='l',las=1,family='serif',mgp=c(1.85,.5,0))
#Figure 3 Panel (a) - Full Y Scale
Wjs=c(0-sd(exdata$Wj),0,0+sd(exdata$Wj))
xlb=mean(exdata$Xc)-sd(exdata$Xc);xub=mean(exdata$Xc)+sd(exdata$Xc)
ylb=1;yub=7
curve(0+1*x,xlb,xub,xlab='LMX',ylab='Individual Empowerment',lwd=2,type='n',
ylim=c(ylb,yub))
for(i in 1:length(Wjs)){
B0j=gammas[1]+gammas[3]*Wjs[i]
B1j=gammas[2]+gammas[4]*Wjs[i]
curve(B0j+B1j*x,xlb,xub,add=T,xlab='LMX',ylab='Individual Empowerment',lwd=2,lty=i)
}
labs=c(expression(W[j]==-1*~~SD),expression(W[j]==0*~~SD),expression(W[j]==1*~~SD))
legend(xlb,5,legend=c("Leadership Climate",labs[1],labs[2],labs[3]),bty='n',lty=c(0:3))
#Figure 3 Panel (b) - Reduced Y Scale
ylb=5;yub=6.5
curve(0+1*x,xlb,xub,xlab='LMX',ylab='Individual Empowerment',lwd=2,type='n',
ylim=c(ylb,yub))
for(i in 1:length(Wjs)){
B0j=gammas[1]+gammas[3]*Wjs[i]
B1j=gammas[2]+gammas[4]*Wjs[i]
curve(B0j+B1j*x,xlb,xub,add=T,xlab='LMX',ylab='Individual Empowerment',lwd=2,lty=i)
}
labs=c(expression(W[j]==-1*~~SD),expression(W[j]==0*~~SD),expression(W[j]==1*~~SD))
legend(xlb,6.5,legend=c("Leadership Climate",labs[1],labs[2],labs[3]),bty='n',lty=c(0:3))
dev.off()
