The daily costs of workaholism

Supplementary material S6: Multilevel regression modeling

Aims and content

The present document includes the analytical steps implemented to model the within-individual relationships between workaholism WHLSM and the focused strain outcomes (i.e., afternoon and evening blood pressure SBP and DBP, emotional exhaustion EE, and sleep disturbances SD) and the within-individual interactions between WHLSM and psychological detachment PD. The analyses are conducted on the daily diary data collected with the Qualtrics platform (Qualtrics, Seattle, WA, USA) from an heterogeneous samples of workers over two weeks, pre-processed as shown in Supplementary Material S3, and aggregated as shown in Supplementary Material S4.

Here, we remove all objects from the R global environment.

# removing all objets from the workspace
rm(list=ls())

The following R packages are used in this document (see References section):

# required packages
packages <- c("lme4","MuMIn","sjPlot","plyr","car","fitdistrplus","ggplot2","gridExtra","influence.ME","mediation","psych","knitr")

# generate packages references
knitr::write_bib(c(.packages(), packages),"packagesMod.bib")

# # run to install missing packages
# xfun::pkg_attach2(packages, message = FALSE); rm(list=ls())

1. Data reading

First, we read daily diary exported from the previous step (see Supplementary Material S4, and the model formulas with the covariates selected from Supplementary Material S5. Moreover, we derive the the preliminary questionnaire dataset prelqs from the diary dataset.

# reading data
load("DATI/diary_aggregated.RData") # daily diary data
load("DATI/mformulas.RData") # selected covariates

# deriving prelqs from diary data (only including variables from the preliminary questionnaire)
prelqs <- diary[!duplicated(diary$ID),c(1,which(colnames(diary)=="gender"):ncol(diary))]

# original sample sizes
cat("diary:",nrow(diary),"responses from",nlevels(diary$ID),"participants")

## diary: 1084 responses from 135 participants

cat("prelqs:",nrow(prelqs),"responses from",nlevels(prelqs$ID),"participants")

## prelqs: 135 responses from 135 participants

2. Data filtering

As we pre-registered here, we filter the data based on participant compliance with the protocol, that is we exclude the participants with less than 3 full days of participation (i.e., with nonmissing response to the afternoon, evening, and morning questionnaire).

# filtering participants with less than 3 observations
clean <- diary[0,]
for(ID in levels(diary$ID)){ 
  if(nrow(diary[diary$ID==ID & diary$aft==1 & diary$eve==1 & diary$mor==1,]) >= 3){ 
    clean <- rbind(clean,diary[diary$ID==ID,]) }}
clean$ID <- as.factor(as.character(clean$ID)) # resetting ID levels
clean_prelqs <- prelqs[prelqs$ID %in% levels(clean$ID),] # filtering prelqs data
clean_prelqs$ID <- as.factor(as.character(clean_prelqs$ID)) # resetting ID levels
cat("diary: Excluded",nlevels(diary$ID)-nlevels(clean$ID),"participants and",nrow(diary)-nrow(clean),"observations")

## diary: Excluded 21 participants and 126 observations

# updating sample sizes
cat("diary:",nrow(clean),"responses from",nlevels(clean$ID),"participants")

## diary: 958 responses from 114 participants

cat("prelqs:",nrow(clean_prelqs),"responses from",nlevels(clean_prelqs$ID),"participants")

## prelqs: 114 responses from 114 participants

2.1. BP data filtering

In addition to participant compliance, we pre-registered the exclusion of all participants reporting taking blood pressure medications bp_drugs or suffering from a cardiovascular dysfunction cv_dysf from the analyses of blood pressure. Here, we exclude such participants showing substantially higher blood pressure (see section 3.5 of Supplementary Material S5). In contrast, deviating from the pre-registration (see Supplementary Material S7), we do not exclude those taking psychoactive psy_drugs or hormonal medication horm_drugs, or those reporting sleep dysfunctions sleep_dysf, as we did not find substantial differences in their blood pressure values (see section 3.5 of Supplementary Material S5). Instead, we exclude these participants as a robustness check (see section 4).

# filtering participants with bp_drugs or cv_dysf
cleanBP <- clean[clean$bp_drugs=="No" & clean$cv_dysf=="No",]
cleanBP$ID <- as.factor(as.character(cleanBP$ID)) # resetting ID levels
cleanBP_prelqs <- prelqs[prelqs$ID %in% levels(cleanBP$ID),] # filtering prelqs data
cleanBP_prelqs$ID <- as.factor(as.character(cleanBP_prelqs$ID)) # resetting ID levels
cat("diary: Excluded further",nlevels(clean$ID)-nlevels(cleanBP$ID),"participants and",
    nrow(clean)-nrow(cleanBP),"observations from BP analyses")

## diary: Excluded further 8 participants and 60 observations from BP analyses

# updating sample sizes
cat("diary (BP):",nrow(cleanBP),"responses from",nlevels(cleanBP$ID),"participants")

## diary (BP): 898 responses from 106 participants

cat("prelqs (BP):",nrow(cleanBP_prelqs),"responses from",nlevels(cleanBP_prelqs$ID),"participants")

## prelqs (BP): 106 responses from 106 participants

2.2. Sleep data filtering

Finally, we also pre-registered the exclusion of all participants reporting sleep dysfunctions from the analyses of sleep disturbances. However, the number of such participants is relatively high, whereas they do not seem to show substantial differences in sleep disturbances comparing to the other participants (see section 3.5 of Supplementary Material S5). Thus, deviating from the pre-registration (see Supplementary Material S7), we do not exclude such participants from the main analyses. Instead, we exclude them as a robustness check (see section 4.5).

# number of participants reporting sleep dysfunctions
summary(clean[!duplicated(clean$ID),"sleep_dysf"])

##  No Yes 
## 100  14

3. Multilevel modeling

Here, we specify, compare, and inspect the results of a series of multilevel models for each of the following pre-registered time-varying outcomes:

• afternoon systolic SBP_aft and dyastolic blood pressure DBP_aft

• evening systolic SBP_eve and dyastolic blood pressure DBP_eve

• emotional exhaustion EE

• sleep disturbances SD

For each outcome, we implement a hierarchical regression with the following steps:

1. m0: null model only including the sample intercept, the level-2 variability around the intercept, and the residual term

2. m1: including the covariates selected in the previous step (see Supplementary Material S5), in addition to psychological detachment PD (i.e., only for evening and next morning outcomes), and the level-2 component of workaholism WHLSM.cm

3. m2: including the level-1 component of workaholism WHLSM.mc

4. m3: including the level-1 interactions between WHLSM.mc and PD (i.e., only for outcomes measured in the evening or the following morning).

The following packages are used to optimize the analyses:

library(lme4); library(MuMIn); library(sjPlot); library(plyr); library(car); library(fitdistrplus)
library(gridExtra); library(influence.ME)

3.1. Blood pressure

First, we analyse systolic and diastolic blood pressure. As a main confirmatory analysis, we evaluate the relationships between daily levels of state workaholism and daily averages of blood pressure (i.e., average of afternoon, evening, and next morning measurements). Moreover, as a supplementary pre-registered analysis (see Supplementary Material S7) we apply the same procedure for each time point, that is we analyze the relationship between workaholism and afternoon, evening, and morning blood pressure measurement, respectively. Finally, based on the results, we conduct pre-registered exploratory mediation analyses on time-specific blood pressure measurements.

Daily averages of systolic SBP and diastolic blood pressure DBP are predicted by concurrent state workaholism WHLSM.mc, in addition to trait workaholism WHLSM.cm, and three covariates selected from the previous step (see Supplementary Material S5), namely gender, age, and body mass index BMI.

mformulas[1] # covariates selected for SBP

## [1] "SBP_aft ~ gender + age + BMI"

mformulas[2] # covariates selected for DBP

## [1] "DBP_aft ~ gender + age + BMI"

3.1.1. Afternoon BP

3.1.1.1. Data preparation

First, we prepare the data for the analyses by removing all cases of missing responses in the dependent variable or any predictor or covariate (list-wise deletion), by centering level-2 continuous predictors on the grand mean (grand-mean-centering), and by centering level-1 continuous predictors on the individual mean (person-mean-centering).

# list-wise deletion
cleanBP_aft <- as.data.frame(na.omit(cleanBP[,c("ID","SBP_aft","DBP_aft", # grouping and dependent variables
                                                "gender","age","BMI","WHLSM", # core predictors
                                                "WE","WC","sleep_dysf","psy_drugs","horm_drugs", # for robustness checks
                                                "confounders_aft","flagBP_aft","flagTime","careless",
                                                "position","children")])) 
cleanBP_aft$ID <- as.factor(as.character(cleanBP_aft$ID)) # resetting participant identifier levels
cat("Considering",nrow(cleanBP_aft),"complete obs from",nlevels(as.factor(as.character(cleanBP_aft$ID))),"participants")

## Considering 787 complete obs from 106 participants

# person-mean-centering lv-1 continuous predictors
wide <- cleanBP_aft[!duplicated(cleanBP_aft$ID),] # wide-form dataset
for(Var in c("WHLSM","WE","WC")){
  wide <- cbind(wide,aggregate(cleanBP_aft[,Var],list(cleanBP_aft$ID),mean)[,2]) # individual means
  colnames(wide)[ncol(wide)] <- paste0(Var,".cm")
  cleanBP_aft <- join(cleanBP_aft,wide[,c("ID",paste0(Var,".cm"))],by="ID",type="left") # joining to long-form df
  cleanBP_aft[,paste0(Var,".mc")] <- cleanBP_aft[,Var] - cleanBP_aft[,paste0(Var,".cm")] } # mean-centered scores

# grand-mean-centering lv-2 continuous predictors
for(Var in c("age","BMI","WHLSM.cm")){ cleanBP_aft[,paste0(Var,".gmc")] <- cleanBP_aft[,Var] - mean(wide[,Var]) }

# showing data
cleanBP_aft[1:3,] # first three rows

3.1.1.2. Model fit

Here, we fit the multilevel models to the selected data using the default restricted maximum likelihood estimator (REML).

# m0: null model
m0_SBP_aft <- lmer(SBP_aft ~ (1|ID), # only fixed and random intercept + residual term
                   data=cleanBP_aft)

# m1: covariates
m1_SBP_aft <- lmer(SBP_aft ~ gender + age.gmc + BMI.gmc + WHLSM.cm.gmc + (1|ID), # covariates
                   data=cleanBP_aft)

# m2: state workaholism
m2_SBP_aft <- lmer(SBP_aft ~ gender + age.gmc + BMI.gmc + WHLSM.cm.gmc + WHLSM.mc + (1|ID),
                   data=cleanBP_aft)

The same models are specified for diastolic blood pressure.

# m0: null model
m0_DBP_aft <- lmer(DBP_aft ~ (1|ID), # only fixed and random intercept + residual term
                   data=cleanBP_aft)

# m1: covariates
m1_DBP_aft <- lmer(DBP_aft ~ gender + age.gmc + BMI.gmc + WHLSM.cm.gmc + (1|ID), # covariates
                   data=cleanBP_aft)

# m2: state workaholism
m2_DBP_aft <- lmer(DBP_aft ~ gender + age.gmc + BMI.gmc + WHLSM.cm.gmc + WHLSM.mc + (1|ID),
                   data=cleanBP_aft)

From the previous chunks, we see that all models converged without problems. Here, we inspect the diagnostics (i.e., normality of residual and random effect distributions, homoscedasticity, and multicollinearity) of the most complex model m2.bis. Influential cases are analyzed in a dedicated section below.

SBP_aft

Model m2.bis shows some deviation from normality especially in the lower tail of the distribution of residuals and both tails of the distribution of random effects. Particularly, participants S082 and S096 are associated with the highest extreme deviations from the distributions of both random effects, and will be removed as a robustness check (see section 4.1). Besides that, we can see that the homoscedsticity assumption holds and that none of the variance inflation factors (VIFs) shows extreme values, ruling out the risk of multicollinearity.

# normality and homoscedasticity
p <- plot_model(m2_SBP_aft,type="diag",dot.size=1)
p[[2]] <- p[[2]]$ID
plot_grid(p,tags=TRUE,margin=c(0,0,0,0))

# participant with highest random effects (i.e., BLUPS)
re <- ranef(m2_SBP_aft)$ID
re[re$WHLSM.mc==max(re$WHLSM.mc)|re$`(Intercept)`==max(re$`(Intercept)`),]

## numeric(0)

# homoscedasticity and multicollinearity
par(mfrow=c(1,2))
for(Var in c("gender")){ boxplot(resid(m2_SBP_aft) ~ cleanBP_aft[,Var],main=paste("Residuals by",Var)) }
barplot(vif(m2_SBP_aft),main="VIF Values",xlim=c(0,10),las=2,horiz=TRUE) # variance inflation factors (VIFs)
abline(v = 5, lwd = 5, lty = 2)

Here, we better inspect the residual distribution and the fit of models specified with alternative family distributions. We can see that none of the alternative families substantially better approximate the distribution of model residuals, with the log-transformed solution (i.e., normal distribution with log-transformed dependent variable) showing the best fit. Yet, it is not so better than the original model. Thus, we initially rely on the normal distribution and then consider the logarithmic transformation as a robustness check (see section 4.1).

# inspecting residual distribution
descdist(resid(m2_SBP_aft)) # unknown best-fit distribution

## summary statistics
## ------
## min:  -36.76048   max:  49.42156 
## median:  -0.4067047 
## mean:  4.281329e-14 
## estimated sd:  9.662342 
## estimated skewness:  0.3987242 
## estimated kurtosis:  4.901672

# fitting model with alternative families
models <- list(
  m2_SBP_aft,
  glmer(formula=formula(m2_SBP_aft),family=Gamma(link="log"),data=cleanBP_aft), # gamma log (doesn't converge)
  glmer(formula=formula(m2_SBP_aft),family=Gamma(link="identity"),data=cleanBP_aft), # gamma id (doesn't converge)
  glmer(formula=formula(m2_SBP_aft),family=gaussian(link="log"),data=cleanBP_aft), # log-normal (singular fit)
  lmer(formula=as.formula(paste("log(SBP_aft) ~",as.character(formula(m2_SBP_aft))[3])),data=cleanBP_aft)) # log transf

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00424267 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00352728 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

# normal Q-Q plot of model residuals
par(mfrow=c(1,5))
for(i in 1:length(models)){ 
  qqnorm(resid(models[[i]]),main=c("Norm","Gamma-log","Gamma-id","log-norm","log-transf")[i]); qqline(resid(models[[i]]))}

DBP_aft

Model m2.bis shows some deviation from normality in both tails of the distribution of residuals and random effects. Particularly, participant S082 is associated with the highest extreme deviation from the distributions of both random effects, and will be removed as a robustness check (see section 4.1). Besides that, we can see that the homoscedsticity assumption holds and that none of the variance inflation factors (VIFs) shows extreme values, ruling out the risk of multicollinearity.

# normality and homoscedasticity
p <- plot_model(m2_DBP_aft,type="diag",dot.size=1)
p[[2]] <- p[[2]]$ID
plot_grid(p,tags=TRUE,margin=c(0,0,0,0))

# participant with highest random effects (i.e., BLUPS)
re <- ranef(m2_DBP_aft)$ID
re[re$WHLSM.mc==max(re$WHLSM.mc)|re$`(Intercept)`==max(re$`(Intercept)`),]

## numeric(0)

# homoscedasticity and multicollinearity
par(mfrow=c(1,2))
for(Var in c("gender")){ boxplot(resid(m2_DBP_aft) ~ cleanBP_aft[,Var],main=paste("Residuals by",Var)) }
barplot(vif(m2_DBP_aft),main="VIF Values",xlim=c(0,10),las=2,horiz=TRUE) # variance inflation factors (VIFs)
abline(v = 5, lwd = 5, lty = 2)

Here, we better inspect the residual distribution and the fit of models specified with alternative family distributions. We can see that some of alternative families slightly better approximate the distribution of model residuals, although not reaching convergence. The log-transformed solution (i.e., normal distribution with log-transformed dependent variable) shows the best fit while reaching convergence. Yet, it is not so better than the original model. Thus, we initially rely on the normal distribution and then consider the logarithmic transformation as a robustness check.

# inspecting residual distribution
descdist(resid(m2_DBP_aft)) # unknown best-fit distribution

## summary statistics
## ------
## min:  -30.20225   max:  36.78257 
## median:  -0.01672168 
## mean:  4.119735e-14 
## estimated sd:  8.413706 
## estimated skewness:  0.1099257 
## estimated kurtosis:  4.379899

# fitting model with alternative families
models <- list(
  m2_DBP_aft,
  glmer(formula=formula(m2_DBP_aft),family=Gamma(link="log"),data=cleanBP_aft), # gamma log (doesn't converge)
  glmer(formula=formula(m2_DBP_aft),family=Gamma(link="identity"),data=cleanBP_aft), # gamma id (doesn't converge)
  glmer(formula=formula(m2_DBP_aft),family=gaussian(link="log"),data=cleanBP_aft), # log-normal (singular fit)
  lmer(formula=as.formula(paste("log(SBP_aft) ~",as.character(formula(m2_DBP_aft))[3])),data=cleanBP_aft)) # log transf

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00264134 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00739431 (tol = 0.002, component 1)

## boundary (singular) fit: see help('isSingular')

# normal Q-Q plot of model residuals
par(mfrow=c(1,5))
for(i in 1:length(models)){ 
  qqnorm(resid(models[[i]]),main=c("Norm","Gamma-log","Gamma-id","log-norm","log-transf")[i]); qqline(resid(models[[i]]))}

3.1.1.3. Results

Here, we compare the specified models based on the Akaike weight and the likelihood ratio test (with type-I error set to p < .05), and we inspect the results of the selected model(s).

SBP_aft

We can see that the inclusion of state WHLSM is associated with stronger evidence and significantly higher likelihood than models including less predictors (Aw = .99, \(\chi^2\)(1) = 15.07, p < .001), with model m2 being selected as the best model.

# Akaike weight adding one model at time
Weights(AIC(m0_SBP_aft,m1_SBP_aft)) # covariates: better

##  model weights 
## [1] 0 1

Weights(AIC(m0_SBP_aft,m1_SBP_aft,m2_SBP_aft)) # state workaholism: better

##  model weights 
## [1] 0.000 0.001 0.999

# Likelihood ratio test
anova(m1_SBP_aft,m2_SBP_aft) # best model is m2

## refitting model(s) with ML (instead of REML)

Here, we inspect the coefficients estimated by the selected model m2.bis and those estimated by more parsimonious models. We can see that state WHLSM.mc is positively related to BP_aft, whereas trait WHLSM.cm.gmc is not. Among the included covariates, gender, age and BMI, but not trait WHLSM.cm.gmc predict higher SBP_aft.

# regression table
tab_model(m1_SBP_aft,m2_SBP_aft,
          dv.labels=c("Covariates","State WHLSM"),
          show.icc=FALSE,show.p=FALSE,show.se=TRUE,show.r2=FALSE,show.ci=FALSE,
          collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Covariates		State WHLSM
Predictors	b (SE)	t	b (SE)	t
(Intercept)	118.25 (1.97)	59.96	118.25 (1.97)	59.96
gender [M]	5.39 (2.92)	1.84	5.39 (2.92)	1.84
age gmc	0.41 (0.12)	3.43	0.41 (0.12)	3.43
BMI gmc	1.10 (0.43)	2.54	1.10 (0.43)	2.54
WHLSM cm gmc	1.33 (1.14)	1.16	1.33 (1.14)	1.16
WHLSM mc			1.56 (0.40)	3.90
Random Effects
σ2	108.96		106.74
τ00	195.27 ID		195.55 ID
N	106 ID		106 ID
Observations	787		787

# plotting main effects from selected model
grid.arrange(plot_model(m2_SBP_aft,type="pred",terms="WHLSM.cm.gmc"), # trait WHLSM
             plot_model(m2_SBP_aft,type="pred",terms="WHLSM.mc"),nrow=1) # state WHLSM

DBP_aft

Results are in line with those found for SBP_aft: the inclusion of state WHLSM is associated with stronger evidence and significantly higher likelihood than models including less predictors (Aw = .99, \(\chi^2\)(1) = 13.84, p < .001), with model m2 being selected as the best model.

# Akaike weight adding one model at time
Weights(AIC(m0_DBP_aft,m1_DBP_aft)) # covariates: better

##  model weights 
## [1] 0 1

Weights(AIC(m0_DBP_aft,m1_DBP_aft,m2_DBP_aft)) # state workaholism: better

##  model weights 
## [1] 0.000 0.003 0.997

# Likelihood ratio test
anova(m1_DBP_aft,m2_DBP_aft) # best model is m2

## refitting model(s) with ML (instead of REML)

Here, we inspect the coefficients estimated by the selected model m2 and those estimated by the other models. In line with the results found for SBP_aft, we can see that state WHLSM.mc is positively related to DBP_aft, whereas trait WHLSM.cm.gmc is not. Among the included covariates, both age and BMI, but not gender predict higher DBP_aft.

# regression table
tab_model(m1_DBP_aft,m2_DBP_aft,
          dv.labels=c("Covariates","State WHLSM"),
          show.icc=FALSE,show.p=FALSE,show.se=TRUE,show.r2=FALSE,show.ci=FALSE,
          collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Covariates		State WHLSM
Predictors	b (SE)	t	b (SE)	t
(Intercept)	77.32 (1.46)	52.80	77.32 (1.46)	52.81
gender [M]	-0.38 (2.17)	-0.18	-0.38 (2.17)	-0.18
age gmc	0.28 (0.09)	3.07	0.28 (0.09)	3.07
BMI gmc	0.96 (0.32)	2.98	0.96 (0.32)	2.98
WHLSM cm gmc	1.11 (0.85)	1.31	1.11 (0.85)	1.31
WHLSM mc			1.30 (0.35)	3.74
Random Effects
σ2	82.14		80.61
τ00	104.48 ID		104.66 ID
N	106 ID		106 ID
Observations	787		787

# plotting main effects from selected model
grid.arrange(plot_model(m2_DBP_aft,type="pred",terms="WHLSM.cm.gmc"), # trait WHLSM
             plot_model(m2_DBP_aft,type="pred",terms="WHLSM.mc"),nrow=1) # state WHLSM

3.1.1.4. Influential cases

Here, we evaluate the presence of influential cases in the selected model m2. Specifically, the Cook’s distance is considered as the main measure of individual-level (i.e., participant) influence on the estimated parameters, and it is recomputed by progressively excluding the most influential participants (i.e., based on the rule-of-thumb of 4/N) until all extreme values are removed.

3.1.1.4.1. Cook’s distance

SBP_aft

We can see that participants S096 and S082 are potentially influential cases. Note that these are the same participants showing the most extreme deviations from the distributions of random effects (see section 3.1.2).

# cook's distance on the whole sample
infl <- influence(m2_SBP_aft,"ID")
plot(infl,which="cook",cutoff=4/nlevels(cleanBP_aft$ID),xlab="Cook distance",ylab="ID",sort=TRUE)

# progressively excluding participants
infl <- list(
  influence(exclude.influence(m2_SBP_aft,"ID","S082"),"ID"),
  influence(exclude.influence(m2_SBP_aft,"ID",c("S082","S096")),"ID"))
for(i in 1:length(infl)){ 
  plot(infl[[i]],which="cook",cutoff=4/(nlevels(cleanBP_aft$ID)-i),xlab="Cook distance",ylab="ID",sort=TRUE) }

DBP_aft

We can see that participant S082 is a potentially influential case. Note that this is the same participants showing the most extreme deviations from the distributions of random effects (see section 3.1.2).

# cook's distance on the whole sample
infl <- influence(m2_DBP_aft,"ID")
plot(infl,which="cook",cutoff=4/nlevels(cleanBP_aft$ID),xlab="Cook distance",ylab="ID",sort=TRUE)

# progressively excluding participants
infl <- list(influence(exclude.influence(m2_DBP_aft,"ID","S082"),"ID"))
for(i in 1:length(infl)){ 
  plot(infl[[i]],which="cook",cutoff=4/(nlevels(cleanBP_aft$ID)-i),xlab="Cook distance",ylab="ID",sort=TRUE) }

3.1.1.4.2. Coefficient change

SBP_aft

Here, we inspect the magnitude of the changes in the estimated coefficients after the removal of potentially influential cases. We can see that the coefficients estimated by the updated model do not substantially differ from those estimated by the original model. Thus, we choose to rely on the results obtained with the full sample.

# refitting model without influential cases
m2_SBP_aft.noInfl <- update(m2_SBP_aft,data=cleanBP_aft[!cleanBP_aft$ID%in%c("S096","S082"),])

# plotting coefficients original vs. updated model
plot_models(m2_SBP_aft,m2_SBP_aft.noInfl)

# showing regression table original vs. updated model
tab_model(m2_SBP_aft,m2_SBP_aft.noInfl,dv.labels=c("Original","Updated"),show.icc=FALSE,show.p=FALSE,show.se=TRUE,
          show.r2=FALSE,collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Original			Updated
Predictors	b (SE)	CI	t	b (SE)	CI	t
(Intercept)	118.25 (1.97)	114.38 – 122.12	59.96	116.39 (1.76)	112.94 – 119.84	66.17
gender [M]	5.39 (2.92)	-0.35 – 11.13	1.84	7.31 (2.59)	2.23 – 12.39	2.83
age gmc	0.41 (0.12)	0.18 – 0.65	3.43	0.42 (0.11)	0.21 – 0.63	3.91
BMI gmc	1.10 (0.43)	0.25 – 1.95	2.54	1.05 (0.38)	0.30 – 1.80	2.74
WHLSM cm gmc	1.33 (1.14)	-0.91 – 3.57	1.16	0.87 (1.01)	-1.11 – 2.85	0.87
WHLSM mc	1.56 (0.40)	0.77 – 2.34	3.90	1.19 (0.39)	0.43 – 1.94	3.08
Random Effects
σ2	106.74			95.98
τ00	195.55 ID			148.12 ID
N	106 ID			104 ID
Observations	787			771

DBP_aft

Here, we inspect the magnitude of the changes in the estimated coefficients after the removal of potentially influential cases. We can see that the removal of the influential participant mainly affects the parameter estimated for gender, with men showing substantially higher DBP_aft than women in the updated model. Thus, since gender is not our main focus, we rely on the results obtained with the full sample.

# refitting model without influential cases
m2_DBP_aft.noInfl <- update(m2_SBP_aft,data=cleanBP_aft[cleanBP_aft$ID!="S082",])

# plotting coefficients original vs. updated model
plot_models(m2_DBP_aft,m2_DBP_aft.noInfl)

# showing regression table original vs. updated model
tab_model(m2_DBP_aft,m2_DBP_aft.noInfl,dv.labels=c("Original","Updated"),show.icc=FALSE,show.p=FALSE,show.se=TRUE,
          show.r2=FALSE,collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Original			Updated
Predictors	b (SE)	CI	t	b (SE)	CI	t
(Intercept)	77.32 (1.46)	74.45 – 80.19	52.81	116.94 (1.82)	113.37 – 120.50	64.42
gender [M]	-0.38 (2.17)	-4.65 – 3.88	-0.18	6.91 (2.68)	1.65 – 12.17	2.58
age gmc	0.28 (0.09)	0.10 – 0.45	3.07	0.47 (0.11)	0.26 – 0.69	4.27
BMI gmc	0.96 (0.32)	0.33 – 1.59	2.98	0.95 (0.39)	0.18 – 1.72	2.41
WHLSM cm gmc	1.11 (0.85)	-0.55 – 2.77	1.31	0.74 (1.05)	-1.32 – 2.79	0.70
WHLSM mc	1.30 (0.35)	0.61 – 1.98	3.74	1.42 (0.40)	0.63 – 2.21	3.55
Random Effects
σ2	80.61			104.74
τ00	104.66 ID			159.34 ID
N	106 ID			105 ID
Observations	787			778

3.1.2. Evening BP

Systolic SBP_eve and diastolic blood pressure DBP_eve recoded in the evening are predicted by concurrent state workaholism WHLSM.mc, by psychological detachment PD, trait workaholism WHLSM.cm, and three covariates selected from the previous step (see Supplementary Material S5), namely gender, age, and BMI. A further potential confounder highlighted from the previous step, namely day (i.e., expressing the linear time trend over the study protocol), is not included at this point to get more parsimonious models and more comparable results across BP models, also considering the lack of linear temporal trends in state WHLSM measures (see below). Particularly, as shown in Supplementary Material S5, such linear trend is possibly due to measurement reactivity (i.e., increased physiological activation during the first - less familiar - blood pressure recordings), and it will be considered as a robustness check in section 4.2.

mformulas[3] # covariates selected for SBP_aft

## [1] "SBP_eve ~ gender + age + BMI"

mformulas[4] # covariates selected for DBP_aft

## [1] "DBP_eve ~ gender + age + BMI"

# plotting WHLSM & SBP by day
par(mfrow=c(1,3)); boxplot(SBP_eve ~ day,data=clean); boxplot(DBP_eve ~ day,data=clean); boxplot(WHLSM ~ day,data=clean)

# model weight and estimated parameter for WHLSM by time
m0 <- lmer(WHLSM ~ (1|ID),data=cleanBP)
m1 <- lmer(WHLSM ~ day + (1|ID),data=cleanBP)
Weights(AIC(m0,m1)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.965 0.035

summary(m1)$coefficients # coefficient: |t| < 2

##               Estimate Std. Error   t value
## (Intercept) 3.26380181  0.1377577 23.692339
## day         0.01819296  0.0116718  1.558711

3.1.2.1. Data preparation

First, we prepare the data for the analyses by removing all cases of missing responses in the dependent variable or any predictor or covariate (list-wise deletion), by centering level-2 continuous predictors on the grand mean (grand-mean-centering), and by centering level-1 continuous predictors on the individual mean (person-mean-centering).

# list-wise deletion
cleanBP_eve <- as.data.frame(na.omit(cleanBP[,c("ID","SBP_eve","DBP_eve", # grouping and dependent variables
                                                "gender","age","BMI","WHLSM","PD", # core predictors
                                                "WE","WC","sleep_dysf","psy_drugs","horm_drugs", # for robustness checks
                                                "confounders_eve","flagBP_eve","flagTime","careless","day",
                                                "position","children")])) 
cleanBP_eve$ID <- as.factor(as.character(cleanBP_eve$ID)) # resetting participant identifier levels
cat("Considering",nrow(cleanBP_eve),"complete obs from",nlevels(as.factor(as.character(cleanBP_eve$ID))),"participants")

## Considering 721 complete obs from 106 participants

# person-mean-centering lv-1 continuous predictors
wide <- cleanBP_eve[!duplicated(cleanBP_eve$ID),] # wide-form dataset
for(Var in c("WHLSM","PD","WE","WC")){
  wide <- cbind(wide,aggregate(cleanBP_eve[,Var],list(cleanBP_eve$ID),mean)[,2]) # individual means
  colnames(wide)[ncol(wide)] <- paste0(Var,".cm")
  cleanBP_eve <- join(cleanBP_eve,wide[,c("ID",paste0(Var,".cm"))],by="ID",type="left") # joining to long-form df
  cleanBP_eve[,paste0(Var,".mc")] <- cleanBP_eve[,Var] - cleanBP_eve[,paste0(Var,".cm")] } # mean-centered scores

# grand-mean-centering lv-2 continuous predictors
for(Var in c("age","BMI","WHLSM.cm")){ cleanBP_eve[,paste0(Var,".gmc")] <- cleanBP_eve[,Var] - mean(wide[,Var]) }

# showing data
cleanBP_eve[1:3,] # first three rows

3.1.2.2. Model fit

Here, we fit the multilevel models to the selected data using the default restricted maximum likelihood estimator (REML).

# m0: null model
m0_SBP_eve <- lmer(SBP_eve ~ (1|ID), # only fixed and random intercept + residual term
                   data=cleanBP_eve)

# m1: covariates
m1_SBP_eve <- lmer(SBP_eve ~ gender + age.gmc + BMI.gmc + PD.mc + WHLSM.cm.gmc + (1|ID), # covariates
                       data=cleanBP_eve)

# m2: state workaholism
m2_SBP_eve <- lmer(SBP_eve ~ gender + age.gmc + BMI.gmc + PD.mc + WHLSM.cm.gmc + WHLSM.mc + (1|ID),
                   data=cleanBP_eve)

# m3: interaction
m3_SBP_eve <- lmer(SBP_eve ~ gender + age.gmc + BMI.gmc + PD.mc + WHLSM.cm.gmc + WHLSM.mc + WHLSM.mc:PD.mc + (1|ID),
                   data=cleanBP_eve)

The same models are specified for diastolic blood pressure.

# m0: null model
m0_DBP_eve <- lmer(DBP_eve ~ (1|ID), # only fixed and random intercept + residual term
                   data=cleanBP_eve)

# m1: covariates
m1_DBP_eve <- lmer(DBP_eve ~ gender + age.gmc + BMI.gmc + PD.mc + WHLSM.cm.gmc + (1|ID),
                   data=cleanBP_eve)

# m2: state workaholism
m2_DBP_eve <- lmer(DBP_eve ~ gender + age.gmc + BMI.gmc + PD.mc + WHLSM.cm.gmc + WHLSM.mc + (1|ID),
                   data=cleanBP_eve)

# m3: interactions
m3_DBP_eve <- lmer(DBP_eve ~ gender + age.gmc + BMI.gmc + PD.mc + WHLSM.cm.gmc + WHLSM.mc + WHLSM.mc:PD.mc + (1|ID),
                   data=cleanBP_eve)

From the previous chunks, we see that all models converged without problems. Here, we inspect the diagnostics (i.e., normality of residual and random effect distributions, homoscedasticity, and multicollinearity) of the most complex model m3. Influential cases are analyzed in a dedicated section below.

SBP_eve

Model m3_SBP_eve shows some deviation from normality especially in the upper tail of the distribution of residuals and both tails of the distribution of random effects. Particularly, participant S082 is associated with the highest extreme deviation from the distributions of random intercepts, and will be removed as a robustness check (see section 4.1). Besides that, we can see that the homoscedsticity assumption holds and that none of the variance inflation factors (VIFs) shows extreme values, ruling out the risk of multicollinearity.

# normality and homoscedasticity
p <- plot_model(m3_SBP_eve,type="diag",dot.size=1)
p[[2]] <- p[[2]]$ID
plot_grid(p,tags=TRUE,margin=c(0,0,0,0))

# participant with highest random effects (i.e., BLUPS)
re <- ranef(m3_SBP_eve)$ID
as.character(cleanBP_prelqs[which(re$`(Intercept)`==max(re$`(Intercept)`)),"ID"])

## [1] "S082"

# homoscedasticity and multicollinearity
par(mfrow=c(1,2))
for(Var in c("gender")){ boxplot(resid(m3_SBP_eve) ~ cleanBP_eve[,Var],main=paste("Residuals by",Var)) }
barplot(vif(m3_SBP_eve),main="VIF Values",xlim=c(0,10),las=2,horiz=TRUE) # variance inflation factors (VIFs)
abline(v = 5, lwd = 5, lty = 2)

Here, we better inspect the residual distribution and the fit of models specified with alternative family distributions. We can see that none of the alternative families substantially better approximate the distribution of model residuals, with the log-transformed solution (i.e., normal distribution with log-transformed dependent variable) showing the best fit. Yet, it is not so better than the original model. Thus, we initially rely on the normal distribution and then consider the logarithmic transformation as a robustness check (see section 4.1).

# inspecting residual distribution
descdist(resid(m3_SBP_eve)) # unknown best-fit distribution

## summary statistics
## ------
## min:  -37.79119   max:  46.31558 
## median:  -0.2969038 
## mean:  6.652117e-14 
## estimated sd:  9.132578 
## estimated skewness:  0.149237 
## estimated kurtosis:  4.770943

# fitting model with alternative families
models <- list(
  m3_SBP_eve,
  glmer(formula=formula(m3_SBP_eve),family=Gamma(link="log"),data=cleanBP_eve), # gamma log (doesn't converge)
  glmer(formula=formula(m3_SBP_eve),family=Gamma(link="identity"),data=cleanBP_eve), # gamma id (doesn't converge)
  glmer(formula=formula(m3_SBP_eve),family=gaussian(link="log"),data=cleanBP_eve), # log-normal (singular fit)
  lmer(formula=as.formula(paste("log(SBP_eve) ~",as.character(formula(m3_SBP_eve))[3])),data=cleanBP_eve)) # log transf

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00308287 (tol = 0.002, component 1)

## boundary (singular) fit: see help('isSingular')

# normal Q-Q plot of model residuals
par(mfrow=c(1,5))
for(i in 1:length(models)){ 
  qqnorm(resid(models[[i]]),main=c("Norm","Gamma-log","Gamma-id","log-norm","log-transf")[i]); qqline(resid(models[[i]]))}

DBP_eve

Model m3_DBP_eve shows some deviation from normality in both tails of the distribution of residuals and random effects. Particularly, participant S082 is associated with the highest extreme deviation from the distributions of both random effects, and will be removed as a robustness check (see section 4.1). Besides that, we can see that the homoscedsticity assumption holds and that none of the variance inflation factors (VIFs) shows extreme values, ruling out the risk of multicollinearity.

# normality and homoscedasticity
p <- plot_model(m3_DBP_eve,type="diag",dot.size=1)
p[[2]] <- p[[2]]$ID
plot_grid(p,tags=TRUE,margin=c(0,0,0,0))

# participant with highest random effects (i.e., BLUPS)
re <- ranef(m3_DBP_eve)$ID
re[re$WHLSM.mc==max(re$WHLSM.mc)|re$`(Intercept)`==max(re$`(Intercept)`),]

## numeric(0)

# homoscedasticity and multicollinearity
par(mfrow=c(1,2))
for(Var in c("gender")){ boxplot(resid(m3_DBP_eve) ~ cleanBP_eve[,Var],main=paste("Residuals by",Var)) }
barplot(vif(m3_DBP_eve),main="VIF Values",xlim=c(0,10),las=2,horiz=TRUE) # variance inflation factors (VIFs)
abline(v = 5, lwd = 5, lty = 2)

Here, we better inspect the residual distribution and the fit of models specified with alternative family distributions. We can see that none of the alternative families better approximates the distribution of model residuals, and none of them reached convergence. Thus, we initially rely on the normal distribution and then consider the gamma log as a robustness check (i.e., the only solution that converges).

# inspecting residual distribution
descdist(resid(m3_DBP_eve)) # unknown best-fit distribution but close to normal

## summary statistics
## ------
## min:  -22.70821   max:  31.60269 
## median:  -0.1333256 
## mean:  4.346043e-14 
## estimated sd:  7.838394 
## estimated skewness:  0.1657032 
## estimated kurtosis:  3.719358

# fitting model with alternative families
models <- list(
  m3_DBP_eve,
  glmer(formula=formula(m3_DBP_eve),family=Gamma(link="log"),data=cleanBP_eve), # gamma log
  glmer(formula=formula(m3_DBP_eve),family=Gamma(link="identity"),data=cleanBP_eve), # gamma id (doesn't converge)
  glmer(formula=formula(m3_DBP_eve),family=gaussian(link="log"),data=cleanBP_eve), # log-normal (singular fit)
  lmer(formula=as.formula(paste("log(SBP_eve) ~",as.character(formula(m2_DBP_eve))[3])),data=cleanBP_eve)) # log tr (singular)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0035261 (tol = 0.002, component 1)

## boundary (singular) fit: see help('isSingular')

# normal Q-Q plot of model residuals
par(mfrow=c(1,5))
for(i in 1:length(models)){ 
  qqnorm(resid(models[[i]]),main=c("Norm","Gamma-log","Gamma-id","log-norm","log-transf")[i]); qqline(resid(models[[i]]))}

3.1.2.3. Results

Here, we compare the specified models based on the Akaike weight and the likelihood ratio test (with type-I error set to p < .05), and we inspect the results of the selected model(s).

SBP_eve

We can see that the inclusion of state WHLSM is not associated with stronger evidence or significant likelihood ratio compared to models including less predictors (Aw = .33, \(\chi^2\)(1) = 0.58, p = .44). Similar, the interaction does not imply stronger evidence or significant likelihood ratio (Aw = .11, \(\chi^2\)(1) = 0.65, p = .42).

# Akaike weight adding one model at time
Weights(AIC(m0_SBP_eve,m1_SBP_eve)) # covariates: better

##  model weights 
## [1] 0 1

Weights(AIC(m0_SBP_eve,m1_SBP_eve,m2_SBP_eve)) # state workaholism: worse

##  model weights 
## [1] 0.000 0.668 0.332

Weights(AIC(m0_SBP_eve,m1_SBP_eve,m2_SBP_eve,m3_SBP_eve)) # interaction: worse

##  model weights 
## [1] 0.000 0.593 0.294 0.113

# Likelihood ratio test with m3.RDet
anova(m1_SBP_eve,m2_SBP_eve,m3_SBP_eve) # best model is m1

## refitting model(s) with ML (instead of REML)

Here, we inspect the coefficients estimated by the selected model m1.bis and those estimated by the target models. We can see that state WHLSM.mc is not substantially related to BP_eve, similar to trait WHLSM.cm.gmc is not. Among the included covariates, age, BMI, and low PD predict higher SBP_eve.

# regression table
tab_model(m1_SBP_eve,m2_SBP_eve,m3_SBP_eve,
          dv.labels=c("Covariates","State WHLSM","Interaction"),
          show.icc=FALSE,show.p=FALSE,show.se=TRUE,show.r2=FALSE,show.ci=FALSE,
          collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Covariates		State WHLSM		Interaction
Predictors	b (SE)	t	b (SE)	t	b (SE)	t
(Intercept)	113.08 (1.87)	60.61	113.08 (1.87)	60.61	113.11 (1.86)	60.66
gender [M]	2.48 (2.76)	0.90	2.48 (2.76)	0.90	2.48 (2.76)	0.90
age gmc	0.33 (0.11)	2.91	0.33 (0.11)	2.91	0.33 (0.11)	2.89
BMI gmc	1.40 (0.41)	3.44	1.40 (0.41)	3.44	1.40 (0.41)	3.45
PD mc	-0.82 (0.26)	-3.11	-0.81 (0.27)	-3.02	-0.80 (0.27)	-2.99
WHLSM cm gmc	1.34 (1.06)	1.26	1.34 (1.06)	1.26	1.35 (1.06)	1.27
WHLSM mc			0.31 (0.40)	0.76	0.31 (0.40)	0.76
PD mc × WHLSM mc					0.24 (0.30)	0.80
Random Effects
σ2	96.62		96.69		96.77
τ00	173.03 ID		173.02 ID		172.72 ID
N	106 ID		106 ID		106 ID
Observations	721		721		721

# plotting main effects from model m2
grid.arrange(plot_model(m2_SBP_eve,type="pred",terms="WHLSM.cm.gmc"), # trait WHLSM
             plot_model(m2_SBP_eve,type="pred",terms="WHLSM.mc"),nrow=1) # state WHLSM

# plotting interactions
sd(cleanBP_eve$PD.mc) # RDet: 1 SD = 1.38

## [1] 1.382712

plot_model(m3_SBP_eve,type="pred",terms=c("WHLSM.mc","PD.mc [-1.38,1.38]"))

DBP_eve

Results are similar to those found for SBP_eve: the inclusion of state WHLSM is associated with stronger evidence (Aw = .68) but not with significant likelihood ratio compared to models including less predictors (\(\chi^2\)(1) = 3.71, p = .05). Similar, the interaction does not imply stronger evidence or significant likelihood ratio (Aw = .14, \(\chi^2\)(1) = 0.01, p = .99).

# Akaike weight adding one model at time
Weights(AIC(m0_DBP_eve,m1_DBP_eve)) # covariates: better

##  model weights 
## [1] 0 1

Weights(AIC(m0_DBP_eve,m1_DBP_eve,m2_DBP_eve)) # state workaholism: better

##  model weights 
## [1] 0.000 0.319 0.680

Weights(AIC(m0_DBP_eve,m1_DBP_eve,m2_DBP_eve,m3_DBP_eve)) # state whlsm by detachment: worse

##  model weights 
## [1] 0.000 0.275 0.585 0.140

# Likelihood ratio test with m3
anova(m1_DBP_eve,m2_DBP_eve,m3_DBP_eve) # best model is m1

## refitting model(s) with ML (instead of REML)

Here, we inspect the coefficients estimated by the models specified for DBP_eve. We can see that state WHLSM.mc is not substantially related to DBP_eve. Among the included covariates, age and BMI predict higher SBP_eve.

# regression table
tab_model(m1_DBP_eve,m2_DBP_eve,m3_DBP_eve,
          dv.labels=c("Covariates","State WHLSM","Interaction"),
          show.icc=FALSE,show.p=FALSE,show.se=TRUE,show.r2=FALSE,show.ci=FALSE,
          collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Covariates		State WHLSM		Interaction
Predictors	b (SE)	t	b (SE)	t	b (SE)	t
(Intercept)	72.72 (1.45)	50.17	72.72 (1.45)	50.17	72.72 (1.45)	50.15
gender [M]	-1.31 (2.15)	-0.61	-1.31 (2.15)	-0.61	-1.31 (2.15)	-0.61
age gmc	0.27 (0.09)	3.05	0.27 (0.09)	3.05	0.27 (0.09)	3.05
BMI gmc	1.16 (0.32)	3.66	1.16 (0.32)	3.66	1.16 (0.32)	3.66
PD mc	-0.23 (0.23)	-0.99	-0.18 (0.23)	-0.80	-0.18 (0.23)	-0.80
WHLSM cm gmc	0.46 (0.83)	0.55	0.46 (0.83)	0.55	0.46 (0.83)	0.55
WHLSM mc			0.67 (0.34)	1.95	0.67 (0.35)	1.95
PD mc × WHLSM mc					0.00 (0.26)	0.01
Random Effects
σ2	71.29		70.97		71.08
τ00	102.38 ID		102.42 ID		102.43 ID
N	106 ID		106 ID		106 ID
Observations	721		721		721

# plotting main effects from model m2
grid.arrange(plot_model(m2_SBP_eve,type="pred",terms="WHLSM.cm.gmc"), # trait WHLSM
             plot_model(m2_SBP_eve,type="pred",terms="WHLSM.mc"),nrow=1) # state WHLSM

# plotting interactions
sd(cleanBP_eve$PD.mc) # RDet: 1 SD = 1.38

## [1] 1.382712

plot_model(m3_DBP_eve,type="pred",terms=c("WHLSM.mc","PD.mc [-1.38,1.38]"))

3.1.2.4. Influential cases

Here, we evaluate the presence of influential cases in the target model m2. Specifically, the Cook’s distance is considered as the main measure of individual-level (i.e., participant) influence on the estimated parameters, and it is recomputed by progressively excluding the most influential participants (i.e., based on the rule-of-thumb of 4/N) until all extreme values are removed.

3.1.2.4.1. Cook’s distance

SBP_eve

We can see that participants S082 and S096 are potentially influential cases. Note that the latter is the same participants showing the most extreme deviations from the distributions of random effects (see section 3.3.2).

# cook's distance on the whole sample
infl <- influence(m2_SBP_eve,"ID")
plot(infl,which="cook",cutoff=4/nlevels(cleanBP_eve$ID),xlab="Cook distance",ylab="ID",sort=TRUE)

# progressively excluding participants
infl <- list(
  influence(exclude.influence(m2_SBP_eve,"ID","S082"),"ID"),
  influence(exclude.influence(m2_SBP_eve,"ID",c("S082","S096")),"ID"))
for(i in 1:length(infl)){ 
  plot(infl[[i]],which="cook",cutoff=4/(nlevels(cleanBP_eve$ID)-i),xlab="Cook distance",ylab="ID",sort=TRUE) }

DBP_eve

We can see that participants S082 and S080 is a potentially influential case. Note that the former is the same participant showing the most extreme deviations from the distributions of random effects (see section 3.3.2).

# cook's distance on the whole sample
infl <- influence(m2_DBP_eve,"ID")
plot(infl,which="cook",cutoff=4/nlevels(cleanBP_eve$ID),xlab="Cook distance",ylab="ID",sort=TRUE)

# progressively excluding participants
infl <- list(influence(exclude.influence(m2_DBP_eve,"ID","S082"),"ID"),
             influence(exclude.influence(m2_DBP_eve,"ID",c("S082","S080")),"ID"))
for(i in 1:length(infl)){ 
  plot(infl[[i]],which="cook",cutoff=4/(nlevels(cleanBP_eve$ID)-i),xlab="Cook distance",ylab="ID",sort=TRUE) }

3.1.2.4.2. Coefficient change

SBP_eve

Here, we inspect the magnitude of the changes in the estimated coefficients after the removal of potentially influential cases. We can see that the coefficients estimated by the updated model do not substantially differ from those estimated by the original model. Thus, we choose to rely on the results obtained with the full sample.

# refitting model without influential cases
m2_SBP_eve.noInfl <- update(m2_SBP_eve,data=cleanBP_eve[!cleanBP_eve$ID%in%c("S082","S096"),])

# plotting coefficients original vs. updated model
plot_models(m2_SBP_eve,m2_SBP_eve.noInfl)

# same thing with interactions
m3_DBP_eve.noInfl <- update(m3_SBP_eve,data=cleanBP_eve[!cleanBP_eve$ID%in%c("S082","S096"),])
plot_models(m3_SBP_eve,m3_DBP_eve.noInfl)

# showing regression table original vs. updated model
tab_model(m2_SBP_eve,m2_SBP_eve.noInfl,dv.labels=c("Original","Updated"),show.icc=FALSE,show.p=FALSE,show.se=TRUE,
          show.r2=FALSE,collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Original			Updated
Predictors	b (SE)	CI	t	b (SE)	CI	t
(Intercept)	113.08 (1.87)	109.41 – 116.74	60.61	111.19 (1.63)	108.00 – 114.39	68.33
gender [M]	2.48 (2.76)	-2.94 – 7.91	0.90	4.43 (2.39)	-0.26 – 9.12	1.86
age gmc	0.33 (0.11)	0.11 – 0.56	2.91	0.34 (0.10)	0.15 – 0.54	3.45
BMI gmc	1.40 (0.41)	0.60 – 2.20	3.44	1.34 (0.35)	0.65 – 2.03	3.81
PD mc	-0.81 (0.27)	-1.33 – -0.28	-3.02	-0.88 (0.27)	-1.41 – -0.35	-3.26
WHLSM cm gmc	1.34 (1.06)	-0.75 – 3.43	1.26	0.87 (0.92)	-0.94 – 2.67	0.94
WHLSM mc	0.31 (0.40)	-0.48 – 1.10	0.76	0.16 (0.40)	-0.64 – 0.95	0.39
Random Effects
σ2	96.69			94.08
τ00	173.02 ID			123.34 ID
N	106 ID			104 ID
Observations	721			709

DBP_eve

Here, we inspect the magnitude of the changes in the estimated coefficients after the removal of potentially influential cases. We can see that the coefficient estimated for state WHLSM.mc by the updated model is strongly reduced. Thus, these results question the generalizability of the estimated relationship between WHLSM.mc and DBP_eve.

# refitting model without influential cases
m2_DBP_eve.noInfl <- update(m2_SBP_eve,data=cleanBP_eve[!cleanBP_eve$ID%in%c("S082","S080"),])

# plotting coefficients original vs. updated model
plot_models(m2_DBP_eve,m2_DBP_eve.noInfl)

# same thing with interaction
m3_DBP_eve.noInfl <- update(m3_SBP_eve,data=cleanBP_eve[!cleanBP_eve$ID%in%c("S082","S080"),])
plot_models(m3_DBP_eve,m3_DBP_eve.noInfl)

# showing regression table original vs. updated model
tab_model(m2_DBP_eve,m2_DBP_eve.noInfl,dv.labels=c("Original","Updated"),show.icc=FALSE,show.p=FALSE,show.se=TRUE,
          show.r2=FALSE,collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Original			Updated
Predictors	b (SE)	CI	t	b (SE)	CI	t
(Intercept)	72.72 (1.45)	69.87 – 75.57	50.17	111.67 (1.70)	108.33 – 115.01	65.63
gender [M]	-1.31 (2.15)	-5.53 – 2.90	-0.61	4.09 (2.50)	-0.81 – 8.99	1.64
age gmc	0.27 (0.09)	0.10 – 0.45	3.05	0.39 (0.10)	0.19 – 0.59	3.78
BMI gmc	1.16 (0.32)	0.54 – 1.78	3.66	1.25 (0.37)	0.53 – 1.97	3.41
PD mc	-0.18 (0.23)	-0.63 – 0.26	-0.80	-0.73 (0.26)	-1.25 – -0.22	-2.79
WHLSM cm gmc	0.46 (0.83)	-1.17 – 2.08	0.55	0.77 (0.96)	-1.11 – 2.66	0.80
WHLSM mc	0.67 (0.34)	-0.00 – 1.35	1.95	0.16 (0.40)	-0.62 – 0.94	0.40
Random Effects
σ2	70.97			89.95
τ00	102.42 ID			136.26 ID
N	106 ID			104 ID
Observations	721			706

3.1.3. Mediation

In light of the relationships found between state workaholism and afternoon blood pressure, and the substantial correlation found between afternoon and evening blood pressure, here we explore the mediating paths between workaholism and later blood pressure measurements by evaluating the potential mediating role of afternoon blood pressure. This is done by considering the same covariates included above, namely gender, age, BMI, PD, and trait WHLSM.cm.

library(mediation) # loding mediation package

# data preparation (list-wise deletion and mean centering)
cleanBP_med_eve <- as.data.frame(na.omit(cleanBP[,c("ID","SBP_aft","DBP_aft","SBP_eve","DBP_eve", # grouping and response vars
                                                "gender","age","BMI","WHLSM","PD", # core predictors
                                                "WE","WC","sleep_dysf","psy_drugs","horm_drugs", # for robustness checks
                                                "confounders_aft","confounders_eve",
                                                "flagBP_aft","flagBP_eve","flagTime","careless","day",
                                                "children","position")])) 
cleanBP_med_eve$ID <- as.factor(as.character(cleanBP_med_eve$ID)) # resetting participant identifier levels
wide <- cleanBP_med_eve[!duplicated(cleanBP_med_eve$ID),] # wide-form dataset
for(Var in c("WHLSM","PD","SBP_aft","DBP_aft","WE","WC")){
  wide <- cbind(wide,aggregate(cleanBP_med_eve[,Var],list(cleanBP_med_eve$ID),mean)[,2]) # individual means
  colnames(wide)[ncol(wide)] <- paste0(Var,".cm")
  cleanBP_med_eve <- join(cleanBP_med_eve,wide[,c("ID",paste0(Var,".cm"))],by="ID",type="left") # joining to long-form df
  cleanBP_med_eve[,paste0(Var,".mc")] <- cleanBP_med_eve[,Var] - cleanBP_med_eve[,paste0(Var,".cm")] } # mean-centered scores
for(Var in c("age","BMI","WHLSM.cm")){ cleanBP_med_eve[,paste0(Var,".gmc")] <- cleanBP_med_eve[,Var] - mean(wide[,Var]) } # gmc
cat("Considering",nrow(cleanBP_med_eve),"complete obs from",nlevels(as.factor(as.character(cleanBP_med_eve$ID))),"participants")

## Considering 717 complete obs from 106 participants

We can see that a significant indirect effect is estimated in both cases, whereas the direct effect is not significant for both systolic and diastolic blood pressure. Results are consistent after the removal of influential cases.

SBP_eve

# output and mediation models
mOut_SBP_eve <- lmer(SBP_eve ~ gender + age.gmc + BMI.gmc + PD.mc + WHLSM.cm.gmc + WHLSM.mc + SBP_aft.mc + 
                       (1|ID),data=cleanBP_med_eve)
mMed_SBP_eve <- lmer(SBP_aft ~ gender + age.gmc + BMI.gmc + PD.mc + WHLSM.cm.gmc + WHLSM.mc + 
                       (1|ID),data=cleanBP_med_eve)

# mediation results
NSIM <- 10000
med_SBP_eve <- mediation::mediate(model.m=mMed_SBP_eve,model.y=mOut_SBP_eve,treat="WHLSM.mc",mediator="SBP_aft.mc",
                       covariates=cleanBP_med_eve[,c("gender","age.gmc","BMI.gmc","PD.mc","WHLSM.cm.gmc")],
                       boot=FALSE,sims=NSIM) # quasi-Bayesian CI
summary(med_SBP_eve) # significant indirect effect

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
## Mediator Groups: ID 
## 
## Outcome Groups: ID 
## 
## Output Based on Overall Averages Across Groups 
## 
##                Estimate 95% CI Lower 95% CI Upper p-value   
## ACME             0.1967       0.0565         0.38  0.0032 **
## ADE              0.1126      -0.6597         0.90  0.7800   
## Total Effect     0.3093      -0.4800         1.10  0.4548   
## Prop. Mediated   0.3392      -5.5738         5.61  0.4548   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 717 
## 
## 
## Simulations: 10000

# mediation results without influential cases
mOut_SBP_eve.noInfl <- update(mOut_SBP_eve,data=cleanBP_med_eve[!cleanBP_med_eve$ID%in%c("S082","S096"),])
mMed_SBP_eve.noInfl <- update(mMed_SBP_eve,data=cleanBP_med_eve[!cleanBP_med_eve$ID%in%c("S082","S096"),])
med_SBP_eve.noInfl <- mediation::mediate(model.m=mMed_SBP_eve.noInfl,model.y=mOut_SBP_eve.noInfl,
                              treat="WHLSM.mc",mediator="SBP_aft.mc",
                              covariates=cleanBP_med_eve[!cleanBP_med_eve$ID%in%c("S082","S096"),
                                                         c("gender","age.gmc","BMI.gmc","PD.mc","WHLSM.cm.gmc")],
                              boot=FALSE,sims=NSIM) # quasi-Bayesian CI
summary(med_SBP_eve.noInfl) # indirect effect stays significant

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
## Mediator Groups: ID 
## 
## Outcome Groups: ID 
## 
## Output Based on Overall Averages Across Groups 
## 
##                Estimate 95% CI Lower 95% CI Upper p-value  
## ACME             0.1379       0.0199         0.30   0.018 *
## ADE              0.0188      -0.7724         0.79   0.961  
## Total Effect     0.1567      -0.6359         0.94   0.700  
## Prop. Mediated   0.2099      -5.3281         4.61   0.703  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 705 
## 
## 
## Simulations: 10000

DBP_eve

# output and mediation models
mOut_DBP_eve <- lmer(DBP_eve ~ gender + age.gmc + BMI.gmc + PD.mc + WHLSM.cm.gmc + WHLSM.mc + DBP_aft.mc + 
                       (1|ID),data=cleanBP_med_eve)
mMed_DBP_eve <- lmer(DBP_aft ~ gender + age.gmc + BMI.gmc + PD.mc + WHLSM.cm.gmc + WHLSM.mc + 
                       (1|ID),data=cleanBP_med_eve)

# mediation results
med_DBP_eve <- mediation::mediate(model.m=mMed_DBP_eve,model.y=mOut_DBP_eve,treat="WHLSM.mc",mediator="DBP_aft.mc",
                       covariates=cleanBP_med_eve[,c("gender","age.gmc","BMI.gmc","PD.mc","WHLSM.cm.gmc")],
                       boot=FALSE,sims=NSIM) # quasi-Bayesian CI
summary(med_DBP_eve) # significant indirect effect

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
## Mediator Groups: ID 
## 
## Outcome Groups: ID 
## 
## Output Based on Overall Averages Across Groups 
## 
##                Estimate 95% CI Lower 95% CI Upper p-value   
## ACME             0.1700       0.0453         0.33  0.0044 **
## ADE              0.5198      -0.1547         1.21  0.1304   
## Total Effect     0.6898       0.0265         1.37  0.0428 * 
## Prop. Mediated   0.2376       0.0164         1.48  0.0468 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 717 
## 
## 
## Simulations: 10000

# mediation results without influential cases
mOut_DBP_eve.noInfl <- update(mOut_DBP_eve,data=cleanBP_med_eve[!cleanBP_med_eve$ID%in%c("S082","S080"),])
mMed_DBP_eve.noInfl <- update(mMed_DBP_eve,data=cleanBP_med_eve[!cleanBP_med_eve$ID%in%c("S082","S080"),])
med_DBP_eve.noInfl <- mediation::mediate(model.m=mMed_DBP_eve.noInfl,model.y=mOut_DBP_eve.noInfl,
                              treat="WHLSM.mc",mediator="DBP_aft.mc",
                              covariates=cleanBP_med_eve[!cleanBP_med_eve$ID%in%c("S082","S080"),
                                                         c("gender","age.gmc","BMI.gmc","PD.mc","WHLSM.cm.gmc")],
                              boot=FALSE,sims=NSIM) # quasi-Bayesian CI
summary(med_DBP_eve.noInfl) # indirect effect stays significant

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
## Mediator Groups: ID 
## 
## Outcome Groups: ID 
## 
## Output Based on Overall Averages Across Groups 
## 
##                Estimate 95% CI Lower 95% CI Upper p-value   
## ACME             0.1847       0.0608         0.35  0.0016 **
## ADE              0.5452      -0.1353         1.23  0.1146   
## Total Effect     0.7299       0.0475         1.41  0.0372 * 
## Prop. Mediated   0.2444       0.0415         1.32  0.0388 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 702 
## 
## 
## Simulations: 10000

3.2. Emotional Exhaustion

Emotional exhaustion EE rated in the evening is predicted by state workaholism WHLSM.mc, in addition to psychological detachment PD, trait workaholism WHLSM.cm, and one covariate selected from the previous step (see Supplementary Material S5), namely gender.

mformulas[5] # covariates selected for EE

## [1] "EE ~ gender"

3.2.1. Data preparation

First, we prepare the data for the analyses by removing all cases of missing responses in the dependent variable or any predictor or covariate (list-wise deletion), by centering level-2 continuous predictors on the grand mean (grand-mean-centering), and by centering level-1 continuous predictors on the individual mean (person-mean-centering).

# list-wise deletion
cleanEE <- as.data.frame(na.omit(clean[,c("ID","EE","gender","PD","WHLSM", # core variables
                                          "WE","WC","flagTime","careless","day", # for robustness checks
                                          "position","children")])) 
cleanEE$ID <- as.factor(as.character(cleanEE$ID)) # resetting participant identifier levels
cat("Considering",nrow(cleanEE),"complete observations from",nlevels(as.factor(as.character(cleanEE$ID))),"participants")

## Considering 779 complete observations from 114 participants

# person-mean-centering lv-1 continuous predictors
wide <- cleanEE[!duplicated(cleanEE$ID),] # wide-form dataset
for(Var in c("PD","WHLSM","WE","WC")){
  wide <- cbind(wide,aggregate(cleanEE[,Var],list(cleanEE$ID),mean)[,2]) # computing individual means
  colnames(wide)[ncol(wide)] <- paste0(Var,".cm")
  cleanEE <- join(cleanEE,wide[,c("ID",paste0(Var,".cm"))],by="ID",type="left") # joining with long-form data
  cleanEE[,paste0(Var,".mc")] <- cleanEE[,Var] - cleanEE[,paste0(Var,".cm")] } # computing mean-centered scores

# grand-mean-centering lv-2 continuous predictors
for(Var in c("WHLSM.cm")){ cleanEE[,paste0(Var,".gmc")] <- cleanEE[,Var] - mean(wide[,Var]) }

# showing data
cleanEE[1:3,] # first three rows

3.2.2. Model fit

Here, we fit the multilevel models to the selected data using the default restricted maximum likelihood estimator (REML).

# m0: null model
m0_EE <- lmer(EE ~ (1|ID), # only fixed and random intercept + residual term
              data=cleanEE)

# m1: covariates
m1_EE <- lmer(EE ~ gender + PD.mc + WHLSM.cm.gmc + (1|ID), # covariates + recovery exp
              data=cleanEE)

# m2: state workaholism
m2_EE <- lmer(EE ~ gender + PD.mc + WHLSM.cm.gmc + WHLSM.mc + (1|ID),
              data=cleanEE)

# m3: interactions
m3_EE <- lmer(EE ~ gender + PD.mc + WHLSM.cm.gmc + WHLSM.mc + WHLSM.mc + WHLSM.mc:PD.mc + (1|ID), 
              data=cleanEE)

From the previous chunk, we see that all models converged without problems. Here, we inspect the diagnostics (i.e., normality of residual and random effect distributions, homoscedasticity, and multicollinearity) of the most complex model m3. Influential cases are analyzed in a dedicated section below.

Model m3 shows some deviation from normality in both tails of the distribution of residuals and both random effects. Besides that, we can see that the homoscedsticity assumption holds and that none of the variance inflation factors (VIFs) shows extreme values, ruling out the risk of multicollinearity.

# normality and homoscedasticity
p <- plot_model(m3_EE,type="diag",dot.size=1) 
p[[2]] <- p[[2]]$ID
plot_grid(p,tags=TRUE,margin=c(0,0,0,0))

# homoscedasticity and multicollinearity
par(mfrow=c(1,2)) 
for(Var in c("gender")){ boxplot(resid(m3_EE) ~ cleanEE[,Var],main=paste("Residuals by",Var)) }
barplot(vif(m3_EE),main="VIF Values",xlim=c(0,10),las=2,horiz=TRUE) # variance inflation factors (VIFs)
abline(v = 5, lwd = 5, lty = 2)

Here, we better inspect the residual distribution and the fit of models specified with alternative family distributions. We can see that alternative families might better approximate the distribution of model residuals, with the Gamma-log solution (i.e., Gamma distribution with logarithmic link function) showing the best fit. Yet, it is not so better than the original model. Thus, we initially rely on the normal distribution and then consider the Gamma-log solution as a robustness check.

# inspecting residual distribution
descdist(resid(m3_EE)) # best fit for lognormal?

## summary statistics
## ------
## min:  -3.710471   max:  3.628516 
## median:  -0.05379407 
## mean:  -1.253042e-15 
## estimated sd:  0.9577171 
## estimated skewness:  0.2331125 
## estimated kurtosis:  3.691947

# fitting model with alternative families
models <- list(
  m3_EE,
  glmer(formula=formula(m3_EE),family=Gamma(link="log"),data=cleanEE), # gamma with log
  glmer(formula=formula(m3_EE),family=Gamma(link="identity"),data=cleanEE), # [FAILS TO CONV, max|grad| = 0.09]
  glmer(formula=formula(m3_EE),family=gaussian(link="log"),data=cleanEE), # log-norm
  lmer(formula=as.formula(paste("log(EE) ~",as.character(formula(m3_EE))[3])),data=cleanEE)) # norm log

# normal Q-Q plot of model residuals
par(mfrow=c(1,5))
for(i in 1:length(models)){ 
  qqnorm(resid(models[[i]]),main=c("Norm","Gamma-log","Log-norm","log-transf")[i]); qqline(resid(models[[i]])) }

3.2.3. Results

Here, we compare the specified models based on the Akaike weight and the likelihood ratio test (with type-I error set to p < .05), and we inspect the results of the selected model(s). We can see that the inclusion of state WHLSM is associated with stronger evidence and significantly higher likelihood than the baseline model (Aw = .99, \(\chi^2\)(1) = 20.46, p < .001). In contrast, the interactive model is not better than the previous ones (Aw = .06, \(\chi^2\)(1) = 1.68, p = .19).

# Akaike weight adding one model at time
Weights(AIC(m0_EE,m1_EE)) # covariates: better

##  model weights 
## [1] 0 1

Weights(AIC(m0_EE,m1_EE,m2_EE)) # state workaholism: better

##  model weights 
## [1] 0.000 0.001 0.999

Weights(AIC(m0_EE,m1_EE,m2_EE,m3_EE)) # interaction: worse

##  model weights 
## [1] 0.000 0.001 0.939 0.061

# Likelihood ratio test with m3
anova(m1_EE,m2_EE,m3_EE) # best model is m2

## refitting model(s) with ML (instead of REML)

Here, we inspect the coefficients estimated by the selected model m2 and those estimated by other models. We can see that state WHLSM.mc is substantially and negatively related to EE, whereas the interaction is not substantial. Among the included covariates, low PD and high trait WHLSM.gmc predict substantially higher EE.

# regression table
tab_model(m1_EE,m2_EE,m3_EE,
          dv.labels=c("Baseline","State WHLSM","Interaction"),
          show.icc=FALSE,show.p=FALSE,show.se=TRUE,show.r2=FALSE,show.ci=FALSE,
          collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Baseline		State WHLSM		Interaction
Predictors	b (SE)	t	b (SE)	t	b (SE)	t
(Intercept)	3.12 (0.14)	21.89	3.12 (0.14)	21.89	3.12 (0.14)	21.84
gender [M]	0.18 (0.21)	0.88	0.18 (0.21)	0.88	0.18 (0.21)	0.88
PD mc	-0.08 (0.03)	-3.05	-0.07 (0.03)	-2.70	-0.07 (0.03)	-2.75
WHLSM cm gmc	0.47 (0.08)	5.91	0.47 (0.08)	5.91	0.47 (0.08)	5.88
WHLSM mc			0.19 (0.04)	4.67	0.19 (0.04)	4.67
PD mc × WHLSM mc					-0.04 (0.03)	-1.27
Random Effects
σ2	1.09		1.05		1.05
τ00	1.02 ID		1.03 ID		1.03 ID
N	114 ID		114 ID		114 ID
Observations	779		779		779

# plotting main effects from selected model
grid.arrange(plot_model(m3_EE,type="pred",terms="WHLSM.cm.gmc"), # trait WHLSM
             plot_model(m3_EE,type="pred",terms="WHLSM.mc"),nrow=1) # state WHLSM

# plotting interactions
sd(cleanEE$PD.mc) # RDet: 1 SD = 1.36

## [1] 1.359362

plot_model(m3_EE,type="pred",terms=c("WHLSM.mc","PD.mc [-1.36,1.36]"))

3.2.4. Influential cases

Here, we evaluate the presence of influential cases in the selected model m2.bis. Specifically, the Cook’s distance is considered as the main measure of individual-level (i.e., participant) influence on the estimated parameters, and it is recomputed by progressively excluding the most influential participants (i.e., based on the rule-of-thumb of 4/N) until all extreme values are removed.

3.2.4.1. Cook’s distance

We can see that participant S049 is a potentially influential case.

# cook's distance on the whole sample
infl <- influence(m2_EE,"ID")
plot(infl,which="cook",cutoff=4/nlevels(cleanEE$ID),xlab="Cook distance",ylab="ID",sort=TRUE)

# progressively excluding participants
infl <- influence(exclude.influence(m2_EE,"ID","S049"),"ID")
plot(infl,which="cook",cutoff=4/(nlevels(cleanEE$ID)-1),xlab="Cook distance",ylab="ID",sort=TRUE)

3.2.4.2. Coefficient change

Here, we inspect the magnitude of the changes in the estimated coefficients after the removal of potentially influential cases. We can see that the updated models do not imply substantial changes in the estimated coefficients. Thus, we choose to rely on the results obtained with the full sample.

# refitting model without influential cases
m2_EE.noInfl <- update(m2_EE,data=cleanEE[!cleanEE$ID%in%c("S049"),])

# plotting coefficients original vs. updated model
plot_models(m2_EE,m2_EE.noInfl) # main effect

plot_models(m3_EE,update(m3_EE,data=cleanEE[!cleanEE$ID%in%c("S049"),])) # interaction

# showing regression table original vs. updated model
tab_model(m2_EE,m2_EE.noInfl,dv.labels=c("Original","Updated"),show.icc=FALSE,show.p=FALSE,show.se=TRUE,
          show.r2=FALSE,collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Original			Updated
Predictors	b (SE)	CI	t	b (SE)	CI	t
(Intercept)	3.12 (0.14)	2.84 – 3.40	21.89	3.18 (0.14)	2.90 – 3.45	22.72
gender [M]	0.18 (0.21)	-0.22 – 0.58	0.88	0.12 (0.20)	-0.27 – 0.52	0.61
PD mc	-0.07 (0.03)	-0.13 – -0.02	-2.70	-0.07 (0.03)	-0.12 – -0.02	-2.54
WHLSM cm gmc	0.47 (0.08)	0.31 – 0.62	5.91	0.53 (0.08)	0.37 – 0.68	6.60
WHLSM mc	0.19 (0.04)	0.11 – 0.27	4.67	0.19 (0.04)	0.11 – 0.27	4.67
Random Effects
σ2	1.05			1.06
τ00	1.03 ID			0.96 ID
N	114 ID			113 ID
Observations	779			772

3.3. Sleep disturbances

Sleep disturbances SD rated in the morning is predicted by previous day state workaholism WHLSM.mc, in addition to psychological detachment PD, trait workaholism WHLSM.cm, and one covariate selected from the previous step (see Supplementary Material S5), namely gender.

mformulas[6] # covariates selected for SD

## [1] "SD ~ gender"

3.3.1. Data preparation

First, we prepare the data for the analyses by removing all cases of missing responses in the dependent variable or any predictor or covariate (list-wise deletion), by centering level-2 continuous predictors on the grand mean (grand-mean-centering), and by centering level-1 continuous predictors on the individual mean (person-mean-centering).

# list-wise deletion
cleanSD <- as.data.frame(na.omit(clean[,c("ID","SD","gender","PD","WHLSM", # core variables
                                          "WE","WC","sleep_dysf","flagTime","careless", # for robustness checks
                                          "position","children")])) 
cleanSD$ID <- as.factor(as.character(cleanSD$ID)) # resetting participant identifier levels
cat("Considering",nrow(cleanSD),"complete observations from",nlevels(as.factor(as.character(cleanSD$ID))),"participants")

## Considering 723 complete observations from 114 participants

# person-mean-centering lv-1 continuous predictors
wide <- cleanSD[!duplicated(cleanSD$ID),] # wide-form dataset
for(Var in c("PD","WHLSM","WE","WC")){
  wide <- cbind(wide,aggregate(cleanSD[,Var],list(cleanSD$ID),mean)[,2]) # computing individual means
  colnames(wide)[ncol(wide)] <- paste0(Var,".cm")
  cleanSD <- join(cleanSD,wide[,c("ID",paste0(Var,".cm"))],by="ID",type="left") # joining with long-form data
  cleanSD[,paste0(Var,".mc")] <- cleanSD[,Var] - cleanSD[,paste0(Var,".cm")] } # computing mean-centered scores

# grand-mean-centering lv-2 continuous predictors
for(Var in c("WHLSM.cm")){ cleanSD[,paste0(Var,".gmc")] <- cleanSD[,Var] - mean(wide[,Var]) }

# showing data
cleanSD[1:3,] # first three rows

3.3.2. Model fit

Here, we fit the multilevel models to the selected data using the default restricted maximum likelihood estimator (REML).

# m0: null model
m0_SD <- lmer(SD ~ (1|ID), # only fixed and random intercept + residual term
              data=cleanSD)

# m1: covariates
m1_SD <- lmer(SD ~ gender + PD.mc + WHLSM.cm.gmc + (1|ID), # covariates
              data=cleanSD)

# m2: state workaholism
m2_SD <- lmer(SD ~ gender + PD.mc + WHLSM.cm.gmc + WHLSM.mc + (1|ID),
              data=cleanSD)

# m3: interaction
m3_SD <- lmer(SD ~ gender + PD.mc + WHLSM.cm.gmc + WHLSM.mc + WHLSM.mc + WHLSM.mc:PD.mc + (1|ID), 
              data=cleanSD)

From the previous chunk, we see that all models converged without problems. Here, we inspect the diagnostics (i.e., normality of residual and random effect distributions, homoscedasticity, and multicollinearity) of the most complex model m3. Influential cases are analyzed in a dedicated section below.

Model m3 shows some deviation from normality especially in the lower tail of the residual distribution and both tails of random effects, although deviation from normality is slightly less marked than in model m2. Besides that, we can see that the homoscedsticity assumption holds and that none of the variance inflation factors (VIFs) shows extreme values, ruling out the risk of multicollinearity.

# normality and homoscedasticity
p <- plot_model(m3_SD,type="diag",dot.size=1) 
p[[2]] <- p[[2]]$ID
plot_grid(p,tags=TRUE,margin=c(0,0,0,0))

# homoscedasticity and multicollinearity
par(mfrow=c(1,2)) 
for(Var in c("gender")){ boxplot(resid(m3_SD) ~ cleanSD[,Var],main=paste("Residuals by",Var)) }
barplot(vif(m3_SD),main="VIF Values",xlim=c(0,10),las=2,horiz=TRUE) # variance inflation factors (VIFs)
abline(v = 5, lwd = 5, lty = 2)

Here, we better inspect the residual distribution and the fit of models specified with alternative family distributions. We can see that alternative families might better approximate the distribution of model residuals, with the log-normal solution (i.e., normal distribution with logarithmic link function) showing the best fit. Yet, it is not so better than the original model. Thus, we initially rely on the normal distribution and then consider the log-normal solution as a robustness check.

# inspecting residual distribution
descdist(resid(m3_SD)) # best fit for Gamma

## summary statistics
## ------
## min:  -2.821223   max:  4.242021 
## median:  -0.1752208 
## mean:  -1.441861e-15 
## estimated sd:  0.9856311 
## estimated skewness:  0.8569118 
## estimated kurtosis:  4.099986

# fitting model with alternative families
models <- list(
  m3_SD,
  glmer(formula=formula(m3_SD),family=Gamma(link="log"),data=cleanSD), # Gamma with log [FAILS TO CONV, max|grad| = 0.05]
  glmer(formula=formula(m3_SD),family=Gamma(link="identity"),data=cleanSD), # Gamma with identity
  glmer(formula=formula(m3_SD),family=gaussian(link="log"),data=cleanSD), # log-norm
  lmer(formula=as.formula(paste("log(SD) ~",as.character(formula(m3_SD))[3])),data=cleanSD)) # norm log

# normal Q-Q plot of model residuals
par(mfrow=c(1,5))
for(i in 1:length(models)){ 
  qqnorm(resid(models[[i]]),main=c("Norm","Gamma-log","Log-norm","log-transf")[i]); qqline(resid(models[[i]])) }

3.3.3. Results

Here, we compare the specified models based on the Akaike weight and the likelihood ratio test (with type-I error set to p < .05), and we inspect the results of the selected model(s). We can see that the inclusion of both state WHLSM (Aw = .51, \(\chi^2\)(1) = 5.67, p = .01), and its interaction with PD (Aw = .21, \(\chi^2\)(1) = 5.68, p = .02) are associated with stronger evidence (not the interaction) and significantly higher likelihood than the baseline model. Model m3 is selected as the best model.

# Akaike weight adding one model at time
Weights(AIC(m0_SD,m1_SD)) # covariates: better

##  model weights 
## [1] 0.001 0.999

Weights(AIC(m0_SD,m1_SD,m2_SD)) # state workaholism: better

##  model weights 
## [1] 0.000 0.449 0.551

Weights(AIC(m0_SD,m1_SD,m2_SD,m3_SD)) # interaction: worse

##  model weights 
## [1] 0.000 0.364 0.447 0.189

# Likelihood ratio test with m3
anova(m1_SD,m2_SD,m3_SD) # best model is m3

## refitting model(s) with ML (instead of REML)

Here, we inspect the coefficients estimated by the selected model m3 and those estimated by more parsimonious models. We can see that both trait WHLSM.cm.gmc and state WHLSM.mc are positively related to SD. In the selected model, a substantial interaction is shown such that WHLSM.mc is positively related to SD in those working days with lower RDet and RRel. Among the included covariates, gender predict substantial differences in SD, with higher sleep disturbances for females.

# regression table
tab_model(m1_SD,m2_SD,m3_SD,
          dv.labels=c("Baseline","State WHLSM","Interaction"),
          show.icc=FALSE,show.p=FALSE,show.se=TRUE,show.r2=FALSE,show.ci=FALSE,
          collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Baseline		State WHLSM		Interaction
Predictors	b (SE)	t	b (SE)	t	b (SE)	t
(Intercept)	2.75 (0.12)	23.77	2.75 (0.12)	23.77	2.75 (0.12)	23.59
gender [M]	-0.40 (0.17)	-2.36	-0.40 (0.17)	-2.36	-0.39 (0.17)	-2.34
PD mc	-0.03 (0.03)	-0.89	-0.02 (0.03)	-0.67	-0.02 (0.03)	-0.75
WHLSM cm gmc	0.34 (0.06)	5.30	0.34 (0.06)	5.30	0.34 (0.06)	5.24
WHLSM mc			0.11 (0.04)	2.62	0.11 (0.04)	2.59
PD mc × WHLSM mc					-0.07 (0.03)	-2.31
Random Effects
σ2	1.13		1.12		1.11
τ00	0.60 ID		0.60 ID		0.61 ID
N	114 ID		114 ID		114 ID
Observations	723		723		723

# plotting main effects from selected model
grid.arrange(plot_model(m3_SD,type="pred",terms="WHLSM.cm.gmc"), # trait WHLSM
             plot_model(m3_SD,type="pred",terms="WHLSM.mc"),nrow=1) # state WHLSM

# plotting interactions
sd(cleanSD$PD.mc) # PD: 1 SD = 1.36

## [1] 1.362052

plot_model(m3_SD,type="pred",terms=c("WHLSM.mc","PD.mc [-1.36,1.36]"))

3.3.4. Influential cases

Here, we evaluate the presence of influential cases in the selected model m3. Specifically, the Cook’s distance is considered as the main measure of individual-level (i.e., participant) influence on the estimated parameters, and it is recomputed by progressively excluding the most influential participants (i.e., based on the rule-of-thumb of 4/N) until all extreme values are removed.

3.3.4.1. Cook’s distance

We can see that participants S132, S049, S079, and S002 are potentially influential cases.

# cook's distance on the whole sample
infl <- influence(m3_SD,"ID")
plot(infl,which="cook",cutoff=4/nlevels(cleanSD$ID),xlab="Cook distance",ylab="ID",sort=TRUE)

# progressively excluding participants
infl <- list(
  influence(exclude.influence(m3_SD,"ID","S132"),"ID"),
  influence(exclude.influence(m3_SD,"ID",c("S132","S049")),"ID"),
  influence(exclude.influence(m3_SD,"ID",c("S132","S049","S079")),"ID"),
  influence(exclude.influence(m3_SD,"ID",c("S132","S049","S079","S002")),"ID"))
for(i in 1:length(infl)){ 
  plot(infl[[i]],which="cook",cutoff=4/(nlevels(cleanSD$ID)-i),xlab="Cook distance",ylab="ID",sort=TRUE) }

3.3.4.2. Coefficient change

Here, we inspect the magnitude of the changes in the estimated coefficients after the removal of potentially influential cases. We can see that the updated model implies an increase in the coefficient estimated for trait WHLSM.cm.gmc and gender, and a decrease in those estimated for the interaction. However, it does not drop far below the cut-off of |t|=2. Thus, we choose to rely on the results obtained with the full sample.

# refitting model without influential cases
m3_SD.noInfl <- update(m3_SD,data=cleanSD[!cleanSD$ID%in%c("S132","S049","S079","S002"),])

# plotting coefficients original vs. updated model
plot_models(m3_SD,m3_SD.noInfl)

# showing regression table original vs. updated model
tab_model(m3_SD,m3_SD.noInfl,dv.labels=c("Original","Updated"),show.icc=FALSE,show.p=FALSE,show.se=TRUE,
          show.r2=FALSE,collapse.se=TRUE,string.est="b (SE)",show.stat=TRUE,string.stat="t")

	Original			Updated
Predictors	b (SE)	CI	t	b (SE)	CI	t
(Intercept)	2.75 (0.12)	2.52 – 2.97	23.59	2.78 (0.11)	2.57 – 2.99	26.01
gender [M]	-0.39 (0.17)	-0.72 – -0.06	-2.34	-0.49 (0.15)	-0.79 – -0.19	-3.21
PD mc	-0.02 (0.03)	-0.08 – 0.03	-0.75	-0.01 (0.03)	-0.07 – 0.05	-0.41
WHLSM cm gmc	0.34 (0.06)	0.21 – 0.47	5.24	0.38 (0.06)	0.26 – 0.50	6.26
WHLSM mc	0.11 (0.04)	0.03 – 0.20	2.59	0.13 (0.04)	0.04 – 0.21	2.85
PD mc × WHLSM mc	-0.07 (0.03)	-0.14 – -0.01	-2.31	-0.07 (0.03)	-0.14 – -0.00	-2.08
Random Effects
σ2	1.11			1.08
τ00	0.61 ID			0.45 ID
N	114 ID			110 ID
Observations	723			693

4. Robusness checks

Here, we conduct a series of robustness checks (or multiverse data analysis; see Steegen et al. 2016) for each selected model by using alternative inclusion criteria, models with additional or less covariates, outlier removal, etc. The following packages and functions are used to optimize the analyses.

library(psych)

glmerAn

#' @title Generalized linear (mixed-effects) regression analysis
#' @param modelType = type of model: GLM, mixed-effects (GLMER), or cumulative link mixedm odel (CLMM)
#' @param data = data.frame of data
#' @param resp = name of the response variable (character)
#' @param fix.eff = character vector of names of the predictor(s)
#' @param REML = argument from the lme4::lmer() function, see ?lmer
#' @param ran.eff = character string indicating the random effect by using the lme4 syntax (defult: "(1|ID)")
#' @param family = character string indicating the name of the GLM(ER) family to be used in the models (default: "normal")
#' @param link =character string indicating the name of the GLM(ER) link function to be used in the models (default: "identity")
#' @param nAGQ = argument from the lme4::glmer() function, see ?glmer
#' @param mComp.baseline = character string indicating the name of the last predictor included in the baseline model to be compared with the subsequent models. If equal to NA (defult), the null model is used as the baseline model for comparison
#' @param p.adjust.method = argument from the stats::p.adjust() function (see ?p.adjust) indicating which method should be used to correct the p-values obtained from the likelihood ratio test (default: NA, for no adjustment)
#' @param key.model = character string indicating the name of the predictor(s) whose model(s) should be considered for the "key.res" output
#' @param key.predictor = character string indicating the name of the predictor to be considered by the "key.res" output
#' @param digits = number of digits for all numeric ouputs
#' @param messages = boolean indicating whether a message should be printed for each operation (defult: FALSE)
glmerAn <- function(data,modelType=c("GLMER"),resp,fix.eff,REML=TRUE,ran.eff="(1|ID)",family="normal",
                    link="identity",nAGQ=1,mComp.baseline=NA,p.adjust.method=NA,coeff.models=NA,transform=NULL,
                    plot.model=NA,plot.pred="all",key.model=NA,key.predictor=NA,digits=3,messages=FALSE){ 
  
  if(messages==TRUE){ cat("Running",modelType,"analysis of",resp,"...") }
  
  # modeling .......................................................................................
  
  # creating model formulas
  formulas <- character()
  if(modelType=="GLM"){ ran.eff <- "1" }
  null.f <- paste(resp,"~",ran.eff) # creating null model formula
  for(i in 1:length(fix.eff)){ # creating other formulas
    if(i==1){ formulas[i] <- paste(resp,"~",fix.eff[1]) } else { formulas[i] <- paste(formulas[i-1],"+",fix.eff[i])  }}
  if(modelType%in%c("GLMER","CLMM")){ if(!is.na(ran.eff)){ formulas <- paste(formulas,"+",ran.eff)
      if(substr(ran.eff,2,2)!="1"){ ranSlope <- paste(fix.eff[which(grepl(ran.eff,fix.eff))])[1]
        null.f <- gsub(ranSlope,"1",null.f)  # removing random slope from models without the related predictor
        for(i in 1:length(formulas)){ 
          if(!(grepl(ranSlope,gsub(paste(ranSlope,"[|]",sep=""),"",formulas[i])))){ 
            formulas[i] <- gsub(paste(ranSlope,"[|]",sep=""),"1|",formulas[i]) }}}
    } else { stop(message="Error: GLMER model type without ran.eff specification") }}
  if(messages==TRUE){ cat("\n\nModel specification:\n - model M0 (null):",null.f)
    for(i in 1:length(formulas)){ cat("\n - model M",i,": ",formulas[i],sep="")}}
  
  # fitting models
  models <- list()
  if(modelType=="GLM"){ if(messages==TRUE){ 
    cat("\n\nFitting GLM models of",resp,"on",nrow(data),"participants \n   using the",
        family,"family with the",link,"link function...") }
    if(family=="normal" & link=="identity"){ null.m <- lm(as.formula(null.f),data=data) # normal family
      for(i in 1:length(formulas)){ models[[i]] <- lm(formula=as.formula(formulas[i]),data=data) }
      } else if (family=="gamma") { null.m <- glm(as.formula(null.f),data=data,family=Gamma(link=link),nAGQ=nAGQ) # gamma
        for(i in 1:length(formulas)){ models[[i]] <- glm(formula=as.formula(formulas[i]),data=data,family=Gamma(link=link)) }
        } else if(family=="normal" & link!="identity"){  
          null.m <- glm(as.formula(null.f),data=data,family=gaussian(link=link)) # normal with other link functions
          for(i in 1:length(formulas)){ 
            models[[i]] <- glm(formula=as.formula(formulas[i]),data=data,family=gaussian(link=link)) }
        } else if(family=="binomial"){ 
          null.m <- glm(as.formula(null.f),data=data,family=binomial(link=link)) # logistic regression
          for(i in 1:length(formulas)){ 
            models[[i]] <- glm(formula=as.formula(formulas[i]),data=data,family=binomial(link=link))}
        } else { stop(message="Error: only normal, gamma, and binomial family are allowed, 
                      with identity, inverse, and log link functions") }
  } else if(modelType=="GLMER"){ suppressMessages(suppressWarnings(require(lme4)))
    if(messages==TRUE){ cat("\n\nFitting",modelType,"models of",resp,"on",nrow(data),"observations from",
                            nlevels(as.factor(as.character(data$ID))),"participants \n   using the",family,
                            "family with the",link,"link function using",ifelse(REML==FALSE,"ML","REML"),"estimator...") }
    if(family=="normal" & link=="identity"){ null.m <- lmer(as.formula(null.f),data=data,REML=REML) # normal  identity
      for(i in 1:length(formulas)){ models[[i]] <- lmer(formula=as.formula(formulas[i]),data=data,REML=REML) }
      } else if (family=="gamma") { null.m <- glmer(as.formula(null.f),data=data,family=Gamma(link=link),nAGQ=nAGQ) # gamma
        for(i in 1:length(formulas)){
          models[[i]]<-glmer(formula=as.formula(formulas[i]),data=data,family=Gamma(link=link),nAGQ=nAGQ) }
        } else if(family=="normal" & link!="identity"){ 
          null.m <- glmer(as.formula(null.f),data=data,family=gaussian(link=link),nAGQ=nAGQ) # normal with other links
          for(i in 1:length(formulas)){ models[[i]] <- glmer(formula=as.formula(formulas[i]),data=data,
                                                             family=gaussian(link=link),nAGQ=nAGQ) }
        } else if(family=="binomial"){ 
          null.m <- glmer(as.formula(null.f),data=data,family=binomial(link=link),nAGQ=nAGQ) # logistic
          for(i in 1:length(formulas)){ models[[i]] <- glmer(formula=as.formula(formulas[i]),data=data,
                                                             family=binomial(link=link),nAGQ=nAGQ)}
        } else { stop(message="Error: only normal, logistic, and gamma family are allowed, 
                               with identity, inverse, and log link functions") }
  } else if(modelType=="CLMM"){ suppressMessages(suppressWarnings(require(ordinal))) # cumulative link mixed models
    if(messages==TRUE){ 
      cat("\n\nFitting",modelType,"models of",resp,"on",nrow(data),"observations from",
          nlevels(as.factor(as.character(data$ID))),"participants \n   using Cumulative Link Mixed Models") }
    data[,resp] <- factor(data[,resp],ordered=TRUE) # response variable as ordered factor
    null.m <- suppressWarnings(clmm(as.formula(gsub("~","~ 1 +",null.f)),data=data)) # suppress formula warning (bugged)
    for(i in 1:length(formulas)){ models[[i]] <- suppressWarnings(clmm(formula=as.formula(formulas[i]),data=data,nAGQ=nAGQ)) }
  } else { stop(message="Error: modelType can only be 'GLM', 'GLMER', or 'CLMM'") }
  
  # outputs...............................................................................................................
  if(messages==TRUE){ cat("\n\nGenerating models outputs...") }
  
  # model comparison
    # likelihood ratio test
    if(messages==TRUE){ cat("\n\n - Running likelihood ratio test:") } 
    suppressMessages(suppressWarnings(require(knitr))); suppressMessages(suppressWarnings(require(MuMIn))) 
    m.num <- 1
    if(is.na(mComp.baseline)){ bsl <- null.m  # selecting baseline model
      } else { m.num <- grep(mComp.baseline,fix.eff)[1] + 1
         bsl <- models[[m.num - 1]] }
    if(modelType!="CLMM"){ lrt <- as.data.frame(anova(bsl,models[[m.num]]))
      if(length(models)>m.num){
        for(i in m.num:(length(models)-1)){ lrt <- rbind(lrt,as.data.frame(anova(models[[i]],models[[i+1]]))[2,]) }}
      } else { lrt <- as.data.frame(ordinal:::anova.clm(bsl,models[[m.num]])) # use anova.clm() to avoid env. issue
        if(length(models)>m.num){
          for(i in m.num:(length(models)-1)){ 
            lrt <- rbind(lrt,as.data.frame(ordinal:::anova.clm(models[[i]],models[[i+1]]))[2,]) }}}
    rownames(lrt) <- c(ifelse(is.na(mComp.baseline),"Null model","Baseline"),
                       fix.eff[m.num:length(fix.eff)])
    if(!is.na(p.adjust.method)){ # p-value corrections for multiple comparison
      if(messages==TRUE){ cat(" (applying",p.adjust.method,"p-values correction)")}
      lrt[!is.na(lrt$`Pr(>Chisq)`),"Pr(>Chisq)"] <- p.adjust(lrt[!is.na(lrt$`Pr(>Chisq)`),"Pr(>Chisq)"],
                                                             method=p.adjust.method)  }
    
    # Akaike weights
    AICs <- lrt[1:2,"AIC"] # Akaike weight
    
    # updating key results
    key <- lrt[which(grepl(key.predictor,row.names(lrt))),] # key results
    if(nrow(key)>1){ key <- key[1,] }
    key.results <- data.frame(sig.LRT=key[,ncol(key)]<0.05, # sig.LRT
                              higher.Aw=key$AIC==min(AICs[1:which(AICs==key$AIC)])) # higher.Aw
  
  # estimated parameters from key.model
  modSummary <- summary(models[[which(fix.eff==key.model)]])
  modSummary <- modSummary$coefficients
  if(modelType=="CLMM"){ modSummary <- modSummary[nlevels(data[,resp]):nrow(modSummary),] }
  key <- round(modSummary[which(grepl(key.predictor,row.names(modSummary))),3][1],2) # taking only first coeff for key.results
  key.results <- cbind(key.results,t.196=abs(key)>1.96,t=key)
  
  # returning key results (sig. LRT, Aw higher than previous model, t > 1.96)
  return(key.results) }

glmerMed

#' @title Generalized linear mixed-effects mediation analysis
#' @param data = data.frame of data
#' @param resp = name of the response variable (character)
#' @param treat = name of the treatment variable (character)
#' @param med = name of the mediator variable (character)
#' @param fix.eff = character vector of names of the predictor(s)
#' @param REML = argument from the lme4::lmer() function, see ?lmer
#' @param ran.eff = character string indicating the random effect by using the lme4 syntax (defult: "(1|ID)")
#' @param family = character string indicating the name of the GLM(ER) family to be used in the models (default: "normal")
#' @param link = character string indicating the name of the GLM(ER) link function to be used in the models (default: "identity")
#' @param sims = number of Monte Carlo draws for quasi-Bayesian CI and p-value (default: 1000)
#' @param alpha.level = significance level (default: 0.05)
#' @param noCov = boolean indicating wheter covariates should be considered in the mediation models or not
#' @param messages = boolean indicating whether a message should be printed for each operation (defult: FALSE)
glmerMed <- function(data,resp,treat,med,fix.eff,REML=TRUE,ran.eff="(1|ID)",family="normal",
                     link="identity",sims=1000,alpha.level=0.05,noCov=FALSE,messages=FALSE){ 
  
  library(lme4); library(mediation)
  
  if(messages==TRUE){ cat("Running",modelType,"analysis of",resp,"...") }
  
  # modeling .......................................................................................
  
  # creating model formulas
  mMed.f <- paste(gsub(".mc","",med),"~",fix.eff[1])
  if(length(fix.eff)>1){
    for(i in 2:length(fix.eff)){ mMed.f <- paste(mMed.f,fix.eff[i],sep=" + ") }} # mediator model 
  mOut.f <- paste(gsub(gsub(".mc","",med),resp,mMed.f),med,sep=" + ") # outcome model
  mMed.f <- paste(mMed.f,ran.eff,sep=" + ") # adding random effects
  mOut.f <- paste(mOut.f,ran.eff,sep=" + ")
  
  # fitting models
  if(messages==TRUE){ cat("\n\nFitting GLMR models of",resp,"and",med,"on",nrow(data),"observations from",
                            nlevels(as.factor(as.character(data$ID))),"participants \n   using the",family,
                            "family with the",link,"link function using",ifelse(REML==FALSE,"ML","REML"),"estimator...") }
  if(family=="normal" & link=="identity"){ 
    mMed <- lmer(as.formula(mMed.f),data=data,REML=REML)
    mOut <- lmer(as.formula(mOut.f),data=data,REML=REML)
  } else if(family=="normal" & link!="identity") { 
      mMed <- glmer(as.formula(mMed.f),data=data,REML=REML,family=gaussian(link=link)) 
      mOut <- glmer(as.formula(mOut.f),data=data,REML=REML,family=gaussian(link=link)) 
  } else if(family=="gamma"){
      mMed <- glmer(as.formula(mMed.f),data=data,family=Gamma(link=link)) 
      mOut <- glmer(as.formula(mOut.f),data=data,family=Gamma(link=link))  
  } else { stop(message="Error: only normal, logistic, and gamma family are allowed, 
                               with identity, inverse, and log link functions") }
  
  # mediation analysis
  if(noCov==TRUE){
     medM <- mediation::mediate(model.m=mMed,model.y=mOut,treat=treat,mediator=med, # mediation & output models, treatment & mediator
                 covariates=NULL,boot=FALSE,sims=sims) # quasi-Bayesian confidence intervals
  } else {
    medM <- mediation::mediate(model.m=mMed,model.y=mOut,treat=treat,mediator=med, # mediation & output models, treatment & mediator
                 covariates=data[,fix.eff[fix.eff!=treat]], # covariates
                 boot=FALSE,sims=sims) } # quasi-Bayesian confidence intervals
  
  # returning estimated quasi-Bayesian p-value
  key.results <- data.frame(indirect.p=medM$d.avg.p,ind.p.sig=medM$d.avg.p<alpha.level,
                            direct.p=medM$z.avg.p,ind.p.sig=medM$z.avg.p<alpha.level)
  return(key.results) }

Moreover, since a further robustness check to be implemented for all outcomes is reproducing the analyses with the full sample, we also reprocess the full diary dataset by removing missing responses and mean-centering predictors.

Show code

# afternoon blood pressure
cleanBP_aft_full <- as.data.frame(na.omit(diary[,c("ID","SBP_aft","DBP_aft","gender","age","BMI","WHLSM")])) # listwise del
                                                   
cleanBP_aft_full$ID <- as.factor(as.character(cleanBP_aft_full$ID)) # resetting participant identifier levels
cat(nrow(cleanBP_aft_full),"complete obs from",nlevels(as.factor(as.character(cleanBP_aft_full$ID))),"participants")

## 907 complete obs from 135 participants

for(Var in c("WHLSM")){ # person-mean-centering lv-1 continuous predictors
  cleanBP_aft_prelqs_full <- cbind(prelqs[prelqs$ID%in%levels(cleanBP_aft_full$ID),],
                             aggregate(cleanBP_aft_full[,Var],list(cleanBP_aft_full$ID),mean)[,2]) # individual means
  colnames(cleanBP_aft_prelqs_full)[ncol(cleanBP_aft_prelqs_full)] <- paste0(Var,".cm")
  cleanBP_aft_full <- join(cleanBP_aft_full,cleanBP_aft_prelqs_full[,c("ID",paste0(Var,".cm"))],by="ID",type="left")
  cleanBP_aft_full[,paste0(Var,".mc")] <- cleanBP_aft_full[,Var] - cleanBP_aft_full[,paste0(Var,".cm")] } # mean-centered
for(Var in c("age","BMI","WHLSM.cm")){ 
  cleanBP_aft_full[,paste0(Var,".gmc")] <- cleanBP_aft_full[,Var] - mean(cleanBP_aft_prelqs_full[,Var]) } # gmc

# evening blood pressure
cleanBP_eve_full <- as.data.frame(na.omit(diary[,c("ID","SBP_eve","DBP_eve","gender","age","BMI", # list-wise deletion
                                                   "WHLSM","PD")]))
cleanBP_eve_full$ID <- as.factor(as.character(cleanBP_eve_full$ID)) # resetting participant identifier levels
cat(nrow(cleanBP_eve_full),"complete obs from",nlevels(as.factor(as.character(cleanBP_eve_full$ID))),"participants")

## 825 complete obs from 134 participants

for(Var in c("WHLSM","PD")){ # person-mean-centering lv-1 continuous predictors
  cleanBP_eve_prelqs_full <- cbind(prelqs[prelqs$ID%in%levels(cleanBP_eve_full$ID),],
                             aggregate(cleanBP_eve_full[,Var],list(cleanBP_eve_full$ID),mean)[,2]) # individual means
  colnames(cleanBP_eve_prelqs_full)[ncol(cleanBP_eve_prelqs_full)] <- paste0(Var,".cm")
  cleanBP_eve_full <- join(cleanBP_eve_full,cleanBP_eve_prelqs_full[,c("ID",paste0(Var,".cm"))],by="ID",type="left")
  cleanBP_eve_full[,paste0(Var,".mc")] <- cleanBP_eve_full[,Var] - cleanBP_eve_full[,paste0(Var,".cm")] } # mean-centered
for(Var in c("age","BMI","WHLSM.cm")){ 
  cleanBP_eve_full[,paste0(Var,".gmc")] <- cleanBP_eve_full[,Var] - mean(cleanBP_eve[!duplicated(cleanBP_eve$ID),Var]) } # gmc

# afternoon-to-evening blood pressure
cleanBP_med_eve_full <- as.data.frame(na.omit(diary[,c("ID","SBP_aft","DBP_aft","SBP_eve","DBP_eve", # list-wise deletion
                                                       "gender","age","BMI","WHLSM","PD")])) 
cleanBP_med_eve_full$ID <- as.factor(as.character(cleanBP_med_eve_full$ID)) 
wide <- cleanBP_med_eve_full[!duplicated(cleanBP_med_eve_full$ID),] 
for(Var in c("WHLSM","PD","SBP_aft","DBP_aft")){
  wide <- cbind(wide,aggregate(cleanBP_med_eve_full[,Var],list(cleanBP_med_eve_full$ID),mean)[,2]) 
  colnames(wide)[ncol(wide)] <- paste0(Var,".cm")
  cleanBP_med_eve_full <- join(cleanBP_med_eve_full,wide[,c("ID",paste0(Var,".cm"))],by="ID",type="left") 
  cleanBP_med_eve_full[,paste0(Var,".mc")] <- cleanBP_med_eve_full[,Var] - cleanBP_med_eve_full[,paste0(Var,".cm")] } 
for(Var in c("age","BMI","WHLSM.cm")){ cleanBP_med_eve_full[,paste0(Var,".gmc")] <- 
  cleanBP_med_eve_full[,Var] - mean(wide[,Var]) } # gmc
cat("Considering",nrow(cleanBP_med_eve_full),"complete obs from",
    nlevels(as.factor(as.character(cleanBP_med_eve_full$ID))),"participants")

## Considering 820 complete obs from 134 participants

# emotional exhaustion
cleanEE_full <- as.data.frame(na.omit(diary[,c("ID","EE","gender","PD","WHLSM")])) # listwise deletion
cleanEE_full$ID <- as.factor(as.character(cleanEE_full$ID)) # resetting participant identifier levels
cat(nrow(cleanEE_full),"complete obs from",nlevels(as.factor(as.character(cleanEE_full$ID))),"participants")

## 827 complete obs from 134 participants

for(Var in c("PD","WHLSM")){ # person-mean-centering lv-1 continuous predictors
  clean_prelqs_full <- cbind(prelqs[prelqs$ID%in%levels(cleanEE_full$ID),],
                             aggregate(cleanEE_full[,Var],list(cleanEE_full$ID),mean)[,2]) # individual means
  colnames(clean_prelqs_full)[ncol(clean_prelqs_full)] <- paste0(Var,".cm")
  cleanEE_full <- join(cleanEE_full,clean_prelqs_full[,c("ID",paste0(Var,".cm"))],by="ID",type="left")
  cleanEE_full[,paste0(Var,".mc")] <- cleanEE_full[,Var] - cleanEE_full[,paste0(Var,".cm")] } # mean-centered scores
for(Var in c("WHLSM.cm")){ cleanEE_full[,paste0(Var,".gmc")] <- cleanEE_full[,Var] - mean(cleanEE_full[,Var]) } # gmc

# sleep disturbances
cleanSD_full <- as.data.frame(na.omit(diary[,c("ID","SD","gender","dailyHassles_eve","PD","WHLSM")])) # listwise del
cleanSD_full$ID <- as.factor(as.character(cleanSD_full$ID)) # resetting participant identifier levels
cat(nrow(cleanSD_full),"complete obs from",nlevels(as.factor(as.character(cleanSD_full$ID))),"participants")

## 753 complete obs from 132 participants

for(Var in c("PD","WHLSM")){ # person-mean-centering lv-1 continuous predictors
  clean_prelqs_full <- cbind(prelqs[prelqs$ID%in%levels(cleanSD_full$ID),],
                             aggregate(cleanSD_full[,Var],list(cleanSD_full$ID),mean)[,2]) # individual means
  colnames(clean_prelqs_full)[ncol(clean_prelqs_full)] <- paste0(Var,".cm")
  cleanSD_full <- join(cleanSD_full,clean_prelqs_full[,c("ID",paste0(Var,".cm"))],by="ID",type="left")
  cleanSD_full[,paste0(Var,".mc")] <- cleanSD_full[,Var] - cleanSD_full[,paste0(Var,".cm")] } # mean-centered scores
for(Var in c("WHLSM.cm")){ cleanSD_full[,paste0(Var,".gmc")] <- cleanSD_full[,Var] - mean(cleanSD_full[,Var]) } # gmc

As a further check, we include position as an additional covariate, which we recode into two levels, namely “Employee/Project” vs. “Managers/Employers”.

cleanBP_aft$position <- as.factor(gsub("Employee","employee/project",
                                       gsub("Project","employee/project",cleanBP_aft$position)))
cleanBP_eve$position <- as.factor(gsub("Employee","employee/project",
                                       gsub("Project","employee/project",cleanBP_eve$position)))
cleanBP_med_eve$position <- as.factor(gsub("Employee","employee/project",
                                       gsub("Project","employee/project",cleanBP_med_eve$position)))
cleanEE$position <- as.factor(gsub("Employee","employee/project",
                                       gsub("Project","employee/project",cleanEE$position)))
cleanSD$position <- as.factor(gsub("Employee","employee/project",
                                       gsub("Project","employee/project",cleanSD$position)))

Finally, we load the raw preliminary questionnaire item scores to compute the composite score at the retrospective version of the DUWAS, used for a robustness check.

# isolating raw item scores at the retrospective DUWAS scale
prelqs.retroWHLSM <- prelqs[,c("ID",paste0("duwas",1:10))]
prelqs.retroWHLSM <- prelqs.retroWHLSM[prelqs.retroWHLSM$ID %in% clean$ID,] # subsampling participants included in the clean dataset
prelqs.retroWHLSM$ID <- as.factor(as.character(prelqs.retroWHLSM$ID))
cat("Included participants =",nlevels(prelqs.retroWHLSM$ID))

## Included participants = 114

# computing Cronbach's alpha and 95% CI
psych::alpha(prelqs.retroWHLSM[,paste0("duwas",1:10)])$feldt

## 
##    95% confidence boundaries (Feldt)
##  lower alpha upper
##   0.74   0.8  0.85

# computing grand-mean-centered composite score and adding it to all datasets used below
prelqs.retroWHLSM$WHLSM.retro <- apply(prelqs.retroWHLSM[,paste0("duwas",1:10)],1,mean)
cleanBP_aft <- plyr::join(cleanBP_aft,prelqs.retroWHLSM[,c("ID","WHLSM.retro")],by="ID",type="left")
cleanBP_aft$WHLSM.retro.gmc <- cleanBP_aft$WHLSM.retro - mean(cleanBP_aft$WHLSM.retro)
cleanBP_eve <- plyr::join(cleanBP_eve,prelqs.retroWHLSM[,c("ID","WHLSM.retro")],by="ID",type="left")
cleanBP_eve$WHLSM.retro.gmc <- cleanBP_eve$WHLSM.retro - mean(cleanBP_eve$WHLSM.retro)
cleanBP_med_eve <- plyr::join(cleanBP_med_eve,prelqs.retroWHLSM[,c("ID","WHLSM.retro")],by="ID",type="left")
cleanBP_med_eve$WHLSM.retro.gmc <- cleanBP_med_eve$WHLSM.retro - mean(cleanBP_med_eve$WHLSM.retro)
cleanEE <- plyr::join(cleanEE,prelqs.retroWHLSM[,c("ID","WHLSM.retro")],by="ID",type="left")
cleanEE$WHLSM.retro.gmc <- cleanEE$WHLSM.retro - mean(cleanEE$WHLSM.retro)
cleanSD <- plyr::join(cleanSD,prelqs.retroWHLSM[,c("ID","WHLSM.retro")],by="ID",type="left")
cleanSD$WHLSM.retro.gmc <- cleanSD$WHLSM.retro - mean(cleanSD$WHLSM.retro)

4.1. Blood pressure

4.1.1. Afternoon BP

For afternoon blood pressure, we implement the following robustness checks:

1. No Infl: we remove influential participants

2. No dysf/drugs: we remove all participants reporting sleep dysfunctions, hormonal or psychoactive medications, in addition to those meeting exclusion criteria for blood pressure

3. No Cov: we remove all covariates, that is we only include WHLSM.mc and its interactions as model predictors

4. All in: we include all complete observations from all participants, including those meeting the exclusion criteria for compliance and blood pressure

5. ML: we refit the models by using the Maximum Likelihood estimator, rather than the Restricted Maximum Likelihood

6. Rand slope: we include the random slope for WHLSM.mc

7. logTransf: we log-transform the response variable values before fitting the models

8. confounders_aft: we include potentially confounding factors for blood pressure (e.g., smoking, physical activity) reported in the afternoon as an additional control variable

9. position: we include job position (Employee/Project vs. Manager/(Self-)Employer) as an additional control variable

10. children: we include the number of children as an additional control variable

11. No flagBP: we exclude all observations that were reprocessed due to extreme BP values (see Supplementary Material S3)

12. No flagTime: we exclude all observations that were flagged due to their associated timing (e.g., morning BP recorded in the afternoon) (see Supplementary Material S3)

13. No careless: we exclude one participant S137 flagged as potentially careless (see Supplementary Material S3).

14. WE: we replace the predictor term for state workaholism with the composite score at the working excessively dimension

15. WC: we replace the predictor term for state workaholism with the composite score at the working excessively dimension

16. WHLSM retro: we replace the predictor term for trait workaholism with the composite score at the retrospective version of the DUWAS included in the preliminary questionnaire.

In all cases, the results are consistent with those reported in the main analyses, showing substantial contribution and main effect of state WHLSM.mc.

SBP_aft

checks <- c("Original","No Infl","No dysf/drugs","No Cov","All in","ML","Rand slope","logTransf",
            "confounders_aft","position","children","No flagBP","No flagTime","No careless","WE","WC","WHLSM retro")

# main effect of WHLSM.mc
predictors <- c("gender","age.gmc","BMI.gmc","WHLSM.cm.gmc","WHLSM.mc") # predictors
r <- "SBP_aft" # response variable
p <- cbind(check=checks,
      rbind(glmerAn(data=cleanBP_aft,resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # original
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft[!cleanBP_aft$ID%in%c("S096","S082"),], # without influential
                    resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft[cleanBP_aft$sleep_dysf=="No" & # without participants meeting exclusion criteria
                                       cleanBP_aft$psy_drugs=="No" & cleanBP_aft$psy_drugs=="No",], 
                    resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=c("WHLSM.mc"), # without covariates
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft_full,resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # full sample
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # maximum likelihood
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",REML=FALSE),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # random slope
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",ran.eff="(WHLSM.mc|ID)"),
            glmerAn(data=cleanBP_aft,resp="log(SBP_aft)",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # log-transf
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=c(predictors[1:3],"confounders_aft", # adding confounders_aft
                                                      predictors[4:length(predictors)]), 
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=c(predictors[1:3],"position", # adding position
                                                      predictors[4:length(predictors)]), 
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=c(predictors[1:3],"children", # adding children
                                                      predictors[4:length(predictors)]), 
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft[cleanBP_aft$flagBP_aft==FALSE,], # without flagged BP
                    resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft[cleanBP_aft$flagTime==FALSE,], # without flagged times
                    resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft[cleanBP_aft$careless==FALSE,], # without careless participant
                    resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=gsub("WHLSM.mc","WE.mc",predictors), 
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WE.mc",key.model="WE.mc"), # working excess
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=gsub("WHLSM.mc","WC.mc",predictors), 
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WC.mc",key.model="WC.mc"), # working compulsively
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors), # WHLSM retro
                    mComp.baseline="WHLSM.retro.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc")))
kable(p)

check	sig.LRT	higher.Aw	t.196	t
Original	TRUE	TRUE	TRUE	3.90
No Infl	TRUE	TRUE	TRUE	3.08
No dysf/drugs	TRUE	TRUE	TRUE	3.33
No Cov	TRUE	TRUE	TRUE	3.90
All in	TRUE	TRUE	TRUE	3.43
ML	TRUE	TRUE	TRUE	3.90
Rand slope	TRUE	TRUE	TRUE	3.61
logTransf	TRUE	TRUE	TRUE	3.68
confounders_aft	TRUE	TRUE	TRUE	3.91
position	TRUE	TRUE	TRUE	3.90
children	TRUE	TRUE	TRUE	3.90
No flagBP	TRUE	TRUE	TRUE	3.80
No flagTime	TRUE	TRUE	TRUE	3.88
No careless	TRUE	TRUE	TRUE	3.93
WE	TRUE	TRUE	TRUE	3.75
WC	TRUE	TRUE	TRUE	3.26
WHLSM retro	TRUE	TRUE	TRUE	3.90

DBP_aft

checks <- c("Original","No Infl","No dysf/drugs","No Cov","All in","ML","Rand slope","logTransf",
            "confounders_aft","position","children","No flagBP","No flagTime","No careless","WE","WC","WHLSM retro")

# main effect of WHLSM.mc
predictors <- c("gender","age.gmc","BMI.gmc","WHLSM.cm.gmc","WHLSM.mc") # predictors
r <- "DBP_aft" # response variable
p <- cbind(check=checks,
      rbind(glmerAn(data=cleanBP_aft,resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # original
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft[!cleanBP_aft$ID%in%c("S082"),], # without influential
                    resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft[cleanBP_aft$sleep_dysf=="No" & # without participants meeting exclusion criteria
                                       cleanBP_aft$psy_drugs=="No" & cleanBP_aft$psy_drugs=="No",], 
                    resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=c("WHLSM.mc"), # without covariates
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft_full,resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # full sample
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # maximum likelihood
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",REML=FALSE),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # random slope
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",ran.eff="(WHLSM.mc|ID)"),
            glmerAn(data=cleanBP_aft,resp="log(SBP_aft)",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # log-transf
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=c(predictors[1:3],"confounders_aft", # adding confounders_aft
                                                      predictors[4:length(predictors)]), 
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=c(predictors[1:3],"position", # adding position
                                                      predictors[4:length(predictors)]), 
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=c(predictors[1:3],"children", # adding children
                                                      predictors[4:length(predictors)]), 
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft[cleanBP_aft$flagBP_aft==FALSE,], # without flagged BP
                    resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft[cleanBP_aft$flagTime==FALSE,], # without flagged times
                    resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft[cleanBP_aft$careless==FALSE,], # without careless
                    resp=r,fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=gsub("WHLSM.mc","WE.mc",predictors), 
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WE.mc",key.model="WE.mc"), # working excess
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=gsub("WHLSM.mc","WC.mc",predictors), 
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WC.mc",key.model="WC.mc"), # working compulsively
            glmerAn(data=cleanBP_aft,resp=r,fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors), 
                    mComp.baseline="WHLSM.retro.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc")))
kable(p)

check	sig.LRT	higher.Aw	t.196	t
Original	TRUE	TRUE	TRUE	3.74
No Infl	TRUE	TRUE	TRUE	3.45
No dysf/drugs	TRUE	TRUE	TRUE	3.96
No Cov	TRUE	TRUE	TRUE	3.74
All in	TRUE	TRUE	TRUE	3.17
ML	TRUE	TRUE	TRUE	3.74
Rand slope	TRUE	TRUE	TRUE	3.56
logTransf	TRUE	TRUE	TRUE	3.68
confounders_aft	TRUE	TRUE	TRUE	3.73
position	TRUE	TRUE	TRUE	3.74
children	TRUE	TRUE	TRUE	3.74
No flagBP	TRUE	TRUE	TRUE	3.20
No flagTime	TRUE	TRUE	TRUE	3.72
No careless	TRUE	TRUE	TRUE	3.79
WE	TRUE	TRUE	TRUE	3.45
WC	TRUE	TRUE	TRUE	3.27
WHLSM retro	TRUE	TRUE	TRUE	3.74

4.1.2. Evening BP

For evening blood pressure, we implement the following robustness checks:

1. No Infl: we remove influential participants

2. No dysf/drugs: we remove all participants reporting sleep dysfunctions, hormonal or psychoactive medications, in addition to those meeting exclusion criteria for blood pressure

3. No Cov: we remove all covariates, that is we only include WHLSM.mc and its interactions as model predictors

4. All in: we include all complete observations from all participants, including those meeting the exclusion criteria for compliance and blood pressure

5. ML: we refit the models by using the Maximum Likelihood estimator, rather than the Restricted Maximum Likelihood

6. Rand slope: we include the random slope for WHLSM.mc

7. logTransf: we log-transform the response variable values before fitting the models

8. confounders_eve: we include potentially confounding factors for blood pressure (e.g., smoking, physical activity) reported in the afternoon as an additional control variable

9. position: we include job position (Employee/Project vs. Manager/(Self-)Employer) as an additional control variable

10. children: we include the number of children as an additional control variable

11. No flagBP: we exclude all observations that were reprocessed due to extreme BP values (see Supplementary Material S3)

12. No flagTime: we exclude all observations that were flagged due to their associated timing (e.g., morning BP recorded in the afternoon) (see Supplementary Material S3)

13. No careless: we exclude one participant S137 flagged as potentially careless (see Supplementary Material S3)

14. day1: we include recording day as a further covariate (i.e., 1 = first day, 2 = any other day). In all cases, the results are consistent with those reported in the main analyses, showing no substantial contribution and main effect of state WHLSM.mc.

15. WE: we replace the predictor term for state workaholism with the composite score at the working excessively dimension

16. WC: we replace the predictor term for state workaholism with the composite score at the working excessively dimension

17. WHLSM retro: we replace the predictor term for trait workaholism with the composite score at the retrospective version of the DUWAS included in the preliminary questionnaire.

# setting robustness checks
checks <- c("Original","No Infl","No dysf/drugs","No Cov","All in","ML","Rand slope","logTransf",
            "confounders_eve","position","children","No flagBP","No flagTime","No careless","day1",
            "WE","WC","WHLSM retro")

# creating categorical variable day1 (i.e., first day vs. all other days)
cleanBP_eve$day1 <- 0 
cleanBP_eve[cleanBP_eve$day==1,"day1"] <- 1
cleanBP_eve$day1 <- as.factor(cleanBP_eve$day1)

In all cases, the results are consistent with those reported in the main analyses, showing no substantial contribution or main effect of state WHLSM.mc (only substantial in a few cases for DBP_eve) and no substantial interactions with psychological detachment.

MAIN EFFECT

SBP_eve

predictors <- c("gender","age.gmc","BMI.gmc","PD.mc","WHLSM.cm.gmc","WHLSM.mc") # predictors
r <- "SBP_eve" # response variable
key <- "WHLSM.mc" # key model and key predictor
bsl <- "WHLSM.cm.gmc" # baseline model
p <- cbind(check=checks,
      rbind(glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # original
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[!cleanBP_eve$ID%in%c("S082","S096"),], # without influential
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$sleep_dysf=="No" & # without participants meeting exclusion criteria
                                       cleanBP_eve$psy_drugs=="No" & cleanBP_eve$psy_drugs=="No",], 
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c("WHLSM.mc"), # without covariates
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve_full,resp=r,fix.eff=predictors,mComp.baseline=bsl, # full sample
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # maximum likelihood
                    key.predictor=key,key.model=key,REML=FALSE),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # random slope (SINGULAR FIT)
                    key.predictor=key,key.model=key,ran.eff="(WHLSM.mc|ID)"),
            glmerAn(data=cleanBP_eve,resp="log(SBP_eve)",fix.eff=predictors,mComp.baseline=bsl, # log-transf
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:5],"confounders_eve", # adding confounders_eve
              predictors[6:length(predictors)]),mComp.baseline="confounders_eve",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:5],"position", # adding position
              predictors[6:length(predictors)]),mComp.baseline="position",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:5],"children", # adding children
              predictors[6:length(predictors)]),mComp.baseline="children",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$flagBP_eve==FALSE,], # without flagged BP
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$flagTime==FALSE,], # without flagged times
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$careless==FALSE,], # without careless
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:5],"day1",predictors[6:length(predictors)]), # adding day
                    mComp.baseline="day1",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.mc","WE.mc",predictors), # working excessively
                    mComp.baseline=bsl,key.predictor="WE.mc",key.model="WE.mc"),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.mc","WC.mc",predictors), # working compulsively
                    mComp.baseline=bsl,key.predictor="WC.mc",key.model="WC.mc"),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors), # WHLSM retro
                    mComp.baseline="WHLSM.retro.gmc",key.predictor=key,key.model=key)))
kable(p)

check	sig.LRT	higher.Aw	t.196	t
Original	FALSE	FALSE	FALSE	0.76
No Infl	FALSE	FALSE	FALSE	0.39
No dysf/drugs	FALSE	FALSE	FALSE	0.48
No Cov	FALSE	FALSE	FALSE	1.05
All in	FALSE	FALSE	FALSE	0.92
ML	FALSE	FALSE	FALSE	0.76
Rand slope	FALSE	FALSE	FALSE	0.81
logTransf	FALSE	FALSE	FALSE	0.51
confounders_eve	FALSE	FALSE	FALSE	0.76
position	FALSE	FALSE	FALSE	0.76
children	FALSE	FALSE	FALSE	0.76
No flagBP	FALSE	FALSE	FALSE	0.52
No flagTime	FALSE	FALSE	FALSE	0.65
No careless	FALSE	FALSE	FALSE	0.84
day1	FALSE	FALSE	FALSE	1.04
WE	FALSE	FALSE	FALSE	0.37
WC	FALSE	FALSE	FALSE	1.02
WHLSM retro	FALSE	FALSE	FALSE	0.76

DBP_eve

r <- "DBP_eve" # response variable
p <- cbind(check=checks,
      rbind(glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key), # original
            glmerAn(data=cleanBP_eve[!cleanBP_eve$ID%in%c("S082","S080"),], # without influential
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$sleep_dysf=="No" & # without participants meeting exclusion criteria
                                       cleanBP_eve$psy_drugs=="No" & cleanBP_eve$psy_drugs=="No",], 
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c("WHLSM.mc"),key.predictor=key,key.model=key), # without covariates
            glmerAn(data=cleanBP_eve_full,resp=r,fix.eff=predictors,mComp.baseline=bsl, # full sample
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # maximum likelihood
                    key.predictor=key,key.model=key,REML=FALSE),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # random slope
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp="log(DBP_eve)",fix.eff=predictors,mComp.baseline=bsl, # Gamma log
                    key.predictor=key,key.model=key,family="gamma",link="log"),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:5],"confounders_eve", # adding confounders_eve
              predictors[6:length(predictors)]),mComp.baseline="confounders_eve",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:5],"position", # adding position
              predictors[6:length(predictors)]),mComp.baseline="position",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:5],"children", # adding children
              predictors[6:length(predictors)]),mComp.baseline="children",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$flagBP_eve==FALSE,], # without flagged BP
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$flagTime==FALSE,], # without flagged times
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$careless==FALSE,], # without careless participants
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:5],"day1",predictors[6:length(predictors)]), # adding day
                    mComp.baseline="day1",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.mc","WE.mc",predictors), # working excessively
                    mComp.baseline=bsl,key.predictor="WE.mc",key.model="WE.mc"),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.mc","WC.mc",predictors), # working compulsively
                    mComp.baseline=bsl,key.predictor="WC.mc",key.model="WC.mc"),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors), 
                    mComp.baseline="WHLSM.retro.gmc",key.predictor=key,key.model=key)))
