1 Introduction

The data set contains four variables:

  • Cohesion (COHES),
  • Leadership Climate (LEAD),
  • Well-Being (WBEING) and
  • Work Hours (HRS).

Each of these variables has two variants - a group mean version that replicates each group mean for every individual in the group, and a within-group version where the group mean is subtracted from each individual response (i.e., a group-mean centered or demeaned variable). The group mean version is designated with a G. (e.g., G.HRS), and the within-group version is designated with a W. (e.g., W.HRS).

2 Data

#從multilevel中載入資料並命名為dat
data(bh1996, package="multilevel")
dta <- bh1996
#查看資料結構
str(dta)
'data.frame':   7382 obs. of  13 variables:
 $ GRP     : num  1 1 1 1 1 1 1 1 1 1 ...
 $ COHES   : num  3.75 3.25 3.38 3.75 4 ...
 $ G.COHES : num  2.93 2.93 2.93 2.93 2.93 ...
 $ W.COHES : num  0.821 0.321 0.446 0.821 1.071 ...
 $ LEAD    : num  3.18 3 3.64 3.36 3.55 ...
 $ G.LEAD  : num  2.93 2.93 2.93 2.93 2.93 ...
 $ W.LEAD  : num  0.2527 0.0709 0.7072 0.4345 0.6163 ...
 $ HRS     : num  12 11 12 9 7 8 9 8 12 11 ...
 $ G.HRS   : num  11 11 11 11 11 ...
 $ W.HRS   : num  1.027 0.027 1.027 -1.973 -3.973 ...
 $ WBEING  : num  2.11 3.33 2.11 4.39 1.72 ...
 $ G.WBEING: num  2.79 2.79 2.79 2.79 2.79 ...
 $ W.WBEING: num  -0.682 0.54 -0.682 1.596 -1.071 ...

3 Tables

# Number of companies and data in each company
length(unique(dta$GRP))
[1] 99
with(dta, table(GRP)) |> quantile()
   0%   25%   50%   75%  100% 
 15.0  43.5  64.0  94.0 226.0 
#使用GRP進行分類
dta <- dta %>% group_by(GRP) %>% mutate(nc=n()) %>% ungroup()
with(dta, table(nc)) |> quantile()
    0%    25%    50%    75%   100% 
 15.00  60.50  94.00 166.25 324.00

我們可以知道，本數據包含99間公司，公司員工中位數為94，每個公司有15到324個觀值，共7382筆觀察值。

4 Visualization

# 以每家公司幸福感分數的中位數進行排序畫圖
#X軸為幸福感分數，Y軸為公司
ggplot(dta, aes(x=reorder(factor(GRP), WBEING, median), WBEING)) +
  geom_boxplot() +
  coord_flip()+
  labs(x="Group ID",
       y="Well-being score")

5 Null model - random intercepts only

#估計公司在幸福感分書的隨機效果
m0 <- lme4::lmer(WBEING ~ (1 | GRP), data=dta)
summary(m0, corr=FALSE)
Linear mixed model fit by REML ['lmerMod']
Formula: WBEING ~ (1 | GRP)
   Data: dta

REML criterion at convergence: 19347.3

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.322 -0.648  0.031  0.718  2.667 

Random effects:
 Groups   Name        Variance Std.Dev.
 GRP      (Intercept) 0.0358   0.189   
 Residual             0.7895   0.889   
Number of obs: 7382, groups:  GRP, 99

Fixed effects:
            Estimate Std. Error t value
(Intercept)   2.7743     0.0222     125
VarCorr(m0)
 Groups   Name        Std.Dev.
 GRP      (Intercept) 0.1892  
 Residual             0.8885  
performance::icc(m0)
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.043
  Unadjusted ICC: 0.043

5.1 結論

1.公司的幸福感平均為2.7743，標準差為0.1892。

2.公司中，員工幸福感相關為0.043。

6 Work hours - individual and group

#以個人及公司的工作時間為變項，估計公司的截距隨機效果
m1 <- lme4::lmer(WBEING ~ HRS + G.HRS + (1 | GRP) , data=dta)
summary(m1, corr=FALSE)
Linear mixed model fit by REML ['lmerMod']
Formula: WBEING ~ HRS + G.HRS + (1 | GRP)
   Data: dta

REML criterion at convergence: 19212.3

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.353 -0.650  0.038  0.713  2.709 

Random effects:
 Groups   Name        Variance Std.Dev.
 GRP      (Intercept) 0.0135   0.116   
 Residual             0.7801   0.883   
Number of obs: 7382, groups:  GRP, 99

Fixed effects:
            Estimate Std. Error t value
(Intercept)  4.74083    0.21367   22.19
HRS         -0.04646    0.00489   -9.51
G.HRS       -0.12693    0.01940   -6.54
VarCorr(m1)
 Groups   Name        Std.Dev.
 GRP      (Intercept) 0.1164  
 Residual             0.8832  
