Bliese, P.D., Maltarich, M.A. & Hendricks, J.L. Back to Basics with Mixed-Effects Models: Nine Take-Away Points. J Bus Psychol 33, 1–23 (2018).

https://doi.org/10.1007/s10869-017-9491-z

1 Simulation 1: Independent Data

##Set up Independent data
set.seed(12532) #使用set.seed函數，()內設定隨機序列數據的起始值
#create vector x and vector y
Y<-rnorm(1000) #用rnormY隨機生成一組1000筆數據Y，Y是常態分配，平均為0，標準差為1

X<-.30*Y+sqrt(1-.30^2)*rnorm(1000) #sqrt是square root
#X和Y為低度相關
cor(X,Y)
## [1] 0.2937485
#建立group ID向量，ID從1到100，重複10次，平均一組有10個成員
G.ID<-rep(1:100,10) 
#隨機分派Group IDs給observations(X,Y)
NO.GRP<-data.frame(G.ID=G.ID,Y=Y,X=X) 
head(NO.GRP)
##   G.ID          Y          X
## 1    1 -0.3562218 -0.2533692
## 2    2  0.3926757  1.9169494
## 3    3 -0.4866951 -0.3444651
## 4    4  0.2367052 -1.6619983
## 5    5  0.5555600  1.3272490
## 6    6  0.1011640  2.1926501
#Calculate group means for X.
TDAT<-aggregate(NO.GRP[,c("G.ID","X")],by=list(NO.GRP$G.ID),mean) #list()取得G.ID變量
#Assign means back to members as X.G
NO.GRP<-merge(NO.GRP,TDAT[,c("G.ID","X")],by="G.ID", #merge()將兩個data.frame 透過by="G.ID"進行合併
suffixes=c("",".G")) #The default is c(".x",".y")

2 Table 1

#These data have no group-level properties.
#NO.GRP[1:15,]
knitr::kable(head(NO.GRP,15))
G.ID Y X X.G
1 -0.3562218 -0.2533692 -0.1805595
1 -0.4222437 -1.1081410 -0.1805595
1 -0.4751256 0.1228866 -0.1805595
1 -0.0754015 -0.2089118 -0.1805595
1 0.0401238 -0.5647489 -0.1805595
1 -0.0425110 0.3182596 -0.1805595
1 1.8467950 0.9231918 -0.1805595
1 0.5551599 -0.5044792 -0.1805595
1 1.6542047 0.7838111 -0.1805595
1 -2.3206311 -1.3140941 -0.1805595
2 -0.8275048 -0.6610696 -0.2058514
2 0.3926757 1.9169494 -0.2058514
2 0.3734971 -0.8080668 -0.2058514
2 -2.4918987 -2.2593381 -0.2058514
2 -1.4049858 1.5874561 -0.2058514 
#install package and library
pacman::p_load(multilevel,lme,TinyTeX,LaTeX)
## Installing package into 'C:/Users/Ching-Fang Wu/Documents/R/win-library/4.0'
## (as 'lib' is unspecified)
## Warning: package 'lme' is not available for this version of R
## 
## A version of this package for your version of R might be available elsewhere,
## see the ideas at
## https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages
## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.0:
##   無法開啟 URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.0/PACKAGES'
## Warning: 'BiocManager' not available.  Could not check Bioconductor.
## 
## Please use `install.packages('BiocManager')` and then retry.
## Warning in p_install(package, character.only = TRUE, ...):
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'lme'
## Installing package into 'C:/Users/Ching-Fang Wu/Documents/R/win-library/4.0'
## (as 'lib' is unspecified)
## Warning: package 'TinyTeX' is not available for this version of R
## 
## A version of this package for your version of R might be available elsewhere,
## see the ideas at
## https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages
## Warning: Perhaps you meant 'tinytex' ?
## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.0:
##   無法開啟 URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.0/PACKAGES'
## Warning: 'BiocManager' not available.  Could not check Bioconductor.
## 
## Please use `install.packages('BiocManager')` and then retry.
## Warning in p_install(package, character.only = TRUE, ...):
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'TinyTeX'
## Installing package into 'C:/Users/Ching-Fang Wu/Documents/R/win-library/4.0'
## (as 'lib' is unspecified)
## Warning: package 'LaTeX' is not available for this version of R
## 
## A version of this package for your version of R might be available elsewhere,
## see the ideas at
## https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages
## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.0:
##   無法開啟 URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.0/PACKAGES'
## Warning: 'BiocManager' not available.  Could not check Bioconductor.
## 
## Please use `install.packages('BiocManager')` and then retry.
## Warning in p_install(package, character.only = TRUE, ...):
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'LaTeX'
## Warning in pacman::p_load(multilevel, lme, TinyTeX, LaTeX): Failed to install/load:
## lme, TinyTeX, LaTeX

有兩種方式可計算ICC(1)

##ICC(1) estimate via mixed-effects
tmod<-lme(fixed=Y~1,random=~1|G.ID,data=NO.GRP)
VarCorr(tmod)
## G.ID = pdLogChol(1) 
##             Variance     StdDev      
## (Intercept) 4.299649e-09 6.557171e-05
## Residual    9.986574e-01 9.993285e-01
#ICC(1)
4.299649e-09/(4.299649e-09+9.986574e-01)
## [1] 4.305429e-09
# ICC(1) estimate via ANOVA
tmod<-aov(Y~as.factor(G.ID),data=NO.GRP)
ICC1(tmod)
## [1] -0.01560787

