The daily costs of workaholism

Supplementary material S5: Descriptive statistics and covariate selection

Aims and content

The present document includes the analytical steps implemented to compute the descriptive statistics from the retrospective and 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. The document also depicts the procedure used to select the covariates to be included in the following multilevel models (see Supplementary Material S6).

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("knitr","Rmisc","Hmisc","lme4","ggplot2","reshape2","MuMIn")

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

## tweaking Hmisc

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

1. Data reading

First, we read daily diary dataset diary and we rebuild the preliminary questionnaire prelqs exported from the previous step (Supplementary Material S4).

# loading data
load("DATI/diary_aggregated.RData") 

# rebuilding the prelqs dataset (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 (days) from",nlevels(diary$ID),"participants")

## diary: 1084 responses (days) 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

prelqs <- prelqs[prelqs$ID %in% levels(clean$ID),] # filtering prelqs data
prelqs$ID <- as.factor(as.character(prelqs$ID)) # resetting ID levels
cat("prelqs:",nrow(prelqs),"responses from",nlevels(prelqs$ID),"participants")

## prelqs: 114 responses from 114 participants

3. Descriptives

First, we describe the sample and compute the descriptive statistics and correlations of the core study variables (i.e., workaholism, WHLSM; afternoon, evening, and morning blood pressure, SBP_aft, DBP_aft, SBP_eve DBP_eve; emotional exhaustion, EE; sleep disturbances, SD, and recovery experiences, RDet and RRel) and a set of potential quantitative covariates (i.e., workhours, workHours; time since waking-up, hhFromAwake; age, age; BMI, BMI, and weekly working hours, weekHours).

The following packages and functions are used to optimize the analysis.

library(knitr); library(ggplot2); library(reshape2)

multidesc()

#' @title Descriptive statistics (mean and SD) and intraclass correlation of multilevel datasets
#' @param long = long-form data.frame including all time-varying variables with one row per observation
#' @param wide = wide-form data.frame including all time-invariant variables with one row per subject
#' @param cluster = character string indicating the name of the column with the subject identifiers
#' @param lv1 = character vector indicating the column names of the time-varying variables in the long dataset
#' @param lv2 = character vector indicating the column names of the time-invariant variables in the wide dataset
#' @param group = character string indicating the name of the column of the grouping variable (defult: NA for no group specification)
#' @param group.labels = character vector indicating the name of the group variable levels (default: NA)
multidesc <- function(long,wide,cluster="ID",lv1,lv2){ require(Rmisc); require(lme4)
  
  # lv1 data
  out <- data.frame(Measure=lv1[1],
                    N=summarySE(long,lv1[1],na.rm=TRUE)[,2],
                    Mean=paste(round(summarySE(long,lv1[1],na.rm=TRUE)[,3],2)," (",
                               round(summarySE(long,lv1[1],na.rm=TRUE)[,4],2),")",sep=""))
  for(i in 2:length(lv1)){
    out <- rbind(out,
                 data.frame(Measure=lv1[i],
                            N=summarySE(long,lv1[i],na.rm=TRUE)[,2],
                            Mean=paste(round(summarySE(long,lv1[i],na.rm=TRUE)[,3],2)," (",
                                       round(summarySE(long,lv1[i],na.rm=TRUE)[,4],2),")",sep="")))}
  
  # lv2 data
  if(!is.na(lv2[1])){
    for(i in 1:length(lv2)){
      out <- rbind(out,
                   data.frame(Measure=lv2[i],
                              N=summarySE(wide,lv2[i],na.rm=TRUE)[,2],
                              Mean=paste(round(summarySE(wide,lv2[i],na.rm=TRUE)[,3],2)," (",
                                         round(summarySE(wide,lv2[i],na.rm=TRUE)[,4],2),")",sep="")))}}
  
  # ICC
  out$ICC <- NA
  for(i in 1:length(lv1)){
    m <- lmer(formula=gsub("var",lv1[i],gsub("ID",cluster,"var~(1|ID)")),data=long) # VAR_between / (VAR_between + VAR_within)
    out[out$Measure==lv1[i],"ICC"] <- round(as.data.frame(VarCorr(m))[1,4]/
                                              (as.data.frame(VarCorr(m))[1,4]+as.data.frame(VarCorr(m))[2,4]),2)}
  rownames(out) <- gsub(".cm","",rownames(out))
  return(out)}

Computes and prints descriptive statistics (mean and SD) and intraclass correlation of a multilevel dataset.

multicorr()

#' @title Correlation matrix of multilevel datasets
#' @param long = long-form data.frame including all time-varying variables with one row per observation
#' @param wide = wide-form data.frame including all time-invariant variables with one row per subject
#' @param cluster = character string indicating the name of the column with the subject identifiers
#' @param lv1 = character vector indicating the column names of the time-varying variables in the long dataset
#' @param lv2 = character vector indicating the column names of the time-invariant variables in the wide dataset
#' @param pvalue = logical value specifying whether p-values should be reported or not (default: FALSE)
#' @param correction = character specifying the type of correction to be applied to the computed p-values (see ?p.adjust), default is "bonferroni"
multicorr <- function(long,wide,lv1,lv2,cluster="ID",pvalue=FALSE,correction="bonferroni"){ 
  require(Rmisc); require(Hmisc); require(plyr)
  
  # renaming cluster as "ID"
  colnames(long)[which(colnames(long)==cluster)] <- "ID"
  
  # person-mean-centering lv-1 variables
  for(Var in lv1){
    wide <- cbind(wide,aggregate(long[,Var],list(long$ID),mean,na.rm=TRUE)[,2]) # computing individual means
    colnames(wide)[ncol(wide)] <- paste0(Var,".cm")
    long <- join(long,wide[,c("ID",paste0(Var,".cm"))],by="ID",type="left") # joining with long-form data
    long[,paste(Var,".mc")] <- long[,Var] - long[,paste0(Var,".cm")] # computing mean-centered scores
    colnames(long)[ncol(long)] <- paste0(Var,".mc") }
  
  # between-subjects correlations (lv2)
  if(!is.na(lv2[1])){
    out <- rcorr(as.matrix(wide[,c(paste0(lv1,".cm"),lv2)]), type = "pearson")
  } else { out <- rcorr(as.matrix(wide[,paste0(lv1,".cm")]), type = "pearson") }
  rb <- round(out$r,2) # corr coefficients
  rb[lower.tri(rb)] <- NA
  rownames(rb) <- gsub(".cm","",rownames(rb))
  colnames(rb) <- gsub(".cm","",colnames(rb))
  rb.p <- out$P # pvalues
  rb.p[lower.tri(rb.p)] <- NA
  
  # within-participant correlations (lv1)
  out <- rcorr(as.matrix(long[,paste(lv1,".mc",sep="")]), type = "pearson")
  rw <- round(out$r,2) # corr coeff.
  rw[upper.tri(rw)] <- NA
  rownames(rw) <- gsub(".mc","",rownames(rw))
  colnames(rw) <- gsub(".mc","",colnames(rw))
  rw.p <- out$P # pvalues
  rw.p[upper.tri(rw.p)] <- NA
  
  # adjusting p-values
  if(pvalue==TRUE){
    adj <- data.frame(original=c(rb.p[!is.na(rb.p)],rw.p[!is.na(rw.p)]),
                      sig=rep(NA,length(c(rb.p[!is.na(rb.p)],rw.p[!is.na(rw.p)]))))
    adj$level <- "w"
    adj[1:length(rb.p[!is.na(rb.p)]),"level"] <- "b"
    adj$adjusted=p.adjust(adj$original,method=correction)
    # significance level
    adj[adj$adjusted<.05,"sig"] <- "*"
    adj[adj$adjusted<.01,"sig"] <- "**"
    adj[adj$adjusted<.001,"sig"] <- "***"
    adj[is.na(adj$sig),"sig"] <- ""
    adjb <- adj[adj$level=="b",]
    adjw <- adj[adj$level=="w",]
    for(i in 1:nrow(adjb)){ rb.p[!is.na(rb.p) & rb.p==adjb[i,"original"]] <- adjb[i,"sig"] }
    for(i in 1:nrow(adjw)){ rw.p[!is.na(rw.p) & rw.p==adjw[i,"original"]] <- adjw[i,"sig"] }
    for(i in 1:nrow(rb.p)){ for(j in 1:ncol(rb.p)){ 
      if(!is.na(rb.p[i,j])){ rb[i,j] <- paste(round(as.numeric(rb[i,j]),2),rb.p[i,j],sep="")}}}
    for(i in 1:nrow(rw.p)){ for(j in 1:ncol(rw.p)){ 
      if(!is.na(rw.p[i,j])){ rw[i,j] <- paste(round(as.numeric(rw[i,j]),2),rw.p[i,j],sep="")}}}}
  
  # filling rb empty cells
  rb[1:length(lv1),1:length(lv1)][lower.tri(rb[1:length(lv1),1:length(lv1)])] <- rw[lower.tri(rw)]
  rb <- t(rb) # transposing the matrix (correlations within-individual are shown above the main diagonal)
  return(rb)}

computes the correlation matrix of a multilevel dataset.

3.1. Sample description

Here, we describe the sample by providing frequencies and other descriptive statistics of sociodemographic and work-related variables collected with the preliminary questionnaire, considering the included subsample.

# gender (freq; %)
summary(prelqs$gender); round(100*summary(prelqs$gender)/nrow(prelqs),1)

##  F  M 
## 59 55

##    F    M 
## 51.8 48.2

# age (mean, SD)
round(mean(prelqs$age),2); round(sd(prelqs$age),2)

## [1] 42.12

## [1] 12.61

# BMI (mean, SD)
round(mean(prelqs$BMI),2); round(sd(prelqs$BMI),2)

## [1] 23.91

## [1] 3.53

# education (freq; %)
summary(prelqs$edu); round(100*summary(prelqs$edu)/nrow(prelqs),1)

## highschool- university+ 
##          45          69

## highschool- university+ 
##        39.5        60.5

# job sector (freq; %)
summary(prelqs$sector); round(100*summary(prelqs$sector)/nrow(prelqs),1)

## Private  Public 
##      91      23

## Private  Public 
##    79.8    20.2

# job position (freq; %)
summary(prelqs$position); round(100*summary(prelqs$position)/nrow(prelqs),1)

##          Employee Manager/Employers           Project 
##                57                53                 4

##          Employee Manager/Employers           Project 
##              50.0              46.5               3.5

# occupational groups (%)
jobs <- round(100*summary(prelqs$job)/nrow(prelqs),1)
jobs[order(jobs,decreasing=TRUE)]

##         Business and administration associate professionals 
##                                                        15.8 
##                   Business and administration professionals 
##                                                        13.2 
##                    Legal, social and cultural professionals 
##                                                        10.5 
##                       Science and engineering professionals 
##                                                        10.5 
##                                      Teaching professionals 
##                                                         9.6 
##             Science and engineering associate professionals 
##                                                         8.8 
##                                        Health professionals 
##                                                         5.3 
##                Production and specialised services managers 
##                                                         3.5 
##          Chief executives, senior officials and legislators 
##                                                         2.6 
##                              Health associate professionals 
##                                                         2.6 
##             Hospitality, retail and other services managers 
##                                                         2.6 
##     Information and communications technology professionals 
##                                                         2.6 
##                      Administrative and commercial managers 
##                                                         1.8 
##                                    Personal service workers 
##                                                         1.8 
##                                               Sales workers 
##                                                         1.8 
##                      Stationary plant and machine operators 
##                                                         1.8 
##                                                  Assemblers 
##                                                         0.9 
##                                    Customer services clerks 
##                                                         0.9 
##                    Electrical and electronic trades workers 
##                                                         0.9 
##                             Handicraft and printing workers 
##                                                         0.9 
## Legal, social, cultural and related associate professionals 
##                                                         0.9 
##                 Metal, machinery and related trades workers 
##                                                         0.9

# weekly working time (mean; SD)
round(mean(prelqs$weekHours),2); round(sd(prelqs$weekHours),2)

## [1] 41.22

## [1] 9.76

# having a partner (freq; %)
summary(prelqs$partner); round(100*summary(prelqs$partner)/nrow(prelqs),1)

##  No Yes 
##  61  53

##   No  Yes 
## 53.5 46.5

# having a children (freq; %)
summary(prelqs$children); round(100*summary(prelqs$children)/nrow(prelqs),1)

##  No Yes 
##  56  58

##   No  Yes 
## 49.1 50.9

# living arrangements (freq; %)
summary(prelqs$home); round(100*summary(prelqs$home)/nrow(prelqs),1)

##    alone  partner children  parents   others 
##       13       28       45       19        9

##    alone  partner children  parents   others 
##     11.4     24.6     39.5     16.7      7.9

# smoking status (freq; %)
summary(prelqs$smoker); round(100*summary(prelqs$smoker)/nrow(prelqs),1)

##  No Yes 
##  83  31

##   No  Yes 
## 72.8 27.2

3.2. Response rate

Here, we report the total number of observations and response rate for each variable, considering the included subsample.

# total number of diary entries, days, and participants
cat(nrow(clean[clean$aft==1,])+nrow(clean[clean$eve==1,])+nrow(clean[clean$mor==1,]),"diary entry,",
    nrow(clean),"daily observations from",nlevels(clean$ID),"participants")

## 2534 diary entry, 958 daily observations from 114 participants

# computing Nweeks, Nobs, and complDays variables
prelqs$Nweeks <- prelqs$complObs <- prelqs$complDays <- 1
for(i in 1:nrow(prelqs)){
  if(max(clean[clean$ID==prelqs[i,"ID"],"day"])>6){ prelqs[i,"Nweeks"] <- 2 } # No. of weeks
  prelqs[i,"complObs"] <- nrow(clean[clean$ID==prelqs[i,"ID"] & clean$aft==1,]) + # No. of complete observations
    nrow(clean[clean$ID==prelqs[i,"ID"] & clean$eve==1,]) + nrow(clean[clean$ID==prelqs[i,"ID"] & clean$mor==1,])
  prelqs[i,"complDays"] <- nrow(clean[clean$ID==prelqs[i,"ID"] & clean$aft==1 & clean$eve==1 & clean$aft==1,]) } # No. full days
cat("Sanity check:",min(prelqs$complDays)==3 & max(prelqs$Nweeks)==2) # min 3 complete days and max 2 weeks

## Sanity check: TRUE

# number and % of participants involved for 1 vs. 2 weeks
summary(as.factor(prelqs$Nweeks))

##  1  2 
## 31 83

round(100*summary(as.factor(prelqs$Nweeks))/nrow(prelqs),1)

##    1    2 
## 27.2 72.8

# computing respRate
prelqs[prelqs$Nweeks==1,"respRate"] <- round(100*prelqs[prelqs$Nweeks==1,"complObs"]/15,2) # max 15 obs when Nweeks = 1
prelqs[prelqs$Nweeks==2,"respRate"] <- round(100*prelqs[prelqs$Nweeks==2,"complObs"]/30,2) # max 30 obs when Nweeks = 2

# computing complDaysRate
prelqs[prelqs$Nweeks==1,"complDaysRate"] <- round(100*prelqs[prelqs$Nweeks==1,"complDays"]/5,2) # max 5 days when Nweeks = 1
prelqs[prelqs$Nweeks==2,"complDaysRate"] <- round(100*prelqs[prelqs$Nweeks==2,"complDays"]/10,2) # max 10 days when Nweeks = 2

# printing info
cat("- No. complete responses per participants:",round(mean(prelqs$complObs),2),", SD =",round(sd(prelqs$complObs),2),
    "\n- response rate per participants:",round(mean(prelqs$respRate),2),"%, SD =",round(sd(prelqs$respRate),2),
    "%\n- No. complete days (Aft + Eve + Mor) per participants",round(mean(prelqs$complDays),2),
    ", SD =",round(sd(prelqs$complDays),2),"\n- complete days rate per participants:",round(mean(prelqs$complDaysRate),2),
    "%, SD =",round(sd(prelqs$complDaysRate),2),"%","\n-",nrow(prelqs[prelqs$Nweeks==2,]),"on",nrow(prelqs),
    "participants involved for two weeks (",round(100*nrow(prelqs[prelqs$Nweeks==2,])/nrow(prelqs)),
    "% )\n- mean No. of complete responses for those participating over 2 weeks =",
    round(mean(prelqs[prelqs$Nweeks==2,"complObs"]),2),
    "\n- mean No. of complete responses for those participating over 1 week =",
    round(mean(prelqs[prelqs$Nweeks==1,"complObs"]),2))

## - No. complete responses per participants: 22.23 , SD = 6.09 
## - response rate per participants: 86.87 %, SD = 12.85 %
## - No. complete days (Aft + Eve + Mor) per participants 6.86 , SD = 2.23 
## - complete days rate per participants: 80.61 %, SD = 18.64 % 
## - 83 on 114 participants involved for two weeks ( 73 % )
## - mean No. of complete responses for those participating over 2 weeks = 25.27 
## - mean No. of complete responses for those participating over 1 week = 14.1

3.3. Blood pressure subsample