performance::icc(m1)
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.017
  Unadjusted ICC: 0.016

6.1 結論

1.所有員工的幸福感平均為4.741，標準差為0.1164。

2.員工每增加一小時的工作時間(HRS)，幸福感會下降0.0464。

3.公司平均增加一小時的工作時間(G.HRS)，幸福感會下降0.1269。

4.同一公司中，員工幸福感相關為0.016。

7 Leadership consideration - individual

#為員工數大於100的公司畫圖
#X軸為領導氛圍，Y軸為幸福感
ggplot(subset(dta, nc > 100), aes(LEAD, WBEING))+
  stat_smooth(method='lm', formula=y~x, se=FALSE,
              size=rel(.5), col=1)+
  geom_point(size=rel(.5), col=8, alpha=.5)+
  facet_wrap(. ~ GRP)+
  labs(x="Leadership",
       y="Well-being")+
  theme_minimal()

由圖中可以看到，領導氛圍的提升，幸福感也跟著上升。

#以個人與公司工作時間、領導氛圍作為變項，估計公司與領導氛圍的隨機效果
m2 <- lme4::lmer(WBEING ~ HRS + LEAD + G.HRS + (LEAD | GRP), data=dta)
summary(m2, corr=FALSE)
Linear mixed model fit by REML ['lmerMod']
Formula: WBEING ~ HRS + LEAD + G.HRS + (LEAD | GRP)
   Data: dta

REML criterion at convergence: 17822.6

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.871 -0.656  0.041  0.697  3.958 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 GRP      (Intercept) 0.1466   0.383         
          LEAD        0.0107   0.103    -0.97
 Residual             0.6413   0.801         
Number of obs: 7382, groups:  GRP, 99

Fixed effects:
            Estimate Std. Error t value
(Intercept)  2.46418    0.20754   11.87
HRS         -0.02848    0.00447   -6.37
LEAD         0.49454    0.01687   29.31
G.HRS       -0.07057    0.01782   -3.96
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00924816 (tol = 0.002, component 1)
#模型無法收斂，因此進行優化
library(optimx)
update(m2, control = lmerControl(optimizer= "optimx",
                                 optCtrl  = list(method="nlminb")))
Linear mixed model fit by REML ['lmerMod']
Formula: WBEING ~ HRS + LEAD + G.HRS + (LEAD | GRP)
   Data: dta
REML criterion at convergence: 17822.6
Random effects:
 Groups   Name        Std.Dev. Corr 
 GRP      (Intercept) 0.383         
          LEAD        0.103    -0.97
 Residual             0.801         
Number of obs: 7382, groups:  GRP, 99
Fixed Effects:
(Intercept)          HRS         LEAD        G.HRS  
     2.4641      -0.0285       0.4945      -0.0706  
m3 <- lme4::lmer(WBEING ~ HRS + LEAD + G.HRS + LEAD:G.HRS + (LEAD | GRP), 
                 data=dta,
                 control = lmerControl(optimizer= "optimx",
                                       optCtrl  = list(method="nlminb")))
summary(m3, corr=FALSE)
Linear mixed model fit by REML ['lmerMod']
Formula: WBEING ~ HRS + LEAD + G.HRS + LEAD:G.HRS + (LEAD | GRP)
   Data: dta
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 17825.9

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.837 -0.660  0.041  0.695  3.953 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 GRP      (Intercept) 0.13598  0.3688        
          LEAD        0.00986  0.0993   -0.97
 Residual             0.64129  0.8008        
Number of obs: 7382, groups:  GRP, 99

Fixed effects:
            Estimate Std. Error t value
(Intercept)  3.64326    0.73255    4.97
HRS         -0.02856    0.00447   -6.39
LEAD         0.12895    0.21881    0.59
G.HRS       -0.17402    0.06415   -2.71
LEAD:G.HRS   0.03217    0.01919    1.68
VarCorr(m3)
 Groups   Name        Std.Dev. Corr 
 GRP      (Intercept) 0.36876       
          LEAD        0.09929  -0.97
 Residual             0.80081       
performance::icc(m3)
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.030
  Unadjusted ICC: 0.024

7.1 結論

1.所有員工的幸福感平均值為2.464，標準差為0.3688。

2.員工每增加一小時的工作時間(HRS)，員工的幸福感會下降0.0286。

3.公司平均增加一小時的工作時間(G.HRS)，員工的幸福感會下降0.1740。

4.領導氛圍每增加一單位，員工的幸福感會上升0.1289，公司間標準差為0.0992。

5.隨機效果及領導氣氛相關為-0.967，因此平均幸福感越低的公司，其越容易受領導氛圍影響而提升幸福感。

6.同一公司中，員工幸福感相關為0.024。

8 References

Bliese, P. D. & Halverson, R. R. (1996). Individual and nomothetic models of job stress: An examination of work hours, cohesion, and well-being. Journal of Applied Social Psychology, 26, 1171-1189.