以上用兩種方式計算組內相關係數ICC(1)皆趨近於0。

3 Table 2

#Table 2 regressing Y on X and X.G in OLS regression in simulated data with near-zero ICC(1)
tmod<-lm(Y~X+X.G, data=NO.GRP)
summary(tmod)$coef
##                Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)  0.06158212 0.03026285  2.034908 4.212427e-02
## X            0.30195145 0.03161021  9.552340 9.483257e-21
## X.G         -0.10980587 0.09814650 -1.118796 2.634968e-01
#Table 2 presents the results from regressing Y on X and X.G in mixed-effects random intercept model in simulated data with near-zero ICC(1)
tmod<-lme(fixed=Y~X+X.G,random=~1|G.ID,data=NO.GRP)
summary(tmod)$tTable
##                   Value  Std.Error  DF   t-value      p-value
## (Intercept)  0.06158212 0.03026285 899  2.034908 4.215308e-02
## X            0.30195145 0.03161021 899  9.552341 1.164921e-20
## X.G         -0.10980587 0.09814650  98 -1.118796 2.659622e-01

Take-away point 1： Mixed-effects models can be used with nested data, even if no group-level effects are expected or evident.

4 Re-run

#Re-run simulation with data that returns a slightly positive ICC(1)
set.seed(125321) #變更起始直
Y<-rnorm(1000)
X<-.30*Y+sqrt(1-.30^2)*rnorm(1000)
cor(X,Y)
## [1] 0.2827705
G.ID<-rep(1:100,10)
NO.GRP<-data.frame(G.ID=G.ID,Y=Y,X=X)
# ICC(1) estimate via mixed-effects
tmod<-lme(fixed=Y~1,random=~1|G.ID,data=NO.GRP)
VarCorr(tmod)
## G.ID = pdLogChol(1) 
##             Variance   StdDev   
## (Intercept) 0.01292908 0.1137061
## Residual    1.00240352 1.0012010
(0.01292908)/(0.01292908+1.00240352)
## [1] 0.01273384
# ICC(1) estimate via ANOVA
tmod<-aov(Y~as.factor(G.ID),data=NO.GRP)
summary(tmod)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(G.ID)  99  112.0   1.132   1.129  0.194
## Residuals       900  902.2   1.002
ICC1(tmod)
## [1] 0.01273379
# Calculate group means for X. Assign means back to members as X.G
TDAT<-aggregate(NO.GRP[,c("G.ID","X")],list(NO.GRP$G.ID),mean)
NO.GRP<-merge(NO.GRP,TDAT[,c("G.ID","X")],by="G.ID",
suffixes=c("",".G"))

5 Table 3

#Table 3 (upper panel) presents the results of regressing Y on X and X.G in OLS in simulated data with small positive ICC(1)
tmod<-lm(Y~X+X.G, data=NO.GRP)
summary(tmod)$coef
##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 0.01482297 0.03076182 0.4818625 6.300094e-01
## X           0.28267888 0.03285108 8.6048573 2.935303e-17
## X.G         0.09472407 0.10806277 0.8765653 3.809340e-01
#Table 3 (lower panel) presents the results of regressing Y on X and X.G in mixed-effects random intercept model in simulated data with small positive ICC(1)
tmod<-lme(fixed=Y~X+X.G, random=~1|G.ID,data=NO.GRP)
summary(tmod)$tTable
##                  Value  Std.Error  DF   t-value      p-value
## (Intercept) 0.01482297 0.03204743 899 0.4625321 6.438117e-01
## X           0.28267888 0.03269794 899 8.6451599 2.437036e-17
## X.G         0.09472407 0.11212448  98 0.8448117 4.002738e-01

Take-away point 2: Failing to use mixed effects to model nested data can increase the risk of type I error (for grouplevel effects) and type II error (for lower-level effects).

6 Group-Mean Centering

# Create a group-mean centered variable (W.X)
NO.GRP$W.X<-NO.GRP$X-NO.GRP$X.G

7 Table 4

#Table 4 (upper panel) Results of regressing YonW.X and X.G in OLS (upper panel) in simulated data with small positive ICC(1)

tmod<-lm(Y~W.X+X.G, data=NO.GRP)
knitr::kable(summary(tmod)$coef)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0148230 0.0307618 0.4818625 0.6300094
W.X 0.2826789 0.0328511 8.6048573 0.0000000
X.G 0.3774030 0.1029484 3.6659436 0.0002594 
#Table 4 (lower panel) Results of regressing YonW.X and X.G in mixed-effects random intercept model in simulated data with small positive ICC(1)

tmod<-lme(fixed=Y~W.X+X.G,random=~1|G.ID,data=NO.GRP)
knitr::kable(summary(tmod)$tTable)
Value Std.Error DF t-value p-value
(Intercept) 0.0148230 0.0320474 899 0.4625321 0.6438117
W.X 0.2826789 0.0326979 899 8.6451599 0.0000000
X.G 0.3774030 0.1072509 98 3.5188807 0.0006589