kable(p)

check	sig.LRT	higher.Aw	t.196	t
Original	FALSE	TRUE	FALSE	1.95
No Infl	TRUE	TRUE	TRUE	2.06
No dysf/drugs	FALSE	TRUE	FALSE	1.44
No Cov	TRUE	TRUE	TRUE	2.04
All in	FALSE	TRUE	FALSE	1.64
ML	FALSE	TRUE	FALSE	1.96
Rand slope	FALSE	TRUE	FALSE	1.95
logTransf	FALSE	TRUE	FALSE	1.72
confounders_eve	FALSE	TRUE	FALSE	1.93
position	FALSE	TRUE	FALSE	1.95
children	FALSE	TRUE	FALSE	1.95
No flagBP	TRUE	TRUE	TRUE	2.04
No flagTime	TRUE	TRUE	TRUE	2.02
No careless	FALSE	TRUE	FALSE	1.93
day1	TRUE	TRUE	TRUE	2.29
WE	FALSE	TRUE	FALSE	1.87
WC	FALSE	TRUE	FALSE	1.65
WHLSM retro	FALSE	TRUE	FALSE	1.95

INTERACTION

SBP_eve

# interaction with RDet.mc
r <- "SBP_eve"
predictors <- c(predictors,"PD.mc:WHLSM.mc") # predictors
key <- "PD.mc:WHLSM.mc" # key model and key predictor
bsl <- "WHLSM.mc" # baseline model
p <- cbind(check=checks,
      rbind(glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # original
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[!cleanBP_eve$ID%in%c("S096","S082"),], # without influential
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve_full,resp=r,fix.eff=predictors,mComp.baseline=bsl, # full sample
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$sleep_dysf=="No" & # without participants meeting exclusion criteria
                                       cleanBP_eve$psy_drugs=="No" & cleanBP_eve$psy_drugs=="No",], 
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c("PD.mc","WHLSM.mc","PD.mc:WHLSM.mc"), # without covariates
                    mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # maximum likelihood
                    key.predictor=key,key.model=key,REML=FALSE),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # random slope (SINGULAR FIT)
                    key.predictor=key,key.model=key,ran.eff="(WHLSM.mc|ID)"),
            glmerAn(data=cleanBP_eve,resp="log(SBP_eve)",fix.eff=predictors,mComp.baseline=bsl, # log-transf
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:6],"confounders_eve", # adding confounders_eve
                                                      predictors[7:length(predictors)]),
                    mComp.baseline="confounders_eve",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:6],"position", # adding position
                                                     predictors[7:length(predictors)]),
                    mComp.baseline="position",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:6],"children", # adding children
                                                     predictors[7:length(predictors)]),
                    mComp.baseline="children",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$flagBP_eve==FALSE,], # without flagged BP
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$flagTime==FALSE,], # without flagged times
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$careless==FALSE,], # without careless participant
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:6],"day1",predictors[7:length(predictors)]), # adding day
                    mComp.baseline="day1",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.mc","WE.mc",predictors), # working excessively
                    mComp.baseline="WE.mc",key.predictor="PD.mc:WE.mc",key.model="PD.mc:WE.mc"),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.mc","WC.mc",predictors), # working compulsively
                    mComp.baseline="WC.mc",key.predictor="PD.mc:WC.mc",key.model="PD.mc:WC.mc"),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors),  # WHLSM.retro
                    mComp.baseline=bsl,key.predictor=key,key.model=key)))
