Wearable and mobile technology to characterize daily patterns of sleep, stress, pre-sleep worry, and mood in adolescent insomnia - Appendix C: Data analysis code and output
Part 2: Characterization of sleep and diary ratings
Luca Menghini, PhD, Dilara Yüksel, PhD, Devin Prouty, PhD, Fiona C Baker PhD, Christopher King, PhD, Massimiliano de Zambotti, PhD
Aims and content
The present document includes the analytical steps implemented to characterize and explore the antecedents of sleep patterns in a sample of older adolescents with (N = 47; 31 girls) and without (N = 46; 28 girls) DSM-5 insomnia.
The data processing pipeline is depicted at this page, whereas the first part of data analysis (data preparation, psychometrics, and descriptives) is depicted at this page. The former script resulted in the following long-form ema and wide-form demos datasets:
# removing all objets from the workspace
rm(list=ls())
# loading datasets exported at the end of the first Data Analysis script
load("emaFINAL.RData")
load("demosFINAL.RData")
# setting system time zone to GMT (for consistent temporal synchronization)
Sys.setenv(tz="GMT")Here, we use generalized linear mixed-effects regression (GLMER) models to examine the following relationships:
s.archit,s.timing,s.variab, ands.autonvariables are predicted by meaningful covariates (Model M1:SO.num,TotalSteps1000,weekday,BMI, andage), theinsomniagroup (M2), participants’sex(M3), and the interaction betweensexandinsomniagroups (M4), in order to characterize sleep between groups.daily diary ratings (
stress,worry, andmood) and the aggregate ratings (PsyDistandfs) are predicted by meaningful covariates (M1:TotalSteps1000andweekday), theinsomniagroup (M2), participants’sex(M3), and the interaction betweensexandinsomniagroups (M4), in order to characterize stress, worry, and mood between groups.sleep measures are predicted by diary ratings. This analyses are depicted in Part 3.
For each step and each outcome variable, models are specified by using alternative distribution families based on the nature of the outcome, inspected to evaluate model fit, compared in terms of likelihood and strength of evidence, and inspected to check the results of the selected model.
Note that all level-1 continuous predictors are centered to the individual mean (mean-centered, mc), whereas all level-2 continuous predictors are centered to the grand average (grand-mean-centered, gmc).
Finally, several robustness checks are performed to account for the arbitrariness of data processing and analytical choices, including (1) the consideration of insomnia.group instead of the insomnia variable, (2) the inclusion of covid19 as a further covariate, (3) the exclusion of participants with extreme missing data (marked with majMiss = 1), (4) the exclusion of the TotalSteps1000 variable from the analyses, (5) the use of alternative family distributions for those variables not showing optimal fit, and (6) the exclusion of diary ratings entered on the following day (marked with lateResp = 1).
The following R packages and functions are used.
# required packages
packages <- c("dplyr","Rmisc","lme4","ordinal","sjPlot","ggplot2","car","knitr","MuMIn","knitr","XML","sjPlot","reshape2")
# generate packages references
knitr::write_bib(c(.packages(), packages),"packagesAn.bib")
# # run to install missing packages
# xfun::pkg_attach2(packages, message = FALSE); rm(list=ls())
# loading mainly used packages
library(lme4); library(ordinal); library(reshape2); library(MuMIn); library(sjPlot); library(ggplot2); library(Rmisc)
glmerAn (from https://github.com/Luca-Menghini/LuMenPsy.Rfunctions )
#' @title Generalized linear (mixed-effects) regression analysis
#' @param modelType = type of model: GLM, mixed-effects (GLMER), or cumulative link mixedm odel (CLMM)
#' @param long = long-form dataset (data.frame) considered if modelType = GLMER or CLMM
#' @param wide = wide-form dataset (data.frame)
#' @param resp = name of the response variable (character)
#' @param fix.eff = character vector of names of the predictor(s)
#' @param cluster.ID = name of the cluster variable considered if modelType = GLMER or CLMM (default: "ID")
#' @param timeVar = name of the variable indexing time, considered if modelType = GLMER (defult: "ActivityDate")
#' @param REML = argument from the lme4::lmer() function, see ?lmer
#' @param ran.eff = character string indicating the random effect by using the lme4 syntax (defult: "(1|ID)")
#' @param family = character string indicating the name of the GLM(ER) family to be used in the models (default: "normal")
#' @param link =character string indicating the name of the GLM(ER) link function to be used in the models (default: "identity")
#' @param nAGQ = argument from the lme4::glmer() function, see ?glmer
#' @param mc.predictors = character vector of names of the predictors to be mean-centered
#' @param gmc.predictors = character vector of names of the predictors to be grand-mean-centered
#' @param outputs = character vector of the desired outcomes: "fit" for model diagnostics, "mComp" for model comparison, "coeff" for the coefficient table, "plotEff" for the effect plot(s), and "key.res" for the key results
#' @param mComp.baseline = character string indicating the name of the last predictor included in the baseline model to be compared with the subsequent models. If equal to NA (defult), the null model is used as the baseline model for comparison
#' @param p.adjust.method = argument from the stats::p.adjust() function (see ?p.adjust) indicating which method should be used to correct the p-values obtained from the likelihood ratio test (default: NA, for no adjustment)
#' @param coeff.models = name of the predictor(s) whose corresponding model(s) should be considered for the "coeff" output
#' @param transform = argument from sjPlot::tab_model(), see ?tab_model
#' @param coeff.models = model(s) to be considered for the "coeff" output
#' @param plot.model = character string indicating the name of the predictor(s) whose model(s) should be considered for the "plotEff" output
#' @param plot.pred = character string indicating the name of the predictor(s) to be plotted
#' @param show.data = argument from sjPlot::plot_model(), see ?plot_model
#' @param dot.size = argument from sjPlot::plot_model(), see ?plot_model
#' @param dodge = argument from sjPlot::plot_model(), see ?plot_model
#' @param dot.size = argument from sjPlot::plot_model(), see ?plot_model
#' @param ci.lv = argument from sjPlot::plot_model(), see ?plot_model
#' @param y.lim = numeric vector of length 2 indicating the limits of the y-axis
#' @param return.plot = booean. Should the plot be returned by the function? (default: FALSE)
#' @param key.model = character string indicating the name of the predictor(s) whose model(s) should be considered for the "key.res" output
#' @param key.predictor = character string indicating the name of the predictor to be considered by the "key.res" output
#' @param digits = number of digits for all numeric ouputs
#' @param messages = boolean indicating whether a message should be printed for each operation (defult: FALSE)
glmerAn <- function(modelType=c("GLM","GLMER"),long,wide,resp,fix.eff,cluster.ID="ID",timeVar="ActivityDate",REML=FALSE,
ran.eff="(1|ID)",family="normal",link="identity",nAGQ=1,mc.predictors=NA,gmc.predictors=NA,
outputs=c("fit","mComp","coeff","plotEff","key.results"),mComp.baseline=NA,p.adjust.method=NA,
coeff.models=NA,transform=NULL,plot.model=NA,plot.pred="all",
show.data=FALSE,dot.size=1,dodge=0,ci.lv=.95,y.lim=NA,return.plot=FALSE,
key.model=NA,key.predictor=NA,digits=3,messages=TRUE){
if(messages==TRUE){ cat("Running",modelType,"analysis of",resp,"...") }
# data preparation ..................................................................................................................
if(messages==TRUE){ cat("\n\nPreparing data...") }
# participants' ID recoding as factor
colnames(long)[which(colnames(long)==cluster.ID)] <- colnames(wide)[which(colnames(wide)==cluster.ID)] <- "ID"
long$ID <- as.factor(as.character(long$ID))
wide$ID <- as.factor(as.character(wide$ID))
wide <- wide[order(wide$ID),] # sorting wide by ID
colnames(long)[which(colnames(long)==timeVar)] <- "time" # time variable
long <- long[order(long$ID,long$time),] # sorting long by ID and time
if(nlevels(long$ID)!=nlevels(wide$ID)){ stop(message="Error: different No. of participants in long and wide datasets") }
# listwise deletion
if(modelType%in%c("GLMER","CLMM")){ memory <- long
for(Var in c(resp,fix.eff[which(grepl(":",fix.eff,fixed=TRUE)==FALSE)])){ colnames(long)[which(colnames(long)==Var)] <- "Var"
long <- long[!is.na(long$Var),] # removing obs. with missing response in any resp or fix.eff variable
colnames(long)[which(colnames(long)=="Var")] <- Var }
long$ID <- as.factor(as.character(long$ID))
wide <- wide[wide$ID%in%levels(long$ID),] # excluding wide IDs not included in long IDs
if(messages==TRUE){ cat("\n - Excluding",nrow(memory)-nrow(long),"incomplete observations (",
round(100*(nrow(memory)-nrow(long))/nrow(memory),1),
"% ) in the response var. or in any of the predictors,\n and",
nlevels(memory$ID)-nlevels(long$ID),"participants with no complete variables") }
} else { memory <- wide
for(Var in c(resp,fix.eff[which(grepl(":",fix.eff,fixed=TRUE)==FALSE)])){ colnames(wide)[which(colnames(wide)==Var)] <- "Var"
wide <- wide[!is.na(wide$Var),]
colnames(wide)[which(colnames(wide)=="Var")] <- Var }
if(messages==TRUE){ cat("\n - Excluding",nrow(memory)-nrow(wide),"participants with no complete variables") }}
# mean-centering predictors
if(!is.na(gmc.predictors[1])){ if(messages==TRUE){ cat("\n - Grand-mean-centering",paste(gmc.predictors,collapse=" , ")) }
for(Var in gmc.predictors){ colnames(wide)[which(colnames(wide)==Var)] <- "Var"
wide$Var.gmc <- wide$Var - mean(wide$Var,na.rm=TRUE) # computing gmc values
colnames(wide) <- gsub("Var",Var,colnames(wide))
long <- plyr::join(long,wide[,c("ID",paste(Var,"gmc",sep="."))],type="left",by="ID") # joining gmc to long dataset
fix.eff <- gsub(Var,paste(Var,"gmc",sep="."),fix.eff) }} # updating labels in fix.eff (from Var to Var.gmc)
if(!is.na(mc.predictors[1])){ if(messages==TRUE){ cat("\n - Mean-centering",paste(mc.predictors,collapse=" , ")) }
suppressMessages(suppressWarnings(require(Rmisc)))
for(Var in mc.predictors){ colnames(long)[which(colnames(long)==Var)] <- "Var"
wide$Var.cm <- suppressWarnings(summarySE(long,measurevar="Var",groupvars="ID",na.rm=TRUE)[,3]) # cluster means
long <- plyr::join(long,wide[,c("ID","Var.cm")],type="left",by="ID") # joining cm to long dataset
long$Var.mc <- long$Var - long$Var.cm # computing mc values
colnames(long) <- gsub("Var",Var,colnames(long))
fix.eff <- gsub(Var,paste(Var,"mc",sep="."),fix.eff) }} # updating labels in fix.eff (from Var to Var.mc)
long <- long[,!duplicated(colnames(long))] # removing duplicated columns (don't know why)
# modeling .........................................................................................................................
# creating model formulas
formulas <- character()
if(modelType=="GLM"){ ran.eff <- "1" }
null.f <- paste(resp,"~",ran.eff) # creating null model formula
for(i in 1:length(fix.eff)){ # creating other formulas
if(i==1){ formulas[i] <- paste(resp,"~",fix.eff[1]) } else { formulas[i] <- paste(formulas[i-1],"+",fix.eff[i]) }}
if(modelType%in%c("GLMER","CLMM")){ if(!is.na(ran.eff)){ formulas <- paste(formulas,"+",ran.eff)
if(substr(ran.eff,2,2)!="1"){ ranSlope <- paste(fix.eff[which(grepl(ran.eff,fix.eff))])[1]
null.f <- gsub(ranSlope,"1",null.f) # removing random slope from models without the related predictor
for(i in 1:length(formulas)){
if(!(grepl(ranSlope,gsub(paste(ranSlope,"[|]",sep=""),"",formulas[i])))){
formulas[i] <- gsub(paste(ranSlope,"[|]",sep=""),"1|",formulas[i]) }}}
} else { stop(message="Error: GLMER model type without ran.eff specification") }}
if(messages==TRUE){ cat("\n\nModel specification:\n - model M0 (null):",null.f)
for(i in 1:length(formulas)){ cat("\n - model M",i,": ",formulas[i],sep="")}}
# fitting models
models <- list()
if(modelType=="GLM"){ if(messages==TRUE){ cat("\n\nFitting GLM models of",resp,"on",nrow(wide),"participants \n using the",
family,"family with the",link,"link function...") }
if(family=="normal" & link=="identity"){ null.m <- lm(as.formula(null.f),data=wide) # normal family
for(i in 1:length(formulas)){ models[[i]] <- lm(formula=as.formula(formulas[i]),data=wide) }
} else if (family=="gamma") { null.m <- glm(as.formula(null.f),data=wide,family=Gamma(link=link),nAGQ=nAGQ) # gamma family
for(i in 1:length(formulas)){ models[[i]] <- glm(formula=as.formula(formulas[i]),data=wide,family=Gamma(link=link)) }
} else if(family=="normal" & link!="identity"){
null.m <- glm(as.formula(null.f),data=wide,family=gaussian(link=link)) # normal with other link functions
for(i in 1:length(formulas)){ models[[i]] <- glm(formula=as.formula(formulas[i]),data=wide,family=gaussian(link=link)) }
} else if(family=="binomial"){
null.m <- glm(as.formula(null.f),data=wide,family=binomial(link=link)) # logistic regression
for(i in 1:length(formulas)){ models[[i]] <- glm(formula=as.formula(formulas[i]),data=wide,family=binomial(link=link))}
} else { stop(message="Error: only normal, gamma, and binomial family are allowed,
with identity, inverse, and log link functions") }
} else if(modelType=="GLMER"){ suppressMessages(suppressWarnings(require(lme4)))
if(messages==TRUE){ cat("\n\nFitting",modelType,"models of",resp,"on",nrow(long),"observations from",
nlevels(as.factor(as.character(long$ID))),"participants \n using the",family,
"family with the",link,"link function using",ifelse(REML==FALSE,"ML","REML"),"estimator...") }
if(family=="normal" & link=="identity"){ null.m <- lmer(as.formula(null.f),data=long,REML=REML) # normal identity
for(i in 1:length(formulas)){ models[[i]] <- lmer(formula=as.formula(formulas[i]),data=long,REML=REML) }
} else if (family=="gamma") { null.m <- glmer(as.formula(null.f),data=long,family=Gamma(link=link),nAGQ=nAGQ) # gamma
for(i in 1:length(formulas)){ models[[i]]<-glmer(formula=as.formula(formulas[i]),data=long,family=Gamma(link=link),nAGQ=nAGQ) }
} else if(family=="normal" & link!="identity"){
null.m <- glmer(as.formula(null.f),data=long,family=gaussian(link=link),nAGQ=nAGQ) # normal with other links
for(i in 1:length(formulas)){ models[[i]] <- glmer(formula=as.formula(formulas[i]),data=long,
family=gaussian(link=link),nAGQ=nAGQ) }
} else if(family=="binomial"){ null.m <- glmer(as.formula(null.f),data=long,family=binomial(link=link),nAGQ=nAGQ) # logistic
for(i in 1:length(formulas)){ models[[i]] <- glmer(formula=as.formula(formulas[i]),data=long,
family=binomial(link=link),nAGQ=nAGQ)}
} else { stop(message="Error: only normal, logistic, and gamma family are allowed,
with identity, inverse, and log link functions") }
} else if(modelType=="CLMM"){ suppressMessages(suppressWarnings(require(ordinal))) # cumulative link mixed models
if(messages==TRUE){ cat("\n\nFitting",modelType,"models of",resp,"on",nrow(long),"observations from",
nlevels(as.factor(as.character(long$ID))),"participants \n using Cumulative Link Mixed Models") }
long[,resp] <- factor(long[,resp],ordered=TRUE) # response variable as ordered factor
null.m <- suppressWarnings(clmm(as.formula(gsub("~","~ 1 +",null.f)),data=long)) # suppress formula warning (bugged)
for(i in 1:length(formulas)){ models[[i]] <- suppressWarnings(clmm(formula=as.formula(formulas[i]),data=long,nAGQ=nAGQ)) }
} else { stop(message="Error: modelType can only be 'GLM', 'GLMER', or 'CLMM'") }
# outputs..........................................................................................................................
if(messages==TRUE){ cat("\n\nGenerating models outputs...") }
# model diagnostics
if("fit"%in%outputs & modelType!="CLMM"){ if(messages==TRUE){ cat("\n\n - Plotting diagnostics of the most complex model:") }
suppressMessages(suppressWarnings(require(sjPlot))); suppressMessages(suppressWarnings(require(ggplot2)))
fit <- models[[length(models)]] # most complex model
if(modelType=="GLMER"){ dat <- long } else { dat <- wide } # fitted dataset
plots <- list()
if(family=="normal" & link=="identity"){ p <- plot_model(fit,type="diag",dot.size=dot.size) # LM(ER) diagnostics (normal)
if(modelType=="GLMER"){ p[[2]] <- p[[2]]$ID }
suppressMessages(plot_grid(p,tags=TRUE,margin=c(0,0,0,0)))
} else { if(family=="binomial"){ # logistic regression
dat$logit <- log(predict(fit,type="response")/(1-predict(fit,type="response"))) # computing logit
for(Var in fix.eff[which(grepl(":",fix.eff,fixed=TRUE)==FALSE)]){ # plotting logit by each continuous predictor
if(class(dat[,Var])=="numeric"){ plots[[length(plots)+1]] <- ggplot(dat,aes_string(x=Var,y="logit")) +
geom_point(size=dot.size) + suppressMessages(stat_smooth(method="loess")) + ggtitle("Logit by",Var) }}
} else { plots[[length(plots)+1]] <- ggplot(data.frame(residuals=residuals(fit,type="deviance")), # deviance residuals QQ plot
aes(sample=residuals)) + stat_qq(size=dot.size) + stat_qq_line() +
labs(x="Theoretical qualtiles (normal distr.)",y="Sample quantiles") + ggtitle("Deviance Residuals' normal QQ plot")
plots[[length(plots)+1]] <- ggplot(fortify.merMod(fit), aes(.fitted, resid(fit,type="deviance"))) + geom_point() +
suppressMessages(stat_smooth(method="loess"))+geom_hline(yintercept=0, col="red", linetype="dashed") +
labs(x="Fitted values",y="Residuals") + ggtitle("Deviance residuals vs. fitted values") }
p <- plot_model(fit,type="diag",dot.size=dot.size) # random effects with GLMER
if(modelType=="GLMER"){ p <- p$ID }
plots[[length(plots)+1]] <- p + ggtitle("Random effects' normal QQ plot") }
if(family!="binomial"){ type <- ifelse(family=="gamma","deviance","pearson") # deviance residuals with gamma
for(Var in fix.eff[which(grepl(":",fix.eff,fixed=TRUE)==FALSE)]){ # homoscedasticity: residuals by levels of categorical pred.
if(class(dat[,Var])=="factor"){ Dat <- cbind(dat,residuals=resid(fit,type=type),Var=dat[,Var],Resp=dat[,resp])
plots[[length(plots)+1]] <- ggplot(Dat,aes_string(x="Var",y="residuals")) +
ggtitle(paste(type,"residuals by",Var,"\n mean =",paste(round(with(Dat, tapply(Resp, Var, mean)),1),collapse=", "),
", SD =",paste(round(with(Dat, tapply(Resp, Var, sd)),1),collapse=", "))) + xlab(Var) +
geom_point(size=dot.size,position=position_jitter(),color="gray") + geom_violin(alpha=0.3) }}}
suppressMessages(plot_grid(plots,tags=TRUE,margin=rep(0,4)))
cat("\n Printing Variance Inflation Factors (VIF):\n")
suppressMessages(suppressWarnings(require(car))) # variance inflation factors
print(vif(fit)) }
# model comparison - likelihood ratio test & Akaike weight
if("mComp"%in%outputs | "key.results"%in%outputs){
# likelihood ratio test
if(messages==TRUE){ cat("\n\n - Running likelihood ratio test:") }
suppressMessages(suppressWarnings(require(knitr))); suppressMessages(suppressWarnings(require(MuMIn)))
m.num <- 1
if(is.na(mComp.baseline)){ bsl <- null.m # selecting baseline model
} else { m.num <- grep(mComp.baseline,fix.eff) + 1
bsl <- models[[m.num - 1]] }
if(modelType!="CLMM"){ lrt <- as.data.frame(anova(bsl,models[[m.num]]))
if(length(models)>m.num){
for(i in m.num:(length(models)-1)){ lrt <- rbind(lrt,as.data.frame(anova(models[[i]],models[[i+1]]))[2,]) }}
} else { lrt <- as.data.frame(ordinal:::anova.clm(bsl,models[[m.num]])) # use anova.clm() to avoid env. issue
if(length(models)>m.num){
for(i in m.num:(length(models)-1)){ lrt <- rbind(lrt,as.data.frame(ordinal:::anova.clm(models[[i]],models[[i+1]]))[2,]) }}}
rownames(lrt) <- c(ifelse(is.na(mComp.baseline),"Null model","Baseline"),
fix.eff[m.num:length(fix.eff)])
if(!is.na(p.adjust.method)){ # p-value corrections for multiple comparison
if(messages==TRUE){ cat(" (applying",p.adjust.method,"p-values correction)")}
lrt[!is.na(lrt$`Pr(>Chisq)`),"Pr(>Chisq)"] <- p.adjust(lrt[!is.na(lrt$`Pr(>Chisq)`),"Pr(>Chisq)"],
method=p.adjust.method) }
if(messages==TRUE){ # printing selected model
if(nrow(lrt[!is.na(lrt$`Pr(>Chisq)`) & lrt$`Pr(>Chisq)` < 0.05,])==0 ){
cat("\n Selected model:",ifelse(is.na(mComp.baseline),"Null model","Baseline"))
} else { for(i in nrow(lrt):2){ if(lrt[i,ncol(lrt)] < .05){ if(messages==TRUE){ cat("\n Selected model:",rownames(lrt)[i]) }
break }}}}
if("mComp"%in%outputs){ print(kable(round(lrt,digits))) } # printing LRT table
# Akaike weights
AICs <- lrt[1:2,"AIC"] # Akaike weight
if(messages==TRUE | ("mComp"%in%outputs & !("key.results"%in%outputs))){
cat("\n\n - Computing Aw coefficients for each model vs. all previous models:",
"\n -",ifelse(is.na(mComp.baseline),"Null model","Baseline"),"vs.",fix.eff[m.num],": =",round(Weights(AICs),digits)[1],
"\n -",fix.eff[m.num],"vs.",ifelse(is.na(mComp.baseline),"Null model","Baseline"),"model =",
round(Weights(AICs),digits)[2]) }
if(nrow(lrt)>2){ for(i in 3:nrow(lrt)){ AICs <- c(AICs,lrt[i,"AIC"])
if(messages==TRUE | ("mComp"%in%outputs & !("key.results"%in%outputs))){
cat("\n -",row.names(lrt)[i],"=",round(Weights(AICs),digits)[length(AICs)]) }}}
if(messages==TRUE | ("mComp"%in%outputs & !("key.results"%in%outputs))){
cat("\n Selected model based on Aw:",row.names(lrt[which(lrt$AIC==min(lrt$AIC)),])) }
if("key.results"%in%outputs){
key <- lrt[which(grepl(key.predictor,row.names(lrt))),] # key results
if(nrow(key)>1){ key <- key[1,] }
key.results <- data.frame(sig.LRT=key[,ncol(key)]<0.05, # sig.LRT
higher.Aw=key$AIC==min(AICs[1:which(AICs==key$AIC)])) }} # higher.Aw
# estimated parameters from key.model
if("key.results"%in%outputs){
modSummary <- summary(models[[which(fix.eff==key.model)]])
modSummary <- modSummary$coefficients
if(modelType=="CLMM"){ modSummary <- modSummary[nlevels(long[,resp]):nrow(modSummary),] }
key <- modSummary[which(grepl(key.predictor,row.names(modSummary))),3][1] # taking only first coeff for key.results
key.results <- cbind(key.results,t.196=abs(key)>1.96,t=key) }
if("coeff"%in%outputs){ if(messages==TRUE){ suppressMessages(suppressWarnings(require(knitr)))
cat("\n\n - Printing estimated coefficients for models",paste(coeff.models,collapse=" , "))}
suppressMessages(suppressWarnings(require(XML))); suppressMessages(suppressWarnings(require(sjPlot)))
out <- data.frame(readHTMLTable(htmlParse(tab_model(models[which(fix.eff%in%coeff.models)],transform=transform,
show.icc=FALSE,show.p=FALSE,show.stat=TRUE,show.re.var=FALSE,
show.ci=FALSE,show.r2=FALSE,show.se=TRUE,collapse.se=TRUE,
string.est="B (SE)",string.stat="t")))[1])
colnames(out) <- out[1,]
out <- out[-1,] # sorting rows and columns
out <- rbind(out[which(!substr(out$Predictors,1,3)%in%c("SD ","Cor")),],
out[which(substr(out$Predictors,1,3)=="SD "),],out[which(substr(out$Predictors,1,4)=="Cor "),])
print(knitr::kable(out)) }
# plotting effects
if("plotEff"%in%outputs){ if(messages==TRUE){ cat("\n\n - Plotting effect(s) estimated by model",which(fix.eff==plot.model)) }
suppressMessages(suppressWarnings(require(sjPlot))); suppressMessages(suppressWarnings(require(ggplot2)))
if(plot.pred[1]=="all"){ fix.eff.plot <- fix.eff } else { fix.eff.plot <- plot.pred }
if(length(fix.eff.plot)==1){ if(grepl(":",plot.pred)){
terms <- c(strsplit(plot.pred,split=":")[[1]][1],strsplit(plot.pred,split=":")[[1]][2])
} else { terms <- fix.eff[which(fix.eff==fix.eff.plot)] }
p <- plot_model(models[[which(fix.eff==fix.eff.plot)]],type="pred",terms=terms,ci.lv=ci.lv,dodge=dodge,
show.data=show.data,dot.size=dot.size,jitter=0.15) + ggtitle(paste(resp,"by",fix.eff.plot,fix.eff.plot))
if(!is.na(y.lim[1])){ p <- p + ylim(y.lim) }
if(return.plot==TRUE){ return(p) } else { print(p) } # returning or printing the plot
} else { plots <- list()
for(i in 1:length(fix.eff.plot)){ if(grepl(":",fix.eff.plot[i],fixed=TRUE)==FALSE){ # predicted effects conditioned on ran.eff.
p <- plot_model(models[[i]],type="pred",terms=fix.eff.plot[i],show.data=show.data,dot.size=dot.size,ci.lv=ci.lv,
colors="gray",jitter=0.15,dodge=dodge)
} else { p <- plot_model(models[[i]],type="pred",terms=strsplit(fix.eff.plot[i],split=":")[[1]],ci.lv=ci.lv,dodge=dodge) }
plots[[i]] <- p + ggtitle(paste(resp,"by",fix.eff.plot[i])) }
plot_grid(plots,tags=TRUE,margin=rep(0,4)) }}
# returning key results (sig. LRT, Aw higher than previous model, t > 1.96)
if("key.results"%in%outputs){ return(key.results) }}Fits GLM(ER) models and prints the core numeric and graphical ouputs considered for the analyses.
key.resPlot (from https://github.com/Luca-Menghini/LuMenPsy.Rfunctions )
key.resPlot <- function(res=NA,robCheck=NA,levels=c(s.archit,s.timing,paste(s.variab,"c",sep=""),s.auton)){
suppressMessages(suppressWarnings(require(ggplot2)))
# setting color scale
res[res$t.196==TRUE & res$sig.LRT==TRUE & res$higher.Aw==TRUE,"Value"] <- 10 # 10 = substantial effect
res[res$t.196==TRUE & ((res$sig.LRT==FALSE & res$higher.Aw==TRUE) |
res$sig.LRT==TRUE & res$higher.Aw==FALSE),"Value"] <- 7.5 # 7.5 = substantial but LRT n.s. OR lower Aw
res[res$t.196==TRUE & res$sig.LRT==FALSE & res$higher.Aw==FALSE,"Value"] <- 5 # 5 = substantial but both LRT n.s. AND lower Aw
res[res$t.196==FALSE & ((res$sig.LRT==TRUE & res$higher.Aw==TRUE) | (res$sig.LRT==TRUE & res$higher.Aw==FALSE) |
res$sig.LRT==FALSE & res$higher.Aw==TRUE),"Value"] <- 2.5 # 2.5 = unsubstantial and LRT n.s. OR lower Aw
res[res$t.196==FALSE & res$sig.LRT==FALSE & res$higher.Aw==FALSE,"Value"] <- 1 # 1 = no evidence at all
res$Measure <- factor(res$measure,levels=levels) # Measures as factor
# plotting
if(is.na(robCheck)){ # single set of analyses
ggplot(res,aes(x=Measure,y=as.factor(1),fill=Value)) + geom_tile() + geom_text(aes(label=round(t,2))) +
scale_fill_gradient2(low="red",high="lightgreen",mid="yellow",midpoint=5,limit = c(1,10), space = "Lab") +
labs(x="",y="") + theme(legend.position = "none",axis.text.y = element_blank(),axis.text.x=element_text(angle=45))
} else { require(reshape2) # multiple robustness checks
colnames(res)[which(colnames(res)==robCheck)] <- "robCheck"
ggplot(res,aes(x=Measure, y=as.factor(robCheck), fill=Value)) + geom_tile() + geom_text(aes(label=round(t,2))) +
scale_fill_gradient2(low="red",high="lightgreen",mid="yellow",midpoint=5,limit = c(1,10), space = "Lab") +
labs(x="",y="") + theme(legend.position = "none",axis.text.x=element_text(angle=45)) }}Visualizes the key results (Aw, LRT, and t value) obtained with the glmerAn function.
1. Sleep characterization
Here, s.archit, s.timing, s.variab, and s.auton variables are predicted by (1) a baseline model including meaningful covariates (SO.num, TotalSteps1000, weekday, BMI, and age), (2) the insomnia group, (3) participants’ sex, and (4) the interaction between sex and insomnia groups, in order to characterize sleep between groups.
predictors <- c("SO.num","TotalSteps1000","weekday.sleep","BMI","age","insomnia","sex","sex:insomnia")1.1. Sleep architecture
s.archit variables are modeled by considering the subset of complete s.archit and TotalSteps1000.
s.archit <- c("TIB","TST","WASO","SE","light.p","deep.p","rem.p")TIB
TIB is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="TIB",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of TIB ...
##
## Preparing data...
## - Excluding 1909 incomplete observations ( 30.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): TIB ~ (1|ID)
## - model M1: TIB ~ SO.num.mc + (1|ID)
## - model M2: TIB ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: TIB ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: TIB ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: TIB ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: TIB ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: TIB ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: TIB ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of TIB on 4310 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.028015 1.013177 1.039421 1.044043
## age.gmc insomnia sex insomnia:sex
## 1.109634 1.736097 2.081566 2.546298
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 48862.27| 48913.22| -24423.14| 48846.27| NA| NA| NA|
## |insomnia | 9| 48864.08| 48921.40| -24423.04| 48846.08| 0.189| 1| 0.664|
## |sex | 10| 48861.30| 48924.99| -24420.65| 48841.30| 4.783| 1| 0.086|
## |sex:insomnia | 11| 48862.93| 48932.98| -24420.46| 48840.93| 0.376| 1| 0.664|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.712
## - insomnia vs. Baseline model = 0.288
## - sex = 0.537
## - sex:insomnia = 0.192
## Selected model based on Aw: sex
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:------------|:------|:------------|:------|:------------|:------|:-------------|:------|
## |2 |(Intercept) |450.87(4.52) |99.74 |448.79(6.57) |68.29 |457.00(7.38) |61.91 |459.23(8.21) |55.93 |
## |3 |SO.num.mc |-37.60(0.90) |-41.63 |-37.60(0.90) |-41.63 |-37.60(0.90) |-41.63 |-37.60(0.90) |-41.63 |
## |4 |TotalSteps1000.mc |-0.44(0.27) |-1.64 |-0.44(0.27) |-1.64 |-0.44(0.27) |-1.64 |-0.44(0.27) |-1.64 |
## |5 |weekday.sleep [weekend] |54.37(2.36) |23.08 |54.36(2.36) |23.08 |54.38(2.36) |23.08 |54.38(2.36) |23.09 |
## |6 |BMI.gmc |-0.58(1.39) |-0.42 |-0.49(1.41) |-0.34 |-0.70(1.37) |-0.51 |-0.79(1.37) |-0.57 |
## |7 |age.gmc |2.73(4.83) |0.57 |3.30(4.99) |0.66 |1.53(4.91) |0.31 |1.34(4.91) |0.27 |
## |10 |insomnia [1] | | |4.09(9.38) |0.44 |2.42(9.15) |0.26 |-1.78(11.41) |-0.16 |
## |11 |sex [M] | | | | |-20.42(9.19) |-2.22 |-26.09(13.02) |-2.00 |
## |12 |insomnia [1] * sex [M] | | | | | | |11.16(18.19) |0.61 |
## |13 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |14 |Observations |4310 |4310 |4310 |4310 |NA |NA |NA |NA |
## |8 |SD (Intercept) |41.23(NA) | |41.16(NA) | |39.89(NA) | |39.82(NA) | |
## |9 |SD (Observations) |8.24(NA) | |8.24(NA) | |8.24(NA) | |8.24(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# keeping key results
res <- cbind(measure="TIB",glmerAn(long=ema,wide=demos,resp="TIB",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="sex",messages=FALSE))Comments:
both (1) residuals and (2) random effects are quite normally distributed (despite some deviations in the tails of the residuals distribution) with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT, while model M3 including
sexshowed stronger evidence in terms of Aw compared to the baseline modela large negative effect is estimated for
SO.num(shorterTSTfor later than usualSO), whereas a large positive effect is estimated forweekday(longerTSTduring weekend), and a small difference is estimated betweensexes(longerTSTfor males than females) across models M2-M4. NeitherinsomnianorTotalSteps1000,age,BMI, or the interaction betweeninsomniaandsexshow a substantial effect
TST
TST is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="TST",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of TST ...
##
## Preparing data...
## - Excluding 1909 incomplete observations ( 30.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): TST ~ (1|ID)
## - model M1: TST ~ SO.num.mc + (1|ID)
## - model M2: TST ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: TST ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: TST ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: TST ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: TST ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: TST ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: TST ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of TST on 4310 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.028013 1.013176 1.039421 1.044328
## age.gmc insomnia sex insomnia:sex
## 1.109695 1.736916 2.082362 2.547996
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: sex
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 47845.45| 47896.40| -23914.73| 47829.45| NA| NA| NA|
## |insomnia | 9| 47847.42| 47904.74| -23914.71| 47829.42| 0.032| 1| 0.858|
## |sex | 10| 47841.99| 47905.67| -23910.99| 47821.99| 7.434| 1| 0.019|
## |sex:insomnia | 11| 47843.50| 47913.56| -23910.75| 47821.50| 0.482| 1| 0.731|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.728
## - insomnia vs. Baseline model = 0.272
## - sex = 0.805
## - sex:insomnia = 0.274
## Selected model based on Aw: sex
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:------------|:------|:------------|:------|:------------|:------|:-------------|:------|
## |2 |(Intercept) |398.65(3.87) |103.11 |397.92(5.63) |70.74 |406.62(6.22) |65.32 |408.74(6.92) |59.06 |
## |3 |SO.num.mc |-32.35(0.80) |-40.27 |-32.35(0.80) |-40.27 |-32.35(0.80) |-40.27 |-32.35(0.80) |-40.27 |
## |4 |TotalSteps1000.mc |-0.26(0.24) |-1.09 |-0.26(0.24) |-1.09 |-0.26(0.24) |-1.10 |-0.26(0.24) |-1.10 |
## |5 |weekday.sleep [weekend] |47.78(2.10) |22.81 |47.78(2.10) |22.81 |47.80(2.10) |22.81 |47.80(2.10) |22.81 |
## |6 |BMI.gmc |-0.14(1.19) |-0.12 |-0.11(1.20) |-0.09 |-0.34(1.15) |-0.30 |-0.43(1.16) |-0.37 |
## |7 |age.gmc |0.16(4.12) |0.04 |0.36(4.27) |0.08 |-1.51(4.14) |-0.37 |-1.70(4.14) |-0.41 |
## |10 |insomnia [1] | | |1.43(8.03) |0.18 |-0.32(7.71) |-0.04 |-4.32(9.62) |-0.45 |
## |11 |sex [M] | | | | |-21.62(7.75) |-2.79 |-27.04(10.97) |-2.46 |
## |12 |insomnia [1] * sex [M] | | | | | | |10.65(15.32) |0.69 |
## |13 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |14 |Observations |4310 |4310 |4310 |4310 |NA |NA |NA |NA |
## |8 |SD (Intercept) |35.13(NA) | |35.11(NA) | |33.48(NA) | |33.40(NA) | |
## |9 |SD (Observations) |7.77(NA) | |7.77(NA) | |7.77(NA) | |7.77(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# keeping key results
res <- rbind(res,
cbind(measure="TST",glmerAn(long=ema,wide=demos,resp="TST",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="sex",messages=FALSE)))Comments:
both (1) residuals and (2) random effects are quite normally distributed (despite some deviations in the tails of the residuals distribution) with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
only model 3 including
sexis associated with a significant LRT and a stronger evidence in terms of Aw compared to the baseline modela large negative effect is estimated for
SO.num(shorterTSTfor later than usualSO), whereas a large positive effect is estimated forweekday(longerTSTduring weekend), and a small but substantial effect is estimated for *sex(longerTSTfor males than females). NeitherinsomnianorTotalSteps1000,age,BMI, or the interaction betweeninsomniaandsex` show a substantial effectoverall,
TSTshows a similar pattern of results thanTIB
WASO
WASO is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="WASO",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of WASO ...
##
## Preparing data...
## - Excluding 1909 incomplete observations ( 30.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): WASO ~ (1|ID)
## - model M1: WASO ~ SO.num.mc + (1|ID)
## - model M2: WASO ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: WASO ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: WASO ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: WASO ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: WASO ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: WASO ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: WASO ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of WASO on 4310 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.028011 1.013175 1.039421 1.044579
## age.gmc insomnia sex insomnia:sex
## 1.109765 1.737660 2.083103 2.549498
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 38242.09| 38293.04| -19113.04| 38226.09| NA| NA| NA|
## |insomnia | 9| 38242.17| 38299.49| -19112.08| 38224.17| 1.922| 1| 0.497|
## |sex | 10| 38244.10| 38307.79| -19112.05| 38224.10| 0.064| 1| 0.980|
## |sex:insomnia | 11| 38246.10| 38316.16| -19112.05| 38224.10| 0.001| 1| 0.980|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.51
## - insomnia vs. Baseline model = 0.49
## - sex = 0.157
## - sex:insomnia = 0.055
## Selected model based on Aw: Baseline
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:------|:-----------|:------|:-----------|:------|:-----------|:------|
## |2 |(Intercept) |46.06(1.18) |39.08 |44.35(1.69) |26.20 |44.10(1.96) |22.55 |44.07(2.18) |20.23 |
## |3 |SO.num.mc |-5.41(0.26) |-20.50 |-5.41(0.26) |-20.50 |-5.41(0.26) |-20.50 |-5.41(0.26) |-20.50 |
## |4 |TotalSteps1000.mc |-0.09(0.08) |-1.19 |-0.09(0.08) |-1.19 |-0.09(0.08) |-1.19 |-0.09(0.08) |-1.19 |
## |5 |weekday.sleep [weekend] |7.11(0.69) |10.33 |7.11(0.69) |10.32 |7.11(0.69) |10.32 |7.11(0.69) |10.32 |
## |6 |BMI.gmc |-0.27(0.36) |-0.74 |-0.19(0.36) |-0.52 |-0.18(0.36) |-0.50 |-0.18(0.36) |-0.49 |
## |7 |age.gmc |2.69(1.25) |2.15 |3.16(1.28) |2.46 |3.21(1.30) |2.47 |3.22(1.30) |2.47 |
## |10 |insomnia [1] | | |3.37(2.41) |1.40 |3.42(2.42) |1.41 |3.47(3.03) |1.15 |
## |11 |sex [M] | | | | |0.62(2.43) |0.25 |0.68(3.45) |0.20 |
## |12 |insomnia [1] * sex [M] | | | | | | |-0.12(4.82) |-0.03 |
## |13 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |14 |Observations |4310 |4310 |4310 |4310 |NA |NA |NA |NA |
## |8 |SD (Intercept) |10.62(NA) | |10.47(NA) | |10.47(NA) | |10.47(NA) | |
## |9 |SD (Observations) |4.46(NA) | |4.46(NA) | |4.46(NA) | |4.46(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# keeping key results
res <- rbind(res,
cbind(measure="WASO",glmerAn(long=ema,wide=demos,resp="WASO",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="sex",messages=FALSE)))Comments:
whereas (1) random effects are quite normally distributed, a (2) marked deviation from normality can be seen in the upper tail of the residuals distribution. Thus, we conduct the main analyses by relying on normality, whereas the gamma will be checked in the robustness check. The residuals show (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT or stronger evidence in terms of Aw compared to the baseline model
a large negative effect is estimated for
SO.num(shorterWASOfor later than usualSO), whereas a large positive effect is estimated forweekday* (longerWASOduring weekend), and a small positive effect is estimated forage(longerWASOfor older adolescents). None of the remaining predictors, includinginsomniaand its interaction withsex, show a substantial effect
SE
SE is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="SE",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex","sex:insomnia"), # showing selected model + models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of SE ...
##
## Preparing data...
## - Excluding 1909 incomplete observations ( 30.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): SE ~ (1|ID)
## - model M1: SE ~ SO.num.mc + (1|ID)
## - model M2: SE ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: SE ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: SE ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: SE ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: SE ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: SE ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: SE ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of SE on 4310 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.028010 1.013175 1.039421 1.044690
## age.gmc insomnia sex insomnia:sex
## 1.109800 1.737997 2.083444 2.550168
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 23217.63| 23268.58| -11600.82| 23201.63| NA| NA| NA|
## |insomnia | 9| 23218.46| 23275.78| -11600.23| 23200.46| 1.167| 1| 0.420|
## |sex | 10| 23217.87| 23281.56| -11598.94| 23197.87| 2.592| 1| 0.322|
## |sex:insomnia | 11| 23219.82| 23289.87| -11598.91| 23197.82| 0.056| 1| 0.813|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.603
## - insomnia vs. Baseline model = 0.397
## - sex = 0.348
## - sex:insomnia = 0.116
## Selected model based on Aw: Baseline
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:------|:-----------|:------|:-----------|:------|:-----------|:------|
## |2 |(Intercept) |88.55(0.20) |433.92 |88.78(0.29) |301.55 |89.05(0.34) |264.83 |89.09(0.37) |237.87 |
## |3 |SO.num.mc |0.15(0.05) |3.33 |0.15(0.05) |3.33 |0.15(0.05) |3.33 |0.15(0.05) |3.33 |
## |4 |TotalSteps1000.mc |0.02(0.01) |1.56 |0.02(0.01) |1.56 |0.02(0.01) |1.56 |0.02(0.01) |1.56 |
## |5 |weekday.sleep [weekend] |-0.07(0.12) |-0.55 |-0.07(0.12) |-0.55 |-0.07(0.12) |-0.54 |-0.07(0.12) |-0.54 |
## |6 |BMI.gmc |0.08(0.06) |1.26 |0.07(0.06) |1.08 |0.06(0.06) |0.97 |0.06(0.06) |0.94 |
## |7 |age.gmc |-0.49(0.22) |-2.25 |-0.55(0.22) |-2.48 |-0.61(0.22) |-2.74 |-0.61(0.22) |-2.75 |
## |10 |insomnia [1] | | |-0.46(0.42) |-1.09 |-0.51(0.42) |-1.23 |-0.59(0.52) |-1.13 |
## |11 |sex [M] | | | | |-0.68(0.42) |-1.62 |-0.78(0.59) |-1.31 |
## |12 |insomnia [1] * sex [M] | | | | | | |0.20(0.83) |0.24 |
## |13 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |14 |Observations |4310 |4310 |4310 |4310 |NA |NA |NA |NA |
## |8 |SD (Intercept) |1.84(NA) | |1.82(NA) | |1.80(NA) | |1.80(NA) | |
## |9 |SD (Observations) |1.86(NA) | |1.86(NA) | |1.86(NA) | |1.86(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# keeping key results
res <- rbind(res,
cbind(measure="SE",glmerAn(long=ema,wide=demos,resp="SE",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) residuals and (2) random effects are quite normally distributed, although a marked deviation from normality can be seen in the lower tail of the residuals distribution (i.e., SE is bounded at 100), with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT or stronger evidence in terms of Aw compared to the baseline model
a substantial positive effect is only estimated for
SO.num(higherSEfor later than usualSO), whereas a negative relationship is estimated betweenSEandage(lowerSEfor older adolescents), whereas none of the remaining predictors shows a substantial effect
light.p
light.p is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="light.p",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","BMI.gmc","insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of light.p ...
##
## Preparing data...
## - Excluding 2490 incomplete observations ( 40 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): light.p ~ (1|ID)
## - model M1: light.p ~ SO.num.mc + (1|ID)
## - model M2: light.p ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: light.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: light.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: light.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: light.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: light.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: light.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of light.p on 3729 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.026944 1.012122 1.037327 1.043152
## age.gmc insomnia sex insomnia:sex
## 1.108687 1.733672 2.066402 2.525245
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 25065.96| 25115.75| -12524.98| 25049.96| NA| NA| NA|
## |insomnia | 9| 25065.59| 25121.60| -12523.79| 25047.59| 2.369| 1| 0.371|
## |sex | 10| 25067.47| 25129.70| -12523.73| 25047.47| 0.121| 1| 0.968|
## |sex:insomnia | 11| 25069.46| 25137.93| -12523.73| 25047.46| 0.002| 1| 0.968|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.454
## - insomnia vs. Baseline model = 0.546
## - sex = 0.176
## - sex:insomnia = 0.061
## Selected model based on Aw: insomnia
##
## - Printing estimated coefficients for models age.gmc , BMI.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:------|:-----------|:------|:-----------|:-----|:-----------|:-----|:-----------|:-----|
## |2 |(Intercept) |58.60(0.49) |119.77 |58.60(0.49) |119.95 |57.82(0.70) |82.81 |57.68(0.81) |71.36 |57.66(0.90) |64.05 |
## |3 |SO.num.mc |-0.08(0.10) |-0.83 |-0.08(0.10) |-0.83 |-0.08(0.10) |-0.83 |-0.08(0.10) |-0.83 |-0.08(0.10) |-0.83 |
## |4 |TotalSteps1000.mc |-0.06(0.03) |-1.91 |-0.06(0.03) |-1.91 |-0.06(0.03) |-1.91 |-0.06(0.03) |-1.91 |-0.06(0.03) |-1.91 |
## |5 |weekday.sleep [weekend] |-0.02(0.25) |-0.08 |-0.02(0.25) |-0.08 |-0.02(0.25) |-0.09 |-0.02(0.25) |-0.09 |-0.02(0.25) |-0.09 |
## |6 |BMI.gmc |-0.34(0.15) |-2.24 |-0.34(0.15) |-2.27 |-0.31(0.15) |-2.04 |-0.30(0.15) |-2.01 |-0.30(0.15) |-2.00 |
## |9 |age.gmc | | |0.27(0.52) |0.52 |0.49(0.53) |0.92 |0.52(0.54) |0.96 |0.52(0.54) |0.96 |
## |10 |insomnia [1] | | | | |1.55(1.00) |1.56 |1.59(1.00) |1.58 |1.62(1.25) |1.29 |
## |11 |sex [M] | | | | | | |0.35(1.01) |0.35 |0.39(1.43) |0.27 |
## |12 |insomnia [1] * sex [M] | | | | | | | | |-0.08(2.00) |-0.04 |
## |13 |N |92 ID |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |NA |
## |14 |Observations |3729 |3729 |3729 |3729 |3729 |NA |NA |NA |NA |NA |
## |7 |SD (Intercept) |4.47(NA) | |4.46(NA) | |4.38(NA) | |4.38(NA) | |4.38(NA) | |
## |8 |SD (Observations) |2.59(NA) | |2.59(NA) | |2.59(NA) | |2.59(NA) | |2.59(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# keeping key results
res <- rbind(res,
cbind(measure="light.p",glmerAn(long=ema,wide=demos,resp="light.p",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) residuals and (2) random effects are quite normally distributed (despite some deviations in the tails of the residuals distribution) with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT, whereas model M2 including
insomniashowed stronger evidence in terms of Aw compared to the baseline modela substantial negative effect is estimated only for
BMI(lowerlight.pfor participants with higherBMI), whereas none of the remaining predictors shows a substantial effect
deep.p
deep.p is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="deep.p",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex","sex:insomnia"), # selected model + models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of deep.p ...
##
## Preparing data...
## - Excluding 2490 incomplete observations ( 40 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): deep.p ~ (1|ID)
## - model M1: deep.p ~ SO.num.mc + (1|ID)
## - model M2: deep.p ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: deep.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: deep.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: deep.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: deep.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: deep.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: deep.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of deep.p on 3729 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.026937 1.012119 1.037330 1.043784
## age.gmc insomnia sex insomnia:sex
## 1.108444 1.735302 2.064357 2.524461
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 22128.94| 22178.74| -11056.47| 22112.94| NA| NA| NA|
## |insomnia | 9| 22130.45| 22186.46| -11056.22| 22112.45| 0.498| 1| 0.813|
## |sex | 10| 22132.21| 22194.45| -11056.11| 22112.21| 0.232| 1| 0.813|
## |sex:insomnia | 11| 22134.16| 22202.62| -11056.08| 22112.16| 0.056| 1| 0.813|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.679
## - insomnia vs. Baseline model = 0.321
## - sex = 0.117
## - sex:insomnia = 0.042
## Selected model based on Aw: Baseline
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:-----|:-----------|:-----|:-----------|:-----|:-----------|:-----|
## |2 |(Intercept) |19.91(0.28) |70.88 |20.12(0.41) |49.61 |20.23(0.47) |43.16 |20.29(0.52) |38.88 |
## |3 |SO.num.mc |0.32(0.07) |4.89 |0.32(0.07) |4.89 |0.32(0.07) |4.89 |0.32(0.07) |4.89 |
## |4 |TotalSteps1000.mc |0.02(0.02) |1.25 |0.02(0.02) |1.25 |0.02(0.02) |1.25 |0.02(0.02) |1.25 |
## |5 |weekday.sleep [weekend] |-0.30(0.17) |-1.74 |-0.30(0.17) |-1.74 |-0.30(0.17) |-1.73 |-0.30(0.17) |-1.73 |
## |6 |BMI.gmc |0.11(0.09) |1.30 |0.10(0.09) |1.18 |0.10(0.09) |1.15 |0.10(0.09) |1.12 |
## |7 |age.gmc |-0.20(0.30) |-0.68 |-0.26(0.31) |-0.84 |-0.28(0.31) |-0.91 |-0.29(0.31) |-0.92 |
## |10 |insomnia [1] | | |-0.41(0.58) |-0.71 |-0.44(0.58) |-0.75 |-0.54(0.73) |-0.74 |
## |11 |sex [M] | | | | |-0.28(0.58) |-0.48 |-0.42(0.83) |-0.51 |
## |12 |insomnia [1] * sex [M] | | | | | | |0.27(1.16) |0.24 |
## |13 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |14 |Observations |3729 |3729 |3729 |3729 |NA |NA |NA |NA |
## |8 |SD (Intercept) |2.52(NA) | |2.51(NA) | |2.51(NA) | |2.50(NA) | |
## |9 |SD (Observations) |2.13(NA) | |2.13(NA) | |2.13(NA) | |2.13(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# keeping key results
res <- rbind(res,
cbind(measure="deep.p",glmerAn(long=ema,wide=demos,resp="deep.p",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) residuals and (2) random effects are quite normally distributed (despite some deviations in the tails of the residuals distribution) with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT or stronger evidence in terms of Aw compared to the baseline model
a substantial positive effect is only estimated for
SO.num(higherdeep.pfor later than usual SO), whereas none of the remaining predictors shows a substantial effect
rem.p
rem.p is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="rem.p",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of rem.p ...
##
## Preparing data...
## - Excluding 2490 incomplete observations ( 40 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): rem.p ~ (1|ID)
## - model M1: rem.p ~ SO.num.mc + (1|ID)
## - model M2: rem.p ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: rem.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: rem.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: rem.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: rem.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: rem.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: rem.p ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of rem.p on 3729 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.026950 1.012125 1.037324 1.042593
## age.gmc insomnia sex insomnia:sex
## 1.109050 1.732412 2.069129 2.526471
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 22477.54| 22527.34| -11230.77| 22461.54| NA| NA| NA|
## |insomnia | 9| 22477.91| 22533.93| -11229.96| 22459.91| 1.633| 1| 0.604|
## |sex | 10| 22479.90| 22542.14| -11229.95| 22459.90| 0.015| 1| 0.937|
## |sex:insomnia | 11| 22481.89| 22550.35| -11229.94| 22459.89| 0.006| 1| 0.937|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.546
## - insomnia vs. Baseline model = 0.454
## - sex = 0.144
## - sex:insomnia = 0.05
## Selected model based on Aw: Baseline
##
## - Printing estimated coefficients for models insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:-----|:-----------|:-----|:-----------|:-----|
## |2 |(Intercept) |22.05(0.60) |36.51 |22.09(0.70) |31.58 |22.07(0.78) |28.32 |
## |3 |SO.num.mc |-0.24(0.07) |-3.53 |-0.24(0.07) |-3.53 |-0.24(0.07) |-3.53 |
## |4 |TotalSteps1000.mc |0.03(0.02) |1.51 |0.03(0.02) |1.51 |0.03(0.02) |1.51 |
## |5 |weekday.sleep [weekend] |0.32(0.18) |1.81 |0.32(0.18) |1.81 |0.32(0.18) |1.81 |
## |6 |BMI.gmc |0.20(0.13) |1.54 |0.20(0.13) |1.53 |0.20(0.13) |1.53 |
## |7 |age.gmc |-0.21(0.46) |-0.46 |-0.22(0.47) |-0.48 |-0.22(0.47) |-0.47 |
## |8 |insomnia [1] |-1.11(0.86) |-1.29 |-1.12(0.87) |-1.29 |-1.07(1.08) |-0.99 |
## |11 |sex [M] | | |-0.11(0.87) |-0.12 |-0.04(1.24) |-0.03 |
## |12 |insomnia [1] * sex [M] | | | | |-0.14(1.73) |-0.08 |
## |13 |N |92 ID |92 ID |92 ID |NA |NA |NA |
## |14 |Observations |3729 |3729 |3729 |NA |NA |NA |
## |9 |SD (Intercept) |3.84(NA) | |3.84(NA) | |3.84(NA) | |
## |10 |SD (Observations) |2.17(NA) | |2.17(NA) | |2.17(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# keeping key results
res <- rbind(res,
cbind(measure="rem.p",glmerAn(long=ema,wide=demos,resp="rem.p",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) residuals and (2) random effects are quite normally distributed (despite some deviations in the tails of the residuals distribution) with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT or stronger evidence in terms of Aw compared to the baseline model
a substantial negative effect is only estimated for
SO.num(higherrem.pfor earlier than usualSO), whereas none of the remaining predictors shows a substantial effectoverall,
rem.pshows a similar pattern of results thandeep.p
1.2. Sleep timing
s.timing variables are modeled by considering the subset of complete s.timing and TotalSteps1000.
s.timing <- c("SO.num","WakeUp.num")SO.num
SO.num is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals. Obviously, SO.num is not included as a covariate in model M1.
glmerAn(long=ema,wide=demos,resp="SO.num",fix.eff=predictors[2:length(predictors)],
modelType=c("GLMER"),mc.predictors="TotalSteps1000",gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of SO.num ...
##
## Preparing data...
## - Excluding 1909 incomplete observations ( 30.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering TotalSteps1000
##
## Model specification:
## - model M0 (null): SO.num ~ (1|ID)
## - model M1: SO.num ~ TotalSteps1000.mc + (1|ID)
## - model M2: SO.num ~ TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M3: SO.num ~ TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M4: SO.num ~ TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M5: SO.num ~ TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M6: SO.num ~ TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M7: SO.num ~ TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of SO.num on 4310 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## TotalSteps1000.mc weekday.sleep BMI.gmc age.gmc
## 1.011562 1.011578 1.042749 1.109673
## insomnia sex insomnia:sex
## 1.732908 2.078852 2.538633
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: sex
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 7| 13835.04| 13879.62| -6910.521| 13821.04| NA| NA| NA|
## |insomnia | 8| 13836.94| 13887.89| -6910.471| 13820.94| 0.098| 1| 0.754|
## |sex | 9| 13832.40| 13889.71| -6907.197| 13814.40| 6.548| 1| 0.032|
## |sex:insomnia | 10| 13832.68| 13896.36| -6906.339| 13812.68| 1.717| 1| 0.285|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.721
## - insomnia vs. Baseline model = 0.279
## - sex = 0.73
## - sex:insomnia = 0.388
## Selected model based on Aw: sex
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:-----|:-----------|:-----|:-----------|:-----|:-----------|:-----|
## |2 |(Intercept) |0.42(0.11) |3.74 |0.46(0.17) |2.79 |0.22(0.18) |1.19 |0.10(0.20) |0.50 |
## |3 |TotalSteps1000.mc |-0.01(0.00) |-2.60 |-0.01(0.00) |-2.60 |-0.01(0.00) |-2.60 |-0.01(0.00) |-2.60 |
## |4 |weekday.sleep [weekend] |0.43(0.04) |10.76 |0.43(0.04) |10.76 |0.43(0.04) |10.75 |0.43(0.04) |10.75 |
## |5 |BMI.gmc |-0.08(0.04) |-2.36 |-0.08(0.04) |-2.39 |-0.08(0.03) |-2.29 |-0.07(0.03) |-2.18 |
## |6 |age.gmc |-0.20(0.12) |-1.63 |-0.21(0.13) |-1.65 |-0.16(0.12) |-1.27 |-0.15(0.12) |-1.19 |
## |9 |insomnia [1] | | |-0.07(0.24) |-0.31 |-0.02(0.23) |-0.10 |0.20(0.28) |0.70 |
## |10 |sex [M] | | | | |0.60(0.23) |2.61 |0.90(0.32) |2.80 |
## |11 |insomnia [1] * sex [M] | | | | | | |-0.59(0.45) |-1.32 |
## |12 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |13 |Observations |4310 |4310 |4310 |4310 |NA |NA |NA |NA |
## |7 |SD (Intercept) |1.06(NA) | |1.06(NA) | |1.02(NA) | |1.01(NA) | |
## |8 |SD (Observations) |1.08(NA) | |1.08(NA) | |1.08(NA) | |1.08(NA) | |
##
##
## - Plotting effect(s) estimated by model 5
# keeping key results
res <- rbind(res,
cbind(measure="SO.num",glmerAn(long=ema,wide=demos,resp="SO.num",fix.eff=predictors[2:length(predictors)],
modelType=c("GLMER"),mc.predictors="TotalSteps1000",gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="sex",messages=FALSE)))Comments:
both (1) residuals and (2) random effects are quite normally distributed, despite some marked deviations in the tails of the residuals distribution, with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
only model M3 including
sexshowed a significant LRT, and higher Aw than the baseline modela relatively large effect is estimated for
weekday(laterSOduring weekend), whereas small but substantial effects are estimated forTotalSteps1000(earlierSOin days with more steps than usual),BMI(earlierSOfor participants with higherBMI), andsex(laterSOfor boys than for girls). Neitherinsomnianor its interaction withsexshow a substantial effect.
WakeUp.num
WakeUp.num is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="WakeUp.num",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors="TotalSteps1000",gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of WakeUp.num ...
##
## Preparing data...
## - Excluding 1909 incomplete observations ( 30.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering TotalSteps1000
##
## Model specification:
## - model M0 (null): WakeUp.num ~ (1|ID)
## - model M1: WakeUp.num ~ SO.num + (1|ID)
## - model M2: WakeUp.num ~ SO.num + TotalSteps1000.mc + (1|ID)
## - model M3: WakeUp.num ~ SO.num + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: WakeUp.num ~ SO.num + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: WakeUp.num ~ SO.num + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: WakeUp.num ~ SO.num + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: WakeUp.num ~ SO.num + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: WakeUp.num ~ SO.num + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of WakeUp.num on 4310 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.037628 1.013078 1.037629 1.047510
## age.gmc insomnia sex insomnia:sex
## 1.110853 1.738586 2.094303 2.553597
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 13578.09| 13629.04| -6781.046| 13562.09| NA| NA| NA|
## |insomnia | 9| 13580.08| 13637.40| -6781.039| 13562.08| 0.015| 1| 0.904|
## |sex | 10| 13582.03| 13645.72| -6781.016| 13562.03| 0.046| 1| 0.904|
## |sex:insomnia | 11| 13583.37| 13653.43| -6780.686| 13561.37| 0.659| 1| 0.904|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.73
## - insomnia vs. Baseline model = 0.27
## - sex = 0.092
## - sex:insomnia = 0.045
## Selected model based on Aw: Baseline
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:------|:-----------|:-----|:-----------|:-----|:-----------|:-----|
## |2 |(Intercept) |7.74(0.07) |118.35 |7.73(0.09) |81.46 |7.72(0.11) |70.56 |7.67(0.12) |63.31 |
## |3 |SO.num |0.39(0.01) |26.34 |0.39(0.01) |26.34 |0.39(0.01) |26.27 |0.39(0.01) |26.24 |
## |4 |TotalSteps1000.mc |-0.01(0.00) |-1.28 |-0.01(0.00) |-1.28 |-0.01(0.00) |-1.28 |-0.01(0.00) |-1.28 |
## |5 |weekday.sleep [weekend] |0.91(0.04) |23.04 |0.91(0.04) |23.04 |0.91(0.04) |23.04 |0.91(0.04) |23.04 |
## |6 |BMI.gmc |-0.06(0.02) |-2.80 |-0.06(0.02) |-2.74 |-0.06(0.02) |-2.72 |-0.05(0.02) |-2.64 |
## |7 |age.gmc |-0.07(0.07) |-1.07 |-0.07(0.07) |-1.00 |-0.07(0.07) |-0.96 |-0.07(0.07) |-0.91 |
## |10 |insomnia [1] | | |0.02(0.13) |0.12 |0.02(0.14) |0.14 |0.10(0.17) |0.60 |
## |11 |sex [M] | | | | |0.03(0.14) |0.22 |0.14(0.19) |0.73 |
## |12 |insomnia [1] * sex [M] | | | | | | |-0.22(0.27) |-0.81 |
## |13 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |14 |Observations |4310 |4310 |4310 |4310 |NA |NA |NA |NA |
## |8 |SD (Intercept) |0.58(NA) | |0.58(NA) | |0.58(NA) | |0.58(NA) | |
## |9 |SD (Observations) |1.07(NA) | |1.07(NA) | |1.07(NA) | |1.07(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# keeping key results
res <- rbind(res,
cbind(measure="WakeUp.num",glmerAn(long=ema,wide=demos,resp="WakeUp.num",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) residuals and (2) random effects are quite normally distributed, despite some marked deviations in the tails of the residuals distribution, with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values
none of the specified models showed a significant LRT or stronger evidence in terms of Aw compared to the baseline model
a large positive effect is estimated for
SO.num(laterWakeUptime for later than usualSO) andweekday(laterWakeUptime during weekend), whereas a small but substantial effect is estimated forBMI(earlierWakeUptime for participants with higherBMI). Neitherinsomnianor its interaction withsexshow a substantial effect.
1.3. Daily variability
s.variab variables are modeled by considering the subset of complete s.timing and TotalSteps1000, and only those days with nonmissing current and previous-day sleep variables.
s.variab <- paste(c(s.archit[1:3],s.timing),"SSD",sep=".")TIB.SSD
All s.variab variables are positively skewed and bounded at zero, thus they are modeled by assuming a gamma distribution with a log link function. A constant is added to all TIB.SSD values to adhere with the gamma assumption of positive-only values (i.e., 10 cases have TIB.SSD = 0).
# adding constant to each value
ema$TIB.SSDc <- ema$TIB.SSD + 0.001
# running analysis
glmerAn(long=ema,wide=demos,resp="TIB.SSDc",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
family="gamma",link="log", # Gamma family with log link function
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex","sex:insomnia"),transform="exp", # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of TIB.SSDc ...
##
## Preparing data...
## - Excluding 2094 incomplete observations ( 33.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): TIB.SSDc ~ (1|ID)
## - model M1: TIB.SSDc ~ SO.num.mc + (1|ID)
## - model M2: TIB.SSDc ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: TIB.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: TIB.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: TIB.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: TIB.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: TIB.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: TIB.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of TIB.SSDc on 4125 observations from 92 participants
## using the gamma family with the log link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.003220 1.026918 1.025246 1.048189
## age.gmc insomnia sex insomnia:sex
## 1.117198 1.755512 2.116309 2.600972
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 83036.86| 83087.46| -41510.43| 83020.86| NA| NA| NA|
## |insomnia | 9| 83034.37| 83091.29| -41508.18| 83016.37| 4.494| 1| 0.102|
## |sex | 10| 83035.96| 83099.21| -41507.98| 83015.96| 0.405| 1| 0.715|
## |sex:insomnia | 11| 83037.83| 83107.40| -41507.92| 83015.83| 0.134| 1| 0.715|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.223
## - insomnia vs. Baseline model = 0.777
## - sex = 0.259
## - sex:insomnia = 0.092
## Selected model based on Aw: insomnia
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:----------------|:------|:----------------|:------|:-----------------|:-----|:-----------------|:-----|
## |2 |(Intercept) |10051.48(609.74) |151.92 |11476.56(979.43) |109.54 |11129.60(1088.80) |95.24 |10936.62(1190.76) |85.41 |
## |3 |SO.num.mc |0.97(0.01) |-1.86 |0.97(0.01) |-1.89 |0.97(0.01) |-1.91 |0.97(0.01) |-1.91 |
## |4 |TotalSteps1000.mc |1.00(0.01) |-0.54 |1.00(0.01) |-0.54 |1.00(0.01) |-0.53 |1.00(0.01) |-0.53 |
## |5 |weekday.sleep [weekend] |1.47(0.08) |7.36 |1.47(0.08) |7.37 |1.47(0.08) |7.37 |1.47(0.08) |7.35 |
## |6 |BMI.gmc |1.06(0.02) |3.02 |1.05(0.02) |2.74 |1.05(0.02) |2.79 |1.05(0.02) |2.82 |
## |7 |age.gmc |1.05(0.07) |0.79 |1.01(0.06) |0.19 |1.02(0.07) |0.29 |1.02(0.07) |0.31 |
## |10 |insomnia [1] | | |0.77(0.09) |-2.15 |0.78(0.09) |-2.11 |0.80(0.12) |-1.46 |
## |11 |sex [M] | | | | |1.08(0.13) |0.64 |1.13(0.20) |0.71 |
## |12 |insomnia [1] * sex [M] | | | | | | |0.92(0.22) |-0.37 |
## |13 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |14 |Observations |4125 |4125 |4125 |4125 |NA |NA |NA |NA |
## |8 |SD (Intercept) |2.21(NA) | |2.15(NA) | |2.14(NA) | |2.14(NA) | |
## |9 |SD (Observations) |3.44(NA) | |3.44(NA) | |3.44(NA) | |3.44(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# plotting effects without observed data (zoom in)
glmerAn(long=ema,wide=demos,resp="TIB.SSDc",fix.eff=predictors,
mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),modelType="GLMER",family="gamma",link="log",
outputs="plotEff",plot.model="insomnia",dot.size=2,message=FALSE)
# keeping key results
res <- rbind(res,
cbind(measure="TIB.SSDc",glmerAn(long=ema,wide=demos,resp="TIB.SSDc",fix.eff=predictors,
modelType=c("GLMER"),family="gamma",link="log",transform="exp",
mComp.baseline="age",p.adjust.method="BH",
mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) deviance residuals and (2) random effects are quite normally distributed (despite some deviations in the lower tails of the residuals distribution) with (3) no substantial linear relationship between residuals and fitted values; although categorical factors do not show substantially different dispersion, (4) the SD is slightly proportional to the mean (as assumed by the gamma distribution); there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT, whereas the inclusion of
insomniais associated with a stronger evidence in terms of Awsubstantial effects are estimated for
weekday.sleep(higherTIB.SSDduring weekend),BMI(higherTIB.SSDin participants with higherBMI), andinsomnia(lowerTST.SSDin theinsomniacompared to the control group. Neithersexnor its interaction withinsomniashow a substantial effect
TST.SSD
All s.variab variables are positively skewed and bounded at zero, thus they are modeled by assuming a gamma distribution with a log link function. A constant is added to all TST.SSD values to adhere with the gamma assumption of positive-only values (i.e., 10 cases have TST.SSD = 0).
# adding constant to each value
ema$TST.SSDc <- ema$TST.SSD + 0.001
# running analysis
glmerAn(long=ema,wide=demos,resp="TST.SSDc",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
family="gamma",link="log", # Gamma family with log link function
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex","sex:insomnia"),transform="exp", # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of TST.SSDc ...
##
## Preparing data...
## - Excluding 2094 incomplete observations ( 33.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): TST.SSDc ~ (1|ID)
## - model M1: TST.SSDc ~ SO.num.mc + (1|ID)
## - model M2: TST.SSDc ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: TST.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: TST.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: TST.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: TST.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: TST.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: TST.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of TST.SSDc on 4125 observations from 92 participants
## using the gamma family with the log link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.003403 1.025041 1.023044 1.048177
## age.gmc insomnia sex insomnia:sex
## 1.118339 1.759746 2.118710 2.608996
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 81388.64| 81439.24| -40686.32| 81372.64| NA| NA| NA|
## |insomnia | 9| 81386.46| 81443.38| -40684.23| 81368.46| 4.182| 1| 0.123|
## |sex | 10| 81388.44| 81451.69| -40684.22| 81368.44| 0.015| 1| 0.904|
## |sex:insomnia | 11| 81390.35| 81459.92| -40684.17| 81368.35| 0.091| 1| 0.904|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.251
## - insomnia vs. Baseline model = 0.749
## - sex = 0.217
## - sex:insomnia = 0.077
## Selected model based on Aw: insomnia
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:---------------|:------|:---------------|:------|:---------------|:-----|:---------------|:-----|
## |2 |(Intercept) |8010.52(479.17) |150.26 |9088.97(767.10) |108.00 |9142.09(886.44) |94.06 |9012.02(973.09) |84.34 |
## |3 |SO.num.mc |0.98(0.01) |-1.63 |0.98(0.01) |-1.65 |0.98(0.01) |-1.65 |0.98(0.01) |-1.65 |
## |4 |TotalSteps1000.mc |1.00(0.01) |-0.40 |1.00(0.01) |-0.40 |1.00(0.01) |-0.40 |1.00(0.01) |-0.40 |
## |5 |weekday.sleep [weekend] |1.46(0.08) |7.45 |1.47(0.08) |7.46 |1.47(0.08) |7.46 |1.47(0.08) |7.45 |
## |6 |BMI.gmc |1.05(0.02) |2.87 |1.05(0.02) |2.59 |1.05(0.02) |2.57 |1.05(0.02) |2.59 |
## |7 |age.gmc |1.03(0.06) |0.48 |0.99(0.06) |-0.10 |0.99(0.06) |-0.12 |0.99(0.06) |-0.10 |
## |10 |insomnia [1] | | |0.78(0.09) |-2.07 |0.78(0.09) |-2.08 |0.80(0.12) |-1.48 |
## |11 |sex [M] | | | | |0.99(0.12) |-0.12 |1.02(0.18) |0.13 |
## |12 |insomnia [1] * sex [M] | | | | | | |0.93(0.22) |-0.30 |
## |13 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |14 |Observations |4125 |4125 |4125 |4125 |NA |NA |NA |NA |
## |8 |SD (Intercept) |2.14(NA) | |2.09(NA) | |2.09(NA) | |2.09(NA) | |
## |9 |SD (Observations) |3.39(NA) | |3.39(NA) | |3.39(NA) | |3.39(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# plotting effects without observed data (zoom in)
glmerAn(long=ema,wide=demos,resp="TST.SSDc",fix.eff=predictors,
mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),modelType="GLMER",family="gamma",link="log",
outputs="plotEff",plot.model="insomnia",dot.size=2,messages=FALSE)
# keeping key results
res <- rbind(res,
cbind(measure="TST.SSDc",glmerAn(long=ema,wide=demos,resp="TST.SSDc",fix.eff=predictors,
modelType=c("GLMER"),family="gamma",link="log",transform="exp",
mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) deviance residuals and (2) random effects are quite normally distributed (despite some deviations in the lower tail of the residuals distribution) with (2) no substantial linear relationship between residuals and fitted values; although categorical factors do not show substantially different dispersion, (4) the SD are proportional to the means (as assumed by the gamma distribution); there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT, whereas the inclusion of
insomniais associated with a stronger evidence in terms of Awsubstantial effects are estimated for
weekday.sleep(higherTST.SSDduring weekend),BMI(higherTST.SSDin participants with higherBMI), andinsomnia(lowerTST.SSDin theinsomniacompared to the control group. Neithersexnor its interaction withinsomniashow a substantial effectresults are highly consistent with those reported for
TIB.SSD
WASO.SSD
All s.variab variables are positively skewed and bounded at zero, thus they are modeled by assuming a gamma distribution with a log link function. A constant is added to all WASO.SSD values to adhere with the gamma assumption of positive-only values (i.e., 10 cases have WASO.SSD = 0).
# adding constant to each value
ema$WASO.SSDc <- ema$WASO.SSD + 0.001
# running analysis
glmerAn(long=ema,wide=demos,resp="WASO.SSDc",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
family="gamma",link="log",nAGQ=0, # Gamma family with log link function, nAGQ = 0 to reach convergence
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex"),transform="exp", # selected model + models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of WASO.SSDc ...
##
## Preparing data...
## - Excluding 2094 incomplete observations ( 33.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): WASO.SSDc ~ (1|ID)
## - model M1: WASO.SSDc ~ SO.num.mc + (1|ID)
## - model M2: WASO.SSDc ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: WASO.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: WASO.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: WASO.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: WASO.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: WASO.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: WASO.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of WASO.SSDc on 4125 observations from 92 participants
## using the gamma family with the log link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.028313 1.013425 1.039840 1.046237
## age.gmc insomnia sex insomnia:sex
## 1.111966 1.749233 2.092896 2.575139
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 59054.60| 59105.20| -29519.30| 59038.60| NA| NA| NA|
## |insomnia | 9| 59056.63| 59113.56| -29519.32| 59038.63| 0.000| 1| 1.000|
## |sex | 10| 59058.63| 59121.88| -29519.31| 59038.63| 0.005| 1| 1.000|
## |sex:insomnia | 11| 59058.94| 59128.52| -29518.47| 59036.94| 1.682| 1| 0.584|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.734
## - insomnia vs. Baseline model = 0.266
## - sex = 0.089
## - sex:insomnia = 0.071
## Selected model based on Aw: Baseline
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-------------|:-----|:-------------|:-----|:--------------|:-----|
## |2 |(Intercept) |583.14(60.26) |61.63 |580.59(87.11) |42.42 |578.84(100.14) |36.77 |
## |3 |SO.num.mc |0.93(0.02) |-2.90 |0.93(0.02) |-2.90 |0.93(0.02) |-2.90 |
## |4 |TotalSteps1000.mc |1.00(0.01) |-0.38 |1.00(0.01) |-0.38 |1.00(0.01) |-0.38 |
## |5 |weekday.sleep [weekend] |1.36(0.09) |4.58 |1.36(0.09) |4.58 |1.36(0.09) |4.58 |
## |6 |BMI.gmc |1.01(0.03) |0.40 |1.01(0.03) |0.40 |1.01(0.03) |0.40 |
## |7 |age.gmc |1.07(0.12) |0.60 |1.07(0.12) |0.59 |1.07(0.12) |0.59 |
## |10 |insomnia [1] | | |1.01(0.22) |0.04 |1.01(0.22) |0.04 |
## |11 |sex [M] | | | | |1.01(0.22) |0.03 |
## |12 |N |92 ID |92 ID |92 ID |NA |NA |NA |
## |13 |Observations |4125 |4125 |4125 |NA |NA |NA |
## |8 |SD (Intercept) |2.49(NA) | |2.49(NA) | |2.49(NA) | |
## |9 |SD (Observations) |3.98(NA) | |3.98(NA) | |3.98(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# plotting effects without observed data (zoom in)
glmerAn(long=ema,wide=demos,resp="WASO.SSDc",fix.eff=predictors,
mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),modelType="GLMER",family="gamma",link="log",nAGQ=0,
outputs="plotEff",plot.model="insomnia",dot.size=2,messages=FALSE)
# keeping key results
res <- rbind(res,
cbind(measure="WASO.SSDc",glmerAn(long=ema,wide=demos,resp="WASO.SSDc",fix.eff=predictors,
modelType=c("GLMER"),family="gamma",link="log",transform="exp",nAGQ=0,
mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) deviance residuals and (2) random effects are quite normally distributed (despite some small deviations in the lower tail and some marked deviations in the upper tail of the residuals distribution) with (3) no substantial linear relationship between residuals and fitted values; although categorical factors do not show substantially different dispersion, the (4) SD are slightly proportional to the means (as assumed by the gamma distribution); there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT or stronger evidence in terms of Aw compared to the baseline model
substantial effects are only estimated for
SO.num(higher variability for earlier than usualSO) andweekday.sleep(higherWASO.SSDduring weekend). Neitherinsomnianor its interaction withsexshow a substantial effect
SO.num.SSD
All s.variab variables are positively skewed and bounded at zero, thus they are initially modeled by assuming a gamma distribution with a log link function. A constant is added to all SO.num.SSD values to adhere with the gamma assumption of positive-only values (i.e., 10 cases have SO.num.SSD = 0). Note that SO.num is not included as a covariate predicting its own variability, and that age is not included due to convergence problems.
# adding constant to each value
ema$SO.num.SSDc <- ema$SO.num.SSD + 0.001
# running analysis
glmerAn(long=ema,wide=demos,resp="SO.num.SSDc",fix.eff=c("TotalSteps1000","weekday.sleep","BMI","insomnia","sex","sex:insomnia"),
modelType=c("GLMER"),mc.predictors="TotalSteps1000",gmc.predictors="BMI",
family="gamma",link="log", # Gamma family with log link function
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="BMI",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex","sex:insomnia"),transform="exp", # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of SO.num.SSDc ...
##
## Preparing data...
## - Excluding 2094 incomplete observations ( 33.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI
## - Mean-centering TotalSteps1000
##
## Model specification:
## - model M0 (null): SO.num.SSDc ~ (1|ID)
## - model M1: SO.num.SSDc ~ TotalSteps1000.mc + (1|ID)
## - model M2: SO.num.SSDc ~ TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M3: SO.num.SSDc ~ TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M4: SO.num.SSDc ~ TotalSteps1000.mc + weekday.sleep + BMI.gmc + insomnia + (1|ID)
## - model M5: SO.num.SSDc ~ TotalSteps1000.mc + weekday.sleep + BMI.gmc + insomnia + sex + (1|ID)
## - model M6: SO.num.SSDc ~ TotalSteps1000.mc + weekday.sleep + BMI.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of SO.num.SSDc on 4125 observations from 92 participants
## using the gamma family with the log link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## TotalSteps1000.mc weekday.sleep BMI.gmc insomnia
## 1.011564 1.011670 1.044345 1.615183
## sex insomnia:sex
## 2.037035 2.546089
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: insomnia
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 6| 10726.17| 10764.11| -5357.083| 10714.17| NA| NA| NA|
## |insomnia | 7| 10721.29| 10765.56| -5353.645| 10707.29| 6.877| 1| 0.026|
## |sex | 8| 10722.06| 10772.66| -5353.029| 10706.06| 1.232| 1| 0.401|
## |sex:insomnia | 9| 10724.05| 10780.98| -5353.027| 10706.05| 0.004| 1| 0.950|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.08
## - insomnia vs. Baseline model = 0.92
## - sex = 0.385
## - sex:insomnia = 0.124
## Selected model based on Aw: insomnia
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:----------|:-----|:----------|:-----|:----------|:-----|
## |2 |(Intercept) |1.94(0.21) |6.20 |1.83(0.22) |5.01 |1.83(0.24) |4.53 |
## |3 |TotalSteps1000.mc |1.00(0.01) |-0.49 |1.00(0.01) |-0.49 |1.00(0.01) |-0.49 |
## |4 |weekday.sleep [weekend] |1.36(0.07) |5.87 |1.36(0.07) |5.86 |1.36(0.07) |5.86 |
## |5 |BMI.gmc |1.03(0.02) |1.42 |1.03(0.02) |1.51 |1.03(0.02) |1.49 |
## |6 |insomnia [1] |0.67(0.10) |-2.66 |0.68(0.10) |-2.65 |0.67(0.12) |-2.15 |
## |9 |sex [M] | | |1.18(0.18) |1.06 |1.17(0.25) |0.71 |
## |10 |insomnia [1] * sex [M] | | | | |1.02(0.31) |0.05 |
## |11 |N |92 ID |92 ID |92 ID |NA |NA |NA |
## |12 |Observations |4125 |4125 |4125 |NA |NA |NA |
## |7 |SD (Intercept) |3.34(NA) | |3.30(NA) | |3.30(NA) | |
## |8 |SD (Observations) |3.85(NA) | |3.85(NA) | |3.85(NA) | |
##
##
## - Plotting effect(s) estimated by model 4
# plotting effects without observed data (zoom in)
glmerAn(long=ema,wide=demos,resp="SO.num.SSDc",fix.eff=c("TotalSteps1000","weekday.sleep","BMI","insomnia","sex","sex:insomnia"),
mc.predictors="TotalSteps1000",gmc.predictors=c("BMI"),modelType="GLMER",family="gamma",link="log",
outputs="plotEff",plot.model="insomnia",dot.size=2,messages=FALSE,show.data=FALSE)
# keeping key results
res <- rbind(res,
cbind(measure="SO.num.SSDc",glmerAn(long=ema,wide=demos,resp="SO.num.SSDc",
fix.eff=c("TotalSteps1000","weekday.sleep","BMI","insomnia","sex","sex:insomnia"),
modelType=c("GLMER"),family="gamma",link="log",transform="exp",
mc.predictors="TotalSteps1000",gmc.predictors=c("BMI"),
mComp.baseline="BMI",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) deviance residuals and (2) random effects are quite normally distributed (despite some deviations in the tails of the residuals distribution, and in the upper tail of the random effects) with (3) no substantial linear relationship between residuals and fitted values; although categorical factors do not show substantially different dispersion, the (4) SD are slightly proportional to the means (as assumed by the gamma distribution); there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
only model M2 including
insomniais associated with a significant LRT and stronger evidence in terms of Aw compared to the baseline modelsubstantial effects are estimated for
weekday.sleep(higherSO.num.SSDduring weekend) andinsomnia(lowerSO.num.SSDin theinsomniacompared to the control group. Neithersexnor its interaction withinsomniashow a substantial effect
WakeUp.num.SSD
All s.variab variables are positively skewed and bounded at zero, thus they are modeled by assuming a gamma distribution with a log link function. A constant is added to all WakeUp.num.SSD values to adhere with the gamma assumption of positive-only values (i.e., 10 cases have WakeUp.num.SSD = 0).
# adding constant to each value
ema$WakeUp.num.SSDc <- ema$WakeUp.num.SSD + 0.001
# running analysis
glmerAn(long=ema,wide=demos,resp="WakeUp.num.SSDc",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
family="gamma",link="log",nAGQ=0, # Gamma family with log link function, nAGQ=0 to reach convergence
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("age.gmc","insomnia","sex"),transform="exp", # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of WakeUp.num.SSDc ...
##
## Preparing data...
## - Excluding 2094 incomplete observations ( 33.7 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): WakeUp.num.SSDc ~ (1|ID)
## - model M1: WakeUp.num.SSDc ~ SO.num.mc + (1|ID)
## - model M2: WakeUp.num.SSDc ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: WakeUp.num.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: WakeUp.num.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: WakeUp.num.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: WakeUp.num.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: WakeUp.num.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: WakeUp.num.SSDc ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of WakeUp.num.SSDc on 4125 observations from 92 participants
## using the gamma family with the log link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.028318 1.013427 1.039837 1.045120
## age.gmc insomnia sex insomnia:sex
## 1.111592 1.745574 2.089028 2.568794
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 11888.30| 11938.90| -5936.149| 11872.30| NA| NA| NA|
## |insomnia | 9| 11889.40| 11946.32| -5935.701| 11871.40| 0.896| 1| 0.344|
## |sex | 10| 11888.03| 11951.28| -5934.016| 11868.03| 3.369| 1| 0.199|
## |sex:insomnia | 11| 11888.81| 11958.39| -5933.406| 11866.81| 1.219| 1| 0.344|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.635
## - insomnia vs. Baseline model = 0.365
## - sex = 0.42
## - sex:insomnia = 0.221
## Selected model based on Aw: sex
##
## - Printing estimated coefficients for models age.gmc , insomnia , sex
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:----------|:-----|:----------|:-----|:----------|:-----|
## |2 |(Intercept) |1.94(0.20) |6.59 |2.07(0.30) |4.99 |1.88(0.31) |3.83 |
## |3 |SO.num.mc |1.00(0.02) |-0.12 |1.00(0.02) |-0.12 |1.00(0.02) |-0.11 |
## |4 |TotalSteps1000.mc |0.99(0.01) |-1.48 |0.99(0.01) |-1.48 |0.99(0.01) |-1.47 |
## |5 |weekday.sleep [weekend] |1.51(0.09) |7.27 |1.51(0.09) |7.26 |1.51(0.09) |7.26 |
## |6 |BMI.gmc |1.06(0.03) |1.80 |1.05(0.03) |1.70 |1.06(0.03) |1.82 |
## |7 |age.gmc |1.01(0.11) |0.12 |1.00(0.11) |-0.04 |1.02(0.11) |0.14 |
## |10 |insomnia [1] | | |0.88(0.18) |-0.59 |0.90(0.18) |-0.51 |
## |11 |sex [M] | | | | |1.27(0.26) |1.17 |
## |12 |N |92 ID |92 ID |92 ID |NA |NA |NA |
## |13 |Observations |4125 |4125 |4125 |NA |NA |NA |
## |8 |SD (Intercept) |2.47(NA) | |2.46(NA) | |2.41(NA) | |
## |9 |SD (Observations) |3.55(NA) | |3.55(NA) | |3.54(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# plotting effects without observed data (zoom in)
glmerAn(long=ema,wide=demos,resp="WakeUp.num.SSDc",fix.eff=predictors,
mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),modelType="GLMER",family="gamma",link="log",nAGQ=0,
outputs="plotEff",plot.model="insomnia",dot.size=2,messages=FALSE,show.data=FALSE)
# keeping key results
res <- rbind(res,
cbind(measure="WakeUp.num.SSDc",glmerAn(long=ema,wide=demos,resp="WakeUp.num.SSDc",fix.eff=predictors,
modelType=c("GLMER"),family="gamma",link="log",transform="exp",nAGQ=0,
mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) deviance residuals and (2) random effects are quite normally distributed (despite some deviations in the lower tails of the residuals distribution) with (3) no substantial linear relationship between residuals and fitted values; although categorical factors do not show substantially different dispersion, the (4) SD are proportional to the means (as assumed by the gamma distribution); there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT or stronger evidence in terms of Aw compared to the baseline model
a substantial effect is only estimated for
weekday.sleep(higherWakeUp.num.SSDduring weekend) whereas neitherinsomnianor its interaction withsexshow a substantial effect
1.4. Sleep autonomic func.
s.auton variables are modeled by considering the subset of complete s.auton and TotalSteps1000.
s.auton <- c("HR_NREM","HR_REM")HR_NREM
HR_NREM is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="HR_NREM",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of HR_NREM ...
##
## Preparing data...
## - Excluding 2544 incomplete observations ( 40.9 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): HR_NREM ~ (1|ID)
## - model M1: HR_NREM ~ SO.num.mc + (1|ID)
## - model M2: HR_NREM ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: HR_NREM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: HR_NREM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: HR_NREM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: HR_NREM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: HR_NREM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: HR_NREM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of HR_NREM on 3675 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.028083 1.013166 1.038783 1.042008
## age.gmc insomnia sex insomnia:sex
## 1.109891 1.731575 2.075195 2.529752
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: sex
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|------:|--:|----------:|
## |Baseline | 8| 20745.84| 20795.51| -10364.92| 20729.84| NA| NA| NA|
## |insomnia | 9| 20746.24| 20802.13| -10364.12| 20728.24| 1.596| 1| 0.206|
## |sex | 10| 20737.07| 20799.17| -10358.54| 20717.07| 11.171| 1| 0.002|
## |sex:insomnia | 11| 20734.72| 20803.02| -10356.36| 20712.72| 4.356| 1| 0.055|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.55
## - insomnia vs. Baseline model = 0.45
## - sex = 0.978
## - sex:insomnia = 0.76
## Selected model based on Aw: sex:insomnia
##
## - Printing estimated coefficients for models insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:-----|:-----------|:-----|:-----------|:-----|
## |2 |(Intercept) |60.84(0.97) |62.58 |62.69(1.06) |59.13 |63.77(1.15) |55.22 |
## |3 |SO.num.mc |-0.00(0.06) |-0.05 |-0.00(0.06) |-0.05 |-0.00(0.06) |-0.05 |
## |4 |TotalSteps1000.mc |0.16(0.02) |9.88 |0.16(0.02) |9.88 |0.16(0.02) |9.88 |
## |5 |weekday.sleep [weekend] |0.34(0.15) |2.37 |0.35(0.15) |2.38 |0.35(0.15) |2.38 |
## |6 |BMI.gmc |0.53(0.21) |2.52 |0.48(0.20) |2.45 |0.45(0.19) |2.30 |
## |7 |age.gmc |1.37(0.74) |1.84 |0.96(0.71) |1.35 |0.87(0.69) |1.25 |
## |8 |insomnia [1] |-1.77(1.39) |-1.27 |-2.18(1.32) |-1.66 |-4.21(1.60) |-2.62 |
## |11 |sex [M] | | |-4.56(1.32) |-3.45 |-7.30(1.83) |-3.98 |
## |12 |insomnia [1] * sex [M] | | | | |5.40(2.56) |2.11 |
## |13 |N |92 ID |92 ID |92 ID |NA |NA |NA |
## |14 |Observations |3675 |3675 |3675 |NA |NA |NA |
## |9 |SD (Intercept) |6.31(NA) | |5.93(NA) | |5.80(NA) | |
## |10 |SD (Observations) |1.96(NA) | |1.96(NA) | |1.96(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# keeping key results
res <- rbind(res,
cbind(measure="HR_NREM",glmerAn(long=ema,wide=demos,resp="HR_NREM",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="sex",messages=FALSE)))Comments:
both (1) residuals and (2) random effects are quite normally distributed, despite a positive deviation in the upper tail of the residuals distribution, with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
only M3 including
sexshows a significant LRT, whereas stronger evidence in terms of Aw is shown by both model M3 and model M4 including the interaction betweensexandinsomnia, compared to the baseline modelsubstantial effects are estimated for
TotalSleep1000(higherstageHR_NREMfor higher than usualTotalSleep1000),weekday.sleep(higherstageHR_NREMduring weekends),BMI(higherstageHR_NREMfor participants with higherBMI),sex(higherstageHR_NREMin girls compared to boys), and its interaction withinsomnia(lowerstageHR_NREMin insomnia girls compared to control girls)
HR_REM
HR_REM is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="HR_REM",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="age",p.adjust.method="BH",
coeff.models=c("insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of HR_REM ...
##
## Preparing data...
## - Excluding 2519 incomplete observations ( 40.5 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Grand-mean-centering BMI , age
## - Mean-centering SO.num , TotalSteps1000
##
## Model specification:
## - model M0 (null): HR_REM ~ (1|ID)
## - model M1: HR_REM ~ SO.num.mc + (1|ID)
## - model M2: HR_REM ~ SO.num.mc + TotalSteps1000.mc + (1|ID)
## - model M3: HR_REM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M4: HR_REM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + (1|ID)
## - model M5: HR_REM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + (1|ID)
## - model M6: HR_REM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + (1|ID)
## - model M7: HR_REM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + (1|ID)
## - model M8: HR_REM ~ SO.num.mc + TotalSteps1000.mc + weekday.sleep + BMI.gmc + age.gmc + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of HR_REM on 3700 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## SO.num.mc TotalSteps1000.mc weekday.sleep BMI.gmc
## 1.026506 1.012604 1.037394 1.042012
## age.gmc insomnia sex insomnia:sex
## 1.109866 1.731580 2.074946 2.529550
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: sex
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 21235.74| 21285.47| -10609.87| 21219.74| NA| NA| NA|
## |insomnia | 9| 21235.68| 21291.62| -10608.84| 21217.68| 2.064| 1| 0.151|
## |sex | 10| 21228.17| 21290.33| -10604.08| 21208.17| 9.511| 1| 0.006|
## |sex:insomnia | 11| 21226.16| 21294.53| -10602.08| 21204.16| 4.011| 1| 0.068|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.492
## - insomnia vs. Baseline model = 0.508
## - sex = 0.956
## - sex:insomnia = 0.723
## Selected model based on Aw: sex:insomnia
##
## - Printing estimated coefficients for models insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:-----|:-----------|:-----|:-----------|:-----|
## |2 |(Intercept) |63.04(0.98) |64.09 |64.77(1.08) |59.86 |65.83(1.18) |55.74 |
## |3 |SO.num.mc |0.13(0.06) |2.16 |0.13(0.06) |2.16 |0.13(0.06) |2.16 |
## |4 |TotalSteps1000.mc |0.14(0.02) |7.97 |0.14(0.02) |7.97 |0.14(0.02) |7.97 |
## |5 |weekday.sleep [weekend] |0.29(0.15) |1.91 |0.29(0.15) |1.92 |0.29(0.15) |1.92 |
## |6 |BMI.gmc |0.41(0.21) |1.92 |0.37(0.20) |1.81 |0.33(0.20) |1.66 |
## |7 |age.gmc |0.84(0.75) |1.12 |0.46(0.72) |0.64 |0.37(0.71) |0.52 |
## |8 |insomnia [1] |-2.04(1.41) |-1.45 |-2.43(1.34) |-1.81 |-4.41(1.64) |-2.69 |
## |11 |sex [M] | | |-4.28(1.35) |-3.17 |-6.96(1.87) |-3.72 |
## |12 |insomnia [1] * sex [M] | | | | |5.30(2.62) |2.02 |
## |13 |N |92 ID |92 ID |92 ID |NA |NA |NA |
## |14 |Observations |3700 |3700 |3700 |NA |NA |NA |
## |9 |SD (Intercept) |6.38(NA) | |6.05(NA) | |5.93(NA) | |
## |10 |SD (Observations) |2.01(NA) | |2.01(NA) | |2.01(NA) | |
##
##
## - Plotting effect(s) estimated by model 6
# keeping key results
res <- rbind(res,
cbind(measure="HR_REM",glmerAn(long=ema,wide=demos,resp="HR_REM",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="sex",messages=FALSE)))Comments:
both (1) residuals and (2) random effects are quite normally distributed, despite a positive deviation in the upper tail of the residuals distribution, with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
only M3 including
sexshows a significant LRT, whereas stronger evidence in terms of Aw is shown by both model M3 and model M4 including the interaction betweensexandinsomnia, compared to the baseline modelsubstantial effects are estimated for
SO.num(higherstageHR_REMfor laterSOthan usual),sex(higherstageHR_REMin girls compared to boys), and its interaction withinsomnia(lowerstageHR_REMin insomnia girls compared to control girls)
1.5. Summary of results
Here, we use the key.resPlot function to graphically summarize the results obtained from the models above. The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering either model M2 (main effect of insomnia) or model M3 (main effect of insomnia and sex) when sex differences were substantial. The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
key.resPlot(res)
Comments:
the analysis of sleep characterization highlighted substantial differences between the insomnia and the control group for
TIB.SSD,TST.SSD(only supported by t-values and Aw, but not by the LRT), andSO.num.SSD(with lower night-to-night variability for the insomnia than the control group)with
HR_NREMandHR_REM, the interactive model was selected and predicted lower HR in female insomnia than in female control, but no differences between boys, with no main effect ofinsomnianone of the remaining sleep outcomes showed significant LRT, higher Aw, and substantial relationships for the insomnia predictor
1.6. Robustness checks
Here, we conduct several robustness checks to account for the arbitrariness of our data processing and analytical choices, including (1) the consideration of insomnia.group instead of the insomnia variable, (2) the inclusion of the covid19 variable as a further covariate, (3) the exclusion of participants with extreme missing data, (4) the exclusion of TotalSteps1000, (5) the use of alternative family distributions for the analyses.
Here, we select the predictors, and the sleep outcomes.
# predictors
PREdictors <- predictors[1:(length(predictors)-2)] # excluding sex and interactive model
# sleep outcomes
sleepVars <- c(s.archit[2:length(s.archit)],s.timing,paste(s.variab,"c",sep=""),s.auton)1.6.1. insomnia.group
Here, we inspect the results obtained by considering the insomnia.group (i.e., 46 controls, 26 DSM.ins and 21 sub.ins) rather than the insomnia variable (i.e., 46 controls vs. 47 insomnia).
The figure below shows the t-values associated with the group differences between full DSM insomnia and controls for each sleep outcome, by considering either model M2 (main effect of insomnia) or model M3 (main effect of insomnia and sex) when sex differences were substantial (i.e., HR_NREM and HR_REM). The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# computing key results for TIB
res2 <- cbind(measure="TIB",glmerAn(long=ema,wide=demos,resp="TIB",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia.group",key.model="insomnia.group",messages=FALSE))
# computing key results for the following sleep outcomes
for(sleepVar in sleepVars){ if(grepl("SSD",sleepVar)){ fam <- c("gamma","log") }else{ fam <- c("normal","identity") } # family
if(grepl("SO.num",sleepVar)){ Predictors <- gsub("insomnia","insomnia.group",PREdictors)[2:length(PREdictors)] # insomnia.group
} else { Predictors <- gsub("insomnia","insomnia.group",PREdictors) }
if(grepl("HR",sleepVar)){ key.model="sex" # adding sex with s.auton variables
Predictors <- c(Predictors,"sex") }else{ key.model <- "insomnia.group" }
if(sleepVar%in%c("WASO.SSDc","WakeUp.num.SSDc")){ nAGQ <- 0 }else{ nAGQ <- 1 } # nAGQ = 0 to reach convergence
res2 <- rbind(res2,
cbind(measure=sleepVar,glmerAn(long=ema,wide=demos,resp=sleepVar,fix.eff=Predictors,nAGQ=nAGQ,
family=fam[1],link=fam[2],key.model=key.model,key.predictor="insomnia.group",
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res$robCheck <- "Original"
res2$robCheck <- "insomnia.group"
RES <- rbind(res,res2)
# plotting robustness check
key.resPlot(RES,robCheck = "robCheck")
Comments:
results are highly consistent with the main analyses for most sleep outcomes
the only differences are highlighted for
TST.SSD, showing a no longer significant LRT (but still higher Aw and substantial t), andHR_NREMandHR_REM, showing a no longer higher Aw (but still substantial t)
Here, we can see that the differences in TIB.SSD, TST.SSD, and SO.num.SSD are substantial, and driven by the full DSM insomnia subgroup, whereas no substantial differences are estimated between the control and the sub-threshold insomnia group.
plot_grid(list(glmerAn(long=ema,wide=demos,resp="TIB.SSDc",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
family="gamma",link="log",modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs=c("mComp","coeff","plotEff"),coeff.models="insomnia.group",
plot.model="insomnia.group",plot.pred="insomnia.group",message=FALSE,return.plot=TRUE),
glmerAn(long=ema,wide=demos,resp="TST.SSDc",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
family="gamma",link="log",modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs=c("mComp","coeff","plotEff"),coeff.models="insomnia.group",
plot.model="insomnia.group",plot.pred="insomnia.group",message=FALSE,return.plot=TRUE),
glmerAn(long=ema,wide=demos,resp="SO.num.SSDc",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
family="gamma",link="log",modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),
gmc.predictors=c("BMI","age"),mComp.baseline="age",p.adjust.method="BH",
outputs=c("mComp","coeff","plotEff"),coeff.models="insomnia.group",
plot.model="insomnia.group",plot.pred="insomnia.group",message=FALSE,return.plot=TRUE)),
tags=TRUE,margin=rep(0,4))##
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:--------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 83036.86| 83087.46| -41510.43| 83020.86| NA| NA| NA|
## |insomnia.group | 10| 83034.82| 83098.07| -41507.41| 83014.82| 6.037| 2| 0.049|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia.group : = 0.265
## - insomnia.group vs. Baseline model = 0.735
## Selected model based on Aw: insomnia.group
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:-----------|:------|
## |2 |(Intercept) |9.35(0.08) |110.29 |
## |3 |SO.num.mc |-0.03(0.02) |-1.92 |
## |4 |TotalSteps1000.mc |-0.00(0.01) |-0.53 |
## |5 |weekday.sleep [weekend] |0.38(0.05) |7.38 |
## |6 |BMI.gmc |0.05(0.02) |2.55 |
## |7 |age.gmc |0.01(0.06) |0.11 |
## |8 |insomnia.group [DSM.ins] |-0.36(0.14) |-2.50 |
## |9 |insomnia.group [sub.ins] |-0.15(0.15) |-1.04 |
## |12 |N ID |92 |NA |
## |13 |Observations |4125 |NA |
## |10 |SD (Intercept) |0.76(NA) | |
## |11 |SD (Observations) |1.23(NA) | |
##
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:--------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 81388.64| 81439.24| -40686.32| 81372.64| NA| NA| NA|
## |insomnia.group | 10| 81387.84| 81451.09| -40683.92| 81367.84| 4.795| 2| 0.091|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia.group : = 0.402
## - insomnia.group vs. Baseline model = 0.598
## Selected model based on Aw: insomnia.group
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:-----------|:------|
## |2 |(Intercept) |9.12(0.08) |108.28 |
## |3 |SO.num.mc |-0.03(0.02) |-1.67 |
## |4 |TotalSteps1000.mc |-0.00(0.01) |-0.39 |
## |5 |weekday.sleep [weekend] |0.38(0.05) |7.47 |
## |6 |BMI.gmc |0.04(0.02) |2.46 |
## |7 |age.gmc |-0.01(0.06) |-0.15 |
## |8 |insomnia.group [DSM.ins] |-0.31(0.14) |-2.18 |
## |9 |insomnia.group [sub.ins] |-0.18(0.15) |-1.23 |
## |12 |N ID |92 |NA |
## |13 |Observations |4125 |NA |
## |10 |SD (Intercept) |0.73(NA) | |
## |11 |SD (Observations) |1.22(NA) | |
##
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:--------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 8| 10726.05| 10776.65| -5355.024| 10710.05| NA| NA| NA|
## |insomnia.group | 10| 10722.66| 10785.91| -5351.329| 10702.66| 7.39| 2| 0.025|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia.group : = 0.155
## - insomnia.group vs. Baseline model = 0.845
## Selected model based on Aw: insomnia.group
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:-----------|:-----|
## |2 |(Intercept) |0.64(0.11) |6.07 |
## |3 |SO.num.mc |-0.01(0.01) |-0.89 |
## |4 |TotalSteps1000.mc |-0.00(0.01) |-0.47 |
## |5 |weekday.sleep [weekend] |0.32(0.05) |5.95 |
## |6 |BMI.gmc |0.03(0.02) |1.18 |
## |7 |age.gmc |0.09(0.08) |1.17 |
## |8 |insomnia.group [DSM.ins] |-0.49(0.18) |-2.77 |
## |9 |insomnia.group [sub.ins] |-0.18(0.19) |-0.96 |
## |12 |N ID |92 |NA |
## |13 |Observations |4125 |NA |
## |10 |SD (Intercept) |1.16(NA) | |
## |11 |SD (Observations) |1.35(NA) | |
1.6.2. COVID-19 emergency
Here, we inspect the results obtained by including the covid19 variable as a further covariate. As shown in section 1.2.8, the covid19 factor discriminates the data collected pre-COVID19 (63%) and the data collected post-COVID19 (37%). The covariate is included at the first step (model M1, before SO.num).
The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering either model M2 (main effect of insomnia) or model M3 (main effect of insomnia and sex) when sex differences were substantial (i.e., HR_NREM and HR_REM. The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# computing key results for TIB
res2 <- cbind(measure="TIB",glmerAn(long=ema,wide=demos,resp="TIB",fix.eff=c("covid19",PREdictors),
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE))
# computing key results for the following sleep outcomes
for(sleepVar in sleepVars){ if(grepl("SSD",sleepVar)){ fam <- c("gamma","log") }else{ fam <- c("normal","identity") } # family
if(grepl("SO.num",sleepVar)){ Predictors <- c("covid19",PREdictors[2:length(PREdictors)]) # predictors + covid19
} else { Predictors <- c("covid19",PREdictors) }
if(grepl("HR",sleepVar)){ key.model="sex" # adding sex with s.auton variables
Predictors <- c("covid19",PREdictors,"sex") } else { key.model <- "insomnia" }
if(sleepVar%in%c("WASO.SSDc","WakeUp.num.SSDc")){ nAGQ <- 0 }else{ nAGQ <- 1 } # nAGQ = 0 to converge
res2 <- rbind(res2,
cbind(measure=sleepVar,glmerAn(long=ema,wide=demos,resp=sleepVar,fix.eff=Predictors,nAGQ=nAGQ,
family=fam[1],link=fam[2],key.model=key.model,key.predictor="insomnia",
mComp.baseline="age",p.adjust.method="BH",
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "covid19"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","covid19"),],robCheck = "robCheck")
Comments:
results are partially inconsistent with the main analyses
the main differences are highlighted for
s.variabands.autonoutcomes, with a no longer substantial main effect onSO.num.SSD,HR_NREM, andHR_REM
Here, we can see that the post-COVID19 data collection period is characterized by longer TIB and TST (but not WASO), slightly lower rem.p, later SO and WakeUp times, lower SO.num.SSD, and lower HR_REM compared to the pre-COVID19 period.
# plotting difference by covid (model M1)
for(sleepVar in c("TIB",sleepVars)){
if(grepl("SSD",sleepVar)){ fam <- c("gamma","log") }else{ fam <- c("normal","identity") }
if(grepl("SO.num",sleepVar)){ Predictors <- c("covid19","TotalSteps1000","weekday.sleep","BMI")
} else { Predictors <- c("covid19","SO.num","TotalSteps1000","weekday.sleep","BMI") }
if(grepl("HR",sleepVar)){ key.model="sex" }else{ key.model <- "insomnia" }
if(sleepVar%in%c("WASO.SSDc","WakeUp.num.SSDc","SO.num.SSDc")){ nAGQ <- 0 }else{ nAGQ <- 1 } # nAGQ = 0 to converge
glmerAn(long=ema,wide=demos,resp=sleepVar,fix.eff=Predictors,
mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),modelType="GLMER",family=fam[1],link=fam[2],nAGQ=nAGQ,
outputs="plotEff",plot.model="insomnia",plot.pred="covid19",messages=FALSE,dot.size=3,dodge=0.5) }
In conclusion, the results of this robustness check were partially inconsistent with those obtained with the main analyses, questioning the generalizability of the results observed for SO.num.SSD and for the main effect of insomnia on both s.auton outcomes.
1.6.3. Low compliance
Here, we inspect the results obtained by excluding seven participants (7.5%) with less than two weeks of valid observations in any data type. In section 3.2, these participants were marked with majMiss = 1.
The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering either model M2 (main effect of insomnia) or model M3 (main effect of insomnia and sex) when sex differences were substantial (i.e., HR_NREM and HR_REM. The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# excluding participants with extreme missing data
EMA <- ema[ema$majMiss!=1,]
DEMOS <- demos[demos$ID%in%levels(as.factor(as.character(EMA$ID))),]
cat("Excluded",nrow(ema)-nrow(EMA),"days from",nrow(demos)-nrow(DEMOS),"participants")## Excluded 465 days from 7 participants# computing key results for TIB
res2 <- cbind(measure="TIB",glmerAn(long=EMA,wide=DEMOS,resp="TIB",fix.eff=PREdictors,
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE))
# computing key results for the other sleep outcomes
for(sleepVar in sleepVars){
if(grepl("SSD",sleepVar)){ fam <- c("gamma","log") }else{ fam <- c("normal","identity") }
if(grepl("SO.num",sleepVar)){ Predictors <- PREdictors[2:length(PREdictors)] } else { Predictors <- PREdictors }
if(grepl("HR",sleepVar)){ Predictors <- c(Predictors,"sex")
key.model="sex" }else{ key.model <- "insomnia" }
if(sleepVar%in%c("WASO.SSDc","WakeUp.num.SSDc","SO.num.SSDc")){ nAGQ <- 0 }else{ nAGQ <- 1 } # nAGQ = 0 to converge
res2 <- rbind(res2,
cbind(measure=sleepVar,glmerAn(long=EMA,wide=DEMOS,resp=sleepVar,fix.eff=Predictors,nAGQ=nAGQ,
family=fam[1],link=fam[2],key.model=key.model,key.predictor="insomnia",
mComp.baseline="age",p.adjust.method="BH",
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "majMiss"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","majMiss"),],robCheck = "robCheck")
Comments:
results are overall inconsistent with the main analyses
the main differences is highlighted for
SO.num.SSD, showing no more substantial group differences. Note, however, thatnAGQwas set to zero to reach convergence inSO.num.SSDmodels, resulting in higher standard errors that might have masked the effect
1.6.4. TotalStep1000
Here, we inspect the results obtained by excluding TotalSteps1000 from the models predictors. This is done in order to replicate the analyses by including cases with invalid daily steps data (thus, with higher sample size).
The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering either model M2 (main effect of insomnia) or model M3 (main effect of insomnia and sex) when sex differences were substantial (i.e., HR_NREM and HR_REM. The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# computing key results for TIB
res2 <- cbind(measure="TIB",glmerAn(long=ema,wide=demos,resp="TIB",fix.eff=PREdictors[PREdictors!="TotalSteps1000"],
modelType=c("GLMER"),mc.predictors="SO.num",gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE))
# computing key results for the following sleep outcomes
for(sleepVar in sleepVars){
if(grepl("SSD",sleepVar)){ fam <- c("gamma","log") }else{ fam <- c("normal","identity") }
if(grepl("SO.num",sleepVar)){ Predictors <- PREdictors[PREdictors!="TotalSteps1000" & PREdictors!="SO.num"]
} else { Predictors <- PREdictors[PREdictors!="TotalSteps1000"] }
if(grepl("HR",sleepVar)){ key.model="sex"
Predictors <- c(Predictors,"sex") }else{ key.model <- "insomnia" }
if(sleepVar%in%c("WASO.SSDc","WakeUp.num.SSDc","TIB.SSDc")){ nAGQ <- 0 }else{ nAGQ <- 1 } # nAGQ = 0 to converge
res2 <- rbind(res2,
cbind(measure=sleepVar,glmerAn(long=ema,wide=demos,resp=sleepVar,fix.eff=Predictors,nAGQ=nAGQ,
family=fam[1],link=fam[2],key.model=key.model,key.predictor="insomnia",
mComp.baseline="age",p.adjust.method="BH",
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "no-TotalSteps"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","no-TotalSteps"),],robCheck = "robCheck")
Comments:
results are partially inconsistent with the main analyses
the main differences are highlighted for
s.variabands.autonoutcomes, with a no longer substantial main effect onTIB.SSD,HR_NREM, andHR_REMnote, however, that
nAGQwas set to zero to reach convergence in TIB.SSD models, resulting in higher standard errors that might have masked the effect
Here, we can see that the differences in HR_NREM and HR_REM are no longer substantial, whereas differences in TIB.SSD are smaller but still quite substantial, as it is substantial the interaction between sex and insomnia on s.auton variables. 
In conclusion, the results of this robustness check were overall consistent with those obtained with the main analyses, although questioning the generalizability of the results observed for TIB.SSD and (slightly less) s.auton outcomes.
1.6.5. Family distribution
Here, we replicate the models predicting those variables that showed poor fit in the model diagnostics above by re-specifying the models with different distribution families and link functions:
WASO: initially modeled by assuming a normal distribution for residuals, despite the highlighted deviation in the upper tail of the distribution. SinceWASOis bounded at zero and positively skewed, we inspect how results change when assuming a gamma distribution. Here, we can see that the fit does not clearly improve when the family distribution is changed.
# WASO: normal vs. gamma (not optimal fit in both cases)
ema$WASOc <- ema$WASO + 0.001
p1 <- lmer(WASOc~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema)
p2 <- glmer(WASOc~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema,family=Gamma("log"),nAGQ=0) # nAGQ = 0 to converge
p3 <- glmer(WASOc~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema,family=Gamma("inverse"),nAGQ=0)
p4 <- glmer(WASOc~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema,family=gaussian("log"),nAGQ=0)
par(mfrow=c(2,2))
qqnorm(resid(p1),main="qqnorm residuals - Normal"); qqline(resid(p1),col="red")
qqnorm(resid(p2),main="qqnorm residuals - Gamma(log)"); qqline(resid(p2),col="red")
qqnorm(resid(p3),main="qqnorm residuals - Gamm(inverse"); qqline(resid(p3),col="red")
qqnorm(resid(p4),main="qqnorm residuals - Normal(log)"); qqline(resid(p4),col="red")
s.variab: modeled by assuming a gamma distribution. Here, we re-specify the models by assuming normality (although it can be seen that the fit is strongly impaired).
# TST.SSDc: normal vs. gamma (clearly better fit for gamma)
p1 <- lmer(TST.SSDc~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema)
p2 <- glmer(TST.SSDc~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema,family=Gamma("log"))
p2 <- glmer(TST.SSDc~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema,family=Gamma("log"))
par(mfrow=c(1,2))
qqnorm(resid(p1),main="qqnorm residuals - Normal"); qqline(resid(p1),col="red")
qqnorm(resid(p2),main="qqnorm residuals - Gamma(log)"); qqline(resid(p2),col="red") 
s.auton: modeled by assuming a normal distribution, but showing marked deviations in the upper tail. Again, these are re-modeled by assuming a gamma distribution. Here, we can see that the fit does not clearly improve when the family distribution is changed.
# HR_NREM: normal vs. gamma (unoptimal fit in both cases)
p1 <- lmer(HR_NREM~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema)
p2 <- glmer(HR_NREM~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema,family=Gamma("log"),nAGQ=0)
p3 <- glmer(HR_NREM~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema,family=Gamma("identity"),nAGQ=0)
p4 <- glmer(HR_NREM~SO.num+TotalSteps1000+weekday.sleep+BMI+insomnia+(1|ID),data=ema,family=Gamma("inverse"),nAGQ=0)
par(mfrow=c(2,2))
qqnorm(resid(p1),main="qqnorm residuals - Normal"); qqline(resid(p1),col="red")
qqnorm(resid(p2,type="deviance"),main="qqnorm residuals - Gamma(log)"); qqline(resid(p2),col="red")
qqnorm(resid(p3,type="deviance"),main="qqnorm residuals - Gamma(identity)"); qqline(resid(p3),col="red")
qqnorm(resid(p4,type="deviance"),main="qqnorm residuals - Gamma(inverse)"); qqline(resid(p4),col="red") 
The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering either model M2 (main effect of insomnia) or model M3 (main effect of insomnia and sex) when sex differences were substantial (i.e., HR_NREM and HR_REM. The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# computing key results for WASO
ema$WASOc <- ema$WASO + 0.001
res2 <- cbind(measure="WASO",glmerAn(long=ema,wide=demos,resp="WASOc",fix.eff=PREdictors,nAGQ=0, # nAGQ = 0 to reach convergence
modelType=c("GLMER"),family="gamma",link="log",mComp.baseline="age",p.adjust.method="BH",
mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE))
# computing key results for the following sleep outcomes
for(sleepVar in c(paste(s.variab,"c",sep=""),s.auton)){
if(grepl("SSD",sleepVar)){ fam <- c("normal","identity") }else{ fam <- c("gamma","log") }
if(grepl("HR",sleepVar)){ key.model="sex"
Predictors <- c(PREdictors,"sex") }else{ key.model <- "insomnia"
Predictors <- PREdictors
nAGQ <- 0 } # nAGQ = 0 to converge
res2 <- rbind(res2,
cbind(measure=sleepVar,glmerAn(long=ema,wide=demos,resp=sleepVar,fix.eff=Predictors,nAGQ=nAGQ,
family=fam[1],link=fam[2],key.model=key.model,key.predictor="insomnia",
modelType=c("GLMER"),mc.predictors=c("SO.num","TotalSteps1000"),gmc.predictors=c("BMI","age"),
mComp.baseline="age",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "familyDistr"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","familyDistr"),],robCheck = "robCheck")
Comments:
results are partially inconsistent with the main analyses
the main differences are highlighted for
s.variabands.autonoutcomes, with a no longer substantial main effect onTST.SSD, and no longer significant LRT (bust still higher Aw and substantial t) forHR_NREMandHR_REM
1.6.6. Summary of results
Here, we summarize the results from all robustness checks. The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering either model M2 (main effect of insomnia) or model M3 (main effect of insomnia and sex) when sex differences were substantial (i.e., HR_NREM and HR_REM. The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
key.resPlot(RES,robCheck = "robCheck")
Comments:
the null results observed from the main analyses (
sarchit,s.timing, ands.auton) are highly consistent across the robustness checksthe lower
s.variabinTIB,TST, andSOhighlighted for theinsomniagroup shows some inconsistencies, especially forSO.num.SSD(whose effect is no longer substantial following the exclusion ofmajMissparticipants, or the inclusion ofcovid19)more consistent results are found for
TST.SSD(only affected by the distribution family, but with a better fit for the main analyses than the robustness checks) andTIB.SSD(whose group differences are no longer substantial whenTotalStepsare excluded from the model, likely due to the lack of convergence with nAGQ = 1)
In conclusion, the results obtained with the main analyses are quite robust, with SO.SSD showing the less consistent relationship with insomnia.
2. Diary ratings
Here, bedtime ratings of stress, worry, and mood are predicted by (1) meaningful covariates (TotalSteps1000 and weekday), (2) the insomnia group, (3) participants’ sex, and (4) the interaction between sex and insomnia groups, in order to characterize the two groups in terms of pre-sleep psychological states.
predictors <- c("TotalSteps1000","weekday.sleep","insomnia","sex","sex:insomnia")In addition to considering the raw item scores, we replicate the models by using the aggregated values (PsyDist) and the factor scores (fs) predicted by the MCFA model specified above.
2.1. Raw item scores
Daily diary item scores are modeled by considering the subset of complete dailyDiary and TotalSteps1000. Due to the ordinal nature of these single-item measures, we do not use normality-based GLMER. Instead, we binary transform each item score (i.e., 0 = score lower than or equal to 3, 1 = score higher than 3), and we use logistic regression.
# recoding stress, worry, and mood in two categories
ema[!is.na(ema$stress),c("stress.cat","worry.cat","mood.cat")] <- 0
ema[!is.na(ema$stress) & ema$stress>3,"stress.cat"] <- 1
ema[!is.na(ema$worry) & ema$worry>3,"worry.cat"] <- 1
ema[!is.na(ema$mood) & ema$mood>3,"mood.cat"] <- 1
# summarizing the categorical recoded variables
summary(as.factor(ema$stress.cat))## 0 1 NA's
## 4325 607 1287summary(as.factor(ema$worry.cat))## 0 1 NA's
## 4311 621 1287summary(as.factor(ema$mood.cat))## 0 1 NA's
## 4219 713 1287Comments:
- in the 12-to-14% of the cases, stress, worry, and mood were rated as 3 or higher
Here, we use the glmerAn function to conduct logistic regression for each item score.
STRESS
stress is a one-item ordinal variable, and it was dicotomized to be modeled by assuming a binomial distribution for the residuals (i.e., logistic regression).
glmerAn(long=ema,wide=demos,resp="stress.cat",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
family="binomial",link="logit",transform="exp", # logistic regression (exponentiated coefficients = odds ratios)
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="weekday.sleep",p.adjust.method="BH",
coeff.models=c("weekday.sleep","insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia")## Running GLMER analysis of stress.cat ...
##
## Preparing data...
## - Excluding 2270 incomplete observations ( 36.5 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Mean-centering TotalSteps1000
##
## Model specification:
## - model M0 (null): stress.cat ~ (1|ID)
## - model M1: stress.cat ~ TotalSteps1000.mc + (1|ID)
## - model M2: stress.cat ~ TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M3: stress.cat ~ TotalSteps1000.mc + weekday.sleep + insomnia + (1|ID)
## - model M4: stress.cat ~ TotalSteps1000.mc + weekday.sleep + insomnia + sex + (1|ID)
## - model M5: stress.cat ~ TotalSteps1000.mc + weekday.sleep + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of stress.cat on 3949 observations from 92 participants
## using the binomial family with the logit link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:##
## Printing Variance Inflation Factors (VIF):
## TotalSteps1000.mc weekday.sleep insomnia sex
## 1.008725 1.008607 1.521476 2.041209
## insomnia:sex
## 2.489489
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 4| 2547.329| 2572.454| -1269.665| 2539.329| NA| NA| NA|
## |insomnia | 5| 2545.831| 2577.237| -1267.915| 2535.831| 3.498| 1| 0.141|
## |sex | 6| 2545.032| 2582.719| -1266.516| 2533.032| 2.799| 1| 0.141|
## |sex:insomnia | 7| 2546.591| 2590.559| -1266.295| 2532.591| 0.441| 1| 0.507|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.321
## - insomnia vs. Baseline model = 0.679
## - sex = 0.503
## - sex:insomnia = 0.187
## Selected model based on Aw: sex
##
## - Printing estimated coefficients for models weekday.sleep , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:----------|:------|:----------|:------|:----------|:------|:----------|:------|
## |2 |(Intercept) |0.10(0.01) |-16.77 |0.08(0.01) |-13.40 |0.09(0.02) |-11.34 |0.09(0.02) |-10.51 |
## |3 |TotalSteps1000.mc |0.95(0.01) |-3.25 |0.95(0.01) |-3.25 |0.95(0.01) |-3.26 |0.95(0.01) |-3.26 |
## |4 |weekday.sleep [weekend] |0.67(0.09) |-3.18 |0.67(0.09) |-3.18 |0.67(0.09) |-3.18 |0.67(0.09) |-3.18 |
## |7 |insomnia [1] | | |1.63(0.42) |1.91 |1.60(0.40) |1.86 |1.80(0.56) |1.89 |
## |8 |sex [M] | | | | |0.64(0.17) |-1.69 |0.76(0.29) |-0.70 |
## |9 |insomnia [1] * sex [M] | | | | | | |0.70(0.37) |-0.67 |
## |10 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |11 |Observations |3949 |3949 |3949 |3949 |NA |NA |NA |NA |
## |5 |SD (Intercept) |2.98(NA) | |2.88(NA) | |2.82(NA) | |2.82(NA) | |
## |6 |SD (Observations) |2.72(NA) | |2.72(NA) | |2.72(NA) | |2.72(NA) | |
##
##
## - Plotting effect(s) estimated by model 3## Data were 'prettified'. Consider using `terms="TotalSteps1000.mc [all]"` to get smooth plots.
# keeping key results
res <- cbind(measure="stress.cat",glmerAn(long=ema,wide=demos,resp="stress.cat",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors="TotalSteps1000",family="binomial",link="logit",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",dot.size=3,messages=FALSE))Comments:
the logistic regression assumption of (1) linearity between logit and continuous predictors generally holds, although the slightly higher logit for weekend days, insomnia, and females; (2) random effects are quite normally distributed, with only a slight deviation from normality in the lower tail of the random intercept distribution; there is (3) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models show a significant LRT, while model M2 and M3 including
insomniaandsexshow stronger evidence in terms of Aw compared to the baseline modelsubstantial effects are only estimated for
TotalSteps1000(lowerstressfor higher than usualTotalSteps), andweekday.sleep(lowerstresson weekend days than on weekdays). Neitherinsomnianorsexor their interaction show a substantial effect (i.e., despite the slightly higherstressin insomnia than in controls, and especially in insomnia girls compared to control girls, these differences are not substantial)
WORRY
worry is a one-item ordinal variable, and it was dicotomized to be modeled by assuming a binomial distribution for the residuals (i.e., logistic regression).
glmerAn(long=ema,wide=demos,resp="worry.cat",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
family="binomial",link="logit",transform="exp", # logistic regression (exponentiated coefficients = odds ratios)
nAGQ=0, # nAGQ = 0 to reach convergence in the interactive model
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="weekday.sleep",p.adjust.method="BH",
coeff.models=c("weekday.sleep","insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia")## Running GLMER analysis of worry.cat ...
##
## Preparing data...
## - Excluding 2270 incomplete observations ( 36.5 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Mean-centering TotalSteps1000
##
## Model specification:
## - model M0 (null): worry.cat ~ (1|ID)
## - model M1: worry.cat ~ TotalSteps1000.mc + (1|ID)
## - model M2: worry.cat ~ TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M3: worry.cat ~ TotalSteps1000.mc + weekday.sleep + insomnia + (1|ID)
## - model M4: worry.cat ~ TotalSteps1000.mc + weekday.sleep + insomnia + sex + (1|ID)
## - model M5: worry.cat ~ TotalSteps1000.mc + weekday.sleep + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of worry.cat on 3949 observations from 92 participants
## using the binomial family with the logit link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:##
## Printing Variance Inflation Factors (VIF):
## TotalSteps1000.mc weekday.sleep insomnia sex
## 1.012794 1.012824 1.522278 2.002355
## insomnia:sex
## 2.435146
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 4| 2469.388| 2494.513| -1230.694| 2461.388| NA| NA| NA|
## |insomnia | 5| 2470.157| 2501.563| -1230.079| 2460.157| 1.231| 1| 0.354|
## |sex | 6| 2469.913| 2507.600| -1228.957| 2457.913| 2.244| 1| 0.354|
## |sex:insomnia | 7| 2471.056| 2515.024| -1228.528| 2457.056| 0.857| 1| 0.354|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.595
## - insomnia vs. Baseline model = 0.405
## - sex = 0.314
## - sex:insomnia = 0.151
## Selected model based on Aw: Baseline
##
## - Printing estimated coefficients for models weekday.sleep , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:----------|:------|:----------|:------|:----------|:-----|:----------|:-----|
## |2 |(Intercept) |0.10(0.02) |-14.37 |0.09(0.02) |-10.84 |0.10(0.03) |-9.00 |0.09(0.03) |-8.50 |
## |3 |TotalSteps1000.mc |0.99(0.01) |-0.89 |0.99(0.01) |-0.89 |0.99(0.01) |-0.89 |0.99(0.01) |-0.89 |
## |4 |weekday.sleep [weekend] |0.62(0.08) |-3.57 |0.62(0.08) |-3.57 |0.62(0.08) |-3.57 |0.62(0.08) |-3.57 |
## |7 |insomnia [1] | | |1.38(0.43) |1.03 |1.35(0.42) |0.97 |1.63(0.62) |1.29 |
## |8 |sex [M] | | | | |0.64(0.21) |-1.38 |0.84(0.39) |-0.37 |
## |9 |insomnia [1] * sex [M] | | | | | | |0.57(0.37) |-0.86 |
## |10 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |11 |Observations |3949 |3949 |3949 |3949 |NA |NA |NA |NA |
## |5 |SD (Intercept) |3.80(NA) | |3.76(NA) | |3.70(NA) | |3.69(NA) | |
## |6 |SD (Observations) |2.72(NA) | |2.72(NA) | |2.72(NA) | |2.72(NA) | |
##
##
## - Plotting effect(s) estimated by model 3## Data were 'prettified'. Consider using `terms="TotalSteps1000.mc [all]"` to get smooth plots.
# keeping key results
res <- rbind(res,cbind(measure="worry.cat",glmerAn(long=ema,wide=demos,resp="worry.cat",fix.eff=predictors[1:(length(predictors)-1)],
modelType=c("GLMER"),mc.predictors="TotalSteps1000",family="binomial",link="logit",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",dot.size=3,messages=FALSE)))Comments:
the logistic regression assumption of (1) linear independence between logit and continuous predictors generally holds, although the slightly higher logit for weekend days, insomnia, and females; random effects are quite normally distributed, with only a slight deviation from normality in the lower tail of the random intercept distribution; there is (3) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT or stronger evidence in terms of Aw compared to the baseline model
a substantial effect is only estimated for
weekday.sleep(lowerworryon weekend days than on weekdays). Neitherinsomnianorsex, or their interaction show a substantial effect (i.e., despite the slightly higherworryin insomnia than in controls, and especially in insomnia girls compared to control girls, these differences are not substantial)
MOOD
mood is a one-item ordinal variable, and it was dicotomized to be modeled by assuming a binomial distribution for the residuals (i.e., logistic regression).
glmerAn(long=ema,wide=demos,resp="mood.cat",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
family="binomial",link="logit",transform="exp", # logistic regression (exponentiated coefficients = odds ratios)
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="weekday.sleep",p.adjust.method="BH",
coeff.models=c("weekday.sleep","insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia")## Running GLMER analysis of mood.cat ...
##
## Preparing data...
## - Excluding 2270 incomplete observations ( 36.5 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Mean-centering TotalSteps1000
##
## Model specification:
## - model M0 (null): mood.cat ~ (1|ID)
## - model M1: mood.cat ~ TotalSteps1000.mc + (1|ID)
## - model M2: mood.cat ~ TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M3: mood.cat ~ TotalSteps1000.mc + weekday.sleep + insomnia + (1|ID)
## - model M4: mood.cat ~ TotalSteps1000.mc + weekday.sleep + insomnia + sex + (1|ID)
## - model M5: mood.cat ~ TotalSteps1000.mc + weekday.sleep + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of mood.cat on 3949 observations from 92 participants
## using the binomial family with the logit link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:##
## Printing Variance Inflation Factors (VIF):
## TotalSteps1000.mc weekday.sleep insomnia sex
## 1.010219 1.010192 1.517389 1.928381
## insomnia:sex
## 2.342281
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 4| 2960.838| 2985.962| -1476.419| 2952.838| NA| NA| NA|
## |insomnia | 5| 2962.826| 2994.232| -1476.413| 2952.826| 0.011| 1| 0.915|
## |sex | 6| 2960.979| 2998.666| -1474.489| 2948.979| 3.848| 1| 0.149|
## |sex:insomnia | 7| 2960.727| 3004.695| -1473.363| 2946.727| 2.252| 1| 0.200|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.73
## - insomnia vs. Baseline model = 0.27
## - sex = 0.405
## - sex:insomnia = 0.315
## Selected model based on Aw: sex:insomnia
##
## - Printing estimated coefficients for models weekday.sleep , insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:----------|:------|:----------|:------|:----------|:-----|:----------|:-----|
## |2 |(Intercept) |0.13(0.02) |-15.61 |0.13(0.02) |-11.11 |0.16(0.03) |-9.00 |0.14(0.03) |-8.81 |
## |3 |TotalSteps1000.mc |0.96(0.01) |-2.66 |0.96(0.01) |-2.66 |0.96(0.01) |-2.67 |0.96(0.01) |-2.66 |
## |4 |weekday.sleep [weekend] |0.74(0.08) |-2.68 |0.74(0.08) |-2.68 |0.74(0.08) |-2.68 |0.74(0.08) |-2.68 |
## |7 |insomnia [1] | | |0.97(0.24) |-0.11 |0.95(0.23) |-0.23 |1.23(0.37) |0.69 |
## |8 |sex [M] | | | | |0.60(0.16) |-1.97 |0.88(0.32) |-0.37 |
## |9 |insomnia [1] * sex [M] | | | | | | |0.46(0.24) |-1.51 |
## |10 |N |92 ID |92 ID |92 ID |92 ID |NA |NA |NA |NA |
## |11 |Observations |3949 |3949 |3949 |3949 |NA |NA |NA |NA |
## |5 |SD (Intercept) |2.86(NA) | |2.86(NA) | |2.81(NA) | |2.79(NA) | |
## |6 |SD (Observations) |2.72(NA) | |2.72(NA) | |2.72(NA) | |2.72(NA) | |
##
##
## - Plotting effect(s) estimated by model 3## Data were 'prettified'. Consider using `terms="TotalSteps1000.mc [all]"` to get smooth plots.
# keeping key results
res <- rbind(res,cbind(measure="mood.cat",glmerAn(long=ema,wide=demos,resp="mood.cat",fix.eff=predictors[1:(length(predictors)-1)],
modelType=c("GLMER"),mc.predictors="TotalSteps1000",family="binomial",link="logit",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",dot.size=3,messages=FALSE)))Comments:
the assumption of (1) linear independence between logit and continuous predictors generally holds; (2) random effects are quite normally distributed, with only a slight deviation from normality in the lower tail of the random intercept distribution; there is (3) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT or stronger evidence in terms of Aw compared to the baseline model
substantial effects are only estimated for
TotalSteps1000(lower badmoodfor higher than usualTotalSteps), andweekday.sleep(lower badmoodon weekend days than on weekdays)
2.2. Aggregate scores
Here, we replicate the analyses by considering the aggregate (i.e., mean) score of stress, worry, and mood (PsyDist), and the factor scores predicted by the MCFA model selected above.
PsyDist
PsyDist is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="PsyDist",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="weekday.sleep",p.adjust.method="BH",
coeff.models=c("insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of PsyDist ...
##
## Preparing data...
## - Excluding 2270 incomplete observations ( 36.5 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Mean-centering TotalSteps1000
##
## Model specification:
## - model M0 (null): PsyDist ~ (1|ID)
## - model M1: PsyDist ~ TotalSteps1000.mc + (1|ID)
## - model M2: PsyDist ~ TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M3: PsyDist ~ TotalSteps1000.mc + weekday.sleep + insomnia + (1|ID)
## - model M4: PsyDist ~ TotalSteps1000.mc + weekday.sleep + insomnia + sex + (1|ID)
## - model M5: PsyDist ~ TotalSteps1000.mc + weekday.sleep + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of PsyDist on 3949 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## TotalSteps1000.mc weekday.sleep insomnia sex
## 1.011793 1.011798 1.562056 1.994809
## insomnia:sex
## 2.500507
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 5| 8395.483| 8426.889| -4192.742| 8385.483| NA| NA| NA|
## |insomnia | 6| 8397.453| 8435.140| -4192.726| 8385.453| 0.030| 1| 0.862|
## |sex | 7| 8396.708| 8440.676| -4191.354| 8382.708| 2.745| 1| 0.162|
## |sex:insomnia | 8| 8396.121| 8446.371| -4190.061| 8380.121| 2.586| 1| 0.162|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.728
## - insomnia vs. Baseline model = 0.272
## - sex = 0.283
## - sex:insomnia = 0.275
## Selected model based on Aw: Baseline
##
## - Printing estimated coefficients for models insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:-----|:-----------|:-----|:-----------|:-----|
## |2 |(Intercept) |2.33(0.07) |31.21 |2.40(0.08) |28.58 |2.34(0.09) |25.48 |
## |3 |TotalSteps1000.mc |-0.01(0.00) |-4.85 |-0.01(0.00) |-4.85 |-0.01(0.00) |-4.85 |
## |4 |weekday.sleep [weekend] |-0.21(0.02) |-8.40 |-0.21(0.02) |-8.40 |-0.21(0.02) |-8.40 |
## |5 |insomnia [1] |0.02(0.10) |0.17 |0.01(0.10) |0.11 |0.13(0.13) |1.06 |
## |8 |sex [M] | | |-0.18(0.11) |-1.67 |-0.01(0.15) |-0.06 |
## |9 |insomnia [1] * sex [M] | | | | |-0.34(0.21) |-1.62 |
## |10 |N |92 ID |92 ID |92 ID |NA |NA |NA |
## |11 |Observations |3949 |3949 |3949 |NA |NA |NA |
## |6 |SD (Intercept) |0.49(NA) | |0.48(NA) | |0.47(NA) | |
## |7 |SD (Observations) |0.82(NA) | |0.82(NA) | |0.82(NA) | |
##
##
## - Plotting effect(s) estimated by model 3
# keeping key results
res <- rbind(res,cbind(measure="PsyDist",glmerAn(long=ema,wide=demos,resp="PsyDist",fix.eff=predictors[1:(length(predictors)-1)],
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
both (1) residuals and (2) random effects are quite normally distributed (despite some deviations in the upper tail of the residuals distribution) with (3) homogeneity between the variances across sex and insomnia levels, and (4) no substantial linear relationship between residuals and fitted values; there is (5) not strong evidence of multicollinearity as all VIFs are 2 or less, with the only exception of the interactive term
none of the specified models showed a significant LRT or stronger evidence in terms of Aw compared to the baseline model
substantial effects are only estimated for
TotalSteps1000(lowerPsyDistfor higher than usualTotalSteps), andweekday.sleep(lowerPsyDiston weekend days than on weekdays). Neitherinsomnianorsex, or their interaction show a substantial effect
fs
fs is quite normally distributed, thus it is modeled by assuming a normal distribution for the residuals.
glmerAn(long=ema,wide=demos,resp="fs",fix.eff=predictors,
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
outputs=c("fit","mComp","coeff","plotEff"),mComp.baseline="weekday.sleep",p.adjust.method="BH",
coeff.models=c("insomnia","sex","sex:insomnia"), # showing all models including insomnia
plot.model="insomnia",show.data=TRUE)## Running GLMER analysis of fs ...
##
## Preparing data...
## - Excluding 2270 incomplete observations ( 36.5 % ) in the response var. or in any of the predictors,
## and 1 participants with no complete variables
## - Mean-centering TotalSteps1000
##
## Model specification:
## - model M0 (null): fs ~ (1|ID)
## - model M1: fs ~ TotalSteps1000.mc + (1|ID)
## - model M2: fs ~ TotalSteps1000.mc + weekday.sleep + (1|ID)
## - model M3: fs ~ TotalSteps1000.mc + weekday.sleep + insomnia + (1|ID)
## - model M4: fs ~ TotalSteps1000.mc + weekday.sleep + insomnia + sex + (1|ID)
## - model M5: fs ~ TotalSteps1000.mc + weekday.sleep + insomnia + sex + sex:insomnia + (1|ID)
##
## Fitting GLMER models of fs on 3949 observations from 92 participants
## using the normal family with the identity link function using ML estimator...
##
## Generating models outputs...
##
## - Plotting diagnostics of the most complex model:
##
## Printing Variance Inflation Factors (VIF):
## TotalSteps1000.mc weekday.sleep insomnia sex
## 1.011797 1.011799 1.562363 1.997162
## insomnia:sex
## 2.500309
##
##
## - Running likelihood ratio test: (applying BH p-values correction)
## Selected model: Baseline
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 5| 9752.824| 9784.231| -4871.412| 9742.824| NA| NA| NA|
## |insomnia | 6| 9753.856| 9791.543| -4870.928| 9741.856| 0.968| 1| 0.325|
## |sex | 7| 9751.857| 9795.826| -4868.929| 9737.857| 3.999| 1| 0.137|
## |sex:insomnia | 8| 9751.532| 9801.782| -4867.766| 9735.532| 2.325| 1| 0.191|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia : = 0.626
## - insomnia vs. Baseline model = 0.374
## - sex = 0.504
## - sex:insomnia = 0.372
## Selected model based on Aw: sex:insomnia
##
## - Printing estimated coefficients for models insomnia , sex , sex:insomnia
##
## | |Predictors |B (SE) |t |B (SE) |t |B (SE) |t |
## |:--|:-----------------------|:-----------|:-----|:-----------|:-----|:-----------|:-----|
## |2 |(Intercept) |-0.04(0.14) |-0.29 |0.11(0.16) |0.72 |-0.00(0.17) |-0.01 |
## |3 |TotalSteps1000.mc |-0.02(0.00) |-4.83 |-0.02(0.00) |-4.83 |-0.02(0.00) |-4.83 |
## |4 |weekday.sleep [weekend] |-0.24(0.03) |-8.41 |-0.24(0.03) |-8.41 |-0.24(0.03) |-8.40 |
## |5 |insomnia [1] |0.20(0.20) |0.99 |0.18(0.20) |0.93 |0.40(0.24) |1.67 |
## |8 |sex [M] | | |-0.41(0.20) |-2.02 |-0.10(0.28) |-0.37 |
## |9 |insomnia [1] * sex [M] | | | | |-0.62(0.40) |-1.54 |
## |10 |N |92 ID |92 ID |92 ID |NA |NA |NA |
## |11 |Observations |3949 |3949 |3949 |NA |NA |NA |
## |6 |SD (Intercept) |0.95(NA) | |0.93(NA) | |0.91(NA) | |
## |7 |SD (Observations) |0.89(NA) | |0.89(NA) | |0.89(NA) | |
##
##
## - Plotting effect(s) estimated by model 3
# keeping key results
res <- rbind(res,cbind(measure="fs",glmerAn(long=ema,wide=demos,resp="fs",fix.eff=predictors[1:(length(predictors)-1)],
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE)))Comments:
- results are highly consistent with those reported for
PsyDist, with the only difference thatsexshows higher Aw and a substantial effect
2.3. Summary of results
Here, we use the key.resPlot function to graphically summarize the results obtained from the models above. The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering model M2 (main effect of insomnia). The figure below uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, and |t| > 1.96 ) are shown in green
those relationships showing a |t| > 1.96 but either a nonsignificant LRT or a lower Aw than the previous model are shown in light green (or lime)
those showing a |t| > 1.96 but both nonsignificant LRT and lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT or a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT and lower Aw than the previous models are shown in red
key.resPlot(res,levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
- the analysis of pre-sleep
stress,worry, andmoodratings characterization did not highlight any substantial difference between the insomnia and the control group, with only a (slightly) higher Aw (but non-significant LRT and unsubstantial effect) for the model includinginsomniafor predictingstress
2.4. Robustness checks
Here, we conduct several robustness checks to account for the arbitrariness of our data processing and analytical choices, including (1) the consideration of insomnia.group instead of the insomnia variable, (2) the inclusion of the covid19 variable as a further covariate, (3) the exclusion of participants with extreme missing data, (4) the exclusion of TotalSteps1000, (5) the use of alternative family distributions for the analyses.
2.4.1. insomnia.group
Here, we inspect the results obtained by considering the insomnia.group (i.e., 46 controls, 26 DSM.ins and 21 sub.ins) rather than the insomnia variable (i.e., 46 controls vs. 47 insomnia).
The figure below shows the t-values associated with the group differences between full DSM insomnia and controls for each sleep outcome, by considering either model M2 (main effect of insomnia) or model M3 (main effect of insomnia and sex) when sex differences were substantial (i.e., HR_NREM and HR_REM. The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# computing key results for stress
PREdictors <- predictors[1:(length(predictors)-2)] # excluding sex and interactive model
res2 <- cbind(measure="stress.cat",glmerAn(long=ema,wide=demos,resp="stress.cat",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
family="binomial",link="logit",mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia.group",key.model="insomnia.group",messages=FALSE))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist","fs")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("normal","identity") }else{ fam <- c("binomial","logit") }
res2 <- rbind(res2,
cbind(measure=diaryVar,glmerAn(long=ema,wide=demos,resp=diaryVar,fix.eff=gsub("insomnia","insomnia.group",PREdictors),
family=fam[1],link=fam[2],key.model="insomnia.group",key.predictor="insomnia.group",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res$robCheck <- "Original"
res2$robCheck <- "insomnia.group"
RES <- rbind(res,res2)
# plotting robustness check
key.resPlot(RES,robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
results are partially inconsistent with the main analyses for most sleep outcomes
the most relevant differences are shown for
stress.cat,worry.cat, andfs, now showing significant LRT, higher Aw and a substantial difference forinsomnia(i.e., higherstress,worry, andfsfor the full DSM insomnia group compared to controls)
Here, we can see that the differences in stress, worry, and fs are substantial, and driven by the full DSM insomnia subgroup, whereas no substantial differences are estimated between the control and the sub-threshold insomnia group.
plot_grid(list(glmerAn(long=ema,wide=demos,resp="stress.cat",p.adjust.method="BH",mComp.baseline="weekday.sleep",
fix.eff=gsub("insomnia","insomnia.group",predictors)[1:(length(predictors)-2)],
mc.predictors="TotalSteps1000",modelType="GLMER",family="binomial",link="logit",return.plot=TRUE,
outputs=c("mComp","coeff","plotEff"),coeff.models="insomnia.group",transform="exp",
plot.model="insomnia.group",plot.pred="insomnia.group",dot.size=3,message=FALSE),
glmerAn(long=ema,wide=demos,resp="worry.cat",p.adjust.method="BH",mComp.baseline="weekday.sleep",
fix.eff=gsub("insomnia","insomnia.group",predictors)[1:(length(predictors)-2)],
mc.predictors="TotalSteps1000",modelType="GLMER",family="binomial",link="logit",return.plot=TRUE,
outputs=c("mComp","coeff","plotEff"),coeff.models="insomnia.group",transform="exp",
plot.model="insomnia.group",plot.pred="insomnia.group",dot.size=3,message=FALSE),
glmerAn(long=ema,wide=demos,resp="mood.cat",p.adjust.method="BH",mComp.baseline="weekday.sleep",
fix.eff=gsub("insomnia","insomnia.group",predictors)[1:(length(predictors)-2)],
mc.predictors="TotalSteps1000",modelType="GLMER",family="binomial",link="logit",return.plot=TRUE,
outputs=c("mComp","coeff","plotEff"),coeff.models="insomnia.group",transform="exp",
plot.model="insomnia.group",plot.pred="insomnia.group",dot.size=3,message=FALSE)))##
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:--------------|----:|--------:|--------:|---------:|--------:|------:|--:|----------:|
## |Baseline | 4| 2547.329| 2572.454| -1269.665| 2539.329| NA| NA| NA|
## |insomnia.group | 6| 2541.326| 2579.013| -1264.663| 2529.326| 10.003| 2| 0.007|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia.group : = 0.047
## - insomnia.group vs. Baseline model = 0.953
## Selected model based on Aw: insomnia.group
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:----------|:------|
## |2 |(Intercept) |0.08(0.01) |-13.79 |
## |3 |TotalSteps1000.mc |0.95(0.01) |-3.25 |
## |4 |weekday.sleep [weekend] |0.67(0.09) |-3.17 |
## |5 |insomnia.group [DSM.ins] |2.35(0.67) |3.00 |
## |6 |insomnia.group [sub.ins] |0.95(0.31) |-0.15 |
## |9 |N ID |92 |NA |
## |10 |Observations |3949 |NA |
## |7 |SD (Intercept) |2.74(NA) | |
## |8 |SD (Observations) |2.72(NA) | |
##
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:--------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 4| 2468.627| 2493.752| -1230.314| 2460.627| NA| NA| NA|
## |insomnia.group | 6| 2464.192| 2501.879| -1226.096| 2452.192| 8.436| 2| 0.015|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia.group : = 0.098
## - insomnia.group vs. Baseline model = 0.902
## Selected model based on Aw: insomnia.group
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:----------|:------|
## |2 |(Intercept) |0.08(0.02) |-11.64 |
## |3 |TotalSteps1000.mc |0.99(0.01) |-0.92 |
## |4 |weekday.sleep [weekend] |0.62(0.08) |-3.67 |
## |5 |insomnia.group [DSM.ins] |2.26(0.78) |2.35 |
## |6 |insomnia.group [sub.ins] |0.70(0.28) |-0.89 |
## |9 |N ID |92 |NA |
## |10 |Observations |3949 |NA |
## |7 |SD (Intercept) |3.54(NA) | |
## |8 |SD (Observations) |2.72(NA) | |
##
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:--------------|----:|--------:|--------:|---------:|--------:|-----:|--:|----------:|
## |Baseline | 4| 2960.838| 2985.962| -1476.419| 2952.838| NA| NA| NA|
## |insomnia.group | 6| 2957.211| 2994.898| -1472.606| 2945.211| 7.626| 2| 0.022|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Baseline vs. insomnia.group : = 0.14
## - insomnia.group vs. Baseline model = 0.86
## Selected model based on Aw: insomnia.group
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:----------|:------|
## |2 |(Intercept) |0.14(0.02) |-11.38 |
## |3 |TotalSteps1000.mc |0.96(0.01) |-2.66 |
## |4 |weekday.sleep [weekend] |0.74(0.08) |-2.67 |
## |5 |insomnia.group [DSM.ins] |1.44(0.40) |1.30 |
## |6 |insomnia.group [sub.ins] |0.55(0.17) |-1.89 |
## |9 |N ID |92 |NA |
## |10 |Observations |3949 |NA |
## |7 |SD (Intercept) |2.75(NA) | |
## |8 |SD (Observations) |2.72(NA) | |
plot_grid(list(glmerAn(long=ema,wide=demos,resp="PsyDist",p.adjust.method="BH",
fix.eff=gsub("insomnia","insomnia.group",predictors)[1:(length(predictors)-2)],
mc.predictors="TotalSteps1000",modelType="GLMER",return.plot=TRUE,
outputs=c("mComp","coeff","plotEff"),coeff.models="insomnia.group",
plot.model="insomnia.group",plot.pred="insomnia.group",show.data=TRUE,message=FALSE),
glmerAn(long=ema,wide=demos,resp="fs",p.adjust.method="BH",
fix.eff=gsub("insomnia","insomnia.group",predictors)[1:(length(predictors)-2)],
mc.predictors="TotalSteps1000",modelType="GLMER",return.plot=TRUE,
outputs=c("mComp","coeff","plotEff"),coeff.models="insomnia.group",
plot.model="insomnia.group",plot.pred="insomnia.group",show.data=TRUE,message=FALSE)))##
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:-----------------|----:|--------:|--------:|---------:|--------:|------:|--:|----------:|
## |Null model | 3| 8494.144| 8512.988| -4244.072| 8488.144| NA| NA| NA|
## |TotalSteps1000.mc | 4| 8463.324| 8488.449| -4227.662| 8455.324| 32.820| 1| 0.000|
## |weekday.sleep | 5| 8395.483| 8426.889| -4192.742| 8385.483| 69.841| 1| 0.000|
## |insomnia.group | 7| 8391.267| 8435.236| -4188.634| 8377.267| 8.216| 2| 0.016|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Null model vs. TotalSteps1000.mc : = 0
## - TotalSteps1000.mc vs. Null model model = 1
## - weekday.sleep = 1
## - insomnia.group = 0.892
## Selected model based on Aw: insomnia.group
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:-----------|:-----|
## |2 |(Intercept) |2.33(0.07) |32.73 |
## |3 |TotalSteps1000.mc |-0.01(0.00) |-4.85 |
## |4 |weekday.sleep [weekend] |-0.21(0.02) |-8.39 |
## |5 |insomnia.group [DSM.ins] |0.20(0.12) |1.70 |
## |6 |insomnia.group [sub.ins] |-0.21(0.13) |-1.68 |
## |9 |N ID |92 |NA |
## |10 |Observations |3949 |NA |
## |7 |SD (Intercept) |0.46(NA) | |
## |8 |SD (Observations) |0.82(NA) | |
##
##
## | | npar| AIC| BIC| logLik| deviance| Chisq| Df| Pr(>Chisq)|
## |:-----------------|----:|--------:|--------:|---------:|--------:|------:|--:|----------:|
## |Null model | 3| 9851.424| 9870.268| -4922.712| 9845.424| NA| NA| NA|
## |TotalSteps1000.mc | 4| 9820.813| 9845.938| -4906.407| 9812.813| 32.611| 1| 0.000|
## |weekday.sleep | 5| 9752.824| 9784.231| -4871.412| 9742.824| 69.989| 1| 0.000|
## |insomnia.group | 7| 9748.633| 9792.602| -4867.317| 9734.633| 8.191| 2| 0.017|
##
##
## - Computing Aw coefficients for each model vs. all previous models:
## - Null model vs. TotalSteps1000.mc : = 0
## - TotalSteps1000.mc vs. Null model model = 1
## - weekday.sleep = 1
## - insomnia.group = 0.89
## Selected model based on Aw: insomnia.group
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:-----------|:-----|
## |2 |(Intercept) |-0.04(0.14) |-0.30 |
## |3 |TotalSteps1000.mc |-0.02(0.00) |-4.83 |
## |4 |weekday.sleep [weekend] |-0.24(0.03) |-8.40 |
## |5 |insomnia.group [DSM.ins] |0.53(0.23) |2.33 |
## |6 |insomnia.group [sub.ins] |-0.22(0.24) |-0.88 |
## |9 |N ID |92 |NA |
## |10 |Observations |3949 |NA |
## |7 |SD (Intercept) |0.91(NA) | |
## |8 |SD (Observations) |0.89(NA) | |
In conclusion, the results of this robustness check were inconsistent with those obtained with the main analyses, showing substantial differences in stress, worry, and fs between the full DSM insomnia and the control group.
2.4.2. COVID-19 emergency
Here, we inspect the results obtained by including the covid19 variable as a further covariate. As shown in section 1.2.8, the covid19 factor discriminates the data collected pre-COVID19 (63%) and the data collected post-COVID19 (37%). The covariate is included at the first step (model M1, before SO.num).
The figure below shows the t-values associated with the group differences between full DSM insomnia and controls for each sleep outcome, by considering either model M2 (main effect of insomnia) or model M3 (main effect of insomnia and sex) when sex differences were substantial (i.e., HR_NREM and HR_REM. The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=ema,wide=demos,resp="stress.cat",fix.eff=c("covid19",PREdictors),
modelType=c("GLMER"),mc.predictors="TotalSteps1000",family="binomial",link="logit",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist","fs")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("normal","identity") }else{ fam <- c("binomial","logit") }
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=ema,wide=demos,resp=diaryVar,fix.eff=c("covid19",PREdictors),
family=fam[1],link=fam[2],key.model="insomnia",key.predictor="insomnia",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "covid19"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","covid19"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
results are highly consistent with the main analyses for most outcomes
the only difference is that with
covid19the model predictingstressbyinsomniais not even associated with higher Aw
Here, we can see that the post-COVID19 data collection period does not show substantial differences compared to the pre-COVID19 period, with only slightly higher stress.
# plotting difference by covid (model M1)
for(diaryVar in c("stress.cat","worry.cat","mood.cat","PsyDist","fs")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("normal","identity") }else{ fam <- c("binomial","logit") }
glmerAn(long=ema,wide=demos,resp=diaryVar,fix.eff=c("covid19",PREdictors[1:2]),
mc.predictors="TotalSteps1000",modelType="GLMER",family=fam[1],link=fam[2],
outputs="plotEff",plot.model="insomnia",plot.pred="covid19",messages=FALSE,dot.size=3) }
In conclusion, the results of this robustness check were highly consistent with those obtained with the main analyses, with no substantial group differences in any variable, and no substantial differences between pre- and post-COVID19.
As a further check, we also inspect the results of the covid19 robustness check by considering insomnia.group.
# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=ema,wide=demos,resp="stress.cat",
fix.eff=c("covid19",gsub("insomnia","insomnia.group",PREdictors)),
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
family="binomial",link="logit",messages=FALSE,
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia.group",key.model="insomnia.group"))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist","fs")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("normal","identity") }else{ fam <- c("binomial","logit") }
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=ema,wide=demos,resp=diaryVar,
fix.eff=c("covid19",gsub("insomnia","insomnia.group",PREdictors)),
family=fam[1],link=fam[2],key.model="insomnia.group",key.predictor="insomnia.group",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "covid19.insomnia.group"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","covid19.insomnia.group"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
the results are in line with our first sanity check (section 5.2.4.1), showing substantial differences in
stress,worry, andfsbetween the full DSM insomnia and the control groupin addition, substantial differences in
PsyDist(but not inmood) are found when consideringcovid19as a further covariate
2.4.3. Low compliance
Here, we inspect the results obtained by excluding seven participants (7.5%) with less than two weeks of valid observations in any data type. In section 3.2, these participants were marked with majMiss = 1.
The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering model M2 (main effect of insomnia). The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# excluding participants with extreme missing data
EMA <- ema[ema$majMiss!=1,]
DEMOS <- demos[demos$ID%in%levels(as.factor(as.character(EMA$ID))),]
cat("Excluded",nrow(ema)-nrow(EMA),"days, ",nrow(demos)-nrow(DEMOS),"participants")## Excluded 465 days, 7 participants# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=EMA,wide=DEMOS,resp="stress.cat",fix.eff=PREdictors,
modelType=c("GLMER"),mc.predictors="TotalSteps1000", family="binomial",link="logit",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE))
# computing key results for the following sleep outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist","fs")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("normal","identity") }else{ fam <- c("binomial","logit") }
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=EMA,wide=DEMOS,resp=diaryVar,fix.eff=PREdictors,
family=fam[1],link=fam[2],key.model="insomnia",key.predictor="insomnia",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "majMiss"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","majMiss"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
results are highly consistent with the main analyses for most outcomes
the only difference is that without
majMissthe model predictingstressbyinsomniais not even associated with higher Aw
As a further check, we also inspect the results of the majMiss robustness check by considering insomnia.group.
# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=EMA,wide=DEMOS,resp="stress.cat",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
family="binomial",link="logit",messages=FALSE,
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia.group",key.model="insomnia.group"))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist","fs")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("normal","identity") }else{ fam <- c("binomial","logit") }
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=EMA,wide=DEMOS,resp=diaryVar,
fix.eff=gsub("insomnia","insomnia.group",PREdictors),
family=fam[1],link=fam[2],key.model="insomnia.group",key.predictor="insomnia.group",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "majMiss.insomnia.group"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","majMiss.insomnia.group"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
- the results are highly consistent with our first sanity check (section 5.2.4.1), showing substantial differences in
stress,worry, andfsbetween the full DSM insomnia and the control group
2.4.4. TotalStep1000
Here, we inspect the results obtained by excluding TotalSteps1000 from the models predictors. This is done in order to replicate the analyses by including cases with invalid daily steps data (thus, with higher sample size).
The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering model M2 (main effect of insomnia). The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=ema,wide=demos,resp="stress.cat",fix.eff=c("weekday.sleep","insomnia","sex"),
modelType=c("GLMER"),mc.predictors="TotalSteps1000",family="binomial",link="logit",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist","fs")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("normal","identity") }else{ fam <- c("binomial","logit") }
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=ema,wide=demos,resp=diaryVar,fix.eff=c("weekday.sleep","insomnia","sex"),
family=fam[1],link=fam[2],key.model="insomnia",key.predictor="insomnia",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "no-TotalSteps"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","no-TotalSteps"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
results are overall consistent with the main analyses, with the only exception of
stress.cat, whose model includinginsomnianow shows significant LRT and a substantial positive effect (i.e., higherstressin insomnia compared to controls)none of the remaining variables show substantial differences, with almost identical results than those obtained with the main analyses
As a further check, we also inspect the results of the majMiss robustness check by considering insomnia.group.
# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=EMA,wide=DEMOS,resp="stress.cat",fix.eff=c("weekday.sleep","insomnia.group","sex"),
modelType=c("GLMER"),mc.predictors="TotalSteps1000",family="binomial",link="logit",messages=FALSE,
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia.group",key.model="insomnia.group"))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist","fs")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("normal","identity") }else{ fam <- c("binomial","logit") }
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=EMA,wide=DEMOS,resp=diaryVar,fix.eff=c("weekday.sleep","insomnia.group","sex"),
family=fam[1],link=fam[2],key.model="insomnia.group",key.predictor="insomnia.group",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "no-TotalSteps.insomnia.group"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","no-TotalSteps.insomnia.group"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
- the results are highly consistent with our first sanity check (section 5.2.4.1), showing substantial differences in
stress,worry, andfsbetween the full DSM insomnia and the control group
2.4.5. Family distribution
Here, we replicate the models predicting those variables that showed poor fit in the model diagnostics above by re-specifying the models with different distribution families and link functions:
NORMAL
stress,worry, andmood: initially dichotomized, and modeled with a logistic regression. Here, we can see that the fit of the normal istribution (i.e., assuming item scores as they were continuous) is quite acceptable and better than that shown by the gamma family
# stress: normal vs. gamma (unoptimal fit in both cases)
p1 <- lmer(stress~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema)
p2 <- glmer(stress~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("log"),nAGQ=0) # nAGQ = 0 to converge
p3 <- glmer(stress~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("identity"),nAGQ=0) # nAGQ = 0 to converge
p4 <- glmer(stress~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("inverse"),nAGQ=0) # nAGQ = 0 to converge
par(mfrow=c(2,2))
qqnorm(resid(p1),main="qqnorm residuals - Normal"); qqline(resid(p1),col="red")
qqnorm(resid(p2,type="deviance"),main="qqnorm residuals - Gamma(log)"); qqline(resid(p2),col="red")
qqnorm(resid(p3,type="deviance"),main="qqnorm residuals - Gamma(identity)"); qqline(resid(p3),col="red")
qqnorm(resid(p4,type="deviance"),main="qqnorm residuals - Gamma(inverse)"); qqline(resid(p4),col="red") 
# worry: normal vs. gamma (unoptimal fit in both cases)
p1 <- lmer(worry~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema)
p2 <- glmer(worry~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("log"),nAGQ=0) # nAGQ = 0 to converge
p3 <- glmer(worry~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("identity"),nAGQ=0) # nAGQ = 0 to converge
p4 <- glmer(worry~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("inverse"),nAGQ=0) # nAGQ = 0 to converge
par(mfrow=c(2,2))
qqnorm(resid(p1),main="qqnorm residuals - Normal"); qqline(resid(p1),col="red")
qqnorm(resid(p2,type="deviance"),main="qqnorm residuals - Gamma(log)"); qqline(resid(p2),col="red")
qqnorm(resid(p3,type="deviance"),main="qqnorm residuals - Gamma(identity)"); qqline(resid(p3),col="red")
qqnorm(resid(p4,type="deviance"),main="qqnorm residuals - Gamma(inverse)"); qqline(resid(p4),col="red") 
# mood: normal vs. gamma (unoptimal fit in both cases)
p1 <- lmer(mood~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema)
p2 <- glmer(mood~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("log"),nAGQ=0) # nAGQ = 0 to converge
p3 <- glmer(mood~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("identity"),nAGQ=0) # nAGQ = 0 to converge
p4 <- glmer(mood~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("inverse"),nAGQ=0) # nAGQ = 0 to converge
par(mfrow=c(2,2))
qqnorm(resid(p1),main="qqnorm residuals - Normal"); qqline(resid(p1),col="red")
qqnorm(resid(p2,type="deviance"),main="qqnorm residuals - Gamma(log)"); qqline(resid(p2),col="red")
qqnorm(resid(p3,type="deviance"),main="qqnorm residuals - Gamma(identity)"); qqline(resid(p3),col="red")
qqnorm(resid(p4,type="deviance"),main="qqnorm residuals - Gamma(inverse)"); qqline(resid(p4),col="red") 
PsyDist: modeled by assuming a normal distribution for the residuals. Here, we re-specify the models by assuming a gamma distributon, which shows a slightly better fit (note thatPsyDistis positively skewed)
# PsyDist: normal vs. gamma (better fit for gamma)
p1 <- lmer(PsyDist~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema)
p2 <- glmer(PsyDist~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("log"),nAGQ=0)
p3 <- glmer(PsyDist~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("identity"))
p4 <- glmer(PsyDist~TotalSteps1000+weekday.sleep+insomnia+(1|ID),data=ema,family=Gamma("inverse"),nAGQ=0)
par(mfrow=c(2,2))
qqnorm(resid(p1),main="qqnorm residuals - Normal"); qqline(resid(p1),col="red")
qqnorm(resid(p2,type="deviance"),main="qqnorm residuals - Gamma(log)"); qqline(resid(p2),col="red")
qqnorm(resid(p3,type="deviance"),main="qqnorm residuals - Gamma(identity)"); qqline(resid(p3),col="red")
qqnorm(resid(p4,type="deviance"),main="qqnorm residuals - Gamma(inverse)"); qqline(resid(p4),col="red") 
The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering either M2 (main effect of insomnia), and using a normal and a gamma distribution for modeling raw item scores and the aggregate score, respectively. The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=ema,wide=demos,resp="stress",fix.eff=PREdictors,
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("gamma","log") }else{ fam <- c("normal","identity") }
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=ema,wide=demos,resp=gsub(".cat","",diaryVar),fix.eff=PREdictors,
family=fam[1],link=fam[2],nAGQ=0, # nAGQ = 0 to reach convergence with PsyDist
key.model="insomnia",key.predictor="insomnia",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "distr.Normal.Gamma"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","distr.Normal.Gamma"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
- results are highly consistent with the main analyses
As a further check, we also inspect the results of the distr.Normal.Gamma robustness check by considering insomnia.group.
# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=ema,wide=demos,resp="stress",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia.group",key.model="insomnia.group",messages=FALSE))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist")){
if(diaryVar=="PsyDist"){ fam <- c("gamma","log") }else{ fam <- c("normal","identity") }
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=ema,wide=demos,resp=gsub(".cat","",diaryVar),
fix.eff=gsub("insomnia","insomnia.group",PREdictors),
family=fam[1],link=fam[2],nAGQ=0, # nAGQ = 0 to reach convergence with PsyDist
key.model="insomnia.group",key.predictor="insomnia.group",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "distr.Normal.Gamma.insomnia.group"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","distr.Normal.Gamma.insomnia.group"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
- the results are highly consistent with our first sanity check (section 5.2.4.1), showing substantial differences in
stress,worry, and (also)PsyDistbetween the full DSM insomnia and the control group
ORDINAL REGRESSION
Here, we account for the ordinal nature of item ratings by replicating the analyses by using cumulative link mixed models (CLMM). The figure below uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=ema,wide=demos,resp="stress",fix.eff=PREdictors,
modelType=c("CLMM"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat")){
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=ema,wide=demos,resp=gsub(".cat","",diaryVar),fix.eff=PREdictors,
mComp.baseline="weekday.sleep",p.adjust.method="BH",
modelType=c("CLMM"),mc.predictors="TotalSteps1000",messages=FALSE,
key.model="insomnia",key.predictor="insomnia",outputs="key.results"))) }
# joining with the main results
res2$robCheck <- "distr.CLMM"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","distr.CLMM"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
- results are highly consistent with the main analyses
As a further check, we also inspect the results of the distr.CLMM robustness check by considering insomnia.group.
# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=ema,wide=demos,resp="stress",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
modelType=c("CLMM"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia.group",key.model="insomnia.group",messages=FALSE))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("gamma","log") }else{ fam <- c("normal","identity") }
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=ema,wide=demos,resp=gsub(".cat","",diaryVar),
fix.eff=gsub("insomnia","insomnia.group",PREdictors),
key.model="insomnia.group",key.predictor="insomnia.group",
modelType=c("CLMM"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "distr.CLMM.insomnia.group"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","distr.CLMM.insomnia.group"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
- the results are highly consistent with our first sanity check (section 5.2.4.1), showing substantial differences in
stressandworrybetween the full DSM insomnia and the control group
Here, we inspect the differences in the insomnia effect estimated by the CLMM models specified with the insomnia and the insomnia.group predictor. The figures below show that, whereas no substantial differences are estimated between the insomnia and the control group in almost any response category, substantial differences are estimated in stress and worry between the full DSM insomnia and the control group, with the former showing lower probability of rating stress and mood as 1, and higher probability of rating them as 3, compared to the latter.
# plotting stress by insomnia (CLMM)
glmerAn(long=ema,wide=demos,resp="stress",fix.eff=PREdictors,modelType=c("CLMM"),mc.predictors="TotalSteps1000",
outputs=c("coeff","plotEff"),plot.model="insomnia",plot.pred="insomnia",coeff.models="insomnia",messages=FALSE)##
##
## | |Predictors |B (SE) |t |
## |:--|:-----------------------|:-----------|:-----|
## |2 |1|2 |-0.96(0.20) |-4.76 |
## |3 |2|3 |0.73(0.20) |3.64 |
## |4 |3|4 |2.53(0.21) |12.25 |
## |5 |4|5 |4.36(0.23) |19.04 |
## |6 |TotalSteps1000.mc |-0.02(0.01) |-2.40 |
## |7 |weekday.sleep [weekend] |-0.48(0.07) |-6.83 |
## |8 |insomnia [1] |0.20(0.28) |0.72 |
## |11 |N ID |92 |NA |
## |12 |Observations |3949 |NA |
## |9 |SD (Intercept) |1.30(NA) | |
## |10 |SD (Observations) |1.47(NA) | |
# plotting stress by insomnia.group (CLMM)
glmerAn(long=ema,wide=demos,resp="stress",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
modelType=c("CLMM"),mc.predictors="TotalSteps1000",messages=FALSE,
outputs=c("coeff","plotEff"),plot.model="insomnia.group",plot.pred="insomnia.group",coeff.models="insomnia.group")##
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:-----------|:-----|
## |2 |1|2 |-0.97(0.19) |-5.00 |
## |3 |2|3 |0.73(0.19) |3.80 |
## |4 |3|4 |2.53(0.20) |12.80 |
## |5 |4|5 |4.36(0.22) |19.73 |
## |6 |TotalSteps1000.mc |-0.02(0.01) |-2.40 |
## |7 |weekday.sleep [weekend] |-0.48(0.07) |-6.82 |
## |8 |insomnia.group [DSM.ins] |0.70(0.31) |2.21 |
## |9 |insomnia.group [sub.ins] |-0.44(0.34) |-1.28 |
## |12 |N ID |92 |NA |
## |13 |Observations |3949 |NA |
## |10 |SD (Intercept) |1.23(NA) | |
## |11 |SD (Observations) |1.47(NA) | |
# plotting worry by insomnia (CLMM)
glmerAn(long=ema,wide=demos,resp="worry",fix.eff=PREdictors,modelType=c("CLMM"),mc.predictors="TotalSteps1000",
outputs=c("coeff","plotEff"),plot.model="insomnia",plot.pred="insomnia",coeff.models="insomnia",messages=FALSE)##
##
## | |Predictors |B (SE) |t |
## |:--|:-----------------------|:-----------|:-----|
## |2 |1|2 |-1.12(0.22) |-5.13 |
## |3 |2|3 |0.54(0.22) |2.49 |
## |4 |3|4 |2.54(0.22) |11.48 |
## |5 |4|5 |4.35(0.24) |18.00 |
## |6 |TotalSteps1000.mc |-0.02(0.01) |-2.77 |
## |7 |weekday.sleep [weekend] |-0.54(0.07) |-7.75 |
## |8 |insomnia [1] |0.17(0.30) |0.57 |
## |11 |N ID |92 |NA |
## |12 |Observations |3949 |NA |
## |9 |SD (Intercept) |1.40(NA) | |
## |10 |SD (Observations) |1.48(NA) | |
# plotting worry by insomnia.group (CLMM)
glmerAn(long=ema,wide=demos,resp="worry",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
modelType=c("CLMM"),mc.predictors="TotalSteps1000",messages=FALSE,
outputs=c("coeff","plotEff"),plot.model="insomnia.group",plot.pred="insomnia.group",coeff.models="insomnia.group")##
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:-----------|:-----|
## |2 |1|2 |-1.12(0.21) |-5.34 |
## |3 |2|3 |0.54(0.21) |2.58 |
## |4 |3|4 |2.54(0.21) |11.92 |
## |5 |4|5 |4.34(0.23) |18.57 |
## |6 |TotalSteps1000.mc |-0.02(0.01) |-2.77 |
## |7 |weekday.sleep [weekend] |-0.54(0.07) |-7.74 |
## |8 |insomnia.group [DSM.ins] |0.67(0.34) |1.97 |
## |9 |insomnia.group [sub.ins] |-0.47(0.37) |-1.26 |
## |12 |N ID |92 |NA |
## |13 |Observations |3949 |NA |
## |10 |SD (Intercept) |1.34(NA) | |
## |11 |SD (Observations) |1.48(NA) | |
# plotting mood by insomnia (CLMM)
glmerAn(long=ema,wide=demos,resp="mood",fix.eff=PREdictors,modelType=c("CLMM"),mc.predictors="TotalSteps1000",
outputs=c("coeff","plotEff"),plot.model="insomnia",plot.pred="insomnia",coeff.models="insomnia",messages=FALSE)##
##
## | |Predictors |B (SE) |t |
## |:--|:-----------------------|:-----------|:------|
## |2 |1|2 |-2.10(0.20) |-10.47 |
## |3 |2|3 |0.09(0.20) |0.44 |
## |4 |3|4 |1.89(0.20) |9.44 |
## |5 |4|5 |3.85(0.22) |17.48 |
## |6 |TotalSteps1000.mc |-0.05(0.01) |-5.64 |
## |7 |weekday.sleep [weekend] |-0.40(0.07) |-5.81 |
## |8 |insomnia [1] |-0.48(0.27) |-1.75 |
## |11 |N ID |92 |NA |
## |12 |Observations |3949 |NA |
## |9 |SD (Intercept) |1.27(NA) | |
## |10 |SD (Observations) |1.52(NA) | |
# plotting mood by insomnia.group (CLMM)
glmerAn(long=ema,wide=demos,resp="mood",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
modelType=c("CLMM"),mc.predictors="TotalSteps1000",messages=FALSE,
outputs=c("coeff","plotEff"),plot.model="insomnia.group",plot.pred="insomnia.group",coeff.models="insomnia.group")##
##
## | |Predictors |B (SE) |t |
## |:--|:------------------------|:-----------|:------|
## |2 |1|2 |-2.10(0.20) |-10.77 |
## |3 |2|3 |0.09(0.19) |0.44 |
## |4 |3|4 |1.89(0.19) |9.70 |
## |5 |4|5 |3.85(0.22) |17.88 |
## |6 |TotalSteps1000.mc |-0.05(0.01) |-5.64 |
## |7 |weekday.sleep [weekend] |-0.40(0.07) |-5.81 |
## |8 |insomnia.group [DSM.ins] |-0.09(0.31) |-0.30 |
## |9 |insomnia.group [sub.ins] |-0.97(0.34) |-2.85 |
## |12 |N ID |92 |NA |
## |13 |Observations |3949 |NA |
## |10 |SD (Intercept) |1.23(NA) | |
## |11 |SD (Observations) |1.52(NA) | |
2.4.6. Late responses
Here, we inspect the results obtained by excluding 1,010 cases (20.5%) in which participants filled the diary on the following day. In section 1.2.6, these cases were marked with diary.nextDay = 1.
The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering model M2 (main effect of insomnia). The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
# excluding participants with extreme missing data
EMA <- ema[ema$diary.nextDay!=1,]
DEMOS <- demos[demos$ID%in%levels(as.factor(as.character(EMA$ID))),]
cat("Excluded",nrow(ema)-nrow(EMA),"days, ",nrow(demos)-nrow(DEMOS),"participants")## Excluded 1010 days, 0 participants# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=EMA,wide=DEMOS,resp="stress.cat",fix.eff=PREdictors,
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
family="binomial",link="logit",mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia",key.model="insomnia",messages=FALSE))
# computing key results for the following sleep outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist","fs")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("normal","identity") }else{ fam <- c("binomial","logit") }
if(diaryVar=="worry.cat"){ nAGQ=0 } else { nAGQ=1 } # nAGQ = 0 for worry.cat to reach convergence
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=EMA,wide=DEMOS,resp=diaryVar,fix.eff=PREdictors,nAGQ=nAGQ,
family=fam[1],link=fam[2],key.model="insomnia",key.predictor="insomnia",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "lateResp"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","lateResp"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
- results are highly consistent with the main analyses for most outcomes
As a further check, we also inspect the results of the lateResp robustness check by considering insomnia.group.
# computing key results for stress
res2 <- cbind(measure="stress.cat",glmerAn(long=EMA,wide=DEMOS,resp="stress.cat",fix.eff=gsub("insomnia","insomnia.group",PREdictors),
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
family="binomial",link="logit",messages=FALSE,
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",key.predictor="insomnia.group",key.model="insomnia.group"))
# computing key results for the following diary outcomes
for(diaryVar in c("worry.cat","mood.cat","PsyDist","fs")){
if(diaryVar%in%c("PsyDist","fs")){ fam <- c("normal","identity") }else{ fam <- c("binomial","logit") }
if(diaryVar=="worry.cat"){ nAGQ=0 } else { nAGQ=1 } # nAGQ = 0 for worry.cat to reach convergence
res2 <- rbind(res2,cbind(measure=diaryVar,glmerAn(long=EMA,wide=DEMOS,resp=diaryVar,
fix.eff=gsub("insomnia","insomnia.group",PREdictors),nAGQ=nAGQ,
family=fam[1],link=fam[2],key.model="insomnia.group",key.predictor="insomnia.group",
modelType=c("GLMER"),mc.predictors="TotalSteps1000",
mComp.baseline="weekday.sleep",p.adjust.method="BH",
outputs="key.results",messages=FALSE))) }
# joining with the main results
res2$robCheck <- "lateResp.insomnia.group"
RES <- rbind(RES,res2)
# plotting robustness check
key.resPlot(RES[RES$robCheck%in%c("Original","lateResp.insomnia.group"),],
robCheck = "robCheck",levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
- the results are highly consistent with our first sanity check (section 5.2.4.1), showing substantial differences in
stress,worry, andfsbetween the full DSM insomnia and the control group
2.4.7. Summary of results
Here, we summarize the results from all robustness checks. The figure below shows the t-values associated with the group differences between insomnia and controls for each sleep outcome, by considering model M2 (main effect of insomnia). The figure uses the following color code:
substantial relationships (i.e., significant LRT, Aw higher than previous models, AND |t| > 1.96 ) are shown in green
those showing a |t| > 1.96 but nonsignificant LRT OR lower Aw than the previous model are shown in light green (lime)
those showing a |t| > 1.96 but both nonsignificant LRT AND lower Aw than the previous models are shown in yellow
those showing a |t| < 1.96 and either a nonsignificant LRT OR a lower Aw than the previous models are shown in orange
those showing a |t| < 1.96 and both nonsignificant LRT AND lower Aw than the previous models are shown in red
key.resPlot(RES[!grepl(".insomnia.group",RES$robCheck),],robCheck = "robCheck",
levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
Comments:
the null results observed from the main analyses for
moodandPsyDistare highly consistent with those obtained from the main analysesin contrast, the inclusion of
insomnia.groupresulted in substantial group differences (i.e., higherstress,worry, andfs) between the full DSM insomnia subgroup and the control group, but not between the latter and the insomnia sub-threshold subgroupgroup differences in
stressratings were substantial also whenTotalSteps1000was excluded form the models
Here, we can see that the substantial differences in stress, worry, and fs are highly consistent across all robustness checks including the insomnia.group term.
key.resPlot(RES[grepl(".insomnia.group",RES$robCheck),],robCheck = "robCheck",
levels=c("stress.cat","worry.cat","mood.cat","PsyDist","fs"))
In conclusion, the results obtained with the main analyses are overall robust, but the insomnia.group variable shows substantial and consistent differences that were not observed for the insomnia variable.