Take-away point 3: Group-mean centering of a level-1 variable fundamentally changes the interpretation of the level-2 parameter estimate for the analogue of the same variable. In group-mean centered models, level-2 parameter estimates represent overall group effects; in raw or grand-mean models level-2, parameter estimates represent differences in slopes between individual-level and group-level relationships.

Take-away point 4: Group-mean centering changes the conceptual meaning of the level-1 construct such that the term now reflects relative position in a group rather than absolute values.

8 Table 5

# create a dataset that contains 100 (X.G, Y.G) pairs instead of 1000 (X, Y) pairs

TDAT<-aggregate(NO.GRP[,c("Y","X")],list(NO.GRP$G.ID),mean)
names(TDAT)<-c("G.ID","Y.G","X.G")
#Table 5 (results from 100 group means) Results of OLS group-mean regression of Y.G on X.G in simulated data with small positive ICC(1)

tmod<-lm(Y.G~X.G, data=TDAT)
summary(tmod)$coef
##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 0.01482297 0.03204743 0.4625321 0.6447251131
## X.G         0.37740295 0.10725085 3.5188807 0.0006589109
#Add another group-level predictor, Z.G, with no level-1 analogue
set.seed(125321)
Y<-rnorm(1000)
X<-.30*Y+sqrt(1-.30^2)*rnorm(1000)
G.ID<-rep(1:100,10)
NO.GRP<-data.frame(G.ID=G.ID,Y=Y,X=X)
NO.GRP<-NO.GRP[order(NO.GRP$G.ID),]
NO.GRP$Z.G<-rep(rnorm(100),each=10) #新增變數Z.G
tmod<-lme(fixed=Y~Z.G, random=~1|G.ID, data=NO.GRP)
summary(tmod)$tTable
##                  Value  Std.Error  DF   t-value   p-value
## (Intercept) 0.02725969 0.03377608 900 0.8070708 0.4198389
## Z.G         0.01746231 0.03664744  98 0.4764947 0.6347827
TDAT<-aggregate(NO.GRP[,c("Y","Z.G")],list(NO.GRP$G.ID),mean)
names(TDAT)<-c("G.ID","Y.G","Z.G")
tmod<-lm(Y.G~Z.G, data=TDAT)
summary(tmod)$coef
##               Estimate Std. Error   t value  Pr(>|t|)
## (Intercept) 0.02725969 0.03377607 0.8070711 0.4215800
## Z.G         0.01746231 0.03664742 0.4764948 0.6347826

Take-away point 5: In mixed-effects models, higher-level variables predict the group means of the lower-level outcome.

9 Simulate 2：pure group-level data

This simulation relies on the sim.icc function in the {multilevel} library.

library(multilevel)

set.seed(632342) 
ALL.GRP<-sim.icc(gsize=10,ngrp=100,icc1=.15,nitems=2) #未設定 item.cor，表示items with no level-1 correlation
names(ALL.GRP)<-c("G.ID","Y","X")

若設定 item.cor=.3，表示有.30 level-1 item correlations

# covariance theorem decomposition results
with(ALL.GRP, waba(x=X, y=Y, grpid=G.ID))
## $Cov.Theorem
##     RawCorr     EtaBx     EtaBy     CorrB     EtaWx     EtaWy       CorrW
## 1 0.1474048 0.4892028 0.4856851 0.6101503 0.8721701 0.8741338 0.003192502
## 
## $n.obs
## [1] 1000
## 
## $n.grps
## [1] 100
# Randomly sort X on a group-by-group basis into X2 (footnote5)
set.seed(467432)
ALL.GRP$X2<-unlist(tapply(ALL.GRP$X,ALL.GRP$G.ID,sample))
ALL.GRP[1:10,] #observe that X2 values are X values in new order
##    G.ID          Y          X         X2
## 1     1 -0.3501393 -0.8712853 -0.2615682
## 2     1 -2.3714188 -0.7405538 -0.7405538
## 3     1 -1.3145565 -0.7956528 -0.6557583
## 4     1 -1.7994279 -0.6557583 -0.1659889
## 5     1  0.5989924 -0.2615682 -2.8472639
## 6     1 -1.2710712 -1.6651974 -0.7956528
## 7     1  0.2234514 -2.1611040 -1.6651974
## 8     1 -1.5340462 -0.1659889 -0.8712853
## 9     1  0.1501798 -1.3895897 -1.3895897
## 10    1 -1.4292006 -2.8472639 -2.1611040
with(ALL.GRP,waba(Y,X2,G.ID))
## $Cov.Theorem
##     RawCorr     EtaBx     EtaBy     CorrB     EtaWx     EtaWy      CorrW
## 1 0.1815365 0.4856851 0.4892028 0.6101503 0.8741338 0.8721701 0.04796176
## 
## $n.obs
## [1] 1000
## 
## $n.grps
## [1] 100
#Estimate the ICC(1) and ICC(2) values for the Y variable using
#the ANOVA method
tmod<-aov(Y~as.factor(G.ID),data=ALL.GRP)
ICC1(tmod)
## [1] 0.1530071
ICC2(tmod)
## [1] 0.643681
# Calculate group means for X. Assign means back to members as X.G
TDAT<-aggregate(ALL.GRP[,c("G.ID","X")],list(ALL.GRP$G.ID),mean)
ALL.GRP<-merge(ALL.GRP,TDAT[,c("G.ID","X")],
by="G.ID",suffixes=c("",".G"))

