## ------------------- SUPPLEMENTARY MATERIAL ------------------- ##
# Cardiorespiratory temporal causal links by sport (or lack thereof)

# R code - TiMINo, exploratory, RFE, FFTrees 
# Marcel Młyńczak, 2018

# --------------------------------------
# REFERENCES
# [1] Peters, J., Janzing, D., and Scholkopf, B. (2013). Causal inference on time series using 
#     restricted structural equation models. In Advances in Neural Information Processing Systems. 154–162

# --------------------------------------
# PREPARATION - Loading functions and libraries
setwd('/Users/marcel/Documents/Nauka/Publikacje/PostDoc/HAB_3')
library(gam)
library(kernlab)
library(ggplot2)
library(dplyr)
source("timino_causality.R")
source("./util/hammingDistance.R")
source("./util/indtestAll.R")
source("./util/indtestHsic.R")
source("./util/indtestPcor.R")
source("./util/TSindtest.R")
source("./util/fitting_ts.R")

# --------------------------------------
# 4. TiMINo analysis - the code [1]

# mat3 - the name of used variable with four inputs (as described below) 
# In order to deal with breathing phases, very small random numbers were added 
# mat3$Phase <- mat3$Phase + rnorm(n=nrow(mat3),mean=0,sd=0.01)
timino_dag(mat3, alpha = 0.05, max_lag = 4, model = traints_gam, indtest = indtestts_crosscov, output = TRUE, confounder_check = 1, instant = 1)


# --------------------------------------
# G-CAUSALITY DATA LOADING

load("Supp_Data.RData")

# Legend:
# Sport - type of sport based on the Mitchell's classification; control group = 0
# TypeDyn - the level of dynamic component during performance
# TypeStat - the level of static component during performance

# The names of remaining columns, e.g. SuM2-23:
# Su - supine; it can be also "St" for standing
# M2 - method 2 (spectral G-causality); it can be also "M1" for pairwise GC, and "M3" for extended GC assuming zero-lags
# 23 - the causality is assessed for the connection: from "2" to "3"

# For the first set of data (signals sampled at 25Hz):
# 1. TV; 2. iRR; 3. iHR (as written in the paper)

# For the second set of data (signals sampled at 2.5Hz):
# 1. TV; 2. iRR; 3. iHR (as written in the paper)

# For the third set of data (beat-by-beat sequences):
# 1. Lengths of consecutive RR intervals; 2. TV-related impedance amplitude at R peaks;
#  3. iRR at R peaks; 4. breathing phase at R peaks (as written in the paper)

# Table 1 preparation
table1 <- data %>% group_by(Sport) %>% 
  summarize(nF = sum(Sex=="Female"), nM = sum(Sex=="Male"),
            minHT = min(HT), meanHT = mean(HT), sdHT = sd(HT), maxHT = max(HT),
            minBM = min(BM), meanBM = mean(BM), sdBM = sd(BM), maxBM = max(BM)) 
xtable::xtable(table1)


# --------------------------------------
# EXPLORATORY ANALYSIS

# FIRST METHOD

# Summary for supine
summary(data[,grep("SuM1-",colnames(data))])

# Summary for standing
summary(data[,grep("StM1-",colnames(data))])

# Statistical tests mentioned in the text - supine
wilcox.test(dplyr::filter(data,data$GR=="Control")$`SuM1-31`,
       dplyr::filter(data,data$GR=="Sport")$`SuM1-31`,paired=FALSE,na.rm=TRUE)
kruskal.test(dplyr::filter(data,GR!="Control")$`SuM1-31`,as.factor(dplyr::filter(data,GR!="Control")$TypeDyn))

# Statistical tests mentioned in the text - standing
wilcox.test(dplyr::filter(data,data$GR=="Control")$`StM1-31`,
       dplyr::filter(data,data$GR=="Sport")$`StM1-31`,paired=FALSE,na.rm=TRUE)
kruskal.test(dplyr::filter(data,GR!="Control")$`StM1-31`,as.factor(dplyr::filter(data,GR!="Control")$TypeStat))

# Figure 1 - supine
ggplot(data, aes(x = TypeDyn, y = exp(`SuM1-31`)*100-100)) + geom_boxplot(fill="grey96") + theme_bw() + 
  xlab("") + ylab("Prediction Improvement [%]") +
  theme(text = element_text(size=18), axis.title = element_text(size=15, face="bold"))

# Figure 2 - standing
ggplot(data, aes(x = TypeStat, y = exp(`StM1-31`)*100-100)) + geom_boxplot(fill="grey96") + theme_bw() + 
  xlab("") + ylab("Prediction Improvement [%]") +
  theme(text = element_text(size=18), axis.title = element_text(size=15, face="bold"))


# SECOND METHOD

# Summary for supine
summary(data[,grep("SuM2-",colnames(data))])

# Summary for standing
summary(data[,grep("StM2-",colnames(data))])


# THIRD METHOD

# Summary for supine
summary(data[,grep("SuM3-",colnames(data))])

# Summary for standing
summary(data[,grep("StM3-",colnames(data))])

# Statistical tests mentioned in the text - supine
wilcox.test(dplyr::filter(data,data$GR=="Control")$`SuM3-12`,
       dplyr::filter(data,data$GR=="Sport")$`SuM3-12`,paired=FALSE,na.rm=TRUE)
kruskal.test(dplyr::filter(data,GR!="Control")$`SuM3-12`,as.factor(dplyr::filter(data,GR!="Control")$TypeStat))

# Statistical tests mentioned in the text - standing
wilcox.test(dplyr::filter(data,data$GR=="Control")$`StM3-12`,
       dplyr::filter(data,data$GR=="Sport")$`StM3-12`,paired=FALSE,na.rm=TRUE)