kable(p)

check	sig.LRT	higher.Aw	t.196	t
Original	FALSE	FALSE	FALSE	0.80
No Infl	FALSE	FALSE	FALSE	0.58
No dysf/drugs	FALSE	FALSE	FALSE	0.95
No Cov	FALSE	TRUE	FALSE	1.46
All in	FALSE	FALSE	FALSE	0.85
ML	FALSE	FALSE	FALSE	0.81
Rand slope	FALSE	FALSE	FALSE	0.69
logTransf	FALSE	FALSE	FALSE	0.91
confounders_eve	FALSE	FALSE	FALSE	0.79
position	FALSE	FALSE	FALSE	0.79
children	FALSE	FALSE	FALSE	0.80
No flagBP	FALSE	FALSE	FALSE	0.51
No flagTime	FALSE	FALSE	FALSE	1.04
No careless	FALSE	FALSE	FALSE	0.82
day1	FALSE	FALSE	FALSE	0.72
WE	FALSE	FALSE	FALSE	1.05
WC	FALSE	FALSE	FALSE	0.38
WHLSM retro	FALSE	FALSE	FALSE	0.79

DBP_eve

r <- "DBP_eve"
p <- cbind(check=checks,
      rbind(glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # original
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[!cleanBP_eve$ID%in%c("S096","S080"),], # without influential
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$sleep_dysf=="No" & # without participants meeting exclusion criteria (SINGULAR FIT)
                                       cleanBP_eve$psy_drugs=="No" & cleanBP_eve$psy_drugs=="No",], 
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c("PD.mc","WHLSM.mc","PD.mc:WHLSM.mc"), # without covariates
                    mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve_full,resp=r,fix.eff=predictors,mComp.baseline=bsl, # full sample
                    key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # maximum likelihood
                    key.predictor=key,key.model=key,REML=FALSE),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=predictors,mComp.baseline=bsl, # random slope
                    key.predictor=key,key.model=key,ran.eff="(WHLSM.mc|ID)"),
            glmerAn(data=cleanBP_eve,resp="DBP_eve",fix.eff=predictors,mComp.baseline=bsl, # Gamma log
                    key.predictor=key,key.model=key,family="gamma",link="log"),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:6],"confounders_eve", # adding confounders_eve
                                                      predictors[7:length(predictors)]),
                    mComp.baseline="confounders_eve",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:6],"position", # adding position
                                                      predictors[7:length(predictors)]),
                    mComp.baseline="position",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:6],"children", # adding children
                                                      predictors[7:length(predictors)]),
                    mComp.baseline="children",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$flagBP_eve==FALSE,], # without flagged BP
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$flagTime==FALSE,], # without flagged times
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve[cleanBP_eve$careless==FALSE,], # without careless
                    resp=r,fix.eff=predictors,mComp.baseline=bsl,key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=c(predictors[1:6],"day1",predictors[7:length(predictors)]), # adding day
                    mComp.baseline="day1",key.predictor=key,key.model=key),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.mc","WE.mc",predictors), # working excessively
                    mComp.baseline="WE.mc",key.predictor="PD.mc:WE.mc",key.model="PD.mc:WE.mc"),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.mc","WC.mc",predictors), # working compulsively
                    mComp.baseline="WC.mc",key.predictor="PD.mc:WC.mc",key.model="PD.mc:WC.mc"),
            glmerAn(data=cleanBP_eve,resp=r,fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors), # WHLSM.retro
                    mComp.baseline=bsl,key.predictor=key,key.model=key)))