10 Table 6

#Table 6 (upper panel) Results of regressing Y on X and X.G in OLS (upper panel) in simulated data with substantial positive ICC(1)
tmod<-lm(Y~X+X.G,data=ALL.GRP)
summary(tmod)$coef
##               Estimate Std. Error    t value     Pr(>|t|)
## (Intercept) 0.05044334 0.03319800 1.51946908 1.289615e-01
## X           0.00329694 0.03573496 0.09226092 9.265093e-01
## X.G         0.62087736 0.07304733 8.49965800 6.863064e-17
#Table 6 (lower panel) Results of regressing Y on X and X.G in mixed-effects random intercept model in simulated data with substantial positive ICC(1)
tmod<-lme(fixed=Y~X+X.G, random=~1|G.ID, data=ALL.GRP)
summary(tmod)$tTable
##                  Value  Std.Error  DF    t-value      p-value
## (Intercept) 0.05044334 0.04266202 899 1.18239449 2.373618e-01
## X           0.00329694 0.03444278 899 0.09572225 9.237625e-01
## X.G         0.62087736 0.08882183  98 6.99014374 3.367965e-10
#Create a group-mean centered variable (W.X)
ALL.GRP$W.X<-ALL.GRP$X-ALL.GRP$X.G
#Estimate total effect of X.G with a group-mean centered X variable

tmod<-lme(fixed=Y~W.X+X.G, random=~1|G.ID, data=ALL.GRP)
summary(tmod)$tTable
##                  Value  Std.Error  DF    t-value      p-value
## (Intercept) 0.05044334 0.04266204 899 1.18239408 2.373620e-01
## W.X         0.00329694 0.03444277 899 0.09572226 9.237625e-01
## X.G         0.62417430 0.08187196  98 7.62378623 1.596094e-11

Take-away point 6. In mixed-effects models, emergent effects identifying different relationships between level-1 variables and their group mean analogs likely represent important changes in the meaning of constructs across levels.

# Estimate ICC(1) and ICC(2) values for both X and Y using
# ANOVA method
tmod<-aov(Y~as.factor(G.ID),data=ALL.GRP)
ICC1(tmod)
## [1] 0.1530071
ICC2(tmod)
## [1] 0.643681
tmod<-aov(X~as.factor(G.ID),data=ALL.GRP)
ICC1(tmod)
## [1] 0.1568375
ICC2(tmod)
## [1] 0.6503632

11 Re-run：group sizes=2

#Demonstrate the effect of alternative group sizes on ICC(2) values, and on the ability to detect emergent effects. Substitute either 2 or 100 for gsize in sim.icc function.

set.seed(632342)
ALL.GRP<-sim.icc(gsize=2,ngrp=100,icc1=.15,nitems=2)
names(ALL.GRP)<-c("G.ID","Y","X")
ALL.GRP[1:15,]
##    G.ID          Y           X
## 1     1 -0.1093418  0.12424801
## 2     1 -2.1306214 -3.36765654
## 3     2 -1.0169082 -1.36755166
## 4     2 -1.5017796 -1.51738016
## 5     3  0.9077974 -0.23947720
## 6     3 -0.9622663  0.58838004
## 7     4  0.5395617  0.27097823
## 8     4 -1.2179360  0.09219963
## 9     5  0.5633326 -0.32990359
## 10    5 -1.0160478 -1.20213434
## 11    6 -1.5917934 -0.67586873
## 12    6 -1.8433899 -0.91484613
## 13    7 -0.2631602 -1.18739589
## 14    7 -0.1872816 -1.13403820
## 15    8 -3.7387863  0.52169297
with(ALL.GRP,waba(x=X, y=Y, grpid=G.ID))
## $Cov.Theorem
##     RawCorr     EtaBx     EtaBy     CorrB     EtaWx    EtaWy       CorrW
## 1 0.1224128 0.7756321 0.7742287 0.2391332 0.6311853 0.632906 -0.05304541
## 
## $n.obs
## [1] 200
## 
## $n.grps
## [1] 100
tmod<-aov(Y~as.factor(G.ID),data=ALL.GRP)
ICC1(tmod)
## [1] 0.2036817
ICC2(tmod)
## [1] 0.3384311
tmod<-aov(X~as.factor(G.ID),data=ALL.GRP)
ICC1(tmod)
## [1] 0.208023
ICC2(tmod)
## [1] 0.3444024
TDAT<-aggregate(ALL.GRP[,c("G.ID","X")],list(ALL.GRP$G.ID),mean)
ALL.GRP<-merge(ALL.GRP,TDAT[,c("G.ID","X")], by="G.ID",suffixes=c("",".G"))
tmod<-lme(fixed=Y~X+X.G, random=~1|G.ID, data=ALL.GRP)
summary(tmod)$tTable
##                  Value  Std.Error DF    t-value    p-value
## (Intercept)  0.1134809 0.08645205 99  1.3126456 0.19233730
## X           -0.0526323 0.09958068 99 -0.5285392 0.59830802
## X.G          0.2888299 0.13893207 98  2.0789288 0.04023588