kruskal.test(dplyr::filter(data,GR!="Control")$`StM3-12`,as.factor(dplyr::filter(data,GR!="Control")$TypeStat))
pairwise.wilcox.test(dplyr::filter(data,GR!="Control")$`StM3-12`,as.factor(dplyr::filter(data,GR!="Control")$TypeStat))

# Figure 3 - supine
ggplot(data, aes(x = TypeStat, y = exp(`SuM3-12`)*100-100)) + geom_boxplot(fill="grey96") + theme_bw() + 
  xlab("") + ylab("Prediction Improvement [%]") +
  theme(text = element_text(size=18), axis.title = element_text(size=15, face="bold"))

# Figure 4 - standing
ggplot(data, aes(x = TypeStat, y = exp(`StM3-12`)*100-100)) + geom_boxplot(fill="grey96") + theme_bw() + 
  xlab("") + ylab("Prediction Improvement [%]") +
  theme(text = element_text(size=18), axis.title = element_text(size=15, face="bold"))


# --------------------------------------
# RECURSIVE FEATURE ELIMINATION (with Random Forests)

dt <- data[,c(8:25,27:30,32,34:35,37:38,40,42:59,61:64,66,68,69,71:72,74)]; dt[is.na(dt)] <- 0

# Between groups
res <- caret::rfe(dt, as.factor(data$GR), sizes=1:56, rfeControl=caret::rfeControl(functions=rfFuncs, method="cv", number=10))
pred_gr <- predictors(res)

# Between moderate and high dynamic components
ind <- which(data$TypeDyn=="Moderate" | data$TypeDyn=="High")
res <- caret::rfe(dt[ind,],as.factor(as.character(data$TypeDyn[ind])), sizes=1:56, rfeControl=caret::rfeControl(functions=rfFuncs, method="cv", number=10))
pred_dyn <- predictors(res)

# Between various levels of static components
ind <- which(data$GR!="Control")
res <- caret::rfe(dt[ind,],as.factor(as.character(data$TypeStat[ind])), sizes=1:56, rfeControl=caret::rfeControl(functions=rfFuncs, method="cv", number=10))
pred_stat <- predictors(res)

rm(ind,res)


# -------
# Other elements, e.g., preparing tables manually etc.

# Differences in groups, sexes, different static and dynamic components; example for approach 1, during supine
for (i in grep("SuM1-",colnames(data))){
  par(mfrow=c(2,2))
  boxplot(data[,i]~data$GR, main = colnames(data)[i], xlab = "GR", outline = FALSE)
  boxplot(data[,i]~data$Sex, main = colnames(data)[i], xlab = "Sex", outline = FALSE)
  boxplot(data[,i]~data$TypeDyn, main = colnames(data)[i], xlab = "TypeDyn", outline = FALSE)
  boxplot(data[,i]~data$TypeStat, main = colnames(data)[i], xlab = "TypeStat", outline = FALSE)
  invisible(readline(prompt="Press [enter] to continue"))
}; rm(i)

# Table 2 preparation
summ <- function(col_no){
  list(mean(100*exp(data[,col_no])-100,na.rm=TRUE),sd(100*exp(data[,col_no])-100,na.rm=TRUE),
       median(100*exp(data[,col_no])-100,na.rm=TRUE),IQR(100*exp(data[,col_no])-100,na.rm=TRUE),
       sum(is.na(data[,col_no])))}  

# Calculating all lines, without vectorization
round(unlist(summ(8)),2)
round(unlist(summ(9)),2)
round(unlist(summ(10)),2)
round(unlist(summ(11)),2)
round(unlist(summ(12)),2)
round(unlist(summ(13)),2)
round(unlist(summ(42)),2)
round(unlist(summ(43)),2)
round(unlist(summ(44)),2)
round(unlist(summ(45)),2)
round(unlist(summ(46)),2)
round(unlist(summ(47)),2)

# Table 3 preparation
summ2 <- function(col_no){
  list(median(data[,col_no],na.rm=TRUE),IQR(data[,col_no],na.rm=TRUE),
       median(data[,col_no+1],na.rm=TRUE),IQR(data[,col_no+1],na.rm=TRUE),
       sum(is.na(data[,col_no])))}  

# Calculating all lines, without vectorization
round(unlist(summ2(14)),2)
round(unlist(summ2(16)),2)
round(unlist(summ2(18)),2)
round(unlist(summ2(20)),2)
round(unlist(summ2(22)),2)
round(unlist(summ2(24)),2)
round(unlist(summ2(48)),2)
round(unlist(summ2(50)),2)
round(unlist(summ2(52)),2)
round(unlist(summ2(54)),2)
round(unlist(summ2(56)),2)
round(unlist(summ2(58)),2)

# Preparing Table 4 - Calculating all lines, without vectorization
round(unlist(summ(30)),2)
round(unlist(summ(34)),2)
round(unlist(summ(38)),2)
round(unlist(summ(35)),2)
round(unlist(summ(39)),2)
round(unlist(summ(40)),2)
round(unlist(summ(27)),2)
round(unlist(summ(28)),2)
round(unlist(summ(29)),2)
round(unlist(summ(32)),2)
round(unlist(summ(33)),2)
round(unlist(summ(37)),2)
round(unlist(summ(64)),2)
round(unlist(summ(68)),2)
round(unlist(summ(72)),2)
round(unlist(summ(69)),2)
round(unlist(summ(73)),2)
round(unlist(summ(74)),2)
round(unlist(summ(61)),2)
round(unlist(summ(62)),2)
round(unlist(summ(63)),2)
round(unlist(summ(66)),2)
round(unlist(summ(67)),2)
round(unlist(summ(71)),2)