In addition to filtering participant based on their 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 describe such excluded participants and the subsample of participants included in blood pressure analyses.

# filtering participants with bp_drugs or cv_dysf
BPout <- prelqs[prelqs$bp_drugs=="Yes" | prelqs$cv_dysf=="Yes",]
BPin <- prelqs[prelqs$bp_drugs=="No" & prelqs$cv_dysf=="No",]

# sample size for blood pressure
cat(nrow(clean[clean$ID %in% levels(as.factor(as.character(BPin$ID))),]),"ovservations (days) from",
    nlevels(as.factor(as.character(BPin$ID))),"participants included in BP analyses")

## 898 ovservations (days) from 106 participants included in BP analyses

# sociodemographics in excluded participants
summary(BPout$gender); 100*summary(BPout$gender)/nrow(BPout) # gender (freq; %)

## F M 
## 3 5

##    F    M 
## 37.5 62.5

mean(BPout$age); sd(BPout$age) # age (mean, SD)

## [1] 56

## [1] 5.529144

mean(BPout$BMI); sd(BPout$BMI) # BMI (mean, SD)

## [1] 26.49418

## [1] 3.793512

# response rate in included participants
cat("- No. complete responses per participants:",round(mean(BPin$complObs),2),", SD =",round(sd(BPin$complObs),2),
    "\n- response rate per participants:",round(mean(BPin$respRate),2),"%, SD =",round(sd(BPin$respRate),2),
    "%\n- No. complete days (Aft + Eve + Mor) per participants",round(mean(BPin$complDays),2),
    ", SD =",round(sd(BPin$complDays),2),"\n- complete days rate per participants:",round(mean(BPin$complDaysRate),2),
    "%, SD =",round(sd(BPin$complDaysRate),2),"%","\n-",nrow(BPin[BPin$Nweeks==2,]),"on",nrow(BPin),
    "participants involved for two weeks (",round(100*nrow(BPin[BPin$Nweeks==2,])/nrow(BPin)),"% )")

## - No. complete responses per participants: 22.26 , SD = 5.98 
## - response rate per participants: 86.13 %, SD = 12.96 %
## - No. complete days (Aft + Eve + Mor) per participants 6.85 , SD = 2.19 
## - complete days rate per participants: 79.72 %, SD = 18.64 % 
## - 79 on 106 participants involved for two weeks ( 75 % )

3.4. Quantitative variables

Here, we use the multidesc() function for computing the descriptive statistics (mean, SD and ICC) of level-1 and level-2 variables. We can see that WHLSM shows an ICC of .60, indicating that most variance is located at the between-individual level, but within-individual variance is still substantial. Similar results are shown by blood pressure (ICC ranging from .61 to .73) and emotional exhaustion (ICC = .56), whereas higher proportions of within-individual variability are estimated for SD, RDet, and RRel.

# selecting level-1 quantitative variables
lv1 <- c("WHLSM","SBP_aft","DBP_aft", # afternoon
         "SBP_eve","DBP_eve","EE","PD", # evening
         "SD") # next morning

# selecting level-2 quantitative variables
lv2 <- c("age","BMI")

# computing descriptives
desc <- multidesc(long=clean,wide=prelqs,lv1=lv1,lv2=lv2)
kable(desc)

Measure	N	Mean	ICC
WHLSM	847	3.43 (1.54)	0.61
SBP_aft	843	122.18 (19.56)	0.70
DBP_aft	843	78.17 (14.6)	0.61
SBP_eve	841	115.58 (18.22)	0.69
DBP_eve	841	72.96 (14.26)	0.64
EE	841	3.24 (1.57)	0.56
PD	841	4.38 (1.85)	0.35
SD	842	2.57 (1.4)	0.39
age	114	42.12 (12.61)	NA
BMI	114	23.91 (3.53)	NA

3.4.1. Survey duration

Here, we inspect the mean, standard deviation, and distribution of diary survey durations. We can see that, after excluding a few outliers taking more than 30 minutes, the average response time to diary surveys was 4.02 minutes (SD = 3.51 min).

# computing durations in minutes
clean$dur_aft <- difftime(clean$end_aft,clean$start_aft,units="mins")
clean$dur_eve <- difftime(clean$end_eve,clean$start_eve,units="mins")
clean$dur_mor <- difftime(clean$end_mor,clean$start_mor,units="mins")
vars <- paste0("dur_",c("aft","eve","mor"))
clean[,vars] <- lapply(clean[,vars],as.numeric)
clean$dur_mean <- apply(clean[,c("dur_aft","dur_eve","dur_mor")],1,mean,na.rm=TRUE)

# plotting durations
par(mfrow=c(1,4)); for(Var in c(vars,"dur_mean")){ 
  hist(clean[,Var],breaks=100,main=paste(Var,"N =",length(na.omit(clean[,Var])))) }

# No. cases > 30 min (8 to 15)
for(Var in c(vars,"dur_mean")){ 
  colnames(clean)[which(colnames(clean)==Var)] <- "Var"
  cat(Var,"> 60 min:",nrow(clean[!is.na(clean$Var) & clean$Var > 30,]),"  ") 
  colnames(clean)[which(colnames(clean)=="Var")] <- Var }

## dur_aft > 60 min: 10   dur_eve > 60 min: 8   dur_mor > 60 min: 15   dur_mean > 60 min: 13

# mean and standard deviation (in minutes) without considering such outliers
round(mean(clean$dur_mean[clean$dur_mean<30],na.rm=TRUE),2)

## [1] 4.02

round(sd(clean$dur_mean[clean$dur_mean<30],na.rm=TRUE),2)

## [1] 3.51

# removing variables
clean[,c(vars,"dur_mean")] <- NULL

3.5. Correlations

Second, we use the multicorr function to compute the between-individual (shown below the main diagonal) and the within-individual correlations (shown above the main diagonal) among the same continuous variables.

cors <- multicorr(long=clean,wide=prelqs,lv1=lv1,lv2=lv2)
kable(cors)

	WHLSM	SBP_aft	DBP_aft	SBP_eve	DBP_eve	EE	PD	SD	age	BMI
WHLSM	1.00	0.14	0.14	0.04	0.08	0.18	-0.08	0.13	NA	NA
SBP_aft	0.12	1.00	0.46	0.18	0.12	0.05	-0.02	0.09	NA	NA
DBP_aft	0.13	0.83	1.00	0.10	0.12	0.08	-0.02	0.06	NA	NA
SBP_eve	0.11	0.90	0.77	1.00	0.56	0.06	-0.10	0.03	NA	NA
DBP_eve	0.04	0.77	0.86	0.87	1.00	0.08	-0.03	0.09	NA	NA
EE	0.49	0.24	0.22	0.21	0.17	1.00	-0.13	0.09	NA	NA
PD	-0.14	-0.01	-0.05	-0.07	-0.08	-0.13	1.00	-0.06	NA	NA
SD	0.43	0.09	0.15	0.12	0.13	0.47	-0.16	1.00	NA	NA
age	-0.14	0.35	0.35	0.34	0.38	-0.05	0.06	0.02	1.00	NA
BMI	0.06	0.33	0.33	0.40	0.40	0.06	0.07	-0.03	0.22	1

Correlations are also graphically visualized below (blue = negative, white = uncorrelated, red = positive).

ggplot(melt(as.matrix(multicorr(long=clean,wide=prelqs,lv1=lv1,lv2=lv2))),aes(x=Var1, y=Var2, fill=value)) + geom_tile() + 
  geom_text(aes(x = Var1, y = Var2, label = round(value,2)),color="black",size=5)+labs(x="",y="") +
  scale_fill_gradient2(low="darkblue",high="#f03b20",mid="white",midpoint=0,limit = c(-1,1), space = "Lab",
                       guide="legend",breaks=round(seq(1,-1,length.out = 11),2))+
  theme(text=element_text(face="bold",size=10),axis.text.x=element_text(angle=30))

3.6. Univariate distributions

Fourth, we visualize the univariate distributions of the focal variables, for each class of variables. We can see that level-1 self-reported predictors show quite skewed distributions. However, person-mean centered scores are approximately normally distributed.

# level-2
par(mfrow=c(3,6)); for(Var in c(lv2,"gender","edu","mStatus","home_partner","children","home_child","position","sector",
                                "smoker","careless","bp_drugs","horm_drugs","psy_drugs","cv_dysf","sleep_dysf")){
  if(class(clean[,Var])%in%c("factor","logical")){ plot(clean[!duplicated(clean$ID),Var],main=Var,xlab="") 
    } else { hist(clean[!duplicated(clean$ID),Var],main=Var,xlab="") }}

# level-1 (composite scores)
par(mfrow=c(3,6)); for(Var in c(lv1,"dailyHassles_eve","lateWorkHours","flagTime",
                                "where_aft","teleWork","confounders_aft","confounders_eve",
                                "flagBP_aft","flagBP_eve")){ 
  if(class(clean[,Var])%in%c("factor","logical")){ plot(clean[!duplicated(clean$ID),Var],main=Var,xlab="") 
    } else { hist(clean[!duplicated(clean$ID),Var],main=Var,xlab="") }}

# level-1 (person-mean-centered scores)
par(mfrow=c(2,4))

for(Var in lv1){ wide <- clean[!duplicated(clean$ID),]
  wide[,paste0(Var,".cm")] <- aggregate(clean[,Var],list(clean$ID),FUN=mean,na.rm=TRUE)[,2]
  clean <- plyr::join(clean,wide[,c("ID",paste0(Var,".cm"))],by="ID",type="left")
  clean[,paste0(Var,".mc")] <- clean[,Var] - clean[,paste0(Var,".cm")]
  hist(clean[,paste0(Var,".mc")],main=Var,xlab="") 
  clean <- clean[,1:(ncol(clean)-2)]}

3.7. Bivariate distributions

Fifth, we visualize the distribution of each focal variable (i.e., workaholism, blood pressure, emotional exhaustion, and sleep disturbances) and potential covariate against the following categorical variables:

1. Demographics: gender, educational level edu, marital status mStatus and living with a partner home_partner, having children children and living with children home_child, job position position and sector

2. Inclusion criteria: blood pressure bp_drugs, hormonal horm_drugs, and psychoactive medication psy_drugs, cardiovascular cv_dysf and sleep dysfunctions sleep_dysf

3. Level-2 confounders: being a smoker smoker or a careless respondent careless

4. Level-1 confounders: daily nonwork hassles in the evening dailyHassles_eve, presence of late working hours in the evening lateWorkHours, atypical response times flagTime

5. Level-1 covariates: survey location in the afternoon where_aft, remote working teleWork

6. BP confounders: aggregated confounders score confounders_aft and confounders_eve; coffee coffe_aft, coffe_eve; smoking smoke_aft, smoke_eve; vigorous physical activity sport_aft, sport_eve, having a meal meal_aft, meal_eve; and flagged BP cases flagBP_aft, flagBP_eve (see Supplementary Material S3)

The bivPlot function is used to optimize the data visualization.

bivPlot()

bivPlot <- function(data,resp,preds){ require(ggplot2)
 for(Var in preds){ 
   print(
     ggplot(na.omit(data[,c(resp,Var)]),aes_string(x=Var,y=resp,fill=Var)) + 
       geom_point(col="gray",position=position_jitter(width=.15)) + # data points
       geom_violin(alpha=.4) + geom_boxplot(width=0.3,alpha=0.2) + # violin & boxplot
       ggtitle(paste(resp,"by",Var)) + theme(legend.position="none") # title and removing legend
     )}}

DEMOS

WHLSM

Higher WHLSM can be observed in cases where children = "No" or home_child = "No" (with larger differences for the latter), where partner = "No", and where edu = "university+".

bivPlot(diary,resp="WHLSM",pred=c("gender","edu","children","home_child","partner","home_partner","position","sector"))

BP

Blood pressure is mainly affected by children, with SBP_aft, SBP_eve and morning blood pressure also showing differences by gender, while morning blood pressure is also affected by partner.

SBP_aft

Higher SBP_aft can be observed in cases where gender = "M", children = "Yes" or home_child = "Yes" (with larger differences for the former), with slightly higher SBP_aft in cases where partner = "Yes", where edu = "highschool-", and where position = "employer_mng".

bivPlot(diary,resp="SBP_aft",pred=c("gender","edu","children","home_child","partner","home_partner","position","sector"))

DBP_aft

Higher DBP_aft can be observed in cases where children = "Yes" or home_child = "Yes" (with larger differences for the former), with slightly higher DBP_aft in cases where partner = "Yes", and where sector = "public".

bivPlot(diary,resp="DBP_aft",pred=c("gender","edu","children","home_child","partner","home_partner","position","sector"))

SBP_eve

Higher SBP_eve can be observed in cases where gender = "M", children = "Yes" or home_child = "Yes" (with larger differences for the former), with slightly higher SBP_aft in cases where partner = "Yes", where edu = "highschool+", and where sector = "Public".

bivPlot(diary,resp="SBP_eve",pred=c("gender","edu","children","home_child","partner","home_partner","position","sector"))

DBP_eve

Slightly higher DBP_eve can be observed in cases where children = "Yes" or home_child = "Yes" (with larger differences for the former), were partner = "Yes", and where sector = "public".

bivPlot(diary,resp="DBP_eve",pred=c("gender","edu","children","home_child","partner","home_partner","position","sector"))

EE

Higher EE can be observed in cases where home_child = "No".

bivPlot(diary,resp="EE",pred=c("gender","edu","children","home_child","partner","home_partner","position","sector"))

SD

Higher SD can be observed in cases where gender = "F", where home_child = "Yes", where home_partner = "No", and where position = "employee_proj".

bivPlot(diary,resp="SD",pred=c("gender","edu","children","home_child","partner","home_partner","position","sector"))

INCLUSION

WHLSM

Higher WHLSM can be observed in cases where bp_drugs = "Yes".

bivPlot(diary,resp="WHLSM",pred=c("bp_drugs", "horm_drugs", "psy_drugs", "cv_dysf", "sleep_dysf"))

BP

Blood pressure is mainly affected by bp_drugs, and cv_dysf, which further corroborate the validity of our blood pressure measurements and the importance of our pre-registered exclusion criteria.

SBP_aft

Higher SBP_aft can be observed in cases where bp_drugs = "Yes" and cv_dysf = "Yes", with only a slightly lower SBP_aft in cases where horm_drugs = "No".

bivPlot(diary,resp="SBP_aft",pred=c("bp_drugs", "horm_drugs", "psy_drugs", "cv_dysf", "sleep_dysf"))

DBP_aft

Higher DBP_aft can be observed in cases where bp_drugs = "Yes" and cv_dysf = "Yes", with only a slightly lower DBP_aft in cases where horm_drugs = "No".

bivPlot(diary,resp="DBP_aft",pred=c("bp_drugs", "horm_drugs", "psy_drugs", "cv_dysf", "sleep_dysf"))

SBP_eve

Higher SBP_eve can be observed in cases where bp_drugs = "Yes" and cv_dysf = "Yes".

bivPlot(diary,resp="SBP_eve",pred=c("bp_drugs", "horm_drugs", "psy_drugs", "cv_dysf", "sleep_dysf"))

DBP_eve

Higher DBP_eve can be observed in cases where bp_drugs = "Yes" and cv_dysf = "Yes".

bivPlot(diary,resp="DBP_eve",pred=c("bp_drugs", "horm_drugs", "psy_drugs", "cv_dysf", "sleep_dysf"))

EE

Higher EE can be observed in cases where bp_drugs = "Yes" or psy_drugs = "Yes".

bivPlot(diary,resp="EE",pred=c("bp_drugs", "horm_drugs", "psy_drugs", "cv_dysf", "sleep_dysf"))

SD

SD does not substantially differ across inclusion criteria levels, with only slightly higher SD for those reporting bp_drugs = "Yes". Critically, SD does not substantially differ across sleep_dysf levels, corroborating our choice to not rely on this inclusion criterion (see Supplementary Material S7).

bivPlot(diary,resp="SD",pred=c("bp_drugs", "horm_drugs", "psy_drugs", "cv_dysf", "sleep_dysf"))

LV-2 CONF

WHLSM

WHLSM does not substantially differ across level-2 confounders levels, with only slightly higher WHLSM in those with smoker = "quit_less.

bivPlot(diary,resp="WHLSM",pred=c("smoker","careless"))

BP

Blood pressure is not substantially affected by level-2 confounders. Although it is lower in potentially careless participants, this might be due to the very small number of such cases.

SBP_aft

SBP_aft does not substantially differ across level-2 confounders levels. Although it is lower in potentially careless participants, this might be due to the very small number of such cases.

bivPlot(diary,resp="SBP_aft",pred=c("smoker","careless"))

DBP_aft

DBP_aft is slightly lower in smoker participants. Although it is lower in potentially careless participants, this might be due to the very small number of such cases.

bivPlot(diary,resp="DBP_aft",pred=c("smoker","careless"))

SBP_eve

SBP_eve does not substantially differ across level-2 confounders levels. Although it is lower in potentially careless participants, this might be due to the very small number of such cases.

bivPlot(diary,resp="SBP_eve",pred=c("smoker","careless"))

DBP_eve