kable(p)

check	sig.LRT	higher.Aw	t.196	t
Original	FALSE	FALSE	FALSE	0.01
No Infl	FALSE	FALSE	FALSE	-0.01
No dysf/drugs	FALSE	FALSE	FALSE	0.52
No Cov	FALSE	FALSE	FALSE	0.10
All in	FALSE	FALSE	FALSE	0.63
ML	FALSE	FALSE	FALSE	0.01
Rand slope	FALSE	FALSE	FALSE	-0.10
logTransf	FALSE	FALSE	FALSE	-0.09
confounders_eve	FALSE	FALSE	FALSE	-0.04
position	FALSE	FALSE	FALSE	-0.01
children	FALSE	FALSE	FALSE	0.02
No flagBP	FALSE	FALSE	FALSE	0.22
No flagTime	FALSE	FALSE	FALSE	0.05
No careless	FALSE	FALSE	FALSE	0.05
day1	FALSE	FALSE	FALSE	-0.09
WE	FALSE	FALSE	FALSE	0.14
WC	FALSE	FALSE	FALSE	-0.13
WHLSM retro	FALSE	FALSE	FALSE	0.02

4.1.3. Mediations

For Afternoon-to-Evening mediations, we implement the following robustness checks:

1. No Infl: we remove all influential participants found for afternoon and evening BP