Take-away point 7: Substantial ICC (2) values are not necessary for identifying emergent group-level effects (but they help).

12 Real data：Organizational data

The dataset (referred to as bh1996 in R) contains responses from 7382 soldiers nested in 99 groups (army companies).

The dependent variable is a measure of well-being (WBEING) and the three independent variables are work hours (HRS), leadership consideration (LEAD), and cohesion (COHES).

In the dataset, all variables have a group-mean centered variant (W.BEING, W.HRS, W.LEAD, W.COHES) and each have a groupmean variant (G.WBEING, G.HRS, G.LEAD, G.COHES).

Data: bh1996{multilevel}

Column 1: GRP=Group ID

Column 2: COHES=cohesion(individual) Column 3: G.COHES=the cohesion of groupmean Column 4: W.COHES=the cohesion of group-mean centered

Column 5：LEAD=leadership consideration Column 6：G.LEAD Column 7：W.LEAD

Column 8：HRS=work hours Column 9：G.HRS Column 10：W.HRS

Column 11：WBEING=well-being Column 12：G.WBEING Column 13：W.WBEING

library(multilevel)
data(bh1996) #bring data from library to work space
#Estimate ICC(1) using mixed-effects model
null<-lme(fixed=WBEING~1,random=~1|GRP,data=bh1996)
VarCorr(null)
## GRP = pdLogChol(1) 
##             Variance   StdDev   
## (Intercept) 0.03580077 0.1892109
## Residual    0.78949727 0.8885366
0.03580077/(0.03580077+0.78949727) #[1] 0.0433792
## [1] 0.0433792
#likelihood test of significance of ICC(1)
null.gls<-gls(WBEING~1,data=bh1996)
anova(null.gls,null)
##          Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## null.gls     1  2 19540.17 19553.98 -9768.084                        
## null         2  3 19353.34 19374.06 -9673.669 1 vs 2 188.8303  <.0001

13 Table 7 (upper left panel): OLS raw variables.

#Table 7 Results of predicting soldier well-being in OLS and mixed-effects random intercept models in bh1996 data

tmod<-lm(WBEING~HRS+LEAD+COHES+G.HRS+G.LEAD+G.COHES, data=bh1996)
summary(tmod)$coef
##                Estimate  Std. Error     t value      Pr(>|t|)
## (Intercept)  3.67579453 0.197369836  18.6238921  1.098176e-75
## HRS         -0.02559700 0.004508821  -5.6770927  1.422046e-08
## LEAD         0.47093160 0.014379444  32.7503342 9.162611e-220
## COHES        0.08022125 0.012222268   6.5635325  5.610775e-11
## G.HRS       -0.12775356 0.012604704 -10.1353873  5.523360e-24
## G.LEAD      -0.23385790 0.048893364  -4.7830192  1.760341e-06
## G.COHES     -0.03083558 0.063010526  -0.4893718  6.245930e-01

14 Table 7 (upper right panel): Mixed-effect raw variables.

tmod<-lme(fixed=WBEING~HRS+LEAD+COHES+G.HRS+G.LEAD+G.COHES, random=~1|GRP, data=bh1996)
summary(tmod)$tTable
##                   Value  Std.Error   DF    t-value       p-value
## (Intercept)  3.53005480 0.30429144 7280 11.6009006  7.606959e-31
## HRS         -0.02559700 0.00447701 7280 -5.7174316  1.124246e-08
## LEAD         0.47093160 0.01427799 7280 32.9830435 1.783470e-222
## COHES        0.08022125 0.01213603 7280  6.6101700  4.109995e-11
## G.HRS       -0.11540866 0.01883676   95 -6.1267797  2.022912e-08
## G.LEAD      -0.22396286 0.06865998   95 -3.2619127  1.537758e-03
## G.COHES     -0.03787917 0.08913739   95 -0.4249527  6.718322e-01

15 Table 7 (lower left panel): OLS group-mean centered

tmod<-lm(WBEING~W.HRS+W.LEAD+W.COHES+G.HRS+G.LEAD+G.COHES,
data=bh1996)
summary(tmod)$coef
##                Estimate  Std. Error     t value      Pr(>|t|)
## (Intercept)  3.67579453 0.197369836  18.6238921  1.098176e-75
## W.HRS       -0.02559700 0.004508821  -5.6770927  1.422046e-08
## W.LEAD       0.47093160 0.014379444  32.7503342 9.162610e-220
## W.COHES      0.08022125 0.012222268   6.5635325  5.610775e-11
## G.HRS       -0.15335056 0.011770688 -13.0281725  2.239126e-38
## G.LEAD       0.23707370 0.046731067   5.0731499  4.008356e-07
## G.COHES      0.04938568 0.061813773   0.7989429  4.243492e-01

#Table 7 (lower right panel): Mixed-effect group-mean centered

