#
# Summarize data (Table 1)
#
data_diarrhea <- all_data[all_data$diarrhea=='Yes',]
data_not_diarrhea <- all_data[all_data$diarrhea=='No',]

# n
n_diarrhea <- nrow(data_diarrhea)
n_not_diarrhea <- nrow(data_not_diarrhea)
n_diarrhea
n_not_diarrhea

# observational time
summary(data_diarrhea[,'time'])
summary(data_not_diarrhea[,'time'])

# age
mean(data_diarrhea[,'age'])
mean(data_not_diarrhea[,'age'])
sd(data_diarrhea[,'age'])
sd(data_not_diarrhea[,'age'])
t.test(data_diarrhea[,'age'], data_not_diarrhea[,'age'])

# sex
n_male_diarrhea <- nrow(data_diarrhea[data_diarrhea$sex=='male',])
n_male_not_diarrhea <- nrow(data_not_diarrhea[data_not_diarrhea$sex=='male',])
n_male_diarrhea
n_male_not_diarrhea
prop.test(c(n_male_diarrhea, n_male_not_diarrhea), c(n_diarrhea, n_not_diarrhea), correct=F)

# rankin scale
mean(data_diarrhea[,'rankin_scale'])
mean(data_not_diarrhea[,'rankin_scale'])
sd(data_diarrhea[,'rankin_scale'])
sd(data_not_diarrhea[,'rankin_scale'])
t.test(data_diarrhea[,'rankin_scale'], data_not_diarrhea[,'rankin_scale'])

# abx
n_abx_diarrhea <- nrow(data_diarrhea[data_diarrhea$abx>0,])
n_abx_diarrhea
n_abx_not_diarrhea <- nrow(data_not_diarrhea[data_not_diarrhea$abx>0,])
n_abx_not_diarrhea
prop.test(c(n_abx_diarrhea, n_abx_not_diarrhea), c(n_diarrhea, n_not_diarrhea), correct=F)
t.test(data_diarrhea[,'abx'], data_not_diarrhea[,'abx'])

summary(data_diarrhea[,'abx'])
summary(data_not_diarrhea[,'abx'])

# tube_feeding
n_tube_feeding_diarrhea <- nrow(data_diarrhea[data_diarrhea$tube_feeding>0,])
n_tube_feeding_diarrhea
n_tube_feeding_not_diarrhea <- nrow(data_not_diarrhea[data_not_diarrhea$tube_feeding>0,])
n_tube_feeding_not_diarrhea
prop.test(c(n_tube_feeding_diarrhea, n_tube_feeding_not_diarrhea), c(n_diarrhea, n_not_diarrhea), correct=F)
t.test(data_diarrhea[,'tube_feeding'], data_not_diarrhea[,'tube_feeding'])

summary(data_diarrhea[,'tube_feeding'])
summary(data_not_diarrhea[,'tube_feeding'])

# ppi
n_ppi_diarrhea <- nrow(data_diarrhea[data_diarrhea$ppi>0,])
n_ppi_diarrhea
n_ppi_not_diarrhea <- nrow(data_not_diarrhea[data_not_diarrhea$ppi>0,])
n_ppi_not_diarrhea
prop.test(c(n_ppi_diarrhea, n_ppi_not_diarrhea), c(n_diarrhea, n_not_diarrhea), correct=F)
t.test(data_diarrhea[,'ppi'], data_not_diarrhea[,'ppi'])

summary(data_diarrhea[,'ppi'])
summary(data_not_diarrhea[,'ppi'])

# h2ra
n_h2ra_diarrhea <- nrow(data_diarrhea[data_diarrhea$h2_blocker>0,])
n_h2ra_diarrhea
n_h2ra_not_diarrhea <- nrow(data_not_diarrhea[data_not_diarrhea$h2_blocker>0,])
n_h2ra_not_diarrhea
prop.test(c(n_h2ra_diarrhea, n_h2ra_not_diarrhea), c(n_diarrhea, n_not_diarrhea), correct=F)
t.test(data_diarrhea[,'h2_blocker'], data_not_diarrhea[,'h2_blocker'])

summary(data_diarrhea[,'h2_blocker'])
summary(data_not_diarrhea[,'h2_blocker'])

#
# Check VIF
#
library(fmsb)
sex_num <- ifelse(all_data$sex == "male", 1, 0)
VIF(lm(age~sex_num+rankin_scale+abx+tube_feeding+ppi+h2_blocker, data=all_data))
VIF(lm(sex_num~age+rankin_scale+abx+tube_feeding+ppi+h2_blocker, data=all_data))
VIF(lm(rankin_scale~age+sex_num+abx+tube_feeding+ppi+h2_blocker, data=all_data))
VIF(lm(abx~age+sex_num+rankin_scale+tube_feeding+ppi+h2_blocker, data=all_data))
VIF(lm(tube_feeding~age+sex_num+rankin_scale+abx+ppi+h2_blocker, data=all_data))
VIF(lm(ppi~age+sex_num+rankin_scale+abx+tube_feeding+h2_blocker, data=all_data))
VIF(lm(h2_blocker~age+sex_num+rankin_scale+abx+tube_feeding+ppi, data=all_data))

#
# Univariate analysis within training dataset
#
data_diarrhea_in_train <- train_data[train_data$diarrhea=='Yes',]
data_not_diarrhea_in_train <- train_data[train_data$diarrhea=='No',]

t.test(data_diarrhea_in_train[,'age'], data_not_diarrhea_in_train[,'age'])
prop.test(c(nrow(data_diarrhea_in_train[data_diarrhea_in_train$sex=='male',]),
            nrow(data_not_diarrhea_in_train[data_not_diarrhea_in_train$sex=='male',])),
          c(nrow(data_diarrhea_in_train),
            nrow(data_not_diarrhea_in_train)),
          correct=F)
t.test(data_diarrhea_in_train[,'rankin_scale'], data_not_diarrhea_in_train[,'rankin_scale'])
t.test(data_diarrhea_in_train[,'abx'], data_not_diarrhea_in_train[,'abx'])
t.test(data_diarrhea_in_train[,'tube_feeding'], data_not_diarrhea_in_train[,'tube_feeding'])
t.test(data_diarrhea_in_train[,'ppi'], data_not_diarrhea_in_train[,'ppi'])
t.test(data_diarrhea_in_train[,'h2_blocker'], data_not_diarrhea_in_train[,'h2_blocker'])