2. No dysf/drugs: we remove all participants reporting sleep dysfunctions, hormonal or psychoactive medications, in addition to those meeting exclusion criteria for blood pressure

3. No Cov: we remove all covariates, that is we only include WHLSM.mc and its interactions as model predictors

4. All in: we include all complete observations from all participants, including those meeting the exclusion criteria for compliance and blood pressure

5. ML: we refit the models by using the Maximum Likelihood estimator, rather than the Restricted Maximum Likelihood

6. Rand slope: we include the random slope for WHLSM.mc

7. log-transf: we log-transform blood pressure before fitting the models

8. confounders: we include potentially confounding factors for blood pressure (e.g., smoking, physical activity) reported in the afternoon or in the evening as an additional control variable

9. position: we include job position (Employee/Project vs. Manager/(Self-)Employer) as an additional control variable

10. children: we include the number of children as an additional control variable

11. No flagBP: we exclude all observations that were reprocessed due to extreme BP values either in the afternoon or in the evening (see Supplementary Material S3); morning flagged cases are not considered to avoid loosing too many observations

12. No flagTime: we exclude all observations that were flagged due to their associated timing (e.g., morning BP recorded in the afternoon) (see Supplementary Material S3)

13. No careless: we exclude one participant S137 flagged as potentially careless (see Supplementary Material S3)