DBP_eve is slightly lower in smoker participants. Although it is lower in potentially careless participants, this might be due to the very small number of such cases.

bivPlot(diary,resp="DBP_eve",pred=c("smoker","careless"))

EE

EE does not substantially differ across level-2 confounders levels.

bivPlot(diary,resp="EE",pred=c("smoker","careless"))

SD

SD is lower in individuals with smoker = "Yes".

bivPlot(diary,resp="SD",pred=c("smoker","careless"))

LV-1 CONF

WHLSM

Slightly higher WHLSM can be observed in cases where lateWorkHours = TRUE.

bivPlot(diary,resp="WHLSM",pred=c("dailyHassles_eve","lateWorkHours","flagTime"))

BP

Blood pressure is not substantially affected by level-1 confounders.

SBP_aft

SBP_aft does not substantially differ across level-1 confounders levels.

bivPlot(diary,resp="SBP_aft",pred=c("dailyHassles_eve","lateWorkHours","flagTime"))

DBP_aft

DBP_aft does not substantially differ across level-1 confounders levels, with only slightly higher DBP_aft in cases where lateWorkHours = TRUE.

bivPlot(diary,resp="DBP_aft",pred=c("dailyHassles_eve","lateWorkHours","flagTime"))

SBP_eve

SBP_eve does not substantially differ across level-1 confounders levels, with only slightly higher SBP_eve in cases where lateWorkHours = TRUE.

bivPlot(diary,resp="SBP_eve",pred=c("dailyHassles_eve","lateWorkHours","flagTime"))

DBP_eve

DBP_eve does not substantially differ across level-1 confounders levels, with only slightly lower DBP_eve in cases where flagTime = TRUE.

bivPlot(diary,resp="DBP_eve",pred=c("dailyHassles_eve","lateWorkHours","flagTime"))

EE

Lower EE can be observed in cases where dailyHassles_eve = "Yes"

bivPlot(diary,resp="EE",pred=c("dailyHassles_eve","lateWorkHours","flagTime"))

SD

Higher SD can be observed in cases where dailyHassles_eve = "Yes".

bivPlot(diary,resp="SD",pred=c("dailyHassles_eve","lateWorkHours","flagTime"))

LV-1 COV

WHLSM

Higher WHLSM can be observed in cases where where_aft = "workplace".

bivPlot(diary,resp="WHLSM",pred=c("where_aft","teleWork"))

BP

Blood pressure is not substantially affected by level-1 confounders, with only slightly lower SBP_aft in cases where teleWork = "Yes".

SBP_aft

SBP_aft is slightly lower in cases where teleWork = "teleWork".

bivPlot(diary,resp="SBP_aft",pred=c("where_aft","teleWork"))

DBP_aft

DBP_aft does not substantially differ across the levels of level-1 covariates. Although slightly higher DBP_aft is shown in cases where teleWork = "both", this might be due to the small number of such cases.

bivPlot(diary,resp="DBP_aft",pred=c("where_aft","teleWork"))

SBP_eve

SBP_eve does not substantially differ across the levels of level-1 covariates. Although slightly higher SBP_eve is shown in cases where teleWork = "both" or where where_aft = "other", this might be due to the small number of such cases.

bivPlot(diary,resp="SBP_eve",pred=c("where_aft","teleWork"))

DBP_eve

DBP_eve does not substantially differ across the levels of level-1 covariates. Although slightly higher DBP_eve is shown in cases where teleWork = "both", this might be due to the small number of such cases.

bivPlot(diary,resp="DBP_eve",pred=c("where_aft","teleWork"))

EE

Slightly higher EE can be observed in cases where teleWork = "both".

bivPlot(diary,resp="EE",pred=c("where_aft","teleWork"))

SD

Lower SD can be observed in cases where teleWork = "both".

bivPlot(diary,resp="SD",pred=c("where_aft","teleWork"))

BP CONF

WHLSM

Higher WHLSM can be observed in all cases where confounders = FALSE, with larger differences by confounders_eve. Differences by specific confounder are in line with those showed by the aggregated scores. Lower WHLSM can be observed in cases where flagBP_eve = TRUE.

bivPlot(diary,resp="WHLSM",pred=c("confounders_aft", "confounders_eve",
                                  "coffee_aft", "coffee_eve", "smoke_aft", "smoke_eve", 
                                  "sport_aft", "sport_eve", "meal_aft", "meal_eve", 
                                  "flagBP_aft", "flagBP_eve"))

BP

Only systolic blood pressure is slightly affected by BP confounders.

SBP_aft

Only afternoon confounders are considered for SBP_aft, which is only slightly higher in cases where confounders_aft = TRUE.

bivPlot(diary,resp="SBP_aft",pred=c("confounders_aft","coffee_aft","smoke_aft","sport_aft","meal_aft","flagBP_aft"))

DBP_aft

Only afternoon confounders are considered for DBP_aft, which is only slightly lower in cases where coffe_aft or smoke_aft = 1, slightly higher in cases where sport_aft or meal_aft = TRUE. Although higher DBP_aft is shown in cases where flagBP_aft = TRUE, this might be due to the small number of such cases.

bivPlot(diary,resp="DBP_aft",pred=c("confounders_aft","coffee_aft","smoke_aft","sport_aft","meal_aft","flagBP_aft"))

SBP_eve

Only evening confounders are considered for SBP_eve, which is slightly higher in cases where confounders_eve = TRUE, and higher in cases where flagBP_eve = TRUE.

bivPlot(diary,resp="SBP_eve",pred=c("confounders_eve","coffee_eve","smoke_eve","sport_eve","meal_eve","flagBP_eve"))

DBP_eve

Only evening confounders are considered for DBP_eve, which does not substantially differ across the levels of confounders_eve, whereas slightly higher DBP_eve can be observed in cases where flagBP_eve = TRUE.

bivPlot(diary,resp="DBP_eve",pred=c("confounders_eve","coffee_eve","smoke_eve","sport_eve","meal_eve","flagBP_eve"))

EE

EE does not substantially differ across blood pressure confounder levels, with only higher EE when sport_eve = "TRUE".

bivPlot(diary,resp="EE",pred=c("confounders_aft", "confounders_eve",
                               "coffee_aft", "coffee_eve", "smoke_aft", "smoke_eve",
                               "sport_aft", "sport_eve", "meal_aft", "meal_eve"))

SD

SD does not substantially differ across blood pressure confounder levels, with only a slightly lower SD when smoke_eve = "TRUE".

bivPlot(diary,resp="SD",pred=c("confounders_aft", "confounders_eve",
                               "coffee_aft", "coffee_eve", "smoke_aft", "smoke_eve",
                               "sport_aft", "sport_eve", "meal_aft", "meal_eve"))

4. Covariate selection

As pre-registered, for maximizing model parsimony while minimizing the risk of non-convergence, the covariates to be included in the following regression analyses are selected based on an empirical approach accounting for the following criteria:

1. The No. of cases in each level, prioritizing categorical covariates with balanced number of cases

2. The bivariate distributions and correlations with the response variable, prioritizing covariates showing some potential covariance trend with the response variable

3. Linear modeling, prioritizing the covariates whose inclusion as predictor of the corresponding response variable shows stronger evidence (i.e., higher Akaike weight) than more parsimonious models and absolute t-value of 2 or higher

However, deviating from the pre-registration (see Supplementary Material S7), we also consider theoretical arguments concerning some important covariates. Particularly, we prioritize level-1 over level-2 covariates as our study is focused on within-individual relationships. Moreover, following Spector (2021), we attempt to avoid risking to rule out variables acting as potential exploratory mechanisms (e.g., mediators) of the focused relationships. Theoretical arguments are detailed below for each response variable.

The following packages are used to optimize the analyses.

library(lme4); library(MuMIn)

4.1. Afternoon BP

For afternoon systolic and diastolic blood pressure, the pre-registered potential covariates to be considered are:

1. gender, in light of the robust literature on sex differences in blood pressure and hypertension (e.g., [Sandberg & Ji, 2012]), with gender having been highlighted as a potential moderator of the straining effects of trait workaholism (e.g., Bladucci et al., 2018; 2021).

2. age, in light of the consolidated relationship between aging and high blood pressure; however, no clear differences have been highlighted on the relationship between age and workaholism (Clark et al., 2016).

3. BMI, in light of the robust literature on the increased risk for hypertension and CVD in individuals with higher BMI (e.g., Gregg et al., 2005); however, no clear differences have been highlighted, to our knowledge, on the relationship between BMI and workaholism.

4. Participants’ smoking status smoker, in light of the widely acknowledged relationship between smoking and hypertension; however, no clear differences have been highlighted, to our knowledge, on the relationship between smoking status and workaholism.

5. The daily number of worked hours workHours, as long working hours is a definitional aspect of workaholism which has been found positively associated with higher blood pressure and other cardiovascular risks (Virtanen, 2012); however, controlling for this variable might risk to artificially deflate the relationship between workaholism and blood pressure since the former might influence the latter through working hours. That is, working hours is a potential explanatory mechanism of the focused relationship, and including it as a control variable might bias rather than strengthening the results (see Spector et al., 2021).

6. Remote work teleWork, although with no clear expectations on its relationship with blood pressure, while we originally hypothesized that remote working might act as a moderator of the straining effect of workaholism.

7. BP confounders_aft identifying behavioral factors (e.g., smoking, physical activity) that might affect blood pressure, although with no clear expectation on their relationship with workaholism.

Further pre-registered level-2 covariates (i.e., retrospective measures of trait workaholism, technostress creators) are not included due to our aim to focus the study on the focused level-1 variables, whereas recovery experiences will be included in all models. Moreover, we do not include the previous BP measure as we focus on the covariance of intra-individual fluctuations between workaholism and blood pressure, rather than the change of the latter over time (see Supplementary material S7). In contrast, additional (i.e., not pre-registered) level-2 covariates potentially moderating the relationship between workaholism and afternoon blood pressure are considered as potentially confounding factors.

Prior to evaluate the inclusion of each covariate, we prepare the data for the analyses. That is, we remove cases of missing responses in the response variable or any covariate (list-wise deletion) and we center level-1 continuous covariates on the individual mean (person-mean-centering). Moreover, we apply our inclusion criteria by removing all cases with bp_drugs or cv_dysf = "Yes".

# inclusion criteria
cleanBP <- clean[clean$bp_drugs=="No" & clean$cv_dysf=="No",]

# list-wise deletion
cleanBP <- as.data.frame(na.omit(cleanBP[,c("ID","day","SBP_aft","DBP_aft",
                                        "gender","age","BMI","smoker","workHours",
                                        "teleWork","confounders_aft",
                                        "home_child","children")]))
cat("Considering",nrow(cleanBP),"complete observations from",nlevels(as.factor(as.character(cleanBP$ID))),"participants")

## Considering 719 complete observations from 106 participants

# mean centering lv-1 continuous covariates
for(Var in c("workHours")){
  wide <- prelqs[prelqs$ID%in%levels(as.factor(as.character(cleanBP$ID))),]
  wide <- cbind(wide,aggregate(cleanBP[,Var],list(cleanBP$ID),mean)[,2]) # computing individual means
  colnames(wide)[ncol(wide)] <- paste0(Var,".cm")
  cleanBP <- join(cleanBP,wide[,c("ID",paste0(Var,".cm"))],by="ID",type="left") # joining with long-form data
  cleanBP[,paste0(Var,".mc")] <- cleanBP[,Var] - cleanBP[,paste0(Var,".cm")] } # computing mean-centered scores

4.1.1. gender

gender is a categorical variable with a roughly balanced number of participants. In section 3.5 we observed some tendency with higher SBP_aft in men than in women. Consistently, its inclusion in the null model shows improved evidence (i.e., higher Akaike weight) compared to the null model, although the estimated t-value is lower than 2, with a particularly low value for DBP_aft. In both cases, also considering its theoretical relevance as a potential moderator of the straining effects of workaholism (e.g., Balducci et al. (2018; 2021)) we include gender as a covariate.

SBP_aft

# number of cases by covariate level
summary(cleanBP[!duplicated(cleanBP$ID),"gender"]) # No. participants

##  F  M 
## 56 50

summary(cleanBP$gender) # No. cases

##   F   M 
## 389 330

# model comparison
m0 <- lmer(SBP_aft ~ (1|ID),data=cleanBP)
m1 <- lmer(SBP_aft ~ gender + (1|ID),data=cleanBP)
Weights(AIC(m0,m1)) # Akaike weights: stronger evidence including gender

##  model weights 
## [1] 0.059 0.941

summary(m1)$coefficients # coefficient: proximal to |t|=2

##               Estimate Std. Error   t value
## (Intercept) 118.249199   2.138441 55.296928
## genderM       5.785312   3.114885  1.857312

DBP_aft

# model comparison
m0d <- lmer(DBP_aft ~ (1|ID),data=cleanBP)
m1d <- lmer(DBP_aft ~ gender + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d)) # Akaike weights: stronger evidence including gender

##  model weights 
## [1] 0.317 0.683

summary(m1d)$coefficients # coefficient: lower than |t|=2

##               Estimate Std. Error   t value
## (Intercept) 76.7945450   1.598385 48.045090
## genderM      0.2649914   2.328735  0.113792

4.1.2. age

age is a level-2 continuous variable that in section 3.2 showed moderate-to-strong positive relationships with blood pressure (r from .42 to .52). Consistently, its inclusion in the model including gender implies stronger evidence than more parsimonious models, with an absolute t-value higher than 2. Thus, despite the uncertainty on its relationship with workaholism, we include age as a covariate for afternoon blood pressure.

SBP_aft

