Job Stress
Lynn
2022-10-07
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
data(bh1996, package="multilevel")
dta <- bh1996str(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
with(dta, table(GRP)) |> quantile() 0% 25% 50% 75% 100%
15.0 43.5 64.0 94.0 226.0 library(dplyr)
dta <- dta %>% group_by(GRP) %>% mutate(nc=n()) %>% ungroup()
# mutate() 在現存變數的函數創造新的變數
# Always ungroup() when you’ve finished with your calculations.with(dta, table(nc)) |> quantile() 0% 25% 50% 75% 100%
15.00 60.50 94.00 166.25 324.00 4 Visualization
ggplot(dta, aes(x=reorder(factor(GRP), WBEING, median), WBEING)) +
geom_boxplot() +
coord_flip()+
labs(x="Group ID",
y="Well-being score")
# coord_flip() Flipped cartesian coordinates so that horizontal becomes vertical, and vertical, horizontal. 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# The estimated WBEING score is 2.77 and the variances are 0.36 and 0.79 at the GRP levels.VarCorr(m0) Groups Name Std.Dev.
GRP (Intercept) 0.1892
Residual 0.8885 # This function calculates the estimated variances, standard deviations, and correlations between the random-effects terms in a linear mixed-effects model, of class "lme", or a nonlinear mixed-effects model, of class "nlme".performance::icc(m0)# Intraclass Correlation Coefficient
Adjusted ICC: 0.043
Unadjusted ICC: 0.043# At the GRP level, the intra-class correlation is 0.43. These values indicate the correlation in WBEING between two people in the same GRP, and between two GRPs on the same people. It means that 4.3% of the variation in WBEING can be attributed to attending the same GRP and 4.3% to the people themselves (which includes being in the same GRP).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# The estimated WBEING score is 4.74 and the variances are 0.01 and 0.78 at the GRP levels.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.0167 Leadership consideration - individual
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)# 嘗試使用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 <- 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# Optimization methods are used in many areas of study to find solutions that maximize or minimize some study parameters
# method="nlminb": Uses the function nlminb for optimisation, so that optimising a Mk2/Mkn likelihood function behaves as similarly as possible to ape's ace function.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.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.