tmod<-lme(fixed=WBEING~W.HRS+W.LEAD+W.COHES+G.HRS+G.LEAD+G.COHES,random=~1|GRP, data= bh1996)
summary(tmod)$tTable
##                   Value  Std.Error   DF    t-value       p-value
## (Intercept)  3.53005475 0.30429151 7280 11.6008979  7.607194e-31
## W.HRS       -0.02559700 0.00447701 7280 -5.7174316  1.124246e-08
## W.LEAD       0.47093160 0.01427799 7280 32.9830436 1.783466e-222
## W.COHES      0.08022125 0.01213603 7280  6.6101700  4.109995e-11
## G.HRS       -0.14100566 0.01829699   95 -7.7064933  1.244371e-11
## G.LEAD       0.24696875 0.06715902   95  3.6773730  3.904120e-04
## G.COHES      0.04234208 0.08830739   95  0.4794851  6.326956e-01
# Estimate ICC values for leadership and well-being
mult.icc(bh1996[,c("LEAD","WBEING")],bh1996$GRP)
##   Variable       ICC1      ICC2
## 1     LEAD 0.14746129 0.9280442
## 2   WBEING 0.04337922 0.7717561
# Estimate a group mean analysis
TDAT<-aggregate(bh1996[,c("WBEING","HRS","LEAD","COHES")],
list(bh1996$GRP),mean,na.rm=TRUE)
summary(lm(WBEING~HRS+LEAD+COHES,data=TDAT))$coef
##                Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)  3.40203759 0.32099635 10.598368 8.753644e-18
## HRS         -0.12960704 0.01923640 -6.737592 1.232628e-09
## LEAD         0.25802876 0.06619315  3.898119 1.806074e-04
## COHES        0.03301303 0.08641862  0.382013 7.033049e-01

16 Table 8：Change ICC(2) values

# Table 8 Raw, within-group, and between-group correlations in the full bh1996 data and a random subset
mult.icc(bh1996[,c("HRS","WBEING")],bh1996$GRP)
##   Variable       ICC1      ICC2
## 1      HRS 0.12923699 0.9171286
## 2   WBEING 0.04337922 0.7717561
round(with(bh1996,waba(x=WBEING,y=HRS,grpid=GRP))$Cov.Theorem,dig=4)
##   RawCorr  EtaBx  EtaBy   CorrB  EtaWx  EtaWy   CorrW
## 1 -0.1632 0.2359 0.3788 -0.7122 0.9718 0.9255 -0.1107
set.seed(271843)
bh1996.small<-bh1996[sample(1:7382, 2000),]
length(unique(bh1996.small$GRP))
## [1] 99
mult.icc(bh1996.small[,c("HRS","WBEING")],bh1996.small$GRP)
##   Variable       ICC1      ICC2
## 1      HRS 0.13631938 0.7612564
## 2   WBEING 0.02984101 0.3832458
round(with(bh1996.small,waba(x=WBEING,y=HRS,grpid=GRP))$Cov.Theorem,dig=4)
##   RawCorr  EtaBx  EtaBy   CorrB  EtaWx  EtaWy   CorrW
## 1 -0.1735 0.2761 0.4227 -0.5706 0.9611 0.9063 -0.1228

17 Figure 1

library(lattice)
# Figure 1
xyplot(WBEING~HRS|as.factor(GRP),data=bh1996[1:1582,],
subset=HRS>5&HRS<19,
type=c("p","g","r"),col="black",col.line="black",
xlab="Work Hours",
ylab="Well-Being")

