library(readxl)
zs <- read_xlsx("D:/2.RESEARCH/0. Projects/Zinc_sepsis/Original_data_and_analysis_results/Zinc_Sepsis_THAU.xlsx")
attach(zs)
zs$Zinc31_pos <- ifelse(zs$Zinc31 > 0,1,0)# define new variable
zs$Zinc31_pos
##   [1] 0 0 0 0 1 1 1 1 1 0 1 1 1 1 0 1 1 1 0 0 1 1 1 1 1 0 1 1 1 1 0 1 0 1 0 0 0
##  [38] 1 1 1 1 0 0 1 1 1 1 0 0 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 0 0
##  [75] 0 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 0 0 0 0 1 1 1 1 1
library(table1)
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
table1(~as.factor(zs$Zinc31_pos)|Mortality.7day, data = zs,render.categorical="FREQ (PCTnoNA%)")
## Warning in table1.formula(~as.factor(zs$Zinc31_pos) | Mortality.7day, data =
## zs, : Terms to the right of '|' in formula 'x' define table columns and are
## expected to be factors with meaningful labels.
0
(N=104)		 1
(N=16)		 Overall
(N=120)
as.factor(zs$Zinc31_pos)
0 24 (23.1%) 8 (50.0%) 32 (26.7%)
1 80 (76.9%) 8 (50.0%) 88 (73.3%) 
tab = table(zs$Zinc31_pos, zs$Mortality.7day)
chisq.test(tab)
## Warning in chisq.test(tab): Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab
## X-squared = 3.8553, df = 1, p-value = 0.04959
table1(~Zinc1+ Zinc3|Mortality.7day, data = zs,render.continuous=c(.="Mean(SD)", .="Median[Q1, Q3]"),digits=4)
## Warning in table1.formula(~Zinc1 + Zinc3 | Mortality.7day, data = zs,
## render.continuous = c(. = "Mean(SD)", : Terms to the right of '|' in formula
## 'x' define table columns and are expected to be factors with meaningful labels.
0
(N= 104)		 1
(N= 16)		 Overall
(N= 120)
Zinc1
Mean(SD) 59.19(27.30) 68.87(38.15) 60.48(28.97)
Median[Q1, Q3] 54.00[37.53, 71.00] 55.85[44.43, 94.20] 54.00[38.28, 73.60]
Zinc3
Mean(SD) 73.02(28.28) 61.61(26.14) 71.50(28.17)
Median[Q1, Q3] 68.65[54.53, 85.95] 56.05[46.38, 68.45] 66.05[52.68, 82.70] 
library(compareGroups)
createTable(compareGroups(Mortality.7day ~ Zinc31, data = zs, method = 1))
## 
## --------Summary descriptives table by 'Mortality.7day'---------
## 
## _________________________________________ 
##             0           1       p.overall 
##           N=104        N=16               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## Zinc31 13.8 (23.6) -7.26 (35.2)   0.033   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
table1(~as.factor(zs$Zinc31over7.65)|Mortality.7day, data = zs,render.categorical="FREQ (PCTnoNA%)")
## Warning in table1.formula(~as.factor(zs$Zinc31over7.65) | Mortality.7day, :
## Terms to the right of '|' in formula 'x' define table columns and are expected
## to be factors with meaningful labels.
0
(N=104)		 1
(N=16)		 Overall
(N=120)
as.factor(zs$Zinc31over7.65)
0 35 (33.7%) 11 (68.8%) 46 (38.3%)
1 69 (66.3%) 5 (31.3%) 74 (61.7%) 
tab = table(zs$Zinc31over7.65, zs$Mortality.7day)
chisq.test(tab)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab
## X-squared = 5.817, df = 1, p-value = 0.01587

#ROC curve

library(pROC)
## Warning: package 'pROC' was built under R version 4.4.3
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
roc_result <- roc(Mortality.7day ~ Zinc31, data = zs)
## Setting levels: control = 0, case = 1
## Setting direction: controls > cases
plot(roc_result, print.auc = TRUE, col = "blue", main = "ROC Curve for zinc31")

roc_result <- roc(Mortality.7day ~ Zinc31, data = zs)
## Setting levels: control = 0, case = 1
## Setting direction: controls > cases
plot(roc_result, print.auc = TRUE, col = "blue", main = "ROC Curve for Zinc31")

coords_result <- coords(roc_result, "best", ret = c("threshold", "sensitivity", "specificity"))
print(coords_result)
##   threshold sensitivity specificity
## 1      7.65      0.6875   0.6634615
roc_zinc  <- roc(zs$Mortality.7day, zs$Zinc31, ci = TRUE)
## Setting levels: control = 0, case = 1
## Setting direction: controls > cases
roc_sofa <- roc(zs$Mortality.7day, zs$SOFA, ci = TRUE)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
model_combined <- glm(Mortality.7day ~ SOFA + Zinc31, data = zs, family = binomial)
zs$pred_combined <- predict(model_combined, type = "response")
roc_combined <- roc(zs$Mortality.7day, zs$pred_combined, ci = TRUE)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc.test(roc_zinc, roc_sofa)
## 
##  Bootstrap test for two correlated ROC curves
## 
## data:  roc_zinc and roc_sofa
## D = -0.94884, boot.n = 2000, boot.stratified = 1, p-value = 0.3427
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2 
##   0.6859976   0.7785457
roc.test(roc_sofa, roc_combined)
## 
##  DeLong's test for two correlated ROC curves
## 
## data:  roc_sofa and roc_combined
## Z = -0.87685, p-value = 0.3806
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
##  -0.10985018  0.04194153
## sample estimates:
## AUC of roc1 AUC of roc2 
##   0.7785457   0.8125000
roc.test(roc_sofa, roc_combined)
## 
##  DeLong's test for two correlated ROC curves
## 
## data:  roc_sofa and roc_combined
## Z = -0.87685, p-value = 0.3806
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
##  -0.10985018  0.04194153
## sample estimates:
## AUC of roc1 AUC of roc2 
##   0.7785457   0.8125000
plot(roc_zinc, col = "blue", main = "ROC Curve", legacy.axes = TRUE)
plot(roc_sofa, add = TRUE, col = "red")
plot(roc_combined, add = TRUE, col = "darkgreen")

legend("bottomright",
       legend = c(
         paste("∆Zinc (AUC =", round(auc(roc_zinc), 2), ")"),
         paste("SOFA (AUC =", round(auc(roc_sofa), 2), ")"),
         paste("SOFA + ∆Zinc (AUC =", round(auc(roc_combined), 2), ")")
       ),
       col = c("blue", "red", "darkgreen"), lwd = 2)

# Để xem trong Plots panel, paste xuống Console