# Installing necessary extenral packages # install.packages("ggplot2") # install.packages("generalCorr") # install.packages("CAM") # install.packages("stringr") # install.packages("rJava") # devtools::install_github("bd2kccd/r-causal") # install.packages("bnlearn") # install.packages("medmod") # install.packages("powerMediation") # Loading necessary external packages library(ggplot2) # Loading data file (supplemented to the manuscript as RData file; the name of the file should be changed first) load('Data.RData') # Normality tests (only for exploratory analysis, for all 10 considered parameters, for 2 body positions) for (i in 1:20){ print(i) print(shapiro.test(data[,i])) }; rm(i) # Exploratory figures (Fig. 1-10 in the paper) Lab = c('Heart Rate [BPM]','lnRMSSD','RMSSD [ms]','Respiratory Rate [BPM]','Coefficient of iRR', 'Coefficient of InsT','Coefficient of ExpT','Coefficient of InsV','Coefficient of ExpV','Breathing Regularity [%]') Norm = c(1,1,0,1,0,0,0,0,0,0) # Established manually from normality tests (lines 16-19) # The loop for all parameters for (i in 1:10){ dt <- data.frame(id = c(rep("Supine",100),rep("Standing",100)),value = c(data[,i],data[,i+10])) if(Norm[i]==1){ test <- t.test(dt$value[dt$id=="Supine"],dt$value[dt$id=="Standing"],paired=TRUE,alternative = "two.sided") } else { test <- wilcox.test(dt$value[dt$id=="Supine"],dt$value[dt$id=="Standing"],paired=TRUE,alternative = "two.sided") } if(test$p.value < 0.001){ Label <- "***" } else if(test$p.value < 0.01){ Label <- "**" } else if(test$p.value < 0.05){ Label <- "*" } else if(test$p.value < 0.1){ Label <- "." } else { Label <- "ns" } ggplot(data = dt, aes(x = id, y = value)) + geom_violin(position = position_dodge(width = 0.3),fill="grey96") + ylab(Lab[i]) + geom_boxplot(width=.1, position = position_dodge(width = 0.4)) + theme_bw() + annotate("segment",x=1,xend=2,y=1.048*max(dt$value),yend=1.048*max(dt$value)) + annotate("text",1.5,1.05*max(dt$value),label=Label, fontface = 2,size=6) + scale_x_discrete(limits=c("Supine","Standing")) + xlab("Body position") + theme(text = element_text(size=17), axis.title = element_text(size=14,face = "bold")) ggsave(paste0(i,".pdf"),dpi=600) }; rm("i","dt","Label","Lab","Norm","test") # Example of the T-test power calculation (here for mean heart rate - exploratory) power.t.test(n = 100, delta = mean(data$SupHR-data$StaHR), sd = sd(data$SupHR-data$StaHR), type = "paired", alternative = "two.sided") # Example of the Wilcoxon rank paired test power calculation (here for RMSSD, using rnorm - other distributions can also be inlcuded) pval <- replicate(1000, wilcox.test(rnorm(100,mean = mean(data$SupRMSSD),sd = sd(data$SupRMSSD)), rnorm(100,mean = mean(data$StaRMSSD),sd = sd(data$StaRMSSD)), paired = TRUE)$p.value); sum(pval < .05) / 1000 # Bayesian Correlations (Table 2 in the manuscript) # for supine body position dt <- data[,1:10] mat = matrix(rep(0,100),ncol=10) for (i in 1:10){ for (j in 1:10){ if(i90){ mat[i,j] <- (bcort$values$effects$`dt[, j]`$median)*sd(dt[,j])/sd(dt[,i]) } } } }; rm("i","j","bcort") mat <- data.frame(mat); colnames(mat) <- gsub("Sup","",colnames(dt)); rownames(mat) <- gsub("Sup","",colnames(dt)); View(mat) # for standing body position dt <- data[,11:20] mat = matrix(rep(0,100),ncol=10) for (i in 1:10){ for (j in 1:10){ if(i90){ mat[i,j] <- (bcort$values$effects$`dt[, j]`$median)*sd(dt[,j])/sd(dt[,i]) } } } }; rm("i","j","bcort") mat <- data.frame(mat); colnames(mat) <- gsub("Sta","",colnames(dt)); rownames(mat) <- gsub("Sta","",colnames(dt)); View(mat) rm("mat","dt") ## MAIN ANALYSIS library(powerMediation) library(bnlearn) # First method - GC generalCorr::allPairs(as.matrix(data[,c(1:10)]),typ = 3) generalCorr::allPairs(as.matrix(data[,c(11:20)]),typ = 3) # Second method - CAM with pruning test <- CAM::CAM(as.matrix(data[,c(1:10)]),pruning = TRUE, pruneMethodPars = list(cutOffPVal = 0.05)) test$Adj test <- CAM::CAM(as.matrix(data[,c(11:20)]),pruning = TRUE, pruneMethodPars = list(cutOffPVal = 0.05)) test$Adj rm("test") # Here, the R session should be restarted in order not to make a conflict with psycho and rstanarm packages .rs.restartR() # Third method - FGES FGES <- rcausal::fges(df = data[,c(1:10)], penaltydiscount = 1, maxDegree = -1, faithfulnessAssumed = TRUE, verbose = FALSE) FGES$edges FGES <- rcausal::fges(df = data[,c(11:20)], penaltydiscount = 1, maxDegree = -1, faithfulnessAssumed = TRUE, verbose = FALSE) FGES$edges rm("FGES") # Fourth method - GFCI GFCI <- rcausal::gfci(df = data[,c(1:10)], penaltydiscount = 1, maxDegree = -1, faithfulnessAssumed = TRUE, verbose = FALSE) GFCI$edges GFCI <- rcausal::gfci(df = data[,c(11:20)], penaltydiscount = 1, maxDegree = -1, faithfulnessAssumed = TRUE, verbose = FALSE) GFCI$edges rm("GFCI") # Here, the R session should be restarted in order not to make a conflict with rcausal package .rs.restartR() # Fifth method - HC res <- bnlearn::hc(data[,1:10]) plot(res) res <- bnlearn::hc(data[,11:20]) plot(res) # Sixth method - Tabu res <- bnlearn::tabu(data[,1:10]) plot(res) res <- bnlearn::tabu(data[,11:20]) plot(res) rm("res") # Here, the R session should be restarted in order not to make a conflict with bnlearn package .rs.restartR() # Mediations results # (1) RMSSD $\rightarrow$ HR $\rightarrow$ cInsT (supine), -> p = 0.073185 # (2) HR $\rightarrow$ cInsT $\rightarrow$ cExpT (supine), -> p = 0.08592151 # (3) HR $\rightarrow$ cInsT $\rightarrow$ cInsV (supine), -> p = 0.08800903 # (4) cInsT $\rightarrow$ ciRR $\rightarrow$ HR (standing), -> p = 0.1046831 # (5) cInsV $\rightarrow$ ciRR $\rightarrow$ HR (standing), -> p = 0.05838449 # Example of the calculation for the first relationship pred <- "SupRMSSD" medi <- "SupHR" post <- "SupRR" resmed <- medmod::med(data,post,medi,pred, paths=TRUE, estPlot = FALSE) resmedDF <- data.frame(resmed$paths) resmed$med powerMediation::testMediation.Sobel(resmedDF$est[1],resmedDF$est[2],resmedDF$se[1],resmedDF$se[2],alpha = 0.05)$pval