# Find group-specific OLS based intercepts and slopes
lmList(WBEING~HRS|GRP,bh1996)
## Call:
##   Model: WBEING ~ HRS | GRP 
##    Data: bh1996 
## 
## Coefficients:
##    (Intercept)           HRS
## 1     3.990076 -0.1090918821
## 2     3.656022 -0.0770116565
## 3     3.138605 -0.0377495125
## 4     4.612616 -0.1291435185
## 5     3.943244 -0.0814441809
## 6     4.509408 -0.1686102511
## 7     4.749507 -0.1807241415
## 8     3.541909 -0.0594313280
## 9     3.440793 -0.0329745811
## 10    4.452032 -0.1332075472
## 11    2.152582  0.0710503038
## 12    3.176950 -0.0453647359
## 13    2.588025  0.0295107381
## 14    3.180603 -0.0307705367
## 15    3.380240 -0.0591766694
## 16    3.519242 -0.0293147553
## 17    3.137004 -0.0159023799
## 18    3.292136 -0.0602867383
## 19    3.478979 -0.0593093093
## 20    4.203891 -0.1302083333
## 21    2.923744 -0.0362324736
## 22    4.752647 -0.1664581161
## 23    3.406885 -0.0775510671
## 24    3.226696 -0.0160713460
## 25    2.134408  0.0782083929
## 26    3.406327 -0.0690744160
## 27    4.056744 -0.1089082690
## 28    2.784956  0.0113244507
## 29    3.313504 -0.0291498415
## 30    2.898597 -0.0090528692
## 31    3.762125 -0.0846267262
## 32    4.880781 -0.1624848275
## 33    4.509819 -0.1292233756
## 34    3.264284 -0.0538783415
## 35    3.804898 -0.0910653056
## 36    3.835641 -0.0769648055
## 37    3.684616 -0.0957464513
## 38    3.924652 -0.0894073483
## 39    3.252348 -0.0634034494
## 40    2.522441  0.0069986006
## 41    2.476471  0.0264834588
## 42    3.380790 -0.0481803620
## 43    2.352236 -0.0019531959
## 44    2.542019  0.0016669529
## 45    2.683700 -0.0235012083
## 46    4.083326 -0.1103583224
## 47    5.290547 -0.2043859377
## 48    3.158070 -0.0439466874
## 49    4.268862 -0.1225687607
## 50    3.952318 -0.1002583979
## 51    2.622373 -0.0163326144
## 52    3.864839 -0.0926433289
## 53    3.901901 -0.1166338044
## 54    3.892217 -0.1184869502
## 55    2.218996  0.0087618862
## 56    2.606305 -0.0073997540
## 57    1.857415  0.0194664374
## 58    3.907074 -0.1174722481
## 59    2.960114 -0.0017443384
## 60    3.403956 -0.0616221812
## 61    2.558724  0.0062400708
## 62    2.449520  0.0011302312
## 63    3.826470 -0.1077130254
## 64    3.258092 -0.0506602181
## 65    3.231386 -0.0642348202
## 66    3.211466 -0.0223633482
## 67    3.534361 -0.0676987447
## 68    2.943184  0.0286947883
## 69    2.100722  0.1002843395
## 70    4.083821 -0.1153427488
## 71    2.789482  0.0094337214
## 72    1.044416  0.1406914618
## 73    2.887272 -0.0019154213
## 74    3.045580 -0.0154613185
## 75    3.540192 -0.0667879457
## 76    3.009243 -0.0217539752
## 77    4.574339 -0.1293929281
## 78    2.225326  0.0407785220
## 79    3.554965 -0.0582365550
## 80    2.136822  0.0232159685
## 81    3.399819 -0.0263932975
## 82    2.380841  0.0294641329
## 83    2.873608  0.0020058422
## 84    3.515166 -0.0546926784
## 85    4.091984 -0.1228267764
## 86    3.801917 -0.0840384899
## 87    3.177453 -0.0248849551
## 88    4.121530 -0.1274934141
## 89    3.467493 -0.0549382135
## 90    3.543867 -0.0688363697
## 91    4.671204 -0.1613910803
## 92    3.856942 -0.1015175131
## 93    2.883329 -0.0001943124
## 94    2.923498  0.0108107835
## 95    3.073227 -0.0552779040
## 96    3.162074 -0.0159898444
## 97    2.950074 -0.0204042401
## 98    4.201389 -0.1111111111
## 99    2.180802  0.0817230274
## 
## Degrees of freedom: 7382 total; 7184 residual
## Residual standard error: 0.8822627

18 Figure 2

# Function to generate correlations on a group by group basis.
# The function is used below to generate data for Figure 2.
tfun<-function(gsize,ngroup,rho){
OUTPUT<-data.frame(GRP=rep(NA,gsize*ngroup),X=rep(NA,gsize*ngroup),
Y=rep(NA,gsize*ngroup))
for(i in 1:ngroup){
Y<-rnorm(gsize)
X<-rho*Y+sqrt(1-rho^2)*rnorm(gsize)
OUTPUT[((i * gsize - gsize) + 1):(i *gsize),]<-data.frame(GRP=i,X=X,Y=Y)
}
return(OUTPUT)
}
# Figure 2: Simulate the bh1996 Work Hour and Well-Being
# relationship for 25 groups each with 75 members.
set.seed(12532)
SIM.DAT<-tfun(75,25,-.16)
xyplot(Y~X|as.factor(GRP),SIM.DAT,
type=c("p","g","r"),col="black",col.line="black")

19 Test whether slopes randomly vary in bh1996 data

tmod<-lme(fixed=WBEING~HRS, random=~1|GRP, data=bh1996)
tmod2<-lme(fixed=WBEING~HRS, random=~HRS|GRP,data=bh1996)
anova(tmod,tmod2)
##       Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## tmod      1  4 19249.79 19277.42 -9620.896                        
## tmod2     2  6 19250.67 19292.11 -9619.337 1 vs 2 3.118139  0.2103
#> anova(tmod,tmod2)
# Model df AIC BIC logLik Test L.Ratio p-value
#tmod 1 4 19249.79 19277.42 -9620.896
#tmod2 2 6 19250.67 19292.11 -9619.337 1 vs 2 3.118139 0.2103

Take-away point 8: It is helpful and informative to test and report whether level-1 slopes randomly vary among groups prior to conducting a test of a cross-level interaction, but the results of such tests should not prevent subsequently testing cross-level interactions.

20 Table 9

