Zinc_Sepsis_with Shock and MV
Thau NS
2025-07-28
library(readxl)
zs <- read_xlsx("D:/2.RESEARCH/0. Projects/Zinc_sepsis/Original_data_and_analysis_results/Zinc_Sepsis_THAU.xlsx")
attach(zs)#Display variables with grouped by 7-day mortality
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-For continous variables
table1(~SOFA+Lactat+Age+WBC+PCT+Zinc31|Mortality.7day, data = zs,render.continuous=c(.="Mean(SD)", .="Median[Q1, Q3]"),digits=4)## Warning in table1.formula(~SOFA + Lactat + Age + WBC + PCT + Zinc31 |
## 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) |
|
|---|---|---|---|
| SOFA | |||
| Mean(SD) | 6.788(3.511) | 10.81(3.902) | 7.325(3.804) |
| Median[Q1, Q3] | 6.000[4.000, 9.000] | 10.00[8.000, 13.50] | 7.000[4.000, 10.00] |
| Lactat | |||
| Mean(SD) | 3.996(3.001) | 5.700(3.605) | 4.246(3.137) |
| Median[Q1, Q3] | 2.900[2.000, 5.000] | 4.600[2.600, 7.475] | 3.100[2.000, 5.400] |
| Missing | 11 (10.6%) | 0 (0%) | 11 (9.2%) |
| Age | |||
| Mean(SD) | 64.75(13.48) | 59.31(12.78) | 64.03(13.47) |
| Median[Q1, Q3] | 65.00[57.00, 74.00] | 57.50[54.75, 63.75] | 64.50[56.75, 73.25] |
| WBC | |||
| Mean(SD) | 16.38(8.615) | 17.81(9.786) | 16.57(8.749) |
| Median[Q1, Q3] | 15.77[10.44, 20.82] | 16.00[11.72, 22.96] | 15.77[10.53, 21.66] |
| PCT | |||
| Mean(SD) | 66.74(87.88) | 56.86(70.81) | 65.37(85.50) |
| Median[Q1, Q3] | 29.74[6.325, 100.0] | 26.42[6.058, 77.04] | 29.74[5.805, 100.0] |
| Missing | 5 (4.8%) | 0 (0%) | 5 (4.2%) |
| Zinc31 | |||
| Mean(SD) | 13.83(23.56) | -7.256(35.23) | 11.02(26.24) |
| Median[Q1, Q3] | 15.40[2.600, 30.10] | 0.8000[-29.00, 14.65] | 13.75[-3.175, 29.85] |
##For catagorial variables
table1(~as.factor(Male)+as.factor(CKD)+as.factor(COPD)+as.factor(CHD)+as.factor(CLD)+as.factor(Diabetes)+as.factor(Shock_admit)+as.factor(MV_admit)+as.factor(Male)+as.factor(Pathogen_GramPos)+as.factor(Source_of_Infection)|Mortality.7day, data = zs,render.categorical="FREQ (PCTnoNA%)") #decimal did not work for catagorial variables; digits =1 was fixed for percentage; trying to adjust digits for percentage but I failed.## Warning in table1.formula(~as.factor(Male) + as.factor(CKD) + as.factor(COPD) +
## : 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(Male) | |||
| 0 | 33 (31.7%) | 7 (43.8%) | 40 (33.3%) |
| 1 | 71 (68.3%) | 9 (56.3%) | 80 (66.7%) |
| as.factor(CKD) | |||
| 0 | 101 (97.1%) | 15 (93.8%) | 116 (96.7%) |
| 1 | 3 (2.9%) | 1 (6.3%) | 4 (3.3%) |
| as.factor(COPD) | |||
| 0 | 96 (92.3%) | 14 (87.5%) | 110 (91.7%) |
| 1 | 8 (7.7%) | 2 (12.5%) | 10 (8.3%) |
| as.factor(CHD) | |||
| 0 | 56 (53.8%) | 9 (56.3%) | 65 (54.2%) |
| 1 | 48 (46.2%) | 7 (43.8%) | 55 (45.8%) |
| as.factor(CLD) | |||
| 0 | 88 (84.6%) | 10 (62.5%) | 98 (81.7%) |
| 1 | 16 (15.4%) | 6 (37.5%) | 22 (18.3%) |
| as.factor(Diabetes) | |||
| 0 | 62 (59.6%) | 9 (56.3%) | 71 (59.2%) |
| 1 | 42 (40.4%) | 7 (43.8%) | 49 (40.8%) |
| as.factor(Shock_admit) | |||
| 0 | 62 (59.6%) | 6 (37.5%) | 68 (56.7%) |
| 1 | 42 (40.4%) | 10 (62.5%) | 52 (43.3%) |
| as.factor(MV_admit) | |||
| 0 | 90 (86.5%) | 10 (62.5%) | 100 (83.3%) |
| 1 | 14 (13.5%) | 6 (37.5%) | 20 (16.7%) |
| as.factor(Pathogen_GramPos) | |||
| 0 | 39 (37.5%) | 6 (37.5%) | 45 (37.5%) |
| 1 | 13 (12.5%) | 3 (18.8%) | 16 (13.3%) |
| 2 | 52 (50.0%) | 7 (43.8%) | 59 (49.2%) |
| as.factor(Source_of_Infection) | |||
| Abdominal | 29 (27.9%) | 6 (37.5%) | 35 (29.2%) |
| Bone | 1 (1.0%) | 0 (0%) | 1 (0.8%) |
| CNS | 8 (7.7%) | 3 (18.8%) | 11 (9.2%) |
| Pneumoniae | 20 (19.2%) | 3 (18.8%) | 23 (19.2%) |
| Sinusitis | 1 (1.0%) | 0 (0%) | 1 (0.8%) |
| SSTI | 14 (13.5%) | 4 (25.0%) | 18 (15.0%) |
| Unknown | 5 (4.8%) | 0 (0%) | 5 (4.2%) |
| UTI | 26 (25.0%) | 0 (0%) | 26 (21.7%) |
#P-value compared two groups ##for continous variable
library(compareGroups) # for p-value for continous variablecreateTable(compareGroups(Mortality.7day ~ Age+WBC, data = zs, method = 1))#method =1 mean normal distribution##
## --------Summary descriptives table by 'Mortality.7day'---------
##
## _____________________________________
## 0 1 p.overall
## N=104 N=16
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## Age 64.8 (13.5) 59.3 (12.8) 0.131
## WBC 16.4 (8.62) 17.8 (9.79) 0.586
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯createTable(compareGroups(Mortality.7day ~PCT+Lactat+SOFA, data = zs, method = 1))#method =2 mean unnormal distribution##
## --------Summary descriptives table by 'Mortality.7day'---------
##
## ________________________________________
## 0 1 p.overall
## N=104 N=16
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## PCT 66.7 (87.9) 56.9 (70.8) 0.622
## Lactat 4.00 (3.00) 5.70 (3.61) 0.090
## SOFA 6.79 (3.51) 10.8 (3.90) 0.001
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯#for cat. variable ## for single variable
tab = table(zs$Shock_admit, zs$Mortality.7day)
tab##
## 0 1
## 0 62 6
## 1 42 10prop.table(tab, 2)##
## 0 1
## 0 0.5961538 0.3750000
## 1 0.4038462 0.6250000prop.table(tab, 1)##
## 0 1
## 0 0.91176471 0.08823529
## 1 0.80769231 0.19230769chisq.test(tab)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tab
## X-squared = 1.9347, df = 1, p-value = 0.1642fisher.test(tab)##
## Fisher's Exact Test for Count Data
##
## data: tab
## p-value = 0.111
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.7381492 8.8332947
## sample estimates:
## odds ratio
## 2.441511##Note: p-value: For many variables, but the name of test was not indicated
#If package was not installes
#install.packages("gtsummary")
#install.packages("dplyr") # Required for %>%library(gtsummary)## Warning: package 'gtsummary' was built under R version 4.4.3library(dplyr)## Warning: package 'dplyr' was built under R version 4.4.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionzs %>%
select(Mortality.7day,Age,Male, Stroke,Cancer, CHD ,CKD, CLD, COPD,Diabetes,Source_of_Infection,Pathogen_GramPos,Shock_admit,MV_admit,SOFA,Lactat,WBC,PCT) %>%
tbl_summary(by = Mortality.7day,
digits = all_continuous() ~ 2) %>%
add_p() | Characteristic | 0 N = 1041 |
1 N = 161 |
p-value2 |
|---|---|---|---|
| Age | 65.00 (57.00, 74.00) | 57.50 (54.50, 66.50) | 0.071 |
| Male | 71 (68%) | 9 (56%) | 0.3 |
| Stroke | 8 (7.7%) | 3 (19%) | 0.2 |
| Cancer | 2 (1.9%) | 0 (0%) | >0.9 |
| CHD | 48 (46%) | 7 (44%) | 0.9 |
| CKD | 3 (2.9%) | 1 (6.3%) | 0.4 |
| CLD | 16 (15%) | 6 (38%) | 0.075 |
| COPD | 8 (7.7%) | 2 (13%) | 0.6 |
| Diabetes | 42 (40%) | 7 (44%) | 0.8 |
| Source_of_Infection | 0.15 | ||
| Abdominal | 29 (28%) | 6 (38%) | |
| Bone | 1 (1.0%) | 0 (0%) | |
| CNS | 8 (7.7%) | 3 (19%) | |
| Pneumoniae | 20 (19%) | 3 (19%) | |
| Sinusitis | 1 (1.0%) | 0 (0%) | |
| SSTI | 14 (13%) | 4 (25%) | |
| Unknown | 5 (4.8%) | 0 (0%) | |
| UTI | 26 (25%) | 0 (0%) | |
| Pathogen_GramPos | 0.7 | ||
| 0 | 39 (38%) | 6 (38%) | |
| 1 | 13 (13%) | 3 (19%) | |
| 2 | 52 (50%) | 7 (44%) | |
| Shock_admit | 42 (40%) | 10 (63%) | 0.10 |
| MV_admit | 14 (13%) | 6 (38%) | 0.027 |
| SOFA | 6.00 (4.00, 9.00) | 10.00 (8.00, 14.00) | <0.001 |
| Lactat | 2.90 (2.00, 5.00) | 4.60 (2.60, 7.55) | 0.027 |
| Unknown | 11 | 0 | |
| WBC | 15.77 (10.43, 21.08) | 16.00 (11.71, 23.28) | 0.6 |
| PCT | 29.74 (4.95, 100.00) | 26.42 (5.81, 84.69) | 0.8 |
| Unknown | 5 | 0 | |
| 1 Median (Q1, Q3); n (%) | |||
| 2 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test | |||
#Phan tich don bien va da bien voi bien Shock_admit and MV_admit
library(epiDisplay)## Loading required package: foreign## Loading required package: survival## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## The following object is masked from 'package:gtsummary':
##
## select## Loading required package: nnet7-day mortality
logrit.uni.zinc31 <- glm(Mortality.7day~ Zinc31,family=binomial, data = zs)
logistic.display(logrit.uni.zinc31)##
## Logistic regression predicting Mortality.7day
##
## OR(95%CI) P(Wald's test) P(LR-test)
## Zinc31 (cont. var.) 0.97 (0.96,0.99) 0.006 0.005
##
## Log-likelihood = -43.1638
## No. of observations = 120
## AIC value = 90.3275logrit.uni.shock <- glm(Mortality.7day~ Shock_admit,family=binomial, data = zs)
logistic.display(logrit.uni.shock)##
## Logistic regression predicting Mortality.7day
##
## OR(95%CI) P(Wald's test) P(LR-test)
## Shock_admit: 1 vs 0 2.46 (0.83,7.28) 0.104 0.098
##
## Log-likelihood = -45.7503
## No. of observations = 120
## AIC value = 95.5007logrit.uni.SOFA <- glm(Mortality.7day~ SOFA,family=binomial, data = zs)
logistic.display(logrit.uni.SOFA)##
## Logistic regression predicting Mortality.7day
##
## OR(95%CI) P(Wald's test) P(LR-test)
## SOFA (cont. var.) 1.32 (1.13,1.55) < 0.001 < 0.001
##
## Log-likelihood = -39.5636
## No. of observations = 120
## AIC value = 83.1271logrit.multi <- glm(Mortality.7day~Age+Male+ SOFA+ Zinc31+Stroke+Cancer+ CHD +CKD+ CLD+ COPD+Pathogen01+Shock_admit+MV_admit,family=binomial, data = zs)
logistic.display(logrit.multi, decimal = 3)##
## Logistic regression predicting Mortality.7day
##
## crude OR(95%CI) adj. OR(95%CI)
## Age (cont. var.) 0.972 (0.936,1.009) 0.918 (0.859,0.981)
##
## Male: 1 vs 0 0.598 (0.205,1.743) 0.183 (0.029,1.175)
##
## SOFA (cont. var.) 1.325 (1.134,1.547) 1.988 (1.331,2.97)
##
## Zinc31 (cont. var.) 0.974 (0.956,0.992) 0.952 (0.92,0.985)
##
## Stroke: 1 vs 0 2.769 (0.651,11.78) 5.077 (0.516,49.959)
##
## Cancer: 1 vs 0 0 (0,Inf) 0 (0,Inf)
##
## CHD: 1 vs 0 0.907 (0.314,2.62) 4.565 (0.593,35.167)
##
## CKD: 1 vs 0 2.244 (0.219,23.005) 15.361 (0.657,359.193)
##
## CLD: 1 vs 0 3.3 (1.051,10.357) 6.696 (0.947,47.357)
##
## COPD: 1 vs 0 1.714 (0.33,8.907) 4.026 (0.42,38.6)
##
## Pathogen01: 1 vs 0 1.286 (0.445,3.711) 1.217 (0.257,5.762)
##
## Shock_admit: 1 vs 0 2.46 (0.831,7.283) 0.083 (0.007,0.914)
##
## MV_admit: 1 vs 0 3.857 (1.211,12.282) 0.137 (0.013,1.481)
##
## P(Wald's test) P(LR-test)
## Age (cont. var.) 0.0118 0.0064
##
## Male: 1 vs 0 0.0735 0.059
##
## SOFA (cont. var.) < 0.001 < 0.001
##
## Zinc31 (cont. var.) 0.0045 0.0011
##
## Stroke: 1 vs 0 0.1637 0.1733
##
## Cancer: 1 vs 0 0.9945 0.8821
##
## CHD: 1 vs 0 0.1449 0.1405
##
## CKD: 1 vs 0 0.0894 0.117
##
## CLD: 1 vs 0 0.0567 0.045
##
## COPD: 1 vs 0 0.2272 0.2421
##
## Pathogen01: 1 vs 0 0.8048 0.8046
##
## Shock_admit: 1 vs 0 0.042 0.0268
##
## MV_admit: 1 vs 0 0.1017 0.0753
##
## Log-likelihood = -26.09609
## No. of observations = 120
## AIC value = 80.1921714-day mortality
logrit.uni.zinc31 <- glm(Mortality.14day~ Zinc31,family=binomial, data = zs)
logistic.display(logrit.uni.zinc31)##
## Logistic regression predicting Mortality.14day
##
## OR(95%CI) P(Wald's test) P(LR-test)
## Zinc31 (cont. var.) 0.98 (0.96,1) 0.014 0.012
##
## Log-likelihood = -58.286
## No. of observations = 120
## AIC value = 120.572logrit.uni.shock <- glm(Mortality.14day~ Shock_admit,family=binomial, data = zs)
logistic.display(logrit.uni.shock)##
## Logistic regression predicting Mortality.14day
##
## OR(95%CI) P(Wald's test) P(LR-test)
## Shock_admit: 1 vs 0 1.27 (0.52,3.07) 0.597 0.598
##
## Log-likelihood = -61.2695
## No. of observations = 120
## AIC value = 126.5389logrit.multi <- glm(Mortality.14day~Age+Male+ SOFA+ Zinc31+Stroke+Cancer+ CHD +CKD+ CLD+ COPD+Pathogen01+Shock_admit+MV_admit,family=binomial, data = zs)
logistic.display(logrit.multi,decimal = 3)##
## Logistic regression predicting Mortality.14day
##
## crude OR(95%CI) adj. OR(95%CI)
## Age (cont. var.) 0.981 (0.95,1.013) 0.972 (0.929,1.016)
##
## Male: 1 vs 0 1.079 (0.421,2.768) 0.981 (0.29,3.318)
##
## SOFA (cont. var.) 1.24 (1.095,1.404) 1.382 (1.088,1.756)
##
## Zinc31 (cont. var.) 0.98 (0.964,0.996) 0.983 (0.961,1.005)
##
## Stroke: 1 vs 0 1.483 (0.363,6.055) 1.325 (0.22,7.981)
##
## Cancer: 1 vs 0 0 (0,Inf) 0 (0,Inf)
##
## CHD: 1 vs 0 0.741 (0.302,1.814) 1.488 (0.378,5.866)
##
## CKD: 1 vs 0 1.278 (0.127,12.838) 2.921 (0.222,38.508)
##
## CLD: 1 vs 0 2.074 (0.739,5.825) 1.434 (0.366,5.628)
##
## COPD: 1 vs 0 0.946 (0.188,4.76) 1.333 (0.205,8.666)
##
## Pathogen01: 1 vs 0 0.867 (0.359,2.093) 1.126 (0.369,3.441)
##
## Shock_admit: 1 vs 0 1.269 (0.524,3.072) 0.203 (0.037,1.103)
##
## MV_admit: 1 vs 0 10.038 (3.451,29.198) 2.908 (0.743,11.383)
##
## P(Wald's test) P(LR-test)
## Age (cont. var.) 0.2063 0.2018
##
## Male: 1 vs 0 0.9758 0.9758
##
## SOFA (cont. var.) 0.0081 0.0058
##
## Zinc31 (cont. var.) 0.1192 0.1148
##
## Stroke: 1 vs 0 0.7588 0.7612
##
## Cancer: 1 vs 0 0.9933 0.6525
##
## CHD: 1 vs 0 0.5698 0.5692
##
## CKD: 1 vs 0 0.4153 0.4439
##
## CLD: 1 vs 0 0.6051 0.6072
##
## COPD: 1 vs 0 0.7633 0.7668
##
## Pathogen01: 1 vs 0 0.8349 0.8347
##
## Shock_admit: 1 vs 0 0.0649 0.0566
##
## MV_admit: 1 vs 0 0.1252 0.1261
##
## Log-likelihood = -45.86274
## No. of observations = 120
## AIC value = 119.7254930 day mortality
logrit.uni.zinc31 <- glm(Mortality.30day~ Zinc31,family=binomial, data = zs)
logistic.display(logrit.uni.zinc31)##
## Logistic regression predicting Mortality.30day
##
## OR(95%CI) P(Wald's test) P(LR-test)
## Zinc31 (cont. var.) 0.98 (0.97,1) 0.021 0.019
##
## Log-likelihood = -63.5977
## No. of observations = 120
## AIC value = 131.1954logrit.uni.shock <- glm(Mortality.30day~ Shock_admit,family=binomial, data = zs)
logistic.display(logrit.uni.shock)##
## Logistic regression predicting Mortality.30day
##
## OR(95%CI) P(Wald's test) P(LR-test)
## Shock_admit: 1 vs 0 1.08 (0.47,2.51) 0.852 0.852
##
## Log-likelihood = -66.3419
## No. of observations = 120
## AIC value = 136.6837logrit.multi <- glm(Mortality.30day~Age+Male+ SOFA+ Zinc31+Stroke+Cancer+ CHD +CKD+ CLD+ COPD+Pathogen01+Shock_admit+MV_admit,family=binomial, data = zs)
logistic.display(logrit.multi,decimal = 3)##
## Logistic regression predicting Mortality.30day
##
## crude OR(95%CI) adj. OR(95%CI)
## Age (cont. var.) 0.994 (0.963,1.025) 0.984 (0.944,1.026)
##
## Male: 1 vs 0 1.148 (0.468,2.819) 1.202 (0.386,3.745)
##
## SOFA (cont. var.) 1.217 (1.082,1.368) 1.391 (1.111,1.741)
##
## Zinc31 (cont. var.) 0.982 (0.966,0.997) 0.985 (0.965,1.006)
##
## Stroke: 1 vs 0 1.92 (0.52,7.095) 1.57 (0.308,7.995)
##
## Cancer: 1 vs 0 0 (0,Inf) 0 (0,Inf)
##
## CHD: 1 vs 0 0.948 (0.409,2.196) 1.602 (0.45,5.709)
##
## CKD: 1 vs 0 1.048 (0.105,10.478) 1.795 (0.147,21.902)
##
## CLD: 1 vs 0 1.612 (0.584,4.448) 1.061 (0.283,3.984)
##
## COPD: 1 vs 0 0.7685 (0.1537,3.8416) 1.0012 (0.1585,6.3228)
##
## Pathogen01: 1 vs 0 1.048 (0.454,2.419) 1.473 (0.512,4.241)
##
## Shock_admit: 1 vs 0 1.083 (0.467,2.513) 0.171 (0.034,0.857)
##
## MV_admit: 1 vs 0 9.75 (3.368,28.228) 2.993 (0.814,11.002)
##
## P(Wald's test) P(LR-test)
## Age (cont. var.) 0.4618 0.4607
##
## Male: 1 vs 0 0.7512 0.7506
##
## SOFA (cont. var.) 0.004 0.0024
##
## Zinc31 (cont. var.) 0.1572 0.1531
##
## Stroke: 1 vs 0 0.5872 0.5908
##
## Cancer: 1 vs 0 0.9931 0.5937
##
## CHD: 1 vs 0 0.4672 0.465
##
## CKD: 1 vs 0 0.6467 0.658
##
## CLD: 1 vs 0 0.93 0.9301
##
## COPD: 1 vs 0 0.999 0.999
##
## Pathogen01: 1 vs 0 0.4727 0.4699
##
## Shock_admit: 1 vs 0 0.0317 0.0246
##
## MV_admit: 1 vs 0 0.0989 0.0978
##
## Log-likelihood = -50.74749
## No. of observations = 120
## AIC value = 129.4949790 day mortality
logrit.uni.zinc31 <- glm(Mortality.90day~ Zinc31,family=binomial, data = zs)
logistic.display(logrit.uni.zinc31)##
## Logistic regression predicting Mortality.90day
##
## OR(95%CI) P(Wald's test) P(LR-test)
## Zinc31 (cont. var.) 0.99 (0.97,1) 0.063 0.061
##
## Log-likelihood = -68.8229
## No. of observations = 120
## AIC value = 141.6458logrit.uni.shock <- glm(Mortality.90day~ Shock_admit,family=binomial, data = zs)
logistic.display(logrit.uni.shock)##
## Logistic regression predicting Mortality.90day
##
## OR(95%CI) P(Wald's test) P(LR-test)
## Shock_admit: 1 vs 0 1.13 (0.5,2.52) 0.773 0.773
##
## Log-likelihood = -70.5386
## No. of observations = 120
## AIC value = 145.0773logrit.multi <- glm(Mortality.90day~Age+Male+ SOFA+ Zinc31+Stroke+Cancer+ CHD +CKD+ CLD+ COPD+Pathogen01+Shock_admit+MV_admit,family=binomial, data = zs)
logistic.display(logrit.multi,decimal = 3)##
## Logistic regression predicting Mortality.90day
##
## crude OR(95%CI) adj. OR(95%CI)
## Age (cont. var.) 0.989 (0.96,1.019) 0.982 (0.944,1.022)
##
## Male: 1 vs 0 1.211 (0.51,2.872) 1.125 (0.389,3.255)
##
## SOFA (cont. var.) 1.179 (1.056,1.316) 1.293 (1.057,1.581)
##
## Zinc31 (cont. var.) 0.986 (0.971,1.001) 0.989 (0.97,1.009)
##
## Stroke: 1 vs 0 1.576 (0.43,5.783) 1.421 (0.289,6.988)
##
## Cancer: 1 vs 0 2.687 (0.163,44.255) 11.691 (0.482,283.771)
##
## CHD: 1 vs 0 0.979 (0.438,2.189) 1.521 (0.469,4.931)
##
## CKD: 1 vs 0 0.875 (0.088,8.723) 1.552 (0.131,18.459)
##
## CLD: 1 vs 0 2.135 (0.812,5.612) 1.625 (0.48,5.51)
##
## COPD: 1 vs 0 1.143 (0.277,4.71) 1.503 (0.291,7.752)
##
## Pathogen01: 1 vs 0 1.038 (0.466,2.315) 1.324 (0.498,3.521)
##
## Shock_admit: 1 vs 0 1.126 (0.503,2.522) 0.242 (0.057,1.037)
##
## MV_admit: 1 vs 0 7.429 (2.622,21.045) 2.997 (0.819,10.961)
##
## P(Wald's test) P(LR-test)
## Age (cont. var.) 0.3697 0.368
##
## Male: 1 vs 0 0.8285 0.8283
##
## SOFA (cont. var.) 0.0123 0.0102
##
## Zinc31 (cont. var.) 0.2857 0.2847
##
## Stroke: 1 vs 0 0.6659 0.6682
##
## Cancer: 1 vs 0 0.1308 0.1485
##
## CHD: 1 vs 0 0.4849 0.4832
##
## CKD: 1 vs 0 0.7277 0.7346
##
## CLD: 1 vs 0 0.4354 0.438
##
## COPD: 1 vs 0 0.6265 0.6318
##
## Pathogen01: 1 vs 0 0.5734 0.5721
##
## Shock_admit: 1 vs 0 0.056 0.0493
##
## MV_admit: 1 vs 0 0.0972 0.0945
##
## Log-likelihood = -57.53794
## No. of observations = 120
## AIC value = 143.07588check correlation
cor=zs[, c("Age","Male","Zinc31","SOFA","Shock_admit","MV_admit")]
#if not installed install.packages("pysch")
library(psych)## Warning: package 'psych' was built under R version 4.4.3##
## Attaching package: 'psych'## The following objects are masked from 'package:epiDisplay':
##
## alpha, cs, lookuppairs.panels(cor)
#Cut-off of Zinc31 for predict Mortality.7day
# if needed install.packages("pROC") # if not already installed
library(pROC)## Warning: package 'pROC' was built under R version 4.4.3## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following object is masked from 'package:epiDisplay':
##
## ci## The following objects are masked from 'package:stats':
##
## cov, smooth, var# Assuming:
# zinc31 is numeric
# Mortality.7day is binary (0 = alive, 1 = dead)
roc_result <- roc(Mortality.7day ~ Zinc31, data = zs)## Setting levels: control = 0, case = 1## Setting direction: controls > cases# Plot the ROC curve
plot(roc_result, print.auc = TRUE, col = "blue", main = "ROC Curve for zinc31")
# Find best cut-off (Youden's Index)
coords_result <- coords(roc_result, "best", ret = c("threshold", "sensitivity", "specificity"))
print(coords_result)## threshold sensitivity specificity
## 1 7.65 0.6875 0.6634615# Assuming:
# zinc31 is numeric
# Mortality.7day is binary (0 = alive, 1 = dead)
roc_result <- roc(Mortality.7day ~ Zinc31, data = zs)## Setting levels: control = 0, case = 1
## Setting direction: controls > cases# Plot the ROC curve
plot(roc_result, print.auc = TRUE, col = "blue", main = "ROC Curve for Zinc31")
# Find best cut-off (Youden's Index)
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 curves
# 1. ROC for zinc31
roc_zinc <- roc(zs$Mortality.7day, zs$Zinc31, ci = TRUE)## Setting levels: control = 0, case = 1## Setting direction: controls > cases# 2. ROC for SOFA
roc_sofa <- roc(zs$Mortality.7day, zs$SOFA, ci = TRUE)## Setting levels: control = 0, case = 1## Setting direction: controls < cases# 3. Combined model: SOFA + zinc31
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 < casesroc.test(roc_zinc, roc_sofa)##
## Bootstrap test for two correlated ROC curves
##
## data: roc_zinc and roc_sofa
## D = -0.9722, boot.n = 2000, boot.stratified = 1, p-value = 0.3309
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.6859976 0.7785457roc.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.8125000plot(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)
#Print AUC, 95%CI AUC, Cut-off, Sen, spec
library(pROC)
# Individual ROC curves
roc_zinc <- roc(zs$Mortality.7day, zs$Zinc31)## Setting levels: control = 0, case = 1## Setting direction: controls > casesroc_sofa <- roc(zs$Mortality.7day, zs$SOFA)## Setting levels: control = 0, case = 1## Setting direction: controls < cases# Combined model
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)## Setting levels: control = 0, case = 1
## Setting direction: controls < cases# Updated summarize_roc function
summarize_roc <- function(roc_obj, name) {
auc_val <- round(auc(roc_obj), 2)
coords_best <- coords(
roc_obj,
x = "best",
best.method = "youden",
ret = c("threshold", "sensitivity", "specificity"),
transpose = FALSE
)
# Convert list elements to numeric if needed
threshold <- as.numeric(coords_best["threshold"])
sensitivity <- as.numeric(coords_best["sensitivity"])
specificity <- as.numeric(coords_best["specificity"])
cat(paste0("==== ", name, " ====\n"))
cat("AUC:", auc_val, "\n")
cat("Cutoff:", round(threshold, 2), "\n")
cat("Sensitivity:", round(sensitivity, 2), "\n")
cat("Specificity:", round(specificity, 2), "\n\n")
}
# Call summaries
summarize_roc(roc_zinc, "Zinc31")## ==== Zinc31 ====
## AUC: 0.69
## Cutoff: 7.65
## Sensitivity: 0.69
## Specificity: 0.66summarize_roc(roc_sofa, "SOFA")## ==== SOFA ====
## AUC: 0.78
## Cutoff: 7.5
## Sensitivity: 0.81
## Specificity: 0.64summarize_roc(roc_combined, "SOFA + Zinc31")## ==== SOFA + Zinc31 ====
## AUC: 0.81
## Cutoff: 0.1
## Sensitivity: 0.88
## Specificity: 0.69### fix
# Function to print AUC, CI, sensitivity, specificity
print_roc_info <- function(roc_obj, label) {
auc_val <- auc(roc_obj)
ci_val <- ci.auc(roc_obj)
coords_val <- coords(
roc_obj,
x = "best",
ret = c("threshold", "sensitivity", "specificity"),
best.method = "youden",
transpose = FALSE
)
# Extract values from list as numeric
threshold <- as.numeric(coords_val["threshold"])
sensitivity <- as.numeric(coords_val["sensitivity"])
specificity <- as.numeric(coords_val["specificity"])
cat(paste0(label, ":\n"))
cat(" AUC:", round(auc_val, 3), " (95% CI:", round(ci_val[1], 3), "-", round(ci_val[3], 3), ")\n")
cat(" Cutoff:", round(threshold, 2), "\n")
cat(" Sensitivity:", round(sensitivity, 4), "\n")
cat(" Specificity:", round(specificity, 4), "\n\n")
}
# Print results
print_roc_info(roc_zinc, "Zinc31")## Zinc31:
## AUC: 0.686 (95% CI: 0.526 - 0.846 )
## Cutoff: 7.65
## Sensitivity: 0.6875
## Specificity: 0.6635print_roc_info(roc_sofa, "SOFA")## SOFA:
## AUC: 0.779 (95% CI: 0.665 - 0.892 )
## Cutoff: 7.5
## Sensitivity: 0.8125
## Specificity: 0.6442print_roc_info(roc_combined, "Combined (SOFA + Zinc31)")## Combined (SOFA + Zinc31):
## AUC: 0.812 (95% CI: 0.708 - 0.917 )
## Cutoff: 0.1
## Sensitivity: 0.875
## Specificity: 0.6923#Kaplan Meier 7 day
library(survival)
library(survminer)## Loading required package: ggplot2## Warning: Can't find generic `as.gtable` in package gtable to register S3 method.
## ℹ This message is only shown to developers using devtools.
## ℹ Do you need to update gtable to the latest version?##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alpha## The following object is masked from 'package:epiDisplay':
##
## alpha## Loading required package: ggpubr##
## Attaching package: 'survminer'## The following object is masked from 'package:survival':
##
## myelomaattach(zs)## The following objects are masked from zs (pos = 17):
##
## ∆Zinc, ∆Zinc x -1, Age, Albumin, Bilirubin_total, Bloodculture, Ca,
## Cancer, CKD, CLD, Combine.∆Zinc x -1andSOFA, COPD, Creatinin, CRRT,
## CHD, Day_disease_BA, Day_free_MV_in14ds, Day_free_MV_in30ds,
## Day_free_Vasopressor_in14d, Day_free_Vasopressor_in30d,
## Day_hospitalization, Day_in_MV, Day_invasopressor_in30d,
## Day_Living_in14days, Day_Living_in30days, Day_Living_in7ds,
## Day_onset_Shock, DBP, Department, Diabetes, Fever, FiO2, GCS, GGT,
## Glucose, HCO3, HCT, HGB, HR, ID, K, Lactat, Male, Mortality.14day,
## Mortality.30day, Mortality.7day, Mortality.90day, MV, MV_admit, Na,
## Neu, Noradrenalin, Oxy_support, Pathogen.all, Pathogen_GramPos,
## Pathogen01, PCO2, PCT, PLT, PO2, Protein TP, PT, PH, RBC, RR, SBP,
## Secondary_infection, SGOT, SGPT, Shivering, Shock_admit,
## Shock_Total, Shock_within_3rd, SO2, SOFA, SOFA_Bil0123,
## Source_of_Infection, Stroke, Ure, WBC, Zinc1, Zinc3,
## Zinc3_1_percent, Zinc31, Zinc31.SD, Zinc31over7.65roc.data <- data.frame(zs$Day_Living_in7ds, zs$Mortality.7day)
#show table
baseline = Surv(zs$Day_Living_in7ds, zs$Mortality.7day)
km = survfit(baseline ~1)
summary(km)## Call: survfit(formula = baseline ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 3 120 3 0.975 0.0143 0.947 1.000
## 4 117 7 0.917 0.0252 0.869 0.967
## 5 110 5 0.875 0.0302 0.818 0.936
## 6 105 1 0.867 0.0310 0.808 0.930ggsurvplot(km, risk.table = TRUE, data=roc.data)
#Create a graph
s1 = survfit(Surv(zs$Day_Living_in7ds, zs$Mortality.7day) ~ zs$Zinc31over7.65, data=roc.data)
ggsurvplot(s1, risk.table = TRUE, data=roc.data)
##Custom the ROC curve
g1 <- ggsurvplot(
s1,
data = roc.data,
size = 1, # change line size
palette =
c("#E78","#2E9FDF"),# custom color palettes
conf.int = TRUE, # Add confidence interval
pval = TRUE, # Add p-value
risk.table = TRUE, # Add risk table
risk.table.col = "strata",# Risk table color by groups
legend.labs =
c("∆Zinc < 7.65","∆Zinc ≥ 7.65" ), # Change legend labels
risk.table.height = 0.25, # Useful to change when you have multiple groups
ggtheme = theme_light() , # Change ggplot2 theme
risk.table.y.text.col = T, # colour risk table text annotations.
risk.table.y.text = F, # show bars instead of names in text annotations
# in legend of risk table
xlab = "Time in days", # customize X axis label.
break.time.by = 1, # break X axis in time intervals by 500.
)
g1# then g1 again in console to have editable graph
#Kaplan Meier 30 day
roc.data <- data.frame(zs$Day_Living_in30days, zs$Mortality.30day)
#show table
baseline = Surv(zs$Day_Living_in30days, zs$Mortality.30day)
km = survfit(baseline ~1)
summary(km)## Call: survfit(formula = baseline ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 3 120 3 0.975 0.0143 0.947 1.000
## 4 117 7 0.917 0.0252 0.869 0.967
## 5 110 5 0.875 0.0302 0.818 0.936
## 6 105 1 0.867 0.0310 0.808 0.930
## 7 104 2 0.850 0.0326 0.788 0.916
## 8 102 1 0.842 0.0333 0.779 0.910
## 9 101 1 0.833 0.0340 0.769 0.903
## 10 100 1 0.825 0.0347 0.760 0.896
## 11 99 2 0.808 0.0359 0.741 0.882
## 12 97 2 0.792 0.0371 0.722 0.868
## 15 95 2 0.775 0.0381 0.704 0.853
## 21 93 1 0.767 0.0386 0.695 0.846
## 24 92 1 0.758 0.0391 0.685 0.839ggsurvplot(km, risk.table = TRUE, data=roc.data)
#Create a graph
s1 = survfit(Surv(zs$Day_Living_in30days, zs$Mortality.30day) ~ zs$Zinc31over7.65, data=roc.data)
ggsurvplot(s1, risk.table = TRUE, data=roc.data)
##Custom the ROC curve
g1 <- ggsurvplot(
s1,
data = roc.data,
size = 1, # change line size
palette =
c("#E78","#2E9FDF"),# custom color palettes
conf.int = TRUE, # Add confidence interval
pval = TRUE, # Add p-value
risk.table = TRUE, # Add risk table
risk.table.col = "strata",# Risk table color by groups
legend.labs =
c("∆Zinc < 7.65","∆Zinc ≥ 7.65" ), # Change legend labels
risk.table.height = 0.25, # Useful to change when you have multiple groups
ggtheme = theme_light() , # Change ggplot2 theme
risk.table.y.text.col = T, # colour risk table text annotations.
risk.table.y.text = F, # show bars instead of names in text annotations
# in legend of risk table
xlab = "Time in days", # customize X axis label.
break.time.by = 2, # break X axis in time intervals by 500.
)
g1# then g1 again in console to have editable graph