14. day1: we include recording day as a further covariate (i.e., 1 = first day, 2 = any other day). In all cases, the results are consistent with those reported in the main analyses, showing no substantial contribution and main effect of state WHLSM.mc.

15. WE: we replace the predictor term for state workaholism with the composite score at the working excessively dimension

16. WC: we replace the predictor term for state workaholism with the composite score at the working excessively dimension

17. WHLSM retro: we replace the predictor term for trait workaholism with the composite score at the retrospective version of the DUWAS included in the preliminary questionnaire.

# setting robustness checks
checks <- c("Original","No Infl","No dysf/drugs","No Cov","All in","ML","Rand slope","log-transf",
            "confounders","position","children","No flagBP","No flagTime","No careless","day1",
            "WE","WC","WHLSM.retro")

# recoding variables
cleanBP_med_eve$conf <- FALSE # summary of afternoon and evening confounders
cleanBP_med_eve[cleanBP_med_eve$confounders_aft==TRUE | cleanBP_med_eve$confounders_eve==TRUE,"conf"] <- TRUE
summary(cleanBP_med_eve$conf)

##    Mode   FALSE    TRUE 
## logical     531     186

cleanBP_med_eve$flagBP <- FALSE # summary of afternoon and evening flagBP
cleanBP_med_eve[cleanBP_med_eve$flagBP_aft==TRUE | cleanBP_med_eve$flagBP_eve==TRUE,"flagBP"] <- TRUE
summary(cleanBP_med_eve$flagBP)

##    Mode   FALSE    TRUE 
## logical     701      16