# model comparison vs null
m1.bis <- lmer(SBP_aft ~ age + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.041 0.959

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                Estimate Std. Error   t value
## (Intercept) 104.2508117  5.2271705 19.944024
## age           0.4070428  0.1218572  3.340327

# model comparison vs gender
m2 <- lmer(SBP_aft ~ gender + age + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.001 0.015 0.984

summary(m2)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 99.6056504  5.4739198 18.196403
## genderM      7.0313855  2.9637550  2.372458
## age          0.4394154  0.1199469  3.663417

DBP_aft

# model comparison vs null
m1.bisd <- lmer(DBP_aft ~ age + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.049 0.951

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                Estimate Std. Error   t value
## (Intercept) 104.2508117  5.2271705 19.944024
## age           0.4070428  0.1218572  3.340327

# model comparison vs gender
m2d <- lmer(DBP_aft ~ gender + age + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.021 0.044 0.935

summary(m2)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 99.6056504  5.4739198 18.196403
## genderM      7.0313855  2.9637550  2.372458
## age          0.4394154  0.1199469  3.663417

4.1.3. BMI

BMI is a level-2 continuous variable that in section 3.2 showed weak-to-moderate positive relationships with blood pressure (r from .28 to .43). Consistently, its inclusion in the model including gender and age implies stronger evidence than more parsimonious models, with an absolute t-value higher than 2. Thus, despite the uncertainty on its relationship with workaholism, we include BMI as a covariate for afternoon blood pressure.

SBP_aft

# model comparison vs null
m1.bis <- lmer(SBP_aft ~ BMI + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.008 0.992

summary(m1.bis)$coefficients # coefficient: |t| > 2

##              Estimate Std. Error  t value
## (Intercept) 85.026929 10.4387588 8.145310
## BMI          1.515511  0.4355548 3.479496

# model comparison vs previous model
m3 <- lmer(SBP_aft ~ gender + age + BMI + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.001 0.080 0.919

summary(m3)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error  t value
## (Intercept) 75.9904976 10.4542585 7.268856
## genderM      5.3399277  2.9528268 1.808412
## age          0.3781251  0.1189131 3.179845
## BMI          1.1352957  0.4326673 2.623946

DBP_aft

# model comparison vs null
m1.bisd <- lmer(DBP_aft ~ BMI + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.007 0.993

summary(m1.bisd)$coefficients # coefficient: |t| > 2

##              Estimate Std. Error  t value
## (Intercept) 49.600097  7.6470113 6.486207
## BMI          1.151826  0.3190604 3.610055

# model comparison vs previous model
m3d <- lmer(DBP_aft ~ gender + age + BMI + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.001 0.001 0.029 0.969

summary(m3d)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error    t value
## (Intercept) 42.7168148 7.78109757  5.4898187
## genderM     -0.3646649 2.19626487 -0.1660387
## age          0.2532776 0.08837419  2.8659682
## BMI          1.0105432 0.32190699  3.1392395

4.1.4. smoker

smoker is a level-2 categorical variable that in section 3.5 did not show substantially relationships with blood pressure. Consistently, both models show an absolute t-value lower than 2, although its inclusion implies stronger evidence than more parsimonious models. Also considering the uncertainty of its relationship with workaholism, we do not include it as a covariate for afternoon blood pressure.

SBP_aft

# model comparison vs null
m1.bis <- lmer(SBP_aft ~ smoker + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.155 0.845

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

##             Estimate Std. Error   t value
## (Intercept) 121.9556   1.844173 66.130213
## smokerYes    -3.5964   3.534608 -1.017482

# model comparison vs previous model
m4 <- lmer(SBP_aft ~ gender + age + BMI + smoker + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.013 0.148 0.838

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

##               Estimate Std. Error   t value
## (Intercept) 78.0843784 10.6162759  7.355157
## genderM      6.2564672  3.0660002  2.040596
## age          0.3664817  0.1192519  3.073173
## BMI          1.0912172  0.4340473  2.514051
## smokerYes   -3.6844451  3.3666439 -1.094397

DBP_aft

# model comparison vs null
m1.bisd <- lmer(DBP_aft ~ smoker + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.05 0.95

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

##             Estimate Std. Error   t value
## (Intercept) 121.9556   1.844173 66.130213
## smokerYes    -3.5964   3.534608 -1.017482

# model comparison vs previous model
m4d <- lmer(DBP_aft ~ gender + age + BMI + smoker + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.003 0.097 0.900

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

##               Estimate Std. Error    t value
## (Intercept) 45.0927640 7.84320044  5.7492811
## genderM      0.6816631 2.26441229  0.3010331
## age          0.2403700 0.08793307  2.7335567
## BMI          0.9596918 0.32054971  2.9938940
## smokerYes   -4.1730935 2.48797620 -1.6773044

4.1.5. WorkHours

workHours is a level-1 continuous variable that in section 3.2 was weakly correlated with afternoon blood pressure (r = .12-.04). Here, its inclusion implies stronger evidence (i.e., higher Akaike weight) than more parsimonious models, with an estimated absolute t-value higher than 2. However, we consider the theoretical importance of this variable, which might be a potential explanatory mechanism (e.g., mediator) through which workaholism might impact on blood pressure, and thus, we do not include it as a covariate for afternoon blood pressure.

SBP_aft

# model comparison vs null
m1.bis <- lmer(SBP_aft ~ workHours.mc + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.015 0.985

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                Estimate Std. Error   t value
## (Intercept)  120.976976  1.5735066 76.883679
## workHours.mc   1.039052  0.3140563  3.308489

# model comparison vs previous model
m4 <- lmer(SBP_aft ~ gender + age + BMI + workHours.mc + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.001 0.015 0.984

summary(m4)$coefficients # coefficient: |t| > 2

##                Estimate Std. Error  t value
## (Intercept)  75.9955836 10.4538297 7.269640
## genderM       5.3388469  2.9527772 1.808076
## age           0.3781262  0.1189145 3.179816
## BMI           1.1351249  0.4326554 2.623623
## workHours.mc  1.0390519  0.3140932 3.308100

DBP_aft

# model comparison vs null
m1.bisd <- lmer(DBP_aft ~ workHours.mc + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.142 0.858

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                Estimate Std. Error   t value
## (Intercept)  120.976976  1.5735066 76.883679
## workHours.mc   1.039052  0.3140563  3.308489

# model comparison vs previous model
m4d <- lmer(DBP_aft ~ gender + age + BMI + workHours.mc + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.004 0.142 0.854

summary(m4d)$coefficients # coefficient: |t| > 2

##                Estimate Std. Error   t value
## (Intercept)  42.7211824 7.78035413  5.490905
## genderM      -0.3652985 2.19609521 -0.166340
## age           0.2532791 0.08836923  2.866146
## BMI           1.0103643 0.32187966  3.138950
## workHours.mc  0.7109242 0.28329330  2.509499

4.1.6. teleWork

teleWork is a level-1 categorical variable that in section 3.5 only showed some slight difference in afternoon blood pressure. Here, whereas its inclusion implies weaker evidence than more parsimonious models, absolute t-values are lower than 2. Thus, also considering the unbalanced number of cases per level and the uncertainty of its relationship with both workaholism and blood pressure, we do not include it as a covariate for afternoon blood pressure.

SBP_aft

# number of cases by covariate level
summary(cleanBP$teleWork) # currently 4 categories

##   office teleWork     both   dayOff 
##      543      133       43        0

cleanBP$teleWork <- as.factor(gsub("both","teleWork",gsub("dayOff","teleWork",cleanBP$teleWork))) # recoding to 2 categories
summary(cleanBP$teleWork)

##   office teleWork 
##      543      176

# model comparison vs null
m1.bis <- lmer(SBP_aft ~ teleWork + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.187 0.813

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

##                    Estimate Std. Error   t value
## (Intercept)      121.565029   1.627667 74.686690
## teleWorkteleWork  -2.330782   1.572030 -1.482657

# model comparison vs previous model
m4 <- lmer(SBP_aft ~ gender + age + BMI + teleWork + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.014 0.165 0.821

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

##                    Estimate Std. Error   t value
## (Intercept)      76.3138065 10.4787050  7.282752
## genderM           5.6385384  2.9652628  1.901531
## age               0.3740369  0.1191998  3.137899
## BMI               1.1489370  0.4336873  2.649229
## teleWorkteleWork -2.4552154  1.5545308 -1.579393

DBP_aft

# model comparison vs null
m1.bisd <- lmer(DBP_aft ~ teleWork + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.424 0.576

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

##                    Estimate Std. Error   t value
## (Intercept)      121.565029   1.627667 74.686690
## teleWorkteleWork  -2.330782   1.572030 -1.482657

# model comparison vs previous model
m4d <- lmer(DBP_aft ~ gender + age + BMI + teleWork + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.001 0.013 0.426 0.560

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

##                    Estimate Std. Error    t value
## (Intercept)      42.7764078 7.80143722  5.4831445
## genderM          -0.3118548 2.20747157 -0.1412724
## age               0.2525514 0.08861229  2.8500713
## BMI               1.0129252 0.32275646  3.1383576
## teleWorkteleWork -0.4412502 1.35200750 -0.3263667

4.1.7. confounders_aft

confounders_aft is a level-1 categorical variable that in section 3.5 only showed some slight difference in afternoon blood pressure. Consistently, its inclusion implies weaker evidence than more parsimonious models. However, absolute t-values are lower than 2. Thus, also considering the unbalanced number of cases per level and the uncertainty of its relationship with workaholism, we do not include it as a covariate for afternoon blood pressure.

SBP_aft

# number of cases by covariate level
summary(cleanBP$confounders_aft)

##    Mode   FALSE    TRUE 
## logical     603     116

# model comparison vs null
m1.bis <- lmer(SBP_aft ~ confounders_aft + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.396 0.604

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

##                        Estimate Std. Error    t value
## (Intercept)         121.0797775   1.599142 75.7154629
## confounders_aftTRUE  -0.5699072   1.543322 -0.3692729

# model comparison vs previous model
m4 <- lmer(SBP_aft ~ gender + age + BMI + confounders_aft + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.001 0.032 0.371 0.596

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

##                       Estimate Std. Error    t value
## (Intercept)         76.0718943 10.4567638  7.2748984
## genderM              5.3736525  2.9539259  1.8191562
## age                  0.3804731  0.1190172  3.1967893
## BMI                  1.1330399  0.4327424  2.6182782
## confounders_aftTRUE -0.7757918  1.5303890 -0.5069246

DBP_aft

# model comparison vs null
m1.bisd <- lmer(DBP_aft ~ confounders_aft + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.415 0.585

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

##                        Estimate Std. Error    t value
## (Intercept)         121.0797775   1.599142 75.7154629
## confounders_aftTRUE  -0.5699072   1.543322 -0.3692729

# model comparison vs previous model
m4d <- lmer(DBP_aft ~ gender + age + BMI + confounders_aft + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.001 0.012 0.391 0.597

summary(m4d)$coefficients # coefficient: |t| > 2

##                       Estimate Std. Error    t value
## (Intercept)         42.8075146 7.77106123  5.5085803
## genderM             -0.3255433 2.19385493 -0.1483887
## age                  0.2558796 0.08833487  2.8967000
## BMI                  1.0080372 0.32146015  3.1358076
## confounders_aftTRUE -0.8695439 1.34340881 -0.6472668

4.1.8. Other covariates

In addition to the pre-registered covariates, we consider a further level-2 covariate that in section 3.5 showed substantial trends in both afternoon blood pressure and workaholism, and whose potential influence on both variables can be theoretically argued, namely having children children or living with children home_child. Indeed, respondents reporting to live with their children reported higher afternoon blood pressure, whereas those reporting to have children reported slightly lower workaholism (see section 3.5). Theoretically, it is plausible although not consistently reported by previous research (see Clark et al., 2016) that workaholism levels and overall investment in the work vs. family domain are influenced by family needs and presence of children in the household. On the other hand, family needs and living with children might negatively impact on blood pressure. Thus, both children and home_child might be a common external variable influencing both workaholism and afternoon, and we consider it as a further potential covariate.

Although both children and home_child have an acceptably balanced number of cases per level and both imply stronger evidence than more parsimonious models, none of them show an absolute t value higher than 2. Thus, also in this case, we do not include these two variables as covariates for afternoon blood pressure.

SBP_aft

# number of cases by covariate level
summary(cleanBP[!duplicated(cleanBP$ID),"children"]) # children: No. participants

##  No Yes 
##  55  51

summary(cleanBP$children) # children: No. cases

##  No Yes 
## 389 330

summary(cleanBP[!duplicated(cleanBP$ID),"home_child"]) # home_child: No. participants

##  No Yes 
##  66  40

summary(cleanBP$home_child) # home_child: No. cases

##  No Yes 
## 454 265

# model comparison children vs null
m1.bis <- lmer(SBP_aft ~ children + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.02 0.98

summary(m1.bis)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 117.406750   2.132033 55.067976
## childrenYes   7.443637   3.079216  2.417381

# model comparison home_child vs null
m1.bis <- lmer(SBP_aft ~ home_child + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.03 0.97

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                 Estimate Std. Error   t value
## (Intercept)   118.321044   1.957441 60.446809
## home_childYes   7.050187   3.190032  2.210068

# model comparison children vs previous model
m4 <- lmer(SBP_aft ~ gender + age + BMI + children + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.017 0.201 0.782

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

##               Estimate Std. Error   t value
## (Intercept) 76.3480979 10.8616067 7.0291717
## genderM      5.3307009  2.9678555 1.7961457
## age          0.3620579  0.1730404 2.0923323
## BMI          1.1373504  0.4350812 2.6141107
## childrenYes  0.5394472  4.2020256 0.1283779

# model comparison home_child vs previous model
m4 <- lmer(SBP_aft ~ gender + age + BMI + home_child + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.018 0.211 0.770

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

##                 Estimate Std. Error   t value
## (Intercept)   76.8276612 10.9164560 7.0377842
## genderM        5.4351193  2.9860037 1.8201984
## age            0.3530459  0.1493131 2.3644662
## BMI            1.1246307  0.4362791 2.5777780
## home_childYes  1.0671470  3.8114921 0.2799814

DBP_aft

# model comparison children vs null
m1.bisd <- lmer(DBP_aft ~ children + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.055 0.945

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                 Estimate Std. Error   t value
## (Intercept)   118.321044   1.957441 60.446809
## home_childYes   7.050187   3.190032  2.210068

# model comparison home_child vs null
m1.bisd <- lmer(DBP_aft ~ home_child + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.094 0.906

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                 Estimate Std. Error   t value
## (Intercept)   118.321044   1.957441 60.446809
## home_childYes   7.050187   3.190032  2.210068

# model comparison children vs previous model
m4d <- lmer(DBP_aft ~ gender + age + BMI + children + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.007 0.242 0.751

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

##               Estimate Std. Error    t value
## (Intercept) 41.9172261  8.0757799  5.1904864
## genderM     -0.3432938  2.2059915 -0.1556188
## age          0.2898331  0.1284984  2.2555391
## BMI          1.0052550  0.3234879  3.1075510
## childrenYes -1.2278268  3.1216864 -0.3933216

# model comparison home_child vs previous model
m4d <- lmer(DBP_aft ~ gender + age + BMI + home_child + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.007 0.238 0.755

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

##                 Estimate Std. Error    t value
## (Intercept)   41.3410570  8.1077717  5.0989419
## genderM       -0.5201802  2.2160967 -0.2347281
## age            0.2949959  0.1109493  2.6588355
## BMI            1.0274123  0.3239097  3.1719098
## home_childYes -1.7699884  2.8327357 -0.6248336

4.1.9. Time trends

As a final, not pre-registered, control, we inspect the linear time trends of the outcome variable. This is done by including the number of days since the beginning of the study as a further covariate. We can see that the inclusion of day does not imply stronger evidence than more parsimonious models, with an absolute t value lower than 2. Thus, we chose to not include day as a covariate for afternoon blood pressure.

SBP_aft

# plotting SBP_aft by day
boxplot(SBP_aft ~ day,data=cleanBP)

# model comparison vs null
m1.bis <- lmer(SBP_aft ~ day + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.492 0.508

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

##                Estimate Std. Error  t value
## (Intercept) 122.2885438  1.6916467 72.28965
## day          -0.2625006  0.1253453 -2.09422

# model comparison vs previous model
m4 <- lmer(SBP_aft ~ gender + age + BMI + day + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.000 0.001 0.042 0.479 0.479

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

##               Estimate Std. Error   t value
## (Intercept) 77.1827708 10.4582311  7.380098
## genderM      5.1907294  2.9503987  1.759332
## age          0.3765780  0.1187832  3.170297
## BMI          1.1456268  0.4322120  2.650613
## day         -0.2605224  0.1252345 -2.080277

DBP_aft

# plotting DBP_aft by day
boxplot(DBP_aft ~ day,data=cleanBP)

# model comparison vs null
m1.bisd <- lmer(DBP_aft ~ day + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.889 0.111

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

##                Estimate Std. Error  t value
## (Intercept) 122.2885438  1.6916467 72.28965
## day          -0.2625006  0.1253453 -2.09422

# model comparison vs previous model
m4d <- lmer(DBP_aft ~ gender + age + BMI + day + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.001 0.001 0.025 0.862 0.111

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

##                Estimate Std. Error    t value
## (Intercept) 43.05205334 7.79302433  5.5244346
## genderM     -0.40633049 2.19572316 -0.1850554
## age          0.25285381 0.08831641  2.8630445
## BMI          1.01347650 0.32172096  3.1501725
## day         -0.07331483 0.11250810 -0.6516405

4.2. Evening BP

For evening systolic and diastolic blood pressure, the pre-registered potential covariates to be considered are:

1. gender, in light of the robust literature on sex differences in blood pressure and hypertension (e.g., [Sandberg & Ji, 2012]), with gender having been highlighted as a potential moderator of the straining effects of trait workaholism (e.g., Bladucci et al., 2018; 2021).

2. age, in light of the consolidated relationship between aging and high blood pressure; however, no clear differences have been highlighted on the relationship between age and workaholism (Clark et al., 2016).

3. BMI, in light of the robust literature on the increased risk for hypertension and CVD in individuals with higher BMI (e.g., Gregg et al., 2005); however, no clear differences have been highlighted, to our knowledge, on the relationship between BMI and workaholism.

4. Participants’ smoking status smoker, in light of the widely acknowledged relationship between smoking and hypertension; however, no clear differences have been highlighted, to our knowledge, on the relationship between smoking status and workaholism.

5. The daily number of worked hours workHours, as long working hours is a definitional aspect of workaholism which has been found positively associated with higher blood pressure and other cardiovascular risks (Virtanen, 2012); however, controlling for this variable might risk to artificially deflate the relationship between workaholism and blood pressure since the former might influence the latter through working hours. That is, working hours is a potential explanatory mechanism of the focused relationship, and including it as a control variable might bias rather than strengthening the results (see Spector et al., 2021).

6. The occurrence of late working hours in the evening lateWorkHours, consistently with what argued above for workHours; also in this case, controlling for this variable might risk to rule out a potential explanatory mechanism of the relationship between workaholism and blood pressure.

7. Nonwork hassles in the evening dailyHassles_eve, as hassles have been related to high blood pressure (e.g., Nyklíček et al. 1999) and nonwork hassles such as relationship tension have been recently related to state workaholism (Clark et al., 2021); also in this case, such an argument imply that nonwork hassles are a potential explanatory mechanism of the focused relationship.

8. Remote work teleWork, although with no clear expectations on its relationship with blood pressure, while we originally hypothesized that remote working might act as a moderator of the straining effect of workaholism.

9. BP confounders_eve identifying behavioral factors (e.g., smoking, physical activity) that might affect blood pressure, although with no clear expectation on their relationship with workaholism.

Further pre-registered level-2 covariates (i.e., retrospective measures of trait workaholism, technostress creators) are not included due to our aim to focus the study on the focused level-1 variables, whereas recovery experiences will be included in all models. Moreover, we do not include the previous BP measure as we focus on the covariance of intra-individual fluctuations between workaholism and blood pressure, rather than the change of the latter over time (see Supplementary material S7). In contrast, additional (i.e., not pre-registered) level-2 covariates potentially moderating the relationship between workaholism and evening blood pressure are considered as potentially confounding factors.

Prior to evaluate the inclusion of each covariate, we prepare the data for the analyses. That is, we remove cases of missing responses in the response variable or any covariate (list-wise deletion) and we center level-1 continuous covariates on the individual mean (person-mean-centering). Moreover, we apply our inclusion criteria by removing all cases with bp_drugs or cv_dysf = "Yes".

# inclusion criteria
cleanBP <- clean[clean$bp_drugs=="No" & diary$cv_dysf=="No",]

## Warning in clean$bp_drugs == "No" & diary$cv_dysf == "No": longer object length
## is not a multiple of shorter object length

# list-wise deletion
cleanBP <- as.data.frame(na.omit(cleanBP[,c("ID","day","SBP_eve","DBP_eve",
                                        "gender","age","BMI","smoker","workHours",
                                        "lateWorkHours","dailyHassles_eve","teleWork","confounders_eve",
                                        "home_child","children")]))
cat("Considering",nrow(cleanBP),"complete observations from",nlevels(as.factor(as.character(cleanBP$ID))),"participants")

## Considering 692 complete observations from 104 participants

# mean centering lv-1 continuous covariates
for(Var in c("workHours")){
  wide <- prelqs[prelqs$ID%in%levels(as.factor(as.character(cleanBP$ID))),]
  wide <- cbind(wide,aggregate(cleanBP[,Var],list(cleanBP$ID),mean)[,2]) # computing individual means
  colnames(wide)[ncol(wide)] <- paste0(Var,".cm")
  cleanBP <- join(cleanBP,wide[,c("ID",paste0(Var,".cm"))],by="ID",type="left") # joining with long-form data
  cleanBP[,paste0(Var,".mc")] <- cleanBP[,Var] - cleanBP[,paste0(Var,".cm")] } # computing mean-centered scores

4.2.1. gender

gender is a categorical variable with a roughly balanced number of participants. In section 3.5 we observed some tendency with higher SBP_aft in men than in women. Consistently, its inclusion in the null model shows improved evidence (i.e., higher Akaike weight) compared to the null model, although the estimated t-value is lower than 2, with a particularly low value for DBP_eve. In both cases, also considering its theoretical relevance as a potential moderator of the straining effects of workaholism (e.g., Balducci et al. (2018; 2021)) we include gender as a covariate for evening blood pressure.

SBP_eve

# number of cases by covariate level
summary(cleanBP[!duplicated(cleanBP$ID),"gender"]) # No. participants

##  F  M 
## 56 48

summary(cleanBP$gender) # No. cases

##   F   M 
## 378 314

# model comparison
m0 <- lmer(SBP_eve ~ (1|ID),data=cleanBP)
m1 <- lmer(SBP_eve ~ gender + (1|ID),data=cleanBP)
Weights(AIC(m0,m1)) # Akaike weights: stronger evidence including gender

##  model weights 
## [1] 0.182 0.818

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

##               Estimate Std. Error    t value
## (Intercept) 112.772030   2.047974 55.0651613
## genderM       2.947148   3.009381  0.9793205

DBP_eve

# model comparison
m0d <- lmer(DBP_eve ~ (1|ID),data=cleanBP)
m1d <- lmer(DBP_eve ~ gender + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d)) # Akaike weights: stronger evidence including gender

##  model weights 
## [1] 0.283 0.717

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

##              Estimate Std. Error    t value
## (Intercept) 72.612532   1.623491 44.7261620
## genderM     -1.276225   2.385144 -0.5350723

4.2.2. age

age is a level-2 continuous variable that in section 3.2 showed moderate-to-strong positive relationships with blood pressure (r from .42 to .52). Consistently, its inclusion in the model including gender implies stronger evidence than more parsimonious models, with an absolute t-value higher than 2. Thus, despite the uncertainty on its relationship with workaholism, we include age as a covariate for evening blood pressure.

SBP_eve

# model comparison vs null
m1.bis <- lmer(SBP_eve ~ age + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.092 0.908

summary(m1.bis)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 99.4140642  5.0313896 19.758769
## age          0.3586221  0.1174414  3.053627

# model comparison vs gender
m2 <- lmer(SBP_eve ~ gender + age + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.013 0.060 0.926

summary(m2)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 96.7850331   5.353289 18.079544
## genderM      4.0251342   2.902976  1.386555
## age          0.3772338   0.117619  3.207252

DBP_eve

# model comparison vs null
m1.bisd <- lmer(DBP_eve ~ age + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.02 0.98

summary(m1.bisd)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 58.3154142 3.90114715 14.948273
## age          0.3339043 0.09107584  3.666222

# model comparison vs gender
m2d <- lmer(DBP_eve ~ gender + age + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.009 0.023 0.968

summary(m2d)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 58.5253156 4.19231174 13.960154
## genderM     -0.3232550 2.27397847 -0.142154
## age          0.3324499 0.09212372  3.608733

4.2.3. BMI

BMI is a level-2 continuous variable that in section 3.2 showed weak-to-moderate positive relationships with blood pressure (r from .28 to .43). Consistently, its inclusion in the model including gender and age implies stronger evidence than more parsimonious models, with an absolute t-value higher than 2. Thus, despite the uncertainty on its relationship with workaholism, we include BMI as a covariate for evening blood pressure.

SBP_eve

# model comparison vs null
m1.bis <- lmer(SBP_eve ~ BMI + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.001 0.999

summary(m1.bis)$coefficients # coefficient: |t| > 2

##              Estimate Std. Error  t value
## (Intercept) 74.637445  9.6936361 7.699634
## BMI          1.666018  0.4047412 4.116254

# model comparison vs previous model
m3 <- lmer(SBP_eve ~ gender + age + BMI + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.008 0.992

summary(m3)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 66.7536942  9.9521132 6.7074895
## genderM      1.8651042  2.8177119 0.6619215
## age          0.3015323  0.1134012 2.6589857
## BMI          1.4399848  0.4105817 3.5071816

DBP_eve

# model comparison vs null
m1.bisd <- lmer(DBP_eve ~ BMI + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.002 0.998

summary(m1.bisd)$coefficients # coefficient: |t| > 2

##              Estimate Std. Error  t value
## (Intercept) 41.647371  7.6992896 5.409249
## BMI          1.281191  0.3214748 3.985353

# model comparison vs previous model
m3d <- lmer(DBP_eve ~ gender + age + BMI + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.005 0.995

summary(m3d)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error    t value
## (Intercept) 33.7373006 7.75857430  4.3483892
## genderM     -2.1071024 2.19458836 -0.9601356
## age          0.2703706 0.08827595  3.0627886
## BMI          1.1879624 0.32009100  3.7113271

4.2.4. smoker

smoker is a level-2 categorical variable that in section 3.5 did not show substantially relationships with blood pressure. Consistently, both models show an absolute t-value lower than 2 when the other covariates are considered, although its inclusion implies stronger evidence than more parsimonious models. Also considering the uncertainty of its relationship with workaholism, we do not include it as a covariate for evening blood pressure.

SBP_eve

# model comparison vs null
m1.bis <- lmer(SBP_eve ~ smoker + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.178 0.822

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

##               Estimate Std. Error    t value
## (Intercept) 114.979132   1.770017 64.9593428
## smokerYes    -3.013209   3.348420 -0.8998896

# model comparison vs previous model
m4 <- lmer(SBP_eve ~ gender + age + BMI + smoker + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.002 0.213 0.786

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

##               Estimate Std. Error    t value
## (Intercept) 67.9810435 10.1460246  6.7002640
## genderM      2.4284192  2.9485341  0.8236022
## age          0.2944551  0.1142132  2.5781170
## BMI          1.4147582  0.4134176  3.4221048
## smokerYes   -2.1451698  3.2071411 -0.6688729

DBP_eve

# model comparison vs null
m1.bisd <- lmer(DBP_eve ~ smoker + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.06 0.94

summary(m1.bisd)$coefficients # coefficient: |t| > 2

##              Estimate Std. Error   t value
## (Intercept) 73.437284   1.377685 53.304850
## smokerYes   -5.072819   2.605867 -1.946692

# model comparison vs previous model
m4d <- lmer(DBP_eve ~ gender + age + BMI + smoker + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.001 0.155 0.845

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

##               Estimate Std. Error    t value
## (Intercept) 35.6083679 7.85886356  4.5309818
## genderM     -1.2444063 2.28168622 -0.5453889
## age          0.2595295 0.08832129  2.9384706
## BMI          1.1495439 0.32023108  3.5897325
## smokerYes   -3.2765622 2.48102605 -1.3206481

4.2.5. WorkHours

workHours is a level-1 continuous variable that in section 3.2 was not substantially correlated with evening blood pressure (r = .03-.06). Here, its inclusion implies weaker evidence (i.e., lower Akaike weight) than more parsimonious models, with an estimated absolute t-value lower than 2. Thus, also considering the theoretical importance of this variable, which might be a potential explanatory mechanism (e.g., mediator) through which workaholism might impact on blood pressure, we do not include it as a covariate for evening blood pressure.

SBP_eve

# model comparison vs null
m1.bis <- lmer(SBP_eve ~ workHours.mc + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.783 0.217

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

##                  Estimate Std. Error     t value
## (Intercept)  114.13682845  1.5009580 76.04265396
## workHours.mc   0.02148807  0.3004189  0.07152703

# model comparison vs previous model
m4 <- lmer(SBP_eve ~ gender + age + BMI + workHours.mc + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.006 0.778 0.216

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

##                 Estimate Std. Error    t value
## (Intercept)  66.75274308  9.9521610 6.70736164
## genderM       1.86520431  2.8177044 0.66195883
## age           0.30152409  0.1134005 2.65893175
## BMI           1.44003076  0.4105838 3.50727617
## workHours.mc  0.02148807  0.3004662 0.07151576

DBP_eve

# model comparison vs null
m1.bisd <- lmer(DBP_eve ~ workHours.mc + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.799 0.201

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

##                 Estimate Std. Error    t value
## (Intercept)  72.02074791  1.1848960 60.7823382
## workHours.mc -0.08479742  0.2585833 -0.3279307

# model comparison vs previous model
m4d <- lmer(DBP_eve ~ gender + age + BMI + workHours.mc + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.004 0.796 0.200

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

##                 Estimate Std. Error    t value
## (Intercept)  33.73656247 7.75869758  4.3482250
## genderM      -2.10699373 2.19460644 -0.9600782
## age           0.27036578 0.08827629  3.0627223
## BMI           1.18799666 0.32009611  3.7113749
## workHours.mc -0.08479742 0.25853889 -0.3279871

4.2.6. lateWorkHours

lateWorkHours is a level-1 categorical variable that in section 3.5 only showed slight differences in evening systolic blood pressure. Consistently, its inclusion implies stronger evidence than more parsimonious models, with an absolute t-value higher than 2 for SBP_eve. However, we consider its role as a potential explanatory mechanism (e.g., mediator) through which workaholism might impact on blood pressure. Thus, also in light of the unbalanced number of cases, we do not include it as a covariate for evening blood pressure.

SBP_eve

# summary of cases
summary(cleanBP$lateWorkHours)

##    Mode   FALSE    TRUE 
## logical     530     162

# model comparison vs null
m1.bis <- lmer(SBP_eve ~ lateWorkHours + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.072 0.928

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                     Estimate Std. Error   t value
## (Intercept)       113.488306   1.522966 74.517927
## lateWorkHoursTRUE   2.607868   1.169246  2.230385

# model comparison vs previous model
m4 <- lmer(SBP_eve ~ gender + age + BMI + lateWorkHours + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.000 0.056 0.943

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

##                     Estimate Std. Error   t value
## (Intercept)       65.6998746  9.8959439 6.6390710
## genderM            1.9966948  2.7994116 0.7132552
## age                0.2998013  0.1126424 2.6615311
## BMI                1.4562411  0.4079060 3.5700411
## lateWorkHoursTRUE  2.7321405  1.1628043 2.3496133

DBP_eve

# model comparison vs null
m1.bisd <- lmer(DBP_eve ~ lateWorkHours + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.459 0.541

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

##                     Estimate Std. Error    t value
## (Intercept)       71.8475051   1.208679 59.4429936
## lateWorkHoursTRUE  0.6972278   1.004637  0.6940097

# model comparison vs previous model
m4d <- lmer(DBP_eve ~ gender + age + BMI + lateWorkHours + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.002 0.438 0.559

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

##                     Estimate Std. Error    t value
## (Intercept)       33.4270451 7.74499888  4.3159522
## genderM           -2.0683839 2.18856218 -0.9450880
## age                0.2698522 0.08801337  3.0660368
## BMI                1.1927294 0.31920098  3.7366095
## lateWorkHoursTRUE  0.8089683 0.99692009  0.8114675

4.2.7. dailyHassles_eve

dailyHassles_eve is a level-1 categorical variable that in section 3.5 predicted slightly higher evening systolic blood pressure. Consistently, its inclusion implies stronger evidence (i.e., higher Akaike weight) than more parsimonious models, with an estimated absolute t-value higher than 2 for both systolic and diastolic evening blood pressure. However, the variable shows a strongly unbalanced number of cases, with only 10% of cases with some evening hassle. Moreover, we consider the theoretical importance of this variable, which might be a potential explanatory mechanism (e.g., mediator) through which workaholism might impact on strain. Thus, we do not include it as a covariate for evening blood pressure.

SBP_eve

# number of cases by covariate level
summary(cleanBP$dailyHassles_eve)

##  No Yes 
## 620  72

# model comparison vs null
m1.bis <- lmer(SBP_eve ~ dailyHassles_eve + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.045 0.955

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                       Estimate Std. Error   t value
## (Intercept)         113.796757   1.512925 75.216380
## dailyHassles_eveYes   3.382891   1.430865  2.364227

# model comparison vs previous model
m4 <- lmer(SBP_eve ~ gender + age + BMI + dailyHassles_eve + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.000 0.041 0.959

summary(m4)$coefficients # coefficient: |t| > 2

##                       Estimate Std. Error   t value
## (Intercept)         66.0268707  9.9807438 6.6154258
## genderM              2.0319644  2.8255272 0.7191452
## age                  0.2990615  0.1136900 2.6305004
## BMI                  1.4570446  0.4116337 3.5396631
## dailyHassles_eveYes  3.4402674  1.4282481 2.4087324

DBP_eve

# model comparison vs null
m1.bisd <- lmer(DBP_eve ~ dailyHassles_eve + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.022 0.978

summary(m1.bisd)$coefficients # coefficient: |t| > 2

##                      Estimate Std. Error   t value
## (Intercept)         71.685385   1.196960 59.889560
## dailyHassles_eveYes  3.330774   1.226915  2.714755

# model comparison vs previous model
m4d <- lmer(DBP_eve ~ gender + age + BMI + dailyHassles_eve + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.000 0.023 0.977

summary(m4d)$coefficients # coefficient: |t| > 2

##                       Estimate Std. Error    t value
## (Intercept)         33.0442605 7.80296420  4.2348343
## genderM             -1.9474501 2.20697544 -0.8824068
## age                  0.2679757 0.08875228  3.0193674
## BMI                  1.2041648 0.32180164  3.7419473
## dailyHassles_eveYes  3.3028539 1.22338067  2.6997761

4.2.8. teleWork

teleWork is a level-1 categorical variable that in section 3.5 only showed some slight difference in evening blood pressure. Here, although its inclusion implies stronger evidence than more parsimonious models, absolute t-values are lower than 2. Thus, also considering the unbalanced number of cases per level and the uncertainty of its relationship with both workaholism and blood pressure, we do not include it as a covariate for evening blood pressure.

SBP_eve

# number of cases by covariate level
summary(cleanBP$teleWork) # currently 4 categories

##   office teleWork     both   dayOff 
##      521      129       42        0

cleanBP$teleWork <- as.factor(gsub("both","teleWork",gsub("dayOff","teleWork",cleanBP$teleWork))) # recoding to 2 categories
summary(cleanBP$teleWork)

##   office teleWork 
##      521      171

# model comparison vs null
m1.bis <- lmer(SBP_eve ~ teleWork + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.323 0.677

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

##                    Estimate Std. Error    t value
## (Intercept)      113.794472   1.544992 73.6537680
## teleWorkteleWork   1.366617   1.508391  0.9060097

# model comparison vs previous model
m4 <- lmer(SBP_eve ~ gender + age + BMI + teleWork + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.003 0.335 0.663

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

##                    Estimate Std. Error   t value
## (Intercept)      66.6300806  9.9413202 6.7023372
## genderM           1.7281196  2.8189077 0.6130458
## age               0.3030699  0.1132787 2.6754370
## BMI               1.4317496  0.4102103 3.4902823
## teleWorkteleWork  1.2733757  1.4912985 0.8538705

DBP_eve

# model comparison vs null
m1.bisd <- lmer(DBP_eve ~ teleWork + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.379 0.621

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

##                   Estimate Std. Error    t value
## (Intercept)      71.760041   1.224678 58.5950051
## teleWorkteleWork  1.039961   1.282581  0.8108345

# model comparison vs previous model
m4d <- lmer(DBP_eve ~ gender + age + BMI + teleWork + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.002 0.361 0.637

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

##                    Estimate Std. Error    t value
## (Intercept)      33.6242929 7.73173940  4.3488653
## genderM          -2.2298695 2.19081519 -1.0178264
## age               0.2717364 0.08796724  3.0890637
## BMI               1.1805730 0.31905459  3.7002226
## teleWorkteleWork  1.1519514 1.26191981  0.9128563

4.2.9. confounders_eve

confounders_eve is a level-1 categorical variable that in section 3.5 showed some slight difference in afternoon blood pressure. Here, although its inclusion implies weaker evidence than more parsimonious models, absolute t-values are lower than 2. Thus, also considering the unbalanced number of cases per level and the uncertainty of its relationship with workaholism, we do not include it as a covariate for evening blood pressure.

SBP_eve

# number of cases by covariate level
summary(cleanBP$confounders_eve)

##    Mode   FALSE    TRUE 
## logical     556     136

# model comparison vs null
m1.bis <- lmer(SBP_eve ~ confounders_eve + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.408 0.592

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

##                        Estimate Std. Error    t value
## (Intercept)         114.2968993   1.531994 74.6066118
## confounders_eveTRUE  -0.7463954   1.353309 -0.5515335

# model comparison vs previous model
m4 <- lmer(SBP_eve ~ gender + age + BMI + confounders_eve + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.003 0.415 0.582

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

##                       Estimate Std. Error    t value
## (Intercept)         66.9646776  9.9856414  6.7060968
## genderM              1.9679221  2.8324620  0.6947744
## age                  0.3032453  0.1137364  2.6662115
## BMI                  1.4321009  0.4118961  3.4768496
## confounders_eveTRUE -0.6633106  1.3466258 -0.4925723

DBP_eve

# model comparison vs null
m1.bisd <- lmer(DBP_eve ~ confounders_eve + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.477 0.523

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

##                       Estimate Std. Error    t value
## (Intercept)         71.9611667   1.210820 59.4317760
## confounders_eveTRUE  0.2791632   1.156894  0.2413041

# model comparison vs previous model
m4d <- lmer(DBP_eve ~ gender + age + BMI + confounders_eve + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.002 0.469 0.528

summary(m4d)$coefficients # coefficient: |t| > 2

##                      Estimate Std. Error    t value
## (Intercept)         33.609966 7.76754741  4.3269728
## genderM             -2.169489 2.20197992 -0.9852447
## age                  0.269359 0.08832794  3.0495338
## BMI                  1.192690 0.32039420  3.7225691
## confounders_eveTRUE  0.401611 1.14812688  0.3497967

4.2.10. Other covariates

In addition to the pre-registered covariates, we consider a further level-2 covariate that in section 3.5 showed substantial trends in both evening blood pressure and workaholism, and whose potential influence on both variables can be theoretically argued, namely having children children or living with children home_child. Indeed, respondents reporting to live with their children reported higher evening blood pressure, whereas those reporting to have children reported slightly lower workaholism (see section 3.5). Theoretically, it is plausible although not consistently reported by previous research (see Clark et al., 2016) that workaholism levels and overall investment in the work vs. family domain are influenced by family needs and presence of children in the household. On the other hand, family needs and living with children might negatively impact on blood pressure. Thus, both children and home_child might be a common external variable influencing both workaholism and evening, and we consider it as a further potential covariate.

Although both children and home_child have an acceptably balanced number of cases per level and both imply stronger evidence than more parsimonious models, none of them show an absolute t value higher than 2. Thus, also in this case, we do not include these two variables as covariates for evening blood pressure.

SBP_eve

# number of cases by covariate level
summary(cleanBP[!duplicated(cleanBP$ID),"children"]) # children: No. participants

##  No Yes 
##  54  50

summary(cleanBP$children) # children: No. cases

##  No Yes 
## 384 308

summary(cleanBP[!duplicated(cleanBP$ID),"home_child"]) # home_child: No. participants

##  No Yes 
##  65  39

summary(cleanBP$home_child) # home_child: No. cases

##  No Yes 
## 446 246

# model comparison children vs null
m1.bis <- lmer(SBP_eve ~ children + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.044 0.956

summary(m1.bis)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 111.253421   2.043111 54.452940
## childrenYes   6.052529   2.961321  2.043861

# model comparison home_child vs null
m1.bis <- lmer(SBP_eve ~ home_child + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.071 0.929

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

##                 Estimate Std. Error   t value
## (Intercept)   112.128619   1.875623 59.782080
## home_childYes   5.403577   3.077756  1.755687

# model comparison children vs previous model
m4 <- lmer(SBP_eve ~ gender + age + BMI + children + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.002 0.214 0.784

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

##               Estimate Std. Error     t value
## (Intercept) 66.6125958 10.3254484  6.45130299
## genderM      1.8685767  2.8328983  0.65959896
## age          0.3083002  0.1636152  1.88430051
## BMI          1.4387451  0.4129912  3.48371830
## childrenYes -0.2271017  3.9672903 -0.05724352

# model comparison home_child vs previous model
m4 <- lmer(SBP_eve ~ gender + age + BMI + home_child + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.002 0.226 0.773

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

##                 Estimate Std. Error    t value
## (Intercept)   66.0912873 10.3746008  6.3704897
## genderM        1.7842699  2.8503310  0.6259869
## age            0.3217752  0.1414442  2.2749274
## BMI            1.4481325  0.4139393  3.4984179
## home_childYes -0.8711060  3.6150390 -0.2409672

DBP_eve

# model comparison children vs null
m1.bisd <- lmer(DBP_eve ~ children + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.032 0.968

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

##              Estimate Std. Error  t value
## (Intercept) 69.461465   1.604318 43.29656
## childrenYes  5.381759   2.327368  2.31238

# model comparison home_child vs null
m1.bisd <- lmer(DBP_eve ~ home_child + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.073 0.927

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

##                Estimate Std. Error   t value
## (Intercept)   70.336214   1.477690 47.598754
## home_childYes  4.539922   2.426653  1.870858

# model comparison children vs previous model
m4d <- lmer(DBP_eve ~ gender + age + BMI + children + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.001 0.255 0.744

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

##               Estimate Std. Error    t value
## (Intercept) 33.2719874  8.0446895  4.1358945
## genderM     -2.0934794  2.2056867 -0.9491282
## age          0.2919544  0.1272977  2.2934770
## BMI          1.1845708  0.3218551  3.6804475
## childrenYes -0.7287210  3.0872911 -0.2360390

# model comparison home_child vs previous model
m4d <- lmer(DBP_eve ~ gender + age + BMI + home_child + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.001 0.228 0.771

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

##                 Estimate Std. Error    t value
## (Intercept)   32.1652186  8.0651917  3.9881530
## genderM       -2.2924693  2.2132622 -1.0357875
## age            0.3181343  0.1099083  2.8945428
## BMI            1.2073609  0.3217968  3.7519356
## home_childYes -2.0562008  2.8095946 -0.7318496

4.2.11. Time trends

As a final, not pre-registered, control, we inspect the linear time trends of the outcome variable. This is done by including the number of days since the beginning of the study as a further covariate. We can see that the inclusion of day implies stronger evidence than more parsimonious models, with an absolute t value higher than 2. However, a better inspection of blood pressure distributions through day levels shows that such effect is likely due to measurement reactivity (i.e., increased physiological activation during the first - less familiar - blood pressure recordings). Consistently the linear time trend is no longer substantial after the removal of observations collected on the first day, whereas the effect is larger when day is included as a categorical variable with day 1 vs. other days. Thus, we do not include day as a covariate** for evening blood pressure, but we will consider it as a robustness check.

SBP_eve

# plotting SBP_eve by day
boxplot(SBP_eve ~ day,data=cleanBP)

# model comparison vs null
m1.bis <- lmer(SBP_eve ~ day + (1|ID),data=cleanBP)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.369 0.631

summary(m1.bis)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 115.623226  1.6314245 70.872558
## day          -0.293606  0.1264509 -2.321897

# removing day 1
Weights(AIC(update(m0,data=cleanBP[cleanBP$day!=1,]),update(m1.bis,data=cleanBP[cleanBP$day!=1,]))) # weaker evidence than null

##  model weights 
## [1] 0.887 0.113

summary(update(m1.bis,data=cleanBP[cleanBP$day!=1,]))$coefficients # coefficient: |t| < 2

##                 Estimate Std. Error    t value
## (Intercept) 113.53298973  1.6680511 68.0632574
## day          -0.02416865  0.1360193 -0.1776855

# model comparison vs previous model
m4 <- lmer(SBP_eve ~ gender + age + BMI + day + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.003 0.343 0.654

summary(m4)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error    t value
## (Intercept) 67.9672248  9.9542768  6.8279420
## genderM      1.7906621  2.8148466  0.6361491
## age          0.2999060  0.1132825  2.6474161
## BMI          1.4569939  0.4101931  3.5519709
## day         -0.2990559  0.1262507 -2.3687457

# removing day 1
Weights(AIC(update(m0,data=cleanBP[cleanBP$day!=1,]),update(m1,data=cleanBP[cleanBP$day!=1,]),
            update(m2,data=cleanBP[cleanBP$day!=1,]),update(m3,data=cleanBP[cleanBP$day!=1,]),
            update(m4,data=cleanBP[cleanBP$day!=1,]))) # weaker evidence than null

##  model weights 
## [1] 0.000 0.001 0.008 0.878 0.112

summary(update(m4,data=cleanBP[cleanBP$day!=1,]))$coefficients # coefficient: |t| < 2

##                Estimate Std. Error    t value
## (Intercept) 67.97278905  9.8950264  6.8693894
## genderM      1.67363662  2.8021972  0.5972587
## age          0.27923260  0.1128915  2.4734595
## BMI          1.40698848  0.4069752  3.4571847
## day         -0.02865996  0.1357925 -0.2110570

# isolating day 1 as a factor
cleanBP$day1 <- 0
cleanBP[cleanBP$day==1,"day1"] <- 1
cleanBP$day1 <- as.factor(cleanBP$day1)
summary(cleanBP$day1) # number of cases

##   0   1 
## 614  78

boxplot(SBP_eve ~ day1,data=cleanBP) # plotting SBP by day

# model comparison considering day1 vs previous model
m4 <- lmer(SBP_eve ~ gender + age + BMI + day1 + (1|ID),data=cleanBP)
Weights(AIC(m0,m1,m2,m3,m4)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0 0 0 0 1

summary(m4)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 65.8641068  9.9030861 6.6508668
## genderM      1.8273418  2.8034531 0.6518182
## age          0.2961376  0.1128389 2.6244291
## BMI          1.4590298  0.4084938 3.5717305
## day11        5.2651024  1.2107433 4.3486530

DBP_eve

# plotting DBP_eve by day
boxplot(DBP_eve ~ day,data=cleanBP)

# model comparison vs null
m1.bisd <- lmer(DBP_eve ~ day + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1.bisd)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.308 0.692

summary(m1.bisd)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 73.3960692  1.3062922 56.186565
## day         -0.2714198  0.1085932 -2.499417

# model comparison vs previous model
m4d <- lmer(DBP_eve ~ gender + age + BMI + day + (1|ID),data=cleanBP)
Weights(AIC(m0d,m1d,m2d,m3d,m4d)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.000 0.000 0.001 0.274 0.725

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

##               Estimate Std. Error    t value
## (Intercept) 34.8733854 7.75960457  4.4942220
## genderM     -2.1772548 2.19155573 -0.9934745
## age          0.2688984 0.08815066  3.0504415
## BMI          1.2034270 0.31967205  3.7645675
## day         -0.2778330 0.10831674 -2.5650056

4.3. Emotional Exhaustion

For emotional exhaustion, the pre-registered potential covariates to be considered are:

1. gender, as previous evidence highlighted higher emotional exhaustion in women compared to men (e.g., Purvanova et al 2010) and gender has been highlighted as a potential moderator of the straining effects of trait workaholism (e.g., Bladucci et al., 2018; 2021).

2. age, as previous evidence highlighted some emotional exhaustion trend over age (e.g., Brewer et al., 2004); however, no clear differences have been highlighted for workaholism over age (Clark et al., 2016).

3. The daily number of worked hours workHours, as long working hours is a definitional aspect of workaholism which has been found associated with emotional exhaustion (Seidler et al, 2014); however, controlling for this variable might risk to artificially deflate the relationship between workaholism and sleep disturbances since the former might influence the latter through working hours. That is, working hours is a potential explanatory mechanism of the focused relationship, and including it as a control variable might bias rather than strengthening the results (see Spector et al., 2021).

4. The occurrence of late working hours in the evening lateWorkHours, consistently with what argued above for workHours; also in this case, controlling for this variable might risk to rule out a potential explanatory mechanism of the relationship between workaholism and sleep disturbances.

5. Nonwork hassles in the evening dailyHassles_eve, as private life hassles have been related to emotional exhaustion (e.g., Klusmann et al., 2021) and nonwork hassles such as relationship tension have been recently related to state workaholism (Clark et al., 2021); also in this case, such an argument imply that nonwork hassles are a potential explanatory mechanism of the focused relationship.

6. Remote work teleWork, as emotional exhaustion might be different during remote work days (e.g., Abdel Hadi et al., 2021), although with no clear expectations on its relationship with workaholism.

Further pre-registered level-2 covariates (i.e., retrospective measures of trait workaholism, technostress creators) are not included due to our aim to focus the study on the focused level-1 variables, whereas recovery experiences will be included in all models (see Supplementary material S7). Blood pressure confounders are not considered as well, due to the lack of clear link with exhaustion and workaholism. In contrast, additional (i.e., not pre-registered) level-2 covariates potentially moderating the relationship between workaholism and emotional exhaustion are considered as potentially confounding factors.

Prior to evaluate the inclusion of each covariate, we prepare the data for the analyses. That is, we remove cases of missing responses in the response variable or any covariate (list-wise deletion) and we center level-1 continuous covariates on the individual mean (person-mean-centering).

# list-wise deletion culo
cleanEE <- as.data.frame(na.omit(clean[,c("ID","day","EE","gender","age","workHours","lateWorkHours",
                                          "dailyHassles_eve","teleWork","home_child","children")]))
cat("Considering",nrow(cleanEE),"complete observations from",nlevels(as.factor(as.character(cleanEE$ID))),"participants")

## Considering 774 complete observations from 114 participants

# mean centering lv-1 continuous covariates
for(Var in c("workHours")){
  prelqs <- cbind(prelqs,aggregate(cleanEE[,Var],list(cleanEE$ID),mean)[,2]) # computing individual means
  colnames(prelqs)[ncol(prelqs)] <- paste0(Var,".cm")
  cleanEE <- join(cleanEE,prelqs[,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

4.3.1. gender

gender is a categorical variable with a roughly balanced number of participants. In section 3.5 we did not observed substantial gender differences in emotional exhaustion. Consistently, its inclusion in the null model shows weaker evidence (i.e., lower Akaike weight), with an estimated t-value lower than 2. However, considering its theoretical relevance as a potential moderator of the straining effects of workaholism (e.g., Balducci et al. (2018)) we include gender as a covariate for sleep disturbances.

# number of cases by covariate level
summary(cleanEE[!duplicated(cleanEE$ID),"gender"]) # No. participants

##  F  M 
## 59 55

summary(cleanEE$gender) # No. cases

##   F   M 
## 411 363

# model comparison
m0 <- lmer(EE ~ (1|ID),data=cleanEE)
m1 <- lmer(EE ~ gender + (1|ID),data=cleanEE)
Weights(AIC(m0,m1)) # Akaike weights: weaker evidence including gender

##  model weights 
## [1] 0.725 0.275

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

##              Estimate Std. Error   t value
## (Intercept) 3.1002605  0.1608989 19.268382
## genderM     0.2485523  0.2318810  1.071896

4.3.2. age

age is a level-2 continuous variable that in section 3.2 was uncorrelated with emotional exhaustion (r = 0.00). Consistently, its inclusion in the model including gender implies weaker evidence than more parsimonious models, with an absolute t-value lower than 2. Thus, also considering the uncertainty on its relationship with workaholism, we do not include it as a covariate for emotional exhaustion.

# model comparison vs null
m1.bis <- lmer(EE ~ age + (1|ID),data=cleanEE)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.99 0.01

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

##                 Estimate  Std. Error    t value
## (Intercept)  3.429317614 0.406974930  8.4263608
## age         -0.004969458 0.009255897 -0.5368964

# model comparison vs gender
m2 <- lmer(EE ~ gender + age + (1|ID),data=cleanEE)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.723 0.275 0.003

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

##                 Estimate  Std. Error    t value
## (Intercept)  3.281307057 0.431683535  7.6011865
## genderM      0.240100182 0.233541082  1.0280854
## age         -0.004200372 0.009285946 -0.4523365

4.3.3. workHours

workHours is a level-1 continuous variable that in section 3.2 was positively correlated with emotional exhaustion (r = 0.22). Consistently, its inclusion implies stronger evidence than more parsimonious models, with an absolute t-value higher than 2. However, considering its role as a potential explanatory mechanism (e.g., mediator) through which workaholism might impact on strain, we do not include it as a covariate for emotional exhaustion.

# model comparison vs null
m1.bis <- lmer(EE ~ workHours.mc + (1|ID),data=cleanEE)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.004 0.996

summary(m1.bis)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept)  3.2199170 0.11590941 27.779599
## workHours.mc 0.1257082 0.02941963  4.272938

# model comparison vs gender
m2 <- lmer(EE ~ gender + workHours.mc + (1|ID),data=cleanEE)
Weights(AIC(m0,m1,m2)) # Akaike weights: stronger evidence than previous model

##  model weights 
## [1] 0.012 0.004 0.984

summary(m2)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept)  3.1001459 0.16091111 19.266201
## genderM      0.2487437 0.23189176  1.072671
## workHours.mc 0.1257082 0.02941813  4.273155

4.3.4. lateWorkHours

lateWorkHours is a level-1 categorical variable that in section 3.5 did not predict substantial differences in emotional exhaustion. Consistently, its inclusion implies weaker evidence than more parsimonious models, with an absolute t-value lower than 2. Thus, also considering its role as a potential explanatory mechanism (e.g., mediator) through which workaholism might impact on strain, we do not include it as a covariate for emotional exhaustion.

# number of cases by covariate level
summary(cleanEE$lateWorkHours)

##    Mode   FALSE    TRUE 
## logical     584     190

# model comparison vs null
m1.bis <- lmer(EE ~ lateWorkHours + (1|ID),data=cleanEE)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.794 0.206

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

##                    Estimate Std. Error   t value
## (Intercept)       3.1805575  0.1190303 26.720560
## lateWorkHoursTRUE 0.1526904  0.1126787  1.355095

# model comparison vs gender
m2 <- lmer(EE ~ gender + lateWorkHours + (1|ID),data=cleanEE)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.675 0.256 0.068

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

##                    Estimate Std. Error   t value
## (Intercept)       3.0580254  0.1631294 18.746014
## genderM           0.2534415  0.2308978  1.097635
## lateWorkHoursTRUE 0.1546659  0.1126817  1.372591

4.3.5. dailyHassles_eve

dailyHassles_eve is a level-1 categorical variable that in section 3.5 predicted higher emotional exhaustion. Consistently, its inclusion implies stronger evidence (i.e., higher Akaike weight) than more parsimonious models, with an estimated absolute t-value higher than 2. However, the variable shows a strongly unbalanced number of cases, with only 10% of cases with some evening hassles. Moreover, we consider the theoretical importance of this variable, which might be a potential explanatory mechanism (e.g., mediator) through which workaholism might impact on strain. Thus, we do not include it as a covariate for emotional exhaustion.

# number of cases by covariate level
summary(cleanEE$dailyHassles_eve)

##  No Yes 
## 696  78

# model comparison dailyHassles_eve vs null
m1.bis <- lmer(EE ~ dailyHassles_eve + (1|ID),data=cleanEE)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.001 0.999

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                      Estimate Std. Error   t value
## (Intercept)         3.1576786  0.1150583 27.444156
## dailyHassles_eveYes 0.6301497  0.1437826  4.382657

# model comparison vs gender
m2 <- lmer(EE ~ gender + dailyHassles_eve + (1|ID),data=cleanEE)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

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

summary(m2)$coefficients # coefficient: |t| > 2

##                      Estimate Std. Error   t value
## (Intercept)         3.0172110  0.1592633 18.944797
## genderM             0.2901592  0.2281258  1.271926
## dailyHassles_eveYes 0.6379189  0.1438860  4.433503

4.3.6. teleWork

teleWork is a level-1 categorical variable that in section 3.5 did not predict substantial differences in emotional exhaustion. Consistently, its inclusion implies weaker evidence than more parsimonious models, with an absolute t-value lower than 2. Thus, also considering the unbalanced number of cases per level and the uncertainty of its relationship with workaholism, we do not include it as a covariate for emotional exhaustion.

# number of cases by covariate level
summary(cleanEE$teleWork) # currently 4 categories

##   office teleWork     both   dayOff 
##      580      143       51        0

cleanEE$teleWork <- as.factor(gsub("both","teleWork",gsub("dayOff","teleWork",cleanEE$teleWork))) # recoding to 2 categories
summary(cleanEE$teleWork)

##   office teleWork 
##      580      194

# model comparison vs null
m1.bis <- lmer(EE ~ teleWork + (1|ID),data=cleanEE)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.818 0.182

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

##                    Estimate Std. Error   t value
## (Intercept)       3.2568184  0.1218023 26.738560
## teleWorkteleWork -0.1465454  0.1457767 -1.005273

# model comparison vs gender + dailyHassles_eve
m2 <- lmer(EE ~ gender + teleWork + (1|ID),data=cleanEE)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.680 0.258 0.062

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

##                    Estimate Std. Error   t value
## (Intercept)       3.1316744  0.1637507 19.124646
## genderM           0.2658847  0.2327925  1.142153
## teleWorkteleWork -0.1579611  0.1460864 -1.081285

4.3.7. Other covariates

In addition to the pre-registered covariates, we consider a further level-2 covariate that in section 3.5 showed substantial trends in both emotional exhaustion and workaholism, and whose potential influence on both variables can be theoretically argued, namely having children or living with children home_child. Indeed, respondents reporting to live with their children showed lower emotional exhaustion, whereas those reporting to have children reported slightly lower workaholism (see section 3.5). Theoretically, it is plausible although not consistently reported by previous research (see Clark et al., 2016) that workaholism levels and overall investment in the work vs. family domain are influenced by family needs and presence of children in the household. Thus, both children and home_child might be a common external variable influencing both workaholism and emotional exhaustion, and we consider it as a further potential covariate.

Although both children and home_child have an acceptably balanced number of cases per level, none of them implied stronger evidence than the more parsimonious models. Thus, also in this case, we do not include these two variables as covariates for emotional exhaustion.

# children...........................................................

# number of cases by covariate level
summary(cleanEE[!duplicated(cleanEE$ID),"children"]) # children: No. participants

##  No Yes 
##  56  58

summary(cleanEE$children) # children: No. cases

##  No Yes 
## 398 376

# model comparison home_child vs null
m1.bis <- lmer(EE ~ home_child + (1|ID),data=cleanEE)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.525 0.475

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

##                 Estimate Std. Error   t value
## (Intercept)    3.3766989  0.1476390 22.871319
## home_childYes -0.3985085  0.2354001 -1.692899

# model comparison children vs gender
m2 <- lmer(EE ~ gender + children + (1|ID),data=cleanEE)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.637 0.242 0.121

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

##               Estimate Std. Error    t value
## (Intercept)  3.2623275  0.2028118 16.0854889
## genderM      0.2293104  0.2317081  0.9896521
## childrenYes -0.3024936  0.2315611 -1.3063232

# home_child...........................................................

summary(cleanEE[!duplicated(cleanEE$ID),"home_child"]) # home_child: No. participants

##  No Yes 
##  69  45

summary(cleanEE$home_child) # home_child: No. cases

##  No Yes 
## 474 300

# model comparison children vs null
m1.bis <- lmer(EE ~ children + (1|ID),data=cleanEE)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.647 0.353

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

##               Estimate Std. Error   t value
## (Intercept)  3.3800736  0.1641959 20.585610
## childrenYes -0.3170249  0.2310272 -1.372241

# model comparison home_child vs gender
m2 <- lmer(EE ~ gender + home_child + (1|ID),data=cleanEE)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.606 0.230 0.164

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

##                 Estimate Std. Error   t value
## (Intercept)    3.2727244  0.1954395 16.745456
## genderM        0.1901306  0.2336362  0.813789
## home_childYes -0.3669013  0.2389737 -1.535321

4.3.8. Time trends

As a final, not pre-registered, control, we inspect the linear time trends of the outcome variable. This is done by including the number of days since the beginning of the study as a further covariate. We can see that the inclusion of day implies stronger evidence than the null model unless the model only including gender is considered, with a t-value higher than 2. However, we cannot see a clear linear trend of EE over time. Thus, we do not include day as a covariate for emotional exhaustion, but we will consider it as a robustness check.

# plotting SD by day
boxplot(EE ~ day,data=cleanEE)

# model comparison day vs null
m1.bis <- lmer(EE ~ day + (1|ID),data=cleanEE)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.084 0.916

summary(m1.bis)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 2.98448068 0.13223588 22.569372
## day         0.04712845 0.01268082  3.716516

# model comparison day vs gender + dailyHassles_eve
m2 <- lmer(EE ~ gender + day + (1|ID),data=cleanEE)
Weights(AIC(m0,m1,m2)) # Akaike weights: same evidence than previous models

##  model weights 
## [1] 0.175 0.066 0.759

summary(m2)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 2.85932351 0.17354282 16.476184
## genderM     0.25831871 0.23214377  1.112753
## day         0.04728344 0.01268074  3.728761

4.4. Sleep disturbances

For sleep disturbances, the pre-registered potential covariates to be considered are:

1. gender, as previous evidence highlighted poorer sleep quality in women compared to men (e.g., Zeng et al., 2022) and gender has been highlighted as a potential moderator of the straining effects of trait workaholism (e.g., Bladucci et al., 2018; 2021).

2. age, as previous evidence highlighted poorer sleep quality in aging (e.g., Floyd et al., 2000); however, no clear differences have been highlighted for workaholism over age (Clark et al., 2016).

3. The time since waking up hhFromAwake, as perceived sleep disturbances might change with the time between awakening and the subsequent survey response; however, we have no clear expectation on its relationship with workaholism, and the two variables are uncorrelated (r = 0.02).

4. The daily number of worked hours workHours, as long working hours is a definitional aspect of workaholism which has been found associated with sleep disturbances (Virtanen, 2009); however, controlling for this variable might risk to artificially deflate the relationship between workaholism and sleep disturbances since the former might influence the latter through working hours. That is, working hours is a potential explanatory mechanism of the focused relationship, and including it as a control variable might bias rather than strengthening the results (see Spector et al., 2021).

5. The occurrence of late working hours in the evening lateWorkHours, consistently with what argued above for workHours; also in this case, controlling for this variable might risk to rule out a potential explanatory mechanism of the relationship between workaholism and sleep disturbances

6. Nonwork hassles in the evening dailyHassles_eve, as hassles have been related to sleep disturbances (e.g., Dahl et al., 2007) and nonwork hassles such as relationship tension have been recently related to state workaholism (Clark et al., 2021); also in this case, such an argument imply that nonwork hassles are a potential explanatory mechanism of the focused relationship.

7. Remote work teleWork, although with no clear expectations on its relationship with sleep disturbances, while we originally hypothesized that remote working might act as a moderator of the straining effect of workaholism.

Further pre-registered level-2 covariates (i.e., retrospective measures of trait workaholism, technostress creators) are not included due to our aim to focus the study on the focused level-1 variables, whereas recovery experiences will be included in all models (see Supplementary material S7). In contrast, additional (i.e., not pre-registered) level-2 covariates potentially moderating the relationship between workaholism and sleep disturbances are considered as potentially confounding factors.

Prior to evaluate the inclusion of each covariate, we prepare the data for the analyses. That is, we remove cases of missing responses in the response variable or any covariate (list-wise deletion) and we center level-1 continuous covariates on the individual mean (person-mean-centering).

# list-wise deletion
cleanSD <- as.data.frame(na.omit(clean[,c("ID","day","SD",
                                        "gender","age","workHours","hhFromAwake","lateWorkHours",
                                        "dailyHassles_eve","teleWork",
                                        "home_child","children")]))
cat("Considering",nrow(cleanSD),"complete observations from",nlevels(as.factor(as.character(cleanSD$ID))),"participants")

## Considering 768 complete observations from 114 participants

# mean centering lv-1 continuous covariates
for(Var in c("hhFromAwake","workHours")){
  prelqs <- cbind(prelqs,aggregate(cleanSD[,Var],list(cleanSD$ID),mean)[,2]) # computing individual means
  colnames(prelqs)[ncol(prelqs)] <- paste0(Var,".cm")
  cleanSD <- join(cleanSD,prelqs[,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

4.4.1. gender

gender is a categorical variable with a roughly balanced number of participants. Although in section 3.5 we only observed some slight tendency, its inclusion in the null model shows improved evidence (i.e., higher Akaike weight) compared to the null model, with an estimated t-value close to 2. Thus, also considering its theoretical relevance as a potential moderator of the straining effects of workaholism (e.g., Balducci et al. (2018)) we include gender as a covariate for sleep disturbances.

# number of cases by covariate level
summary(cleanSD[!duplicated(cleanSD$ID),"gender"]) # No. participants

##  F  M 
## 59 55

summary(cleanSD$gender) # No. cases

##   F   M 
## 410 358

# model comparison
m0 <- lmer(SD ~ (1|ID),data=cleanSD)
m1 <- lmer(SD ~ gender + (1|ID),data=cleanSD)
Weights(AIC(m0,m1)) # Akaike weights: stronger evidence including gender

##  model weights 
## [1] 0.468 0.532

summary(m1)$coefficients # coefficient: proximal to |t|=2

##               Estimate Std. Error   t value
## (Intercept)  2.7377768  0.1267734 21.595831
## genderM     -0.3588516  0.1829553 -1.961417

4.4.2. age

age is a level-2 continuous variable that in section 3.2 was unrelated with sleep disturbances (r = 0.02). Consistently, its inclusion in the model including gender implies weaker evidence than more parsimonious models, with an absolute t-value lower than 2. Thus, also considering the uncertainty on its relationship with workaholism, we do not include it as a covariate for sleep disturbances.

# model comparison vs null
m1.bis <- lmer(SD ~ age + (1|ID),data=cleanSD)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.992 0.008

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

##                Estimate  Std. Error   t value
## (Intercept) 2.402643598 0.324103465 7.4131993
## age         0.003861789 0.007368504 0.5240941

# model comparison vs gender
m2 <- lmer(SD ~ gender + age + (1|ID),data=cleanSD)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.466 0.530 0.004

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

##                 Estimate  Std. Error    t value
## (Intercept)  2.618859240 0.339967408  7.7032656
## genderM     -0.353330171 0.184161326 -1.9185905
## age          0.002758294 0.007311898  0.3772336

4.4.3. hhFromAwake

hhFromAwake is a level-1 continuous variable that in section 3.2 was unrelated to sleep disturbances (r = 0.02). Consistently, its inclusion implies weaker evidence than more parsimonious models, with an absolute t-value lower than 2. Thus, also considering the uncertainty of its relationship with both workaholism and sleep disturbances, we do not include it as a covariate for sleep disturbances.

# model comparison vs null
m1.bis <- lmer(SD ~ hhFromAwake.mc + (1|ID),data=cleanSD)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.957 0.043

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

##                   Estimate Std. Error     t value
## (Intercept)     2.56542218 0.09247725 27.74111643
## hhFromAwake.mc -0.00449394 0.04905488 -0.09161045

# model comparison vs gender
m2 <- lmer(SD ~ gender + hhFromAwake.mc + (1|ID),data=cleanSD)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.457 0.520 0.024

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

##                   Estimate Std. Error     t value
## (Intercept)     2.73777344  0.1267739 21.59572270
## genderM        -0.35883930  0.1829568 -1.96133389
## hhFromAwake.mc -0.00449394  0.0490475 -0.09162423

4.4.4. workHours

workHours is a level-1 continuous variable that in section 3.2 was uncorrelated with sleep disturbances (r = 0.02). Consistently, its inclusion implies weaker evidence than more parsimonious models, with an absolute t-value lower than 2. Thus, also considering its role as a potential explanatory mechanism (e.g., mediator) through which workaholism might impact on strain, we do not include it as a covariate for sleep disturbances.

# model comparison vs null
m1.bis <- lmer(SD ~ workHours.mc + (1|ID),data=cleanSD)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.961 0.039

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

##                 Estimate Std. Error    t value
## (Intercept)   2.56563093 0.09248917 27.7397989
## workHours.mc -0.02634331 0.03074465 -0.8568419

# model comparison vs gender
m2 <- lmer(SD ~ gender + workHours.mc + (1|ID),data=cleanSD)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.458 0.521 0.021

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

##                 Estimate Std. Error    t value
## (Intercept)   2.73777677 0.12679557 21.5920533
## genderM      -0.35841232 0.18298806 -1.9586651
## workHours.mc -0.02616138 0.03073998 -0.8510541

4.4.5. lateWorkHours

lateWorkHours is a level-1 categorical variable that in section 3.5 did not predict substantial differences in sleep disturbances. Consistently, its inclusion implies weaker evidence than more parsimonious models, with an absolute t-value lower than 2. Thus, also considering its role as a potential explanatory mechanism (e.g., mediator) through which workaholism might impact on strain, we do not include it as a covariate for sleep disturbances.

# number of cases by covariate level
summary(cleanSD$lateWorkHours)

##    Mode   FALSE    TRUE 
## logical     580     188

# model comparison vs null
m1.bis <- lmer(SD ~ lateWorkHours + (1|ID),data=cleanSD)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.699 0.301

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

##                   Estimate Std. Error   t value
## (Intercept)       2.516734 0.09642726 26.099821
## lateWorkHoursTRUE 0.190358 0.11271379  1.688861

# model comparison vs gender
m2 <- lmer(SD ~ gender + lateWorkHours + (1|ID),data=cleanSD)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.383 0.436 0.180

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

##                     Estimate Std. Error   t value
## (Intercept)        2.6873902  0.1297561 20.711094
## genderM           -0.3537025  0.1821281 -1.942053
## lateWorkHoursTRUE  0.1873724  0.1125491  1.664806

4.4.6. dailyHassles_eve

dailyHassles_eve is a level-1 categorical variable that in section 3.5 predicted slightly higher sleep disturbances. Consistently, its inclusion implies stronger evidence (i.e., higher Akaike weight) than more parsimonious models, with an estimated absolute t-value higher than 2. However, the variable shows a strongly unbalanced number of cases, with only 10% of cases with some evening hassle. Moreover, we consider the theoretical importance of this variable, which might be a potential explanatory mechanism (e.g., mediator) through which workaholism might impact on strain. Thus, we do not include it as a covariate for sleep disturbances.

# number of cases by covariate level
summary(cleanSD$dailyHassles_eve)

##  No Yes 
## 690  78

# model comparison dailyHassles_eve vs null
m1.bis <- lmer(SD ~ dailyHassles_eve + (1|ID),data=cleanSD)
Weights(AIC(m0,m1.bis)) # Akaike weights: stronger evidence than null model

##  model weights 
## [1] 0.096 0.904

summary(m1.bis)$coefficients # coefficient: |t| > 2

##                      Estimate Std. Error   t value
## (Intercept)         2.5225367 0.09252243 27.264057
## dailyHassles_eveYes 0.4270638 0.14607846  2.923523

# model comparison vs gender
m2 <- lmer(SD ~ gender + dailyHassles_eve + (1|ID),data=cleanSD)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.095 0.109 0.796

summary(m2)$coefficients # coefficient: |t| > 2

##                       Estimate Std. Error   t value
## (Intercept)          2.6833050  0.1269841 21.131025
## genderM             -0.3320883  0.1814007 -1.830689
## dailyHassles_eveYes  0.4142470  0.1461602  2.834199

4.4.7. teleWork

teleWork is a level-1 categorical variable that in section 3.5 did not predict substantial differences in sleep disturbances. Consistently, its inclusion implies weaker evidence than more parsimonious models, with an absolute t-value lower than 2. Thus, also considering the unbalanced number of cases per level and the uncertainty of its relationship with both workaholism and sleep disturbances, we do not include it as a covariate for sleep disturbances.

# number of cases by covariate level
summary(cleanSD$teleWork) # currently 4 categories

##   office teleWork     both   dayOff 
##      575      143       50        0

cleanSD$teleWork <- as.factor(gsub("both","teleWork",gsub("dayOff","teleWork",cleanSD$teleWork))) # recoding to 2 categories
summary(cleanSD$teleWork)

##   office teleWork 
##      575      193

# model comparison vs null
m1.bis <- lmer(SD ~ teleWork + (1|ID),data=cleanSD)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.875 0.125

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

##                    Estimate Std. Error    t value
## (Intercept)       2.5814746 0.09907303 26.0562796
## teleWorkteleWork -0.0637545 0.14035665 -0.4542322

# model comparison vs gender + dailyHassles_eve
m2 <- lmer(SD ~ gender + teleWork + (1|ID),data=cleanSD)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.437 0.496 0.067

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

##                     Estimate Std. Error    t value
## (Intercept)       2.74618253  0.1300248 21.1204479
## genderM          -0.35438961  0.1838099 -1.9280224
## teleWorkteleWork -0.04188835  0.1403313 -0.2984961

4.4.8. Other covariates

In addition to the pre-registered covariates, we consider a further level-2 covariate that in section 3.5 showed substantial trends in both sleep disturbances and workaholism, and whose potential influence on both variables can be theoretically argued, namely having children children or living with children home_child. Indeed, respondents reporting to live with their children reported slightly higher sleep disturbances, whereas those reporting to have children reported slightly lower workaholism (see section 3.5). Theoretically, it is plausible although not consistently reported by previous research (see Clark et al., 2016) that workaholism levels and overall investment in the work vs. family domain are influenced by family needs and presence of children in the household. On the other hand, family needs and living with children have been shown to negatively impact on sleep disturbances (e.g., Gay et al., 20040). Thus, both children and home_child might be a common external variable influencing both workaholism and sleep disturbances, and we consider it as a further potential covariate.

Although both children and home_child have an acceptably balanced number of cases per level, none of them implies stronger evidence than the more parsimonious models. Thus, also in this case, we do not include these two variables as covariates for sleep disturbances.

# children.....................................................

# number of cases by covariate level
summary(cleanSD[!duplicated(cleanSD$ID),"children"]) # children: No. participants

##  No Yes 
##  56  58

summary(cleanSD$children) # children: No. cases

##  No Yes 
## 395 373

# model comparison children vs null
m1.bis <- lmer(SD ~ children + (1|ID),data=cleanSD)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.829 0.171

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

##              Estimate Std. Error    t value
## (Intercept) 2.5078091  0.1314603 19.0765451
## childrenYes 0.1145479  0.1853899  0.6178756

# model comparison home_child vs null
m1.bis <- lmer(SD ~ home_child + (1|ID),data=cleanSD)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.637 0.363

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

##               Estimate Std. Error   t value
## (Intercept)   2.451325  0.1179382 20.784826
## home_childYes 0.290675  0.1882206  1.544332

# home_child.....................................................

summary(cleanSD[!duplicated(cleanSD$ID),"home_child"]) # home_child: No. participants

##  No Yes 
##  69  45

summary(cleanSD$home_child) # home_child: No. cases

##  No Yes 
## 470 298

# model comparison home_child vs gender
m2 <- lmer(SD ~ gender + home_child + (1|ID),data=cleanSD)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.388 0.442 0.170

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

##                 Estimate Std. Error   t value
## (Intercept)    2.6264209  0.1542313 17.029100
## genderM       -0.3216865  0.1847922 -1.740801
## home_childYes  0.2382180  0.1890599  1.260013

# model comparison children vs gender
m2 <- lmer(SD ~ gender + children + (1|ID),data=cleanSD)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous model

##  model weights 
## [1] 0.424 0.483 0.093

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

##                Estimate Std. Error   t value
## (Intercept)  2.68790172  0.1602551 16.772642
## genderM     -0.35332513  0.1838018 -1.922316
## childrenYes  0.09391617  0.1836596  0.511360

4.4.9. Time trends

As a final, not pre-registered, control, we inspect the linear time trends of the outcome variable. This is done by including the number of days since the beginning of the study as a further covariate. We can see that the inclusion of day does not imply stronger evidence than more parsimonious models, although a negative time trend (|t| > 2) is estimated for sleep disturbances. Thus, we chose to not include day as a covariate for sleep disturbances.

# plotting SD by day
boxplot(SD ~ day,data=cleanSD)

# model comparison day vs null
m1.bis <- lmer(SD ~ day + (1|ID),data=cleanSD)
Weights(AIC(m0,m1.bis)) # Akaike weights: weaker evidence than null model

##  model weights 
## [1] 0.875 0.125

summary(m1.bis)$coefficients # coefficient: |t| > 2

##               Estimate Std. Error   t value
## (Intercept) 2.42071158 0.11263873 21.490935
## day         0.02891444 0.01296795  2.229684

# model comparison day vs gender + dailyHassles_eve
m2 <- lmer(SD ~ gender + day + (1|ID),data=cleanSD)
Weights(AIC(m0,m1,m2)) # Akaike weights: weaker evidence than previous models

##  model weights 
## [1] 0.437 0.497 0.067

summary(m2)$coefficients # coefficient: |t| > 2

##                Estimate Std. Error   t value
## (Intercept)  2.59157797 0.14266096 18.165992
## genderM     -0.35197859 0.18226085 -1.931180
## day          0.02854292 0.01296185  2.202072

4.5. Summary of results

Here, we summarize the results of our covariate selection procedure:

• for SBP_aft and DBP_aft, we selected gender, age, and BMI as covariates;

• for SBP_eve and DBP_eve, we selected gender, age, and BMI as covariates;

• for EE, we selected gender as a covariate;

• for SD, we selected gender as a covariate.

5. Data export

Here, we export the formulas of the models including the covariates selected for each response variable (see section 4).

# specifying model formulas
mformulas <- c("SBP_aft ~ gender + age + BMI",
               "DBP_aft ~ gender + age + BMI",
               
               "SBP_eve ~ gender + age + BMI",
               "DBP_eve ~ gender + age + BMI",
               
               "EE ~ gender",
               
               "SD ~ gender")

# exporting model formula
save(mformulas,file="DATI/mformulas.RData")

6. Outputs

Finally, we generate and export the outputs to be included in the manuscript, namely the table with the descriptive statistics and correlations that only includes the focused variables and the selected covariates.

# selecting level-1 quantitative variables
lv1 <- c("WHLSM","SBP_aft","DBP_aft", # afternoon
         "SBP_eve","DBP_eve","EE","PD", # evening
         "SD") # next morning
lv2 <- c("age","BMI")

# computing descriptives
tab1 <- cbind(multidesc(long=clean,wide=prelqs,lv1=lv1,lv2=lv2),
              multicorr(long=clean,wide=prelqs,lv1=lv1,lv2=lv2,pvalue=TRUE,correction="BH"))
write.csv(tab1,"RESULTS/Table1.csv",row.names=FALSE)
kable(tab1)

Measure	N	Mean	ICC	WHLSM	SBP_aft	DBP_aft	SBP_eve	DBP_eve	EE	PD	SD	age	BMI
WHLSM	847	3.43 (1.54)	0.61	1	0.14***	0.14***	0.04	0.08	0.18***	-0.08*	0.13**	NA	NA
SBP_aft	843	122.18 (19.56)	0.70	0.12	1	0.46***	0.18***	0.12**	0.05	-0.02	0.09*	NA	NA
DBP_aft	843	78.17 (14.6)	0.61	0.13	0.83***	1	0.1*	0.12**	0.08*	-0.02	0.06	NA	NA
SBP_eve	841	115.58 (18.22)	0.69	0.11	0.9***	0.77***	1	0.56***	0.06	-0.1*	0.03	NA	NA
DBP_eve	841	72.96 (14.26)	0.64	0.04	0.77***	0.86***	0.87***	1	0.08*	-0.03	0.09*	NA	NA
EE	841	3.24 (1.57)	0.56	0.49***	0.24*	0.22*	0.21	0.17	1	-0.13***	0.09*	NA	NA
PD	841	4.38 (1.85)	0.35	-0.14	-0.01	-0.05	-0.07	-0.08	-0.13	1	-0.06	NA	NA
SD	842	2.57 (1.4)	0.39	0.43***	0.09	0.15	0.12	0.13	0.47***	-0.16	1	NA	NA
age	114	42.12 (12.61)	NA	-0.14	0.35***	0.35***	0.34***	0.38***	-0.05	0.06	0.02	1	NA
BMI	114	23.91 (3.53)	NA	0.06	0.33***	0.33**	0.4***	0.4***	0.06	0.07	-0.03	0.22*	1

References

• Abdel Hadi, S., Bakker, A. B., & Häusser, J. A. (2021). The role of leisure crafting for emotional exhaustion in telework during the COVID-19 pandemic. Anxiety, Stress, & Coping, 34(5), 530-544. https://doi.org/10.1080/10615806.2021.1903447

• Balducci, C., Avanzi, L., & Fraccaroli, F. (2018). The Individual “Costs” of Workaholism: An Analysis Based on Multisource and Prospective Data. Journal of Management, 44(7), 2961–2986. https://doi.org/10.1177/0149206316658348

• Brewer, E. W., & Shapard, L. (2004). Employee burnout: A meta-analysis of the relationship between age or years of experience. Human resource development review, 3(2), 102-123. https://doi.org/10.1177/1534484304263335

• Clark, M. A., Hunter, E. M., & Carlson, D. S. (2021). Hidden costs of anticipated workload for individuals and partners: Exploring the role of daily fluctuations in workaholism. Journal of Occupational Health Psychology, 26(5), 393–404. https://doi.org/10.1037/ocp0000284

• Clark, M. A., Michel, J. S., Zhdanova, L., Pui, S. Y., & Baltes, B. B. (2016). All Work and No Play? A Meta-Analytic Examination of the Correlates and Outcomes of Workaholism. Journal of Management, 42(7), 1836–1873. https://doi.org/10.1177/0149206314522301

• Dahl, R. E., & El-Sheikh, M. (2007). Considering sleep in a family context: introduction to the special issue. Journal of Family Psychology, 21(1), 1. https://psycnet.apa.org/doi/10.1037/0893-3200.21.1.1

• Floyd, J. A., Medler, S. M., Ager, J. W., & Janisse, J. J. (2000). Age‐related changes in initiation and maintenance of sleep: a meta‐analysis. Research in nursing & health, 23(2), 106-117. https://doi.org/10.1002/(SICI)1098-240X(200004)23:2%3C106::AID-NUR3%3E3.0.CO;2-A

• Gay, C. L., Lee, K. A., & Lee, S. Y. (2004). Sleep patterns and fatigue in new mothers and fathers. Biological research for nursing, 5(4), 311-318. https://doi.org/10.1177/1099800403262142

• Gregg, E. W., Cheng, Y. J., Cadwell, B. L., Imperatore, G., Williams, D. E., Flegal, K. M., … & Williamson, D. F. (2005). Secular trends in cardiovascular disease risk factors according to body mass index in US adults. Jama, 293(15), 1868-1874. https://doi.org/10.1001/jama.293.15.1868

• Klusmann, U., Aldrup, K., Schmidt, J., & Lüdtke, O. (2021). Is emotional exhaustion only the result of work experiences? A diary study on daily hassles and uplifts in different life domains. Anxiety, Stress, & Coping, 34(2), 173-190.

• Landsbergis, P. A., Dobson, M., Koutsouras, G., & Schnall, P. (2013). Job strain and ambulatory blood pressure: a meta-analysis and systematic review. American journal of public health, 103(3), e61-e71. https://doi.org/10.2105/AJPH.2012.301153

• McQuillan, M. E., Bates, J. E., Staples, A. D., & Deater-Deckard, K. (2022). A 1-year longitudinal study of the stress, sleep, and parenting of mothers of toddlers. Sleep Health, 8(1), 47-53. https://doi.org/10.1016/j.sleh.2021.08.006

• Nyklíček, I., Vingerhoets, A. J., & Van Heck, G. L. (1999). Elevated blood pressure and self-reported symptom complaints, daily hassles, and defensiveness. International journal of behavioral medicine, 6, 177-189. https://doi.org/10.1207/s15327558ijbm0602_5

• Purvanova, R. K., & Muros, J. P. (2010). Gender differences in burnout: A meta-analysis. Journal of vocational behavior, 77(2), 168-185. https://doi.org/10.1016/j.jvb.2010.04.006

• Sandberg, K., & Ji, H. (2012). Sex differences in primary hypertension. Biology of sex differences, 3(1), 1-21. https://doi.org/10.1186/2042-6410-3-7

• Seidler, A., Thinschmidt, M., Deckert, S., Then, F., Hegewald, J., Nieuwenhuijsen, K., & Riedel-Heller, S. G. (2014). The role of psychosocial working conditions on burnout and its core component emotional exhaustion–a systematic review. Journal of occupational medicine and toxicology, 9(1), 1-13.

• Spector, P. E. (2021). Mastering the use of control variables: The hierarchical iterative control (HIC) approach. Journal of Business and Psychology, 36(5), 737-750. https://doi.org/10.1007/s10869-020-09709-0

• Virtanen, M., Ferrie, J. E., Gimeno, D., Vahtera, J., Elovainio, M., Singh-Manoux, A., … & Kivimäki, M. (2009). Long working hours and sleep disturbances: the Whitehall II prospective cohort study. Sleep, 32(6), 737-745. https://doi.org/10.1093/sleep/32.6.737

• Virtanen, M., Heikkilä, K., Jokela, M., Ferrie, J. E., Batty, G. D., Vahtera, J., & Kivimäki, M. (2012). Long working hours and coronary heart disease: a systematic review and meta-analysis. American journal of epidemiology, 176(7), 586-596. https://doi.org/10.1093/aje/kws139

• Zeng, L. N., Zong, Q. Q., Yang, Y., Zhang, L., Xiang, Y. F., Ng, C. H., … & Xiang, Y. T. (2020). Gender difference in the prevalence of insomnia: a meta-analysis of observational studies. Frontiers in Psychiatry, 11, 577429. https://doi.org/10.3389/fpsyt.2020.577429

R packages

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/.

Harrell, Frank E, Jr. 2023. Hmisc: Harrell Miscellaneous. https://hbiostat.org/R/Hmisc/.

Hope, Ryan M. 2022. Rmisc: Ryan Miscellaneous. https://CRAN.R-project.org/package=Rmisc.

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

Wickham, Hadley. 2007. “Reshaping Data with the reshape Package.” Journal of Statistical Software 21 (12): 1–20. http://www.jstatsoft.org/v21/i12/.

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

———. 2020. Reshape2: Flexibly Reshape Data: A Reboot of the Reshape Package. https://github.com/hadley/reshape.

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/.