Job Stress
Your NAME
2022-10-08
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的package，並使用其中的bh1996的date
data(bh1996, package="multilevel")
dta <- bh1996# str() 查看資料
# 發現資料型態是data frame，並且詳細列出每一個變數的資訊
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 schools and students in each school
# quantile() 用來取百分比
# 發現GRP的0,25,50,75,100(%)分別為15,43.5,64,94,226.0
with(dta, table(GRP)) |> quantile() 0% 25% 50% 75% 100%
15.0 43.5 64.0 94.0 226.0 # 將dta導到更細的分類
dta <- dta %>% group_by(GRP) %>% mutate(nc=n()) %>% ungroup()# 發現nc的0,25,50,75,100(%)分別為15,60.5,94,166.25,324
with(dta, table(nc)) |> quantile() 0% 25% 50% 75% 100%
15.00 60.50 94.00 166.25 324.00 4 Visualization
# reorder用在繪圖，如ggplot繪製條形圖，可使x軸按y軸數值大小排序
# 根據中位數，GRP按照WBEING數值大小排序
# coord_flip為橫縱軸位置互換，沒有特殊參數
# 所以X軸變Well-being score，Y軸變Group ID
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
# 查看隨機模型的效果
# Y = WBEING，X = GRP的截距項(平均值)
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() 用來計算變異估計
VarCorr(m0) Groups Name Std.Dev.
GRP (Intercept) 0.1892
Residual 0.8885 # 觀察m0的內部一致性程度
# 內部一致性分數僅0.043，非常低
performance::icc(m0)# Intraclass Correlation Coefficient
Adjusted ICC: 0.043
Unadjusted ICC: 0.0436 Work hours - individual and group
# m1包含固定效果與隨機效果
# 效果並不顯著
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.21368 22.19
HRS -0.04646 0.00489 -9.51
G.HRS -0.12693 0.01940 -6.54VarCorr(m1) Groups Name Std.Dev.
GRP (Intercept) 0.1164
Residual 0.8832 # m1內部一致性
# 內部一致性僅0.017(0.016)，非常低
performance::icc(m1)# Intraclass Correlation Coefficient
Adjusted ICC: 0.017
Unadjusted ICC: 0.0167 Leadership consideration - individual
# x=LEAD, y=WBEING, 進行畫圖
# facet_wrap是基於一個因子進行設置
# facets表示形式爲：~變量（~單元格）
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與m1不同，加入LEAD
# 並在隨機效果更改為(LEAD | GRP)，看GRP中LEAD所造成的隨機影響
# 結果發現，效果不顯著
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.00923328 (tol = 0.002, component 1)# 使用optimx，將資料優化
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加入交互作用的效果，LEAD:G.HRS
# 結果仍不顯著
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.68VarCorr(m3) Groups Name Std.Dev. Corr
GRP (Intercept) 0.36876
LEAD 0.09929 -0.97
Residual 0.80081 # ICC僅0.030、0.024，非常低
performance::icc(m3)# Intraclass Correlation Coefficient
Adjusted ICC: 0.030
Unadjusted ICC: 0.0248 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.