cleanBP_med_eve$day1 <- 0 # creating categorical variable day1 (i.e., first day vs. all other days)
cleanBP_med_eve[cleanBP_med_eve$day==1,"day1"] <- 1
cleanBP_med_eve$day1 <- as.factor(cleanBP_med_eve$day1)
summary(cleanBP_med_eve$day1)

##   0   1 
## 634  83

In all cases, the results for SBP are consistent with those reported in the main analyses, showing significant indirect but not direct relationship between state WHLSM.mc and SBP_eve. Most results for DBP are consistent as well, but the indirect relationship is reduced in one case, i.e., with the removal of participants reporting sleep dysfunctions, hormonal, or psychoactive medications. Similarly, the direct relationship becomes significant with the inclusion of day1 as a covariate. Whereas these findings question the generalizability of such relationship, the high number of consistent robustness checks provides some evidence that the relationship is observable in our sample.

SBP_eve

predictors <- c("gender","age.gmc","BMI.gmc","PD.mc","WHLSM.cm.gmc","WHLSM.mc") # predictors
r <- "SBP_eve" # response variable
t <- "WHLSM.mc" # key model and key predictor
m <- "SBP_aft.mc"
p <- cbind(check=checks,
      rbind(glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=predictors), # original
            glmerMed(data=cleanBP_med_eve[!cleanBP_med_eve$ID%in%c("S082","S096"),], # without influential cases
                    resp=r,treat=t,med=m,fix.eff=predictors),
            glmerMed(data=cleanBP_med_eve[cleanBP_med_eve$sleep_dysf=="No" & # without participants meeting exclusion criteria
                                       cleanBP_med_eve$psy_drugs=="No" & cleanBP_med_eve$psy_drugs=="No",], 
                    resp=r,treat=t,med=m,fix.eff=predictors),
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=c("WHLSM.mc"),noCov=TRUE), # w/o covariates
            glmerMed(data=cleanBP_med_eve_full,resp=r,treat=t,med=m,fix.eff=predictors), # full sample
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=predictors,REML=FALSE), # maximum likelihood
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=predictors, # random slope (singular fit)
                     ran.eff="(WHLSM.mc|ID)"), 
            glmerMed(data=cleanBP_med_eve,resp="log(SBP_eve)",treat=t,med=m,fix.eff=predictors), # log transformation
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=c(predictors[1:5],"conf", # adding confounders
              predictors[6:length(predictors)])),
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=c(predictors[1:5],"position", # adding position
              predictors[6:length(predictors)])),
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=c(predictors[1:5],"children", # adding children
              predictors[6:length(predictors)])),
            glmerMed(data=cleanBP_med_eve[cleanBP_med_eve$flagBP==FALSE,], # without flagged BP cases
                     resp=r,treat=t,med=m,fix.eff=predictors),
            glmerMed(data=cleanBP_med_eve[cleanBP_med_eve$flagTime==FALSE,], # without flagged times
                     resp=r,treat=t,med=m,fix.eff=predictors),
            glmerMed(data=cleanBP_med_eve[cleanBP_med_eve$careless==FALSE,], # without careless participants
                     resp=r,treat=t,med=m,fix.eff=predictors),
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=c(predictors[1:5],"day1", # adding day
              predictors[6:length(predictors)])),
            glmerMed(data=cleanBP_med_eve,resp=r,treat="WE.mc",med=m,
                     fix.eff=gsub("WHLSM.mc","WE.mc",predictors)), # working excessively
            glmerMed(data=cleanBP_med_eve,resp=r,treat="WC.mc",med=m,
                     fix.eff=gsub("WHLSM.mc","WC.mc",predictors)), # working compulsively
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,
                     fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro",predictors)))) # WHLSM retro
kable(p)

check	indirect.p	ind.p.sig	direct.p	ind.p.sig.1
Original	0.002	TRUE	0.748	FALSE
No Infl	0.026	TRUE	0.966	FALSE
No dysf/drugs	0.010	TRUE	0.926	FALSE
No Cov	0.002	TRUE	0.634	FALSE
All in	0.000	TRUE	0.630	FALSE
ML	0.000	TRUE	0.762	FALSE
Rand slope	0.010	TRUE	0.684	FALSE
log-transf	0.000	TRUE	0.970	FALSE
confounders	0.000	TRUE	0.762	FALSE
position	0.006	TRUE	0.826	FALSE
children	0.000	TRUE	0.784	FALSE
No flagBP	0.004	TRUE	0.942	FALSE
No flagTime	0.002	TRUE	0.910	FALSE
No careless	0.000	TRUE	0.744	FALSE
day1	0.002	TRUE	0.598	FALSE
WE	0.006	TRUE	0.922	FALSE
WC	0.012	TRUE	0.534	FALSE
WHLSM.retro	0.002	TRUE	0.796	FALSE

DBP_eve

r <- "DBP_eve" # response variable
p <- cbind(check=checks,
      rbind(glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=predictors), # original
            glmerMed(data=cleanBP_med_eve[!cleanBP_med_eve$ID%in%c("S082","S096"),], # without influential cases
                    resp=r,treat=t,med=m,fix.eff=predictors),
            glmerMed(data=cleanBP_med_eve[cleanBP_med_eve$sleep_dysf=="No" & # without participants meeting exclusion criteria
                                       cleanBP_med_eve$psy_drugs=="No" & cleanBP_med_eve$psy_drugs=="No",], 
                    resp=r,treat=t,med=m,fix.eff=predictors),
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=c("WHLSM.mc"),noCov=TRUE), # w/o covariates
            glmerMed(data=cleanBP_med_eve_full,resp=r,treat=t,med=m,fix.eff=predictors), # full sample
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=predictors,REML=FALSE), # maximum likelihood
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=predictors, # random slope (singular fit)
                     ran.eff="(WHLSM.mc|ID)"), 
            glmerMed(data=cleanBP_med_eve,resp="log(DBP_eve)",treat=t,med=m,fix.eff=predictors), # log-transformation
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=c(predictors[1:5],"conf", # adding confounders
              predictors[6:length(predictors)])),
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=c(predictors[1:5],"position", # adding position
              predictors[6:length(predictors)])),
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=c(predictors[1:5],"children", # adding children
              predictors[6:length(predictors)])),
            glmerMed(data=cleanBP_med_eve[cleanBP_med_eve$flagBP==FALSE,], # without flagged BP cases
                     resp=r,treat=t,med=m,fix.eff=predictors),
            glmerMed(data=cleanBP_med_eve[cleanBP_med_eve$flagTime==FALSE,], # without flagged times
                     resp=r,treat=t,med=m,fix.eff=predictors),
            glmerMed(data=cleanBP_med_eve[cleanBP_med_eve$careless==FALSE,], # without careless participants
                     resp=r,treat=t,med=m,fix.eff=predictors),
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,fix.eff=c(predictors[1:5],"day1", # adding day
              predictors[6:length(predictors)])),
            glmerMed(data=cleanBP_med_eve,resp=r,treat="WE.mc",med=m,
                     fix.eff=gsub("WHLSM.mc","WE.mc",predictors)), # working excessively
            glmerMed(data=cleanBP_med_eve,resp=r,treat="WC.mc",med=m,
                     fix.eff=gsub("WHLSM.mc","WC.mc",predictors)), # working compulsively
            glmerMed(data=cleanBP_med_eve,resp=r,treat=t,med=m,
                     fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors)))) # WHLSM retro
kable(p)

check	indirect.p	ind.p.sig	direct.p	ind.p.sig.1
Original	0.012	TRUE	0.080	FALSE
No Infl	0.024	TRUE	0.064	FALSE
No dysf/drugs	0.142	FALSE	0.214	FALSE
No Cov	0.018	TRUE	0.082	FALSE
All in	0.006	TRUE	0.198	FALSE
ML	0.012	TRUE	0.082	FALSE
Rand slope	0.034	TRUE	0.112	FALSE
log-transf	0.020	TRUE	0.138	FALSE
confounders	0.020	TRUE	0.108	FALSE
position	0.016	TRUE	0.070	FALSE
children	0.014	TRUE	0.094	FALSE
No flagBP	0.034	TRUE	0.068	FALSE
No flagTime	0.026	TRUE	0.076	FALSE
No careless	0.016	TRUE	0.082	FALSE
day1	0.028	TRUE	0.042	TRUE
WE	0.016	TRUE	0.100	FALSE
WC	0.020	TRUE	0.176	FALSE
WHLSM.retro	0.012	TRUE	0.084	FALSE

4.2. Emotional Exhaustion

For emotional exhaustion, we implement the following robustness checks:

1. No Infl: we remove influential participants

2. No Cov: we remove all covariates, that is we only include WHLSM.mc and its interactions as model predictors

3. All in: we include all complete observations from all participants, including those meeting the exclusion criteria for compliance and blood pressure

4. ML: we refit the models by using the Maximum Likelihood estimator, rather than the Restricted Maximum Likelihood

5. Rand slope: we include the random slope for WHLSM.mc

6. Gamma-log: we refit the models using the Gamma family with the logarithmic link function

7. logNorm: we refit the models with using log-normal GLM

8. position: we include job position (Employee/Project vs. Manager/(Self-)Employer) as an additional control variable

9. children: we include the number of children as an additional control variable

10. No flagTime: we exclude all observations that were flagged due to their associated timing (e.g., morning BP recorded in the afternoon) (see Supplementary Material S3)

11. No careless: we exclude one participant S137 flagged as potentially careless (see Supplementary Material S3).

12. day: we include day as a further continuous covariate (i.e., 1 = first day, 2 = second day, etc.)

13. WE: we replace the predictor term for state workaholism with the composite score at the working excessively dimension

14. WC: we replace the predictor term for state workaholism with the composite score at the working excessively dimension

15. WHLSM retro: we replace the predictor term for trait workaholism with the composite score at the retrospective version of the DUWAS included in the preliminary questionnaire.

checks <- c("Original","No Infl","No Cov","All in","ML","Rand slope","Gamma-log","logNorm",
            "position","children","No flagTime","No careless","day","WE","WC","WHLSM retro")

In all cases, the results are consistent with those reported in the main analyses, showing substantial contribution and main effect of state WHLSM.mc, but with no substantial interaction with PD.mc (only significant when using the log-normal distribution). Thus, we interpret these findings as a sign of the consistency of the estimated relationship between state WHLSM.mc and EE.

MAIN EFFECT

