______________________________________________________________________________________________

 

MAIN PAGE

______________________________________________________________________________________________

 

Annotated R Code for Aguinis, Gottfredson, and Culpepper (2013, Journal of Management)

 

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. [pdf version

 

·       Click here to download the data file

·       Click here to download the power calculator included in the article’s Appendix B

 

 

#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()
 

 

______________________________________________________________________________________________

 

MAIN PAGE

______________________________________________________________________________________________

 

<