# Table 9(upper panel): Cross-level interactions and group-mean
# Centering variables
tmod<-lme(fixed=WBEING~HRS*G.COHES, random=~1|GRP, data=bh1996)
summary(tmod)$tTable
##                   Value  Std.Error   DF   t-value     p-value
## (Intercept)  0.75015461 0.82107675 7281  0.913623 0.360945223
## HRS          0.09812558 0.06828664 7281  1.436966 0.150770608
## G.COHES      0.85753355 0.26829047   97  3.196288 0.001879387
## HRS:G.COHES -0.04940481 0.02234020 7281 -2.211476 0.027033848
# Table 9 (Lower panel):
tmod<-lme(fixed=WBEING~W.HRS*G.COHES, random=~1|GRP, data=bh1996)
summary(tmod)$tTable
##                     Value  Std.Error   DF   t-value      p-value
## (Intercept)    1.81207694 0.29827606 7281  6.075167 1.300937e-09
## W.HRS          0.06311464 0.07011506 7281  0.900158 3.680659e-01
## G.COHES        0.31422163 0.09717056   97  3.233712 1.670736e-03
## W.HRS:G.COHES -0.03597709 0.02296508 7281 -1.566600 1.172518e-01

Take-away point 9: When testing cross-level interactions with raw or grand-mean centered variables, it is important to use group-mean centering to verify model results.

Set up a simulation to demonstrate how an interaction involving group means can erroneously be identified as a cross-level interaction and how group-mean centering then corrects the problem. Conceptually, while the code has numerous steps, it functionally just creates 100 group means where 50 group means have a positive X,Y relationship and 50 group means have a negative X,Y relationship. The 100 group means are transformed into individual variables by replicating each group mean 10 times and adding random #error. Thus there are no individual-level relationships except those associated with the group mean. Z.G is a group-level variable where a value of 0 is associated with group-means having a positive relationship, and a value of 1 is associated with the group means having a negative relationship.

set.seed(2500852)
# Create 50 X and Y values to use as group means with
# a POSITIVE correlation (rho) of .50
Y.1<-rnorm(50)
X.1<-.50*Y.1+sqrt(1-.50^2)*rnorm(50)
cor(X.1,Y.1)
## [1] 0.5703528
# Create 50 X and Y values to use as group means with
# a NEGATIVE correlation (rho) of -.50
Y.2<-rnorm(50)
X.2<-(-.50)*Y.2+sqrt(1-(-.50)^2)*rnorm(50)
cor(X.2,Y.2)
## [1] -0.4984011
# Combine the 50 values into a single vector with
# 100 values (half of which are positive and half of
# which are negative)
Y<-c(Y.1,Y.2)
X<-c(X.1,X.2)
# Replicate the 100 values 10 times each for a vector of
# length 1000
Y<-rep(Y,each=10)
X<-rep(X,each=10)
# Add some random error to each of the 1000 values
Y<-Y+rnorm(1000,sd=.1)
X<-X+rnorm(1000,sd=.1)
# Combine the X and Y vectors into a data.frame add Z.G
# as a group-level variable. When Z.G=0 the X and Y group mean
# values were positively correlated; when Z.G=1 X and Yvalues were
# negatively correlated
GMEAN.INT<-data.frame(Y=Y,X=X,Z.G=rep(0:1,each=500),
GRP=rep(1:100,each=10))
# Estimate a cross-level interaction model omitting the random
# slope for X because no slope variability exists and
# include X in random=~X|GRP causes convergence problems.
tmod<-lme(Y~X*Z.G,random=~1|GRP,data=GMEAN.INT)
summary(tmod)$tTable
##                   Value  Std.Error  DF   t-value    p-value
## (Intercept)  0.18754879 0.12945947 898  1.448707 0.14776860
## X            0.07204876 0.04291962 898  1.678690 0.09356029
## Z.G         -0.37171554 0.18307255  98 -2.030428 0.04502244
## X:Z.G       -0.13244156 0.06026903 898 -2.197506 0.02823901
# Create a group-mean variable for X to create a group-mean
# centered variableW.X
TDAT<-aggregate(GMEAN.INT$X,list(GMEAN.INT$GRP),mean)
names(TDAT)<-c("GRP","X.G")
GMEAN.INT<-merge(GMEAN.INT,TDAT,by="GRP")
GMEAN.INT$W.X<-GMEAN.INT$X-GMEAN.INT$X.G
# Estimate a cross-level interaction model using
# Group-mean centered X (W.X). Notice the cross-level interaction
# is no longer significant.
tmod<-lme(Y~W.X*Z.G,random=~1|GRP,data=GMEAN.INT)
summary(tmod)$tTable
##                   Value  Std.Error  DF    t-value    p-value
## (Intercept)  0.19422628 0.13398879 898  1.4495712 0.14752727
## W.X          0.01461219 0.04506368 898  0.3242565 0.74581938
## Z.G         -0.38331127 0.18948876  98 -2.0228708 0.04581048
## W.X:Z.G     -0.02955170 0.06357815 898 -0.4648091 0.64218083
# Show how a fixed-effect OLS model can also recover the cross
# level interaction. Z.G cannot be included as a main-effect because
# the variance is captured by the dummy codes for group, but
# the interaction term X:Z.G can be included.
tmod<-lm(Y~X+X:Z.G+as.factor(GRP),data=GMEAN.INT)
summary(tmod)$coef[c(1,2,102),]
##                Estimate Std. Error     t value      Pr(>|t|)
## (Intercept) -1.47130325 0.04233753 -34.7517536 2.306946e-168
## X            0.01461219 0.04506368   0.3242565  7.458194e-01
## X:Z.G       -0.02955170 0.06357815  -0.4648091  6.421808e-01