# main effect of WHLSM.mc
predictors <- c("gender","PD.mc","WHLSM.cm.gmc","WHLSM.mc")
p <- cbind(check=checks,
      rbind(glmerAn(data=cleanEE,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # original
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanEE[!cleanEE$ID%in%c("S049"),], # without influential
                    resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=c("WHLSM.mc"), # without covariates
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanEE_full,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # full sample
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # maximum likelihood
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",REML=FALSE),
            glmerAn(data=cleanEE,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # random slope
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",ran.eff="(WHLSM.mc|ID)"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # Gamma-log family
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",family="gamma",link="log"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # log-normal family
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",link="log"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=c(predictors[1],"position", # adding position
                                                      predictors[2:length(predictors)]),
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=c(predictors[1],"children", # adding children
                                                      predictors[2:length(predictors)]),
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanEE[cleanEE$flagTime==FALSE,],resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", 
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"), # without flagged times
            glmerAn(data=cleanEE[cleanEE$careless==FALSE,],resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", 
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"), # without careless participants
            glmerAn(data=cleanEE,resp="EE",fix.eff=c(predictors[1:3],"day",predictors[4:length(predictors)]),
                    mComp.baseline="day",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=gsub("WHLSM.mc","WE.mc",predictors),
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WE.mc",key.model="WE.mc"), # working excessively
            glmerAn(data=cleanEE,resp="EE",fix.eff=gsub("WHLSM.mc","WC.mc",predictors),
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WC.mc",key.model="WC.mc"), # working compulsively
            glmerAn(data=cleanEE,resp="EE",fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors),
                    mComp.baseline="WHLSM.retro.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"))) # WHLSM retro
kable(p)

check	sig.LRT	higher.Aw	t.196	t
Original	TRUE	TRUE	TRUE	4.67
No Infl	TRUE	TRUE	TRUE	4.67
No Cov	TRUE	TRUE	TRUE	4.88
All in	TRUE	TRUE	TRUE	4.68
ML	TRUE	TRUE	TRUE	4.68
Rand slope	TRUE	TRUE	TRUE	3.88
Gamma-log	TRUE	TRUE	TRUE	4.96
logNorm	TRUE	TRUE	TRUE	4.41
position	TRUE	TRUE	TRUE	4.67
children	TRUE	TRUE	TRUE	4.67
No flagTime	TRUE	TRUE	TRUE	4.60
No careless	TRUE	TRUE	TRUE	4.54
day	TRUE	TRUE	TRUE	4.45
WE	TRUE	TRUE	TRUE	4.37
WC	TRUE	TRUE	TRUE	4.04
WHLSM retro	TRUE	TRUE	TRUE	4.67

INTERACTION

predictors <- c(predictors,"PD.mc:WHLSM.mc")
p <- cbind(check=checks,
      rbind(glmerAn(data=cleanEE,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.mc", # original
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanEE[!cleanEE$ID%in%c("S049"),], # without influential cases
                    resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.mc",
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=c("PD.mc","WHLSM.mc","PD.mc:WHLSM.mc"), # without covariates
                    mComp.baseline="WHLSM.mc", key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanEE_full,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.mc", # full sample
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.mc", # maximum likelihood
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc",REML=FALSE),
            glmerAn(data=cleanEE,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.mc", # random slope
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc",ran.eff="(WHLSM.mc|ID)"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.mc", # gamma-log family
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc",family="gamma",link="log"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.mc", # log-normal family
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc",link="log"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=c(predictors[1],"position", # adding position
                                                      predictors[2:length(predictors)]),
                    mComp.baseline="WHLSM.mc",key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanEE,resp="EE",fix.eff=c(predictors[1],"children", # adding children
                                                      predictors[2:length(predictors)]),
                    mComp.baseline="WHLSM.mc",key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanEE[cleanEE$flagTime==FALSE,],resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.mc",
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"), # withoug flagged times
            glmerAn(data=cleanEE[cleanEE$careless==FALSE,],resp="EE",fix.eff=predictors,mComp.baseline="WHLSM.mc",
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"), # without careless participant
            glmerAn(data=cleanEE,resp="EE",fix.eff=c(predictors[1:3],"day",predictors[4:length(predictors)]),
                    mComp.baseline="WHLSM.mc",key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"), # day as covariate
            glmerAn(data=cleanEE,resp="EE",fix.eff=gsub("WHLSM.mc","WE.mc",predictors),
                    mComp.baseline="WE.mc",key.predictor="PD.mc:WE.mc",key.model="PD.mc:WE.mc"), # working excessively
            glmerAn(data=cleanEE,resp="EE",fix.eff=gsub("WHLSM.mc","WC.mc",predictors),
                    mComp.baseline="WC.mc",key.predictor="PD.mc:WC.mc",key.model="PD.mc:WC.mc"), # working compulsively
            glmerAn(data=cleanEE,resp="EE",fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors),
                    mComp.baseline="WHLSM.mc",key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"))) # WHLSM retro
kable(p)

check	sig.LRT	higher.Aw	t.196	t
Original	FALSE	FALSE	FALSE	-1.27
No Infl	FALSE	FALSE	FALSE	-1.24
No Cov	FALSE	FALSE	FALSE	-1.35
All in	FALSE	FALSE	FALSE	-1.20
ML	FALSE	FALSE	FALSE	-1.27
Rand slope	FALSE	FALSE	FALSE	-1.34
Gamma-log	FALSE	FALSE	FALSE	-0.62
logNorm	FALSE	TRUE	FALSE	-1.51
position	FALSE	FALSE	FALSE	-1.28
children	FALSE	FALSE	FALSE	-1.25
No flagTime	FALSE	FALSE	FALSE	-1.19
No careless	FALSE	FALSE	FALSE	-1.31
day	FALSE	FALSE	FALSE	-1.24
WE	FALSE	FALSE	FALSE	-1.25
WC	FALSE	FALSE	FALSE	-1.05
WHLSM retro	FALSE	FALSE	FALSE	-1.23

4.3. Sleep disturbances

For sleep disturbances, we implement the following robustness checks:

1. No Infl: we remove influential participants

2. No sleepDysf: we remove all participants reporting sleep dysfunctions

3. No Cov: we remove all covariates, that is we only include WHLSM.mc and its interactions as model predictors

4. All in: we include all complete observations from all participants, including those meeting the exclusion criteria for compliance and blood pressure

5. ML: we refit the models by using the Maximum Likelihood estimator, rather than the Restricted Maximum Likelihood

6. Rand slope: we include the random slope for WHLSM.mc

7. logNorm: we refit the models with using log-normal GLM

8. position: we include job position (Employee/Project vs. Manager/(Self-)Employer) as an additional control variable

9. children: we include the number of children as an additional control variable

10. No flagTime: we exclude all observations that were flagged due to their associated timing (e.g., morning BP recorded in the afternoon) (see Supplementary Material S3)

11. No careless: we exclude one participant S137 flagged as potentially careless (see Supplementary Material S3)

12. WE: we replace the predictor term for state workaholism with the composite score at the working excessively dimension

13. WC: we replace the predictor term for state workaholism with the composite score at the working excessively dimension

14. WHLSM retro: we replace the predictor term for trait workaholism with the composite score at the retrospective version of the DUWAS included in the preliminary questionnaire

checks <- c("Original","No Infl","No sleepDysf","No Cov","All in","ML","Rand slope","logNorm",
            "position","children","No flagTime","No careless","WE","WC","WHLSM retro")

In all but two cases (i.e., only the interaction but not the main effect of state workaholism is substantial when including the random slope and when using the working compulsively dimension rather than the total state workaholism score), the results are consistent with those reported in the main analyses, showing substantial main effect of state WHLSM.mc and interaction. We interpret these findings as a sign of the consistency of the estimated interactions between WHLSM.mc and PD.mc for SD.

MAIN EFFECT

# main effect of WHLSM.mc
predictors <- c("gender","PD.mc","WHLSM.cm.gmc","WHLSM.mc")
p <- cbind(check=checks,
      rbind(glmerAn(data=cleanSD,resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # original
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanSD[!cleanSD$ID%in%c("S132","S049","S079","S002"),], # without influential
                    resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanSD[clean$sleep_dysf!="Yes",], # without participants with sleep dysf
                    resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanSD,resp="SD",fix.eff=c("WHLSM.mc"), # without covariates
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanSD_full,resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # full sample
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanSD,resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # maximum likelihood
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",REML=FALSE),
            glmerAn(data=cleanSD,resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # random slope
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",ran.eff="(WHLSM.mc|ID)"),
            glmerAn(data=cleanSD,resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", # log-normal family
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc",link="log"),
            glmerAn(data=cleanSD,resp="SD",fix.eff=c(predictors[1],"position", # adding position
                                                      predictors[2:length(predictors)]),
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanSD,resp="SD",fix.eff=c(predictors[1],"children", # adding children
                                                      predictors[2:length(predictors)]),
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanSD[cleanSD$flagTime==FALSE,],resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.cm.gmc", 
                    key.predictor="WHLSM.mc",key.model="WHLSM.mc"), # without flagged times
            glmerAn(data=cleanSD[cleanSD$careless==FALSE,],resp="SD",fix.eff=predictors,  # without careless
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"),
            glmerAn(data=cleanSD,resp="SD",fix.eff=gsub("WHLSM.mc","WE.mc",predictors),
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WE.mc",key.model="WE.mc"), # working excessively
            glmerAn(data=cleanSD,resp="SD",fix.eff=gsub("WHLSM.mc","WC.mc",predictors),
                    mComp.baseline="WHLSM.cm.gmc",key.predictor="WC.mc",key.model="WC.mc"), # working compulsively
            glmerAn(data=cleanSD,resp="SD",fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors),
                    mComp.baseline="WHLSM.retro.gmc",key.predictor="WHLSM.mc",key.model="WHLSM.mc"))) # WHLSM retro
kable(p)

check	sig.LRT	higher.Aw	t.196	t
Original	TRUE	TRUE	TRUE	2.62
No Infl	TRUE	TRUE	TRUE	2.87
No sleepDysf	TRUE	TRUE	TRUE	2.61
No Cov	TRUE	TRUE	TRUE	2.68
All in	TRUE	TRUE	TRUE	2.79
ML	TRUE	TRUE	TRUE	2.62
Rand slope	TRUE	TRUE	FALSE	1.21
logNorm	TRUE	TRUE	TRUE	3.32
position	TRUE	TRUE	TRUE	2.62
children	TRUE	TRUE	TRUE	2.62
No flagTime	TRUE	TRUE	TRUE	2.77
No careless	TRUE	TRUE	TRUE	2.65
WE	TRUE	TRUE	TRUE	2.77
WC	FALSE	TRUE	FALSE	1.94
WHLSM retro	TRUE	TRUE	TRUE	2.62

INTERACTION

predictors <- c(predictors,"PD.mc:WHLSM.mc")
p <- cbind(check=checks,
      rbind(glmerAn(data=cleanSD,resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.mc", # original
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanSD[!cleanSD$ID%in%c("S132","S049","S079","S002"),], # without influential cases
                    resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.mc",
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanSD[clean$sleep_dysf!="Yes",], # without participants with sleep dysfunctions
                    resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.mc",
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"), 
            glmerAn(data=cleanSD,resp="SD",fix.eff=c("PD.mc","WHLSM.mc","PD.mc:WHLSM.mc"), # without covariates
                    mComp.baseline="WHLSM.mc", key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanSD_full,resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.mc", # full sample
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanSD,resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.mc", # maximum likelihood
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc",REML=FALSE),
            glmerAn(data=cleanSD,resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.mc", # random slope
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc",ran.eff="(WHLSM.mc|ID)"),
            glmerAn(data=cleanSD,resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.mc", # log-normal family
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc",link="log"),
            glmerAn(data=cleanSD,resp="SD",fix.eff=c(predictors[1],"position", # adding position
                                                      predictors[2:length(predictors)]),
                    mComp.baseline="WHLSM.mc",key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanSD,resp="SD",fix.eff=c(predictors[1],"children", # adding children
                                                      predictors[2:length(predictors)]),
                    mComp.baseline="WHLSM.mc",key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"),
            glmerAn(data=cleanSD[cleanSD$flagTime==FALSE,],resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.mc",
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"), # withoug flagged times
            glmerAn(data=cleanSD[cleanSD$careless==FALSE,],resp="SD",fix.eff=predictors,mComp.baseline="WHLSM.mc",
                    key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"), # without careless participants
            
            glmerAn(data=cleanSD,resp="SD",fix.eff=gsub("WHLSM.mc","WE.mc",predictors),
                    mComp.baseline="WE.mc",key.predictor="PD.mc:WE.mc",key.model="PD.mc:WE.mc"), # working excessively
            glmerAn(data=cleanSD,resp="SD",fix.eff=gsub("WHLSM.mc","WC.mc",predictors),
                    mComp.baseline="WC.mc",key.predictor="PD.mc:WC.mc",key.model="PD.mc:WC.mc"), # working compulsively
            
            glmerAn(data=cleanSD,resp="SD",fix.eff=gsub("WHLSM.cm.gmc","WHLSM.retro.gmc",predictors),
                    mComp.baseline="WHLSM.mc",key.predictor="PD.mc:WHLSM.mc",key.model="PD.mc:WHLSM.mc"))) # WHLSM retro
kable(p)

check	sig.LRT	higher.Aw	t.196	t
Original	TRUE	TRUE	TRUE	-2.31
No Infl	TRUE	TRUE	TRUE	-2.08
No sleepDysf	TRUE	TRUE	TRUE	-2.14
No Cov	TRUE	TRUE	TRUE	-2.44
All in	TRUE	TRUE	TRUE	-2.33
ML	TRUE	TRUE	TRUE	-2.31
Rand slope	TRUE	TRUE	TRUE	-2.76
logNorm	TRUE	TRUE	TRUE	-2.56
position	TRUE	TRUE	TRUE	-2.33
children	TRUE	TRUE	TRUE	-2.39
No flagTime	TRUE	TRUE	TRUE	-2.37
No careless	TRUE	TRUE	TRUE	-2.30
WE	TRUE	TRUE	TRUE	-1.99
WC	TRUE	TRUE	TRUE	-2.23
WHLSM retro	TRUE	TRUE	TRUE	-2.29

5. Outputs

Here, we generate and save the regression tables reporting the results estimated by the selected models. For each model, we visualize the unstandardized coefficients (b), the standard error (SE), and the 95% bootstrap confidence intervals computed with 10,000 iterations.

NSIM = 10000

# coefficients afternoon BP
tab_model(m1_SBP_aft,m2_SBP_aft,m1_DBP_aft,m2_DBP_aft,
          dv.labels=paste0(rep(c("SBP_aft_","DBP_aft_"),each=2),c("baseline","WHLSM.mc")),
          show.icc=FALSE,show.p=FALSE,show.se=TRUE,show.r2=FALSE,collapse.se=TRUE,string.est="b (SE)")

	SBP_aft_baseline		SBP_aft_WHLSM.mc		DBP_aft_baseline		DBP_aft_WHLSM.mc
Predictors	b (SE)	CI	b (SE)	CI	b (SE)	CI	b (SE)	CI
(Intercept)	118.25 (1.97)	114.38 – 122.12	118.25 (1.97)	114.38 – 122.12	77.32 (1.46)	74.45 – 80.20	77.32 (1.46)	74.45 – 80.19
gender [M]	5.39 (2.92)	-0.35 – 11.13	5.39 (2.92)	-0.35 – 11.13	-0.38 (2.17)	-4.65 – 3.88	-0.38 (2.17)	-4.65 – 3.88
age gmc	0.41 (0.12)	0.18 – 0.65	0.41 (0.12)	0.18 – 0.65	0.28 (0.09)	0.10 – 0.45	0.28 (0.09)	0.10 – 0.45
BMI gmc	1.10 (0.43)	0.25 – 1.95	1.10 (0.43)	0.25 – 1.95	0.96 (0.32)	0.33 – 1.59	0.96 (0.32)	0.33 – 1.59
WHLSM cm gmc	1.33 (1.14)	-0.91 – 3.57	1.33 (1.14)	-0.91 – 3.57	1.11 (0.85)	-0.55 – 2.77	1.11 (0.85)	-0.55 – 2.77
WHLSM mc			1.56 (0.40)	0.77 – 2.34			1.30 (0.35)	0.61 – 1.98
Random Effects
σ2	108.96		106.74		82.14		80.61
τ00	195.27 ID		195.55 ID		104.48 ID		104.66 ID
N	106 ID		106 ID		106 ID		106 ID
Observations	787		787		787		787

# bootstrap CI SBP_aft
confint.merMod(m1_SBP_aft,parm=3:(length(fixef(m1_SBP_aft))+2),method="boot",nsim=NSIM) # M1

##                    2.5 %      97.5 %
## (Intercept)  114.4396311 121.9985679
## genderM       -0.4483034  11.1027967
## age.gmc        0.1813017   0.6594479
## BMI.gmc        0.2431582   1.9373830
## WHLSM.cm.gmc  -0.9433809   3.5699117

confint.merMod(m2_SBP_aft,parm=3:(length(fixef(m2_SBP_aft))+2),method="boot",nsim=NSIM) # M2

##                    2.5 %      97.5 %
## (Intercept)  114.2751153 122.1679379
## genderM       -0.3637160  11.0986313
## age.gmc        0.1788889   0.6521283
## BMI.gmc        0.2525861   1.9644597
## WHLSM.cm.gmc  -0.8950650   3.5760164
## WHLSM.mc       0.7611169   2.3512745

# bootstrap CI DBP_aft
confint.merMod(m1_DBP_aft,parm=5:(length(fixef(m1_SBP_aft))+2),method="boot",nsim=NSIM) # M1

##                    2.5 %    97.5 %
## age.gmc       0.09865134 0.4548729
## BMI.gmc       0.32664247 1.6065864
## WHLSM.cm.gmc -0.55033958 2.7842322

confint.merMod(m2_DBP_aft,parm=5:(length(fixef(m2_SBP_aft))+2),method="boot",nsim=NSIM) # M2

##                    2.5 %    97.5 %
## age.gmc       0.09754627 0.4569716
## BMI.gmc       0.33492992 1.5929573
## WHLSM.cm.gmc -0.60166092 2.7686367
## WHLSM.mc      0.61872084 1.9673084

# coefficients evening BP
tab_model(m1_SBP_eve,m2_SBP_eve,m3_SBP_eve,m1_DBP_eve,m2_DBP_eve,m3_DBP_eve,
          dv.labels=paste0(rep(c("SBP_eve_","DBP_eve_"),each=3),c("baseline","WHLSM.mc","interaction")),
          show.icc=FALSE,show.p=FALSE,show.se=TRUE,show.r2=FALSE,collapse.se=TRUE,string.est="b (SE)")

	SBP_eve_baseline		SBP_eve_WHLSM.mc		SBP_eve_interaction		DBP_eve_baseline		DBP_eve_WHLSM.mc		DBP_eve_interaction
Predictors	b (SE)	CI	b (SE)	CI	b (SE)	CI	b (SE)	CI	b (SE)	CI	b (SE)	CI
(Intercept)	113.08 (1.87)	109.41 – 116.74	113.08 (1.87)	109.41 – 116.74	113.11 (1.86)	109.45 – 116.77	72.72 (1.45)	69.87 – 75.57	72.72 (1.45)	69.87 – 75.57	72.72 (1.45)	69.87 – 75.57
gender [M]	2.48 (2.76)	-2.94 – 7.91	2.48 (2.76)	-2.94 – 7.91	2.48 (2.76)	-2.94 – 7.90	-1.31 (2.15)	-5.53 – 2.90	-1.31 (2.15)	-5.53 – 2.90	-1.31 (2.15)	-5.53 – 2.90
age gmc	0.33 (0.11)	0.11 – 0.56	0.33 (0.11)	0.11 – 0.56	0.33 (0.11)	0.11 – 0.56	0.27 (0.09)	0.10 – 0.45	0.27 (0.09)	0.10 – 0.45	0.27 (0.09)	0.10 – 0.45
BMI gmc	1.40 (0.41)	0.60 – 2.20	1.40 (0.41)	0.60 – 2.20	1.40 (0.41)	0.60 – 2.20	1.16 (0.32)	0.54 – 1.78	1.16 (0.32)	0.54 – 1.78	1.16 (0.32)	0.54 – 1.78
PD mc	-0.82 (0.26)	-1.34 – -0.30	-0.81 (0.27)	-1.33 – -0.28	-0.80 (0.27)	-1.32 – -0.27	-0.23 (0.23)	-0.67 – 0.22	-0.18 (0.23)	-0.63 – 0.26	-0.18 (0.23)	-0.63 – 0.27
WHLSM cm gmc	1.34 (1.06)	-0.75 – 3.43	1.34 (1.06)	-0.75 – 3.43	1.35 (1.06)	-0.74 – 3.44	0.46 (0.83)	-1.17 – 2.08	0.46 (0.83)	-1.17 – 2.08	0.46 (0.83)	-1.17 – 2.08
WHLSM mc			0.31 (0.40)	-0.48 – 1.10	0.31 (0.40)	-0.48 – 1.10			0.67 (0.34)	-0.00 – 1.35	0.67 (0.35)	-0.00 – 1.35
PD mc × WHLSM mc					0.24 (0.30)	-0.35 – 0.84					0.00 (0.26)	-0.51 – 0.51
Random Effects
σ2	96.62		96.69		96.77		71.29		70.97		71.08
τ00	173.03 ID		173.02 ID		172.72 ID		102.38 ID		102.42 ID		102.43 ID
N	106 ID		106 ID		106 ID		106 ID		106 ID		106 ID
Observations	721		721		721		721		721		721

# bootstrap CI SBP_eve
confint.merMod(m1_SBP_eve,parm=3:(length(fixef(m1_SBP_eve))+2),method="boot",nsim=NSIM) # M1

##                    2.5 %      97.5 %
## (Intercept)  109.3702872 116.7234842
## genderM       -2.9454330   7.8981746
## age.gmc        0.1097913   0.5571745
## BMI.gmc        0.5999890   2.2064248
## PD.mc         -1.3406626  -0.2902520
## WHLSM.cm.gmc  -0.7206930   3.4018219

confint.merMod(m2_SBP_eve,parm=3:(length(fixef(m2_SBP_eve))+2),method="boot",nsim=NSIM) # M2

##                    2.5 %      97.5 %
## (Intercept)  109.4134076 116.7476503
## genderM       -3.0751138   7.8695452
## age.gmc        0.1075076   0.5581521
## BMI.gmc        0.6171424   2.2154539
## PD.mc         -1.3317888  -0.2849928
## WHLSM.cm.gmc  -0.7500476   3.4051204
## WHLSM.mc      -0.4803054   1.1170363

confint.merMod(m3_SBP_eve,parm=3:(length(fixef(m3_SBP_eve))+2),method="boot",nsim=NSIM) # M3

##                      2.5 %      97.5 %
## (Intercept)    109.4566260 116.7756009
## genderM         -2.9553584   7.9669292
## age.gmc          0.1057149   0.5509935
## BMI.gmc          0.6006390   2.2148860
## PD.mc           -1.3108674  -0.2817490
## WHLSM.cm.gmc    -0.7116954   3.4393428
## WHLSM.mc        -0.4875876   1.0955268
## PD.mc:WHLSM.mc  -0.3694828   0.8242395

# bootstrap CI DBP_eve
confint.merMod(m1_DBP_eve,parm=3:(length(fixef(m1_DBP_eve))+2),method="boot",nsim=NSIM) # M1

##                    2.5 %     97.5 %
## (Intercept)  69.91707201 75.5146361
## genderM      -5.50157954  2.7837334
## age.gmc       0.09630789  0.4471138
## BMI.gmc       0.52973629  1.7700847
## PD.mc        -0.67887995  0.2252417
## WHLSM.cm.gmc -1.14713308  2.1283198

confint.merMod(m2_DBP_eve,parm=3:(length(fixef(m2_DBP_eve))+2),method="boot",nsim=NSIM) # M2

##                     2.5 %     97.5 %
## (Intercept)  69.899178002 75.5317279
## genderM      -5.468736411  2.9722360
## age.gmc       0.093351817  0.4415851
## BMI.gmc       0.544397214  1.7757014
## PD.mc        -0.629106630  0.2678128
## WHLSM.cm.gmc -1.144785651  2.0667171
## WHLSM.mc      0.005499441  1.3510498

confint.merMod(m3_DBP_eve,parm=3:(length(fixef(m3_DBP_eve))+2),method="boot",nsim=NSIM) # M3

##                       2.5 %     97.5 %
## (Intercept)    69.920946995 75.5180192
## genderM        -5.574241061  2.8625886
## age.gmc         0.091355526  0.4513621
## BMI.gmc         0.535751402  1.7784338
## PD.mc          -0.635938890  0.2680905
## WHLSM.cm.gmc   -1.160080851  2.1024836
## WHLSM.mc        0.004045224  1.3491670
## PD.mc:WHLSM.mc -0.501603579  0.5085815

# coefficients EE and SD
tab_model(m1_EE,m2_EE,m3_EE,m1_SD,m2_SD,m3_SD,
          dv.labels=paste0(rep(c("EE_","SD_"),each=3),c("baseline","WHLSM.mc","interaction")),
          show.icc=FALSE,show.p=FALSE,show.se=TRUE,show.r2=FALSE,collapse.se=TRUE,string.est="b (SE)")

	EE_baseline		EE_WHLSM.mc		EE_interaction		SD_baseline		SD_WHLSM.mc		SD_interaction
Predictors	b (SE)	CI	b (SE)	CI	b (SE)	CI	b (SE)	CI	b (SE)	CI	b (SE)	CI
(Intercept)	3.12 (0.14)	2.84 – 3.40	3.12 (0.14)	2.84 – 3.40	3.12 (0.14)	2.84 – 3.40	2.75 (0.12)	2.53 – 2.98	2.75 (0.12)	2.53 – 2.98	2.75 (0.12)	2.52 – 2.97
gender [M]	0.18 (0.21)	-0.22 – 0.58	0.18 (0.21)	-0.22 – 0.58	0.18 (0.21)	-0.22 – 0.59	-0.40 (0.17)	-0.72 – -0.07	-0.40 (0.17)	-0.72 – -0.07	-0.39 (0.17)	-0.72 – -0.06
PD mc	-0.08 (0.03)	-0.14 – -0.03	-0.07 (0.03)	-0.13 – -0.02	-0.07 (0.03)	-0.13 – -0.02	-0.03 (0.03)	-0.08 – 0.03	-0.02 (0.03)	-0.08 – 0.04	-0.02 (0.03)	-0.08 – 0.03
WHLSM cm gmc	0.47 (0.08)	0.31 – 0.62	0.47 (0.08)	0.31 – 0.62	0.47 (0.08)	0.31 – 0.62	0.34 (0.06)	0.22 – 0.47	0.34 (0.06)	0.22 – 0.47	0.34 (0.06)	0.21 – 0.47
WHLSM mc			0.19 (0.04)	0.11 – 0.27	0.19 (0.04)	0.11 – 0.27			0.11 (0.04)	0.03 – 0.20	0.11 (0.04)	0.03 – 0.20
PD mc × WHLSM mc					-0.04 (0.03)	-0.10 – 0.02					-0.07 (0.03)	-0.14 – -0.01
Random Effects
σ2	1.09		1.05		1.05		1.13		1.12		1.11
τ00	1.02 ID		1.03 ID		1.03 ID		0.60 ID		0.60 ID		0.61 ID
N	114 ID		114 ID		114 ID		114 ID		114 ID		114 ID
Observations	779		779		779		723		723		723

 # bootstrap CI EE 
confint.merMod(m1_EE,parm=3:(length(fixef(m1_EE))+2),method="boot",nsim=NSIM) # M1

##                   2.5 %      97.5 %
## (Intercept)   2.8422610  3.39941438
## genderM      -0.2250781  0.58615456
## PD.mc        -0.1378132 -0.02884978
## WHLSM.cm.gmc  0.3184584  0.61927262

confint.merMod(m2_EE,parm=3:(length(fixef(m2_EE))+2),method="boot",nsim=NSIM) # M2

##                   2.5 %      97.5 %
## (Intercept)   2.8507944  3.39588003
## genderM      -0.2193533  0.57305636
## PD.mc        -0.1266091 -0.02031197
## WHLSM.cm.gmc  0.3144118  0.62093079
## WHLSM.mc      0.1106096  0.27035420

confint.merMod(m3_EE,parm=3:(length(fixef(m3_EE))+2),method="boot",nsim=NSIM) # M3

##                     2.5 %      97.5 %
## (Intercept)     2.8379755  3.39360531
## genderM        -0.2182741  0.58687388
## PD.mc          -0.1279269 -0.02116087
## WHLSM.cm.gmc    0.3112911  0.62112684
## WHLSM.mc        0.1107532  0.27214145
## PD.mc:WHLSM.mc -0.1006700  0.02035031

# boostrap CI SD
confint.merMod(m1_SD,parm=3:(length(fixef(m1_SD))+2),method="boot",nsim=NSIM) # M1

##                    2.5 %      97.5 %
## (Intercept)   2.52810605  2.97864335
## genderM      -0.71807474 -0.06310848
## PD.mc        -0.08342478  0.03167821
## WHLSM.cm.gmc  0.21660335  0.47118734

confint.merMod(m2_SD,parm=3:(length(fixef(m2_SD))+2),method="boot",nsim=NSIM) # M2

##                    2.5 %      97.5 %
## (Intercept)   2.52515465  2.98227089
## genderM      -0.71688570 -0.06488114
## PD.mc        -0.07701252  0.03833502
## WHLSM.cm.gmc  0.21358565  0.47199962
## WHLSM.mc      0.02946134  0.19999122

confint.merMod(m3_SD,parm=3:(length(fixef(m3_SD))+2),method="boot",nsim=NSIM) # M3

##                      2.5 %      97.5 %
## (Intercept)     2.52022788  2.97743736
## genderM        -0.71942205 -0.06867787
## PD.mc          -0.07861727  0.03433758
## WHLSM.cm.gmc    0.21444106  0.46769987
## WHLSM.mc        0.02671827  0.19504759
## PD.mc:WHLSM.mc -0.13730658 -0.01152252

# plotting interaction
library(ggplot2); library(gridExtra)
sd(cleanSD$PD.mc) # RDet: 1 SD = 1.36

## [1] 1.362052

p <- plot_model(m3_SD,type="pred",terms=c("WHLSM.mc","PD.mc [-1.36,1.36]"),colors="bw",
           alpha=0.4,legend.title="Psychological\ndetachment",axis.title=c("State workaholism","Sleep disturbances")) +
  scale_color_manual(labels=c("-1 SD","+1 SD"),values=c("black","#666666")) +
  scale_linetype_manual(labels=c("-1 SD","+1 SD"),values=c("solid","dashed")) +
  scale_fill_manual(labels=c("-1 SD","+1 SD"),values=c("black","#666666")) + ggtitle("") +
                    theme(text=element_text(size=15))
p

ggsave("RESULTS/Figure3.tiff",plot=p,dpi=300,width=5,height=4)

References

• Steegen, S., Tuerlinckx, F., Gelman, A., & Vanpaemel, W. (2016). Increasing transparency through a multiverse analysis. Perspectives on Psychological Science, 11(5), 702-712. https://doi.org/10.1177/1745691616658637

R packages

Auguie, Baptiste. 2017. gridExtra: Miscellaneous Functions for "Grid" Graphics. https://CRAN.R-project.org/package=gridExtra.

Bartoń, Kamil. 2023. MuMIn: Multi-Model Inference. https://CRAN.R-project.org/package=MuMIn.

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.

Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2023. Lme4: Linear Mixed-Effects Models Using Eigen and S4. https://github.com/lme4/lme4/.

Delignette-Muller, Marie Laure, and Christophe Dutang. 2015. “fitdistrplus: An R Package for Fitting Distributions.” Journal of Statistical Software 64 (4): 1–34. https://doi.org/10.18637/jss.v064.i04.

Delignette-Muller, Marie-Laure, Christophe Dutang, and Aurelie Siberchicot. 2023. Fitdistrplus: Help to Fit of a Parametric Distribution to Non-Censored or Censored Data. https://CRAN.R-project.org/package=fitdistrplus.

Fox, John, and Sanford Weisberg. 2019. An R Companion to Applied Regression. Third. Thousand Oaks CA: Sage. https://socialsciences.mcmaster.ca/jfox/Books/Companion/.

Fox, John, Sanford Weisberg, and Brad Price. 2023. Car: Companion to Applied Regression. https://CRAN.R-project.org/package=car.

Imai, Kosuke, Luke Keele, and Dustin Tingley. 2010. “A General Approach to Causal Mediation Analysis.” Psychological Methods 15 (4): 309–34. http://imai.princeton.edu/research/BaronKenny.html.

Imai, Kosuke, Luke Keele, Dustin Tingley, and Teppei Yamamoto. 2010. “Causal Mediation Analysis Using r.” In Advances in Social Science Research Using r, edited by H. D. Vinod. New York: Springer-Verlag.

———. 2011. “Unpacking the Black Box of Causality: Learning about Causal Mechanisms from Experimental and Observational Studies.” American Political Science Review 105 (4): 765–89. http://imai.princeton.edu/research/mediationP.html.

Imai, Kosuke, Luke Keele, and Teppei Yamamoto. 2010. “Identification, Inference, and Sensitivity Analysis for Causal Mediation Effects.” Statistical Science 25 (1): 51–71. http://imai.princeton.edu/research/mediation.html.

Imai, Kosuke, and Teppei Yamamoto. 2013. “Identification and Sensitivity Analysis for Multiple Causal Mechanisms: Revisiting Evidence from Framing Experiments.” Political Analysis 21 (2): 141–71. http://imai.princeton.edu/research/medsens.html.

Lüdecke, Daniel. 2023. sjPlot: Data Visualization for Statistics in Social Science. https://strengejacke.github.io/sjPlot/.

Nieuwenhuis, Rense, Ben Pelzer, and Manfred te Grotenhuis. 2017. Influence.ME: Tools for Detecting Influential Data in Mixed Effects Models. http://www.rensenieuwenhuis.nl/r-project/influenceme/.

Nieuwenhuis, Rense, Manfred Te Grotenhuis, and Ben Pelzer. 2012. “Influence.ME: Tools for Detecting Influential Data in Mixed Effects Models.” R Journal 4 (2): 38–47.

R Core Team. 2023. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Revelle, William. 2023. Psych: Procedures for Psychological, Psychometric, and Personality Research. https://personality-project.org/r/psych/ https://personality-project.org/r/psych-manual.pdf.

Tingley, Dustin, Teppei Yamamoto, Kentaro Hirose, Luke Keele, and Kosuke Imai. 2014. “mediation: R Package for Causal Mediation Analysis.” Journal of Statistical Software 59 (5): 1–38. http://www.jstatsoft.org/v59/i05/.

Tingley, Dustin, Teppei Yamamoto, Kentaro Hirose, Luke Keele, Kosuke Imai, Minh Trinh, and Weihuang Wong. 2019. Mediation: Causal Mediation Analysis. https://imai.princeton.edu/projects/mechanisms.html.

Wickham, Hadley. 2011. “The Split-Apply-Combine Strategy for Data Analysis.” Journal of Statistical Software 40 (1): 1–29. https://www.jstatsoft.org/v40/i01/.

———. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

———. 2022. Plyr: Tools for Splitting, Applying and Combining Data. https://CRAN.R-project.org/package=plyr.

Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, Hiroaki Yutani, and Dewey Dunnington. 2023. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.

Xie, Yihui. 2014. “Knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2023. Knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.