Kaplan-Meier survival plots
library(openxlsx)
library(survminer)
## Loading required package: ggplot2
## Loading required package: ggpubr
library(survival)
##
## Attaching package: 'survival'
## The following object is masked from 'package:survminer':
##
## myeloma
data <- read.xlsx("/Users/zhangmengmeng/Desktop/Rmodel/Newmodelcoxx.xlsx", sheet = 5)
Newmodelcoxx <- data
colnames(data)
## [1] "Number" "outcome" "genders" "age" "BMI"
## [6] "typeofdrugs" "injurytype" "latency" "WBC" "PLT"
## [11] "ALTonset" "ASTonset" "ALPonset" "GGTonset" "TBILonset"
## [16] "ALB" "TBA" "Cr" "INR" "Timetoevents"
## [21] "Death28d" "Death90d" "Death6m" "riskscore" "riskgroup"
## [26] "LP" "lntbil" "lninr" "crmgdl" "lncr"
## [31] "meld" "riskgroupRB" "riskgroupnhl"
fit <- survfit(Surv(Timetoevents, outcome) ~ riskgroup, data = Newmodelcoxx)
##The risk groups were stratified according to the previously established linear formula: LP = 0.05 × TBIL(onset) [mg/dL] - 0.11 × albumin [g/L] + 0.45 × INR + 3.22. Using a cutoff value of 1.05, patients with LP > 1.05 were classified into the high-risk group, while those with LP ≤ 1.05 were assigned to the low-risk group.
ggsurvplot(fit,
data = Newmodelcoxx,
xlim = c(0, 180),
break.time.by = 30, # One mark every 30 days, the duration can be adjusted to 100 days if the follow-up period is extended.
# Change legends: title & labels
legend.title = "Testing set",
legend.labs = c("High", "Low"),
palette = c("#C72228", "#0C4E9B"),
# Add p-value and tervals
pval = TRUE,
conf.int = TRUE,
pval.size = 3, # Reduce the font size of the p-value
pval.coord = c(50, 0.25), # Adjust the position of the P-value
# Add risk table
risk.table = TRUE,
tables.height = 0.25,
risk.table.col = "strata",
risk.table.y.text = TRUE,
risk.table.y.text.col = TRUE,
risk.table.height = 0.25,
risk.table.fontsize = 3
,
ncensor.plot = TRUE,
ncensor.plot.height = 0.2,
font.ncensor = 8,
# Risk Table Theme Settings
tables.theme = theme(
text = element_text(size = 4), # Further reduce the font size of the risk table
axis.title = element_text(size = 6),
axis.text = element_text(size = 6),
legend.position = "none" ,
plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm")
),
xlab = "Time (days)",
ylab = "Survival Probability",
ggtheme = theme_bw(base_size = 8)
)
Time-dependent roc curve
library(timeROC)
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
library(survival)
roc_res <- timeROC::timeROC(
T = Newmodelcoxx$Timetoevents, # Timetoevents just Survival time
delta = Newmodelcoxx$outcome, # 1=LT/Death, 0=no LT/Death
marker = Newmodelcoxx$riskscore, # LP score as before
cause = 1,
weighting = "cox",
times = c(31, 90, 180), # This is displayed in this manner due to constraints imposed by our verification queue. time point should be 30, 90, 180, 270
iid = FALSE
)
print(roc_res)
## Time-dependent-Roc curve estimated using IPCW (n=741, without competing risks).
## Cases Survivors Censored AUC (%)
## t=31 18 723 0 93.68
## t=90 18 695 28 95.39
## t=180 46 686 9 94.73
##
## Method used for estimating IPCW:cox
##
## Total computation time : 0.31 secs.
plot(roc_res, time = 31, col = "#C72228", lwd = 2.5,title = FALSE)
plot(roc_res, time = 90, col = "#0C4E9B", lwd = 2.5,add = TRUE)
plot(roc_res, time = 180, col = "#F1b656", lwd = 2.5,add = TRUE)
legend("bottomright", legend = c("30 days AUC: 0.937", "90 days AUC: 0.953", "180 days AUC: 0.947"),
col = c("#C72228", "#0C4E9B", "#F1b656"), lty = 1, lwd = 2.5,cex = 0.8)
Calibration curves
library(rms)
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
library(survival)
library(openxlsx)
dd_validation <- read.xlsx("/Users/zhangmengmeng/Desktop/Rmodel/Newmodelcoxx.xlsx", sheet = 5)
Newmodelcoxx <- dd_validation
validation_data <- Newmodelcoxx
dd_validation_dist <- datadist(validation_data)
options(datadist = "dd_validation_dist")
fit_validation <- cph(Surv(Timetoevents, outcome) ~ LP,
data = validation_data, x = TRUE, y = TRUE, surv = TRUE)
cal_validation <- calibrate(fit_validation,
cmethod = 'KM',
method = "boot",
u = 30, # The assessment time point may be 30, 90, 180,270 days.
m = floor(nrow(validation_data)/3),
B = 500)
## Using Cox survival estimates at 30 Days
plot(cal_validation,
lwd=2,
lty=1,
errbar.col=c("#0C4E9B"),
cex.lab=1.2,cex.axis=1,cex.main=1.2,cex.sub=0.1
)
lines(cal_validation[,c('mean.predicted',"KM")],
type= 'b', # 连线的类型，可以是"p","b","o"
lwd = 3, # 连线的粗细
pch = 16, #点的形状，可以是0-20
col = "#C72228") # 连线的颜色
box(lwd = 2) # 边框粗细
abline(0,1,lty=3,#对角线为虚线
lwd=2,#对角线的粗细
col="grey70"#对角线的颜色
)
DCA
library(rmda)
Newmodelcoxx$event_occurred <- ifelse(Newmodelcoxx$outcome == 1, 1, 0)
Newmodelcoxx$risk_prob <- 1 / (1 + exp(-Newmodelcoxx$LP))
dca_overall <- decision_curve(
event_occurred ~ risk_prob,
data = Newmodelcoxx,
family = binomial,
thresholds = seq(0, 0.5, by = 0.01),
study.design = "cohort",
bootstraps = 100
)
plot_decision_curve(
dca_overall,
curve.names = "Overall Risk Model",
col = c("#0C4E9B"),
lwd = 2,
confidence.intervals = FALSE,
legend.position = "bottomright",
standardize = FALSE,
cost.benefit.axis = FALSE,
xlim = c(0, 0.5)
)
Comparative Performance with Existing Prognostic Models (AUC and DCA
at Total time)
rm(list=ls())
library(pROC)
data <- read.xlsx("/Users/zhangmengmeng/Desktop/Rmodel/Newmodelcoxx.xlsx", sheet = 5)
Newmodelcoxx <- data
colnames(data)
## [1] "Number" "outcome" "genders" "age" "BMI"
## [6] "typeofdrugs" "injurytype" "latency" "WBC" "PLT"
## [11] "ALTonset" "ASTonset" "ALPonset" "GGTonset" "TBILonset"
## [16] "ALB" "TBA" "Cr" "INR" "Timetoevents"
## [21] "Death28d" "Death90d" "Death6m" "riskscore" "riskgroup"
## [26] "LP" "lntbil" "lninr" "crmgdl" "lncr"
## [31] "meld" "riskgroupRB" "riskgroupnhl"
roc1<- roc(Newmodelcoxx$outcome, Newmodelcoxx$LP) ## LP corresponding legend new model
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc2<- roc(Newmodelcoxx$outcome, Newmodelcoxx$riskgroupnhl) ##riskgroupnhl corresponding legend nHy's Law
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc3<- roc(Newmodelcoxx$outcome, Newmodelcoxx$riskgroupRB) ##riskgroupRB corresponding legend Robles-Diaz
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc4<- roc(Newmodelcoxx$outcome, Newmodelcoxx$meld) ##meld corresponding legend MELD
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc1, col = "#C72228", lwd = 2.5,title = FALSE)
plot(roc2, col = "#0C4E9B", lwd = 2.5,add = TRUE)
plot(roc3, col = "#F1b656", lwd = 2.5,add = TRUE)
plot(roc4, col = "#037F77", lwd = 2.5,add = TRUE)
print(roc1)
##
## Call:
## roc.default(response = Newmodelcoxx$outcome, predictor = Newmodelcoxx$LP)
##
## Data: Newmodelcoxx$LP in 686 controls (Newmodelcoxx$outcome 0) < 55 cases (Newmodelcoxx$outcome 1).
## Area under the curve: 0.9195
print(roc2)
##
## Call:
## roc.default(response = Newmodelcoxx$outcome, predictor = Newmodelcoxx$riskgroupnhl)
##
## Data: Newmodelcoxx$riskgroupnhl in 686 controls (Newmodelcoxx$outcome 0) < 55 cases (Newmodelcoxx$outcome 1).
## Area under the curve: 0.6572
print(roc3)
##
## Call:
## roc.default(response = Newmodelcoxx$outcome, predictor = Newmodelcoxx$riskgroupRB)
##
## Data: Newmodelcoxx$riskgroupRB in 686 controls (Newmodelcoxx$outcome 0) < 55 cases (Newmodelcoxx$outcome 1).
## Area under the curve: 0.6121
print(roc4)
##
## Call:
## roc.default(response = Newmodelcoxx$outcome, predictor = Newmodelcoxx$meld)
##
## Data: Newmodelcoxx$meld in 686 controls (Newmodelcoxx$outcome 0) < 55 cases (Newmodelcoxx$outcome 1).
## Area under the curve: 0.9214
legend(x = 0.646, y = 0.29, #"bottomright"
legend = c("new model AUC: 0.920 ", ##new model refers to the LP score
"nHy's Law AUC: 0.657 ", ##TBIL > 2 × ULN, nR ≥ 5, where nR = (ALT – maximum value) or (AST – maximum value) / ULN divided by ALP / ULN
"Robles-Diaz AUC: 0.612 ", ##AST> 17.3×ULN and TBIL> 6.6×ULN；or AST≤17.3×ULN+AST/ALT > 1.5
"MELD AUC: 0.921 " ##medl=10×[(0.957×ln (Cr))+(0.378×ln (TBIL))+(1.12×ln (INR))]+6.43
),
col = c("#C72228", "#0C4E9B", "#F1b656", "#037F77"),
lty = 1, lwd = 2.5,cex = 0.8,
y.intersp = 1)
delong_test <- roc.test(roc1, roc2, method = "delong")
delong_test ##If possible,please further compare the DeLong test.
##
## DeLong's test for two correlated ROC curves
##
## data: roc1 and roc2
## Z = 7.317, p-value = 2.536e-13
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
## 0.1920571 0.3325917
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.9195070 0.6571826
##DCA（Following on from the above）
library(rmda)
dca_newmodel <- decision_curve(outcome ~ LP, data = Newmodelcoxx, family = binomial)
dca_nHysLaw <- decision_curve(outcome ~ riskgroupnhl, data = Newmodelcoxx, family = binomial)
dca_RB <- decision_curve(outcome ~ riskgroupRB, data = Newmodelcoxx, family = binomial)
dca_meld <- decision_curve(outcome ~ meld, data = Newmodelcoxx, family = binomial)
plot_decision_curve(
list(dca_newmodel, dca_nHysLaw, dca_RB, dca_meld),
curve.names = c("newmodel","nHy's Law","Robles-Diaz","MELD"),
col = c("#C72228", "#0C4E9B", "#F1b656","#037F77"),
lty = c(1,1,1,1),
lwd = c(3,3,3,3),
confidence.intervals = FALSE,
legend.position = "none",
standardize = FALSE,
cost.benefit.axis = FALSE
)
## Note: When multiple decision curves are plotted, decision curves for 'All' are calculated using the prevalence from the first DecisionCurve object in the list provided.
legend("bottomright",
legend = c("newmodel","nHy's Law","Robles-Diaz","MELD"),
col = c("#C72228", "#0C4E9B", "#F1b656", "#037F77"),
lty = 1,
lwd = 3,
cex = 0.8)
Comparative Performance with Existing Prognostic Models( AUC and DCA
at 30 days )
rm(list=ls())
library(pROC)
data <- read.xlsx("/Users/zhangmengmeng/Desktop/Rmodel/Newmodelcoxx.xlsx", sheet = 5)
Newmodelcoxx <- data
colnames(data)
## [1] "Number" "outcome" "genders" "age" "BMI"
## [6] "typeofdrugs" "injurytype" "latency" "WBC" "PLT"
## [11] "ALTonset" "ASTonset" "ALPonset" "GGTonset" "TBILonset"
## [16] "ALB" "TBA" "Cr" "INR" "Timetoevents"
## [21] "Death28d" "Death90d" "Death6m" "riskscore" "riskgroup"
## [26] "LP" "lntbil" "lninr" "crmgdl" "lncr"
## [31] "meld" "riskgroupRB" "riskgroupnhl"
roc1<- roc(Newmodelcoxx$Death28d, Newmodelcoxx$LP)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc2<- roc(Newmodelcoxx$Death28d, Newmodelcoxx$riskgroupnhl)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc3<- roc(Newmodelcoxx$Death28d, Newmodelcoxx$riskgroupRB)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc4<- roc(Newmodelcoxx$Death28d, Newmodelcoxx$meld)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc1, col = "#C72228", lwd = 2.5,title = FALSE)
plot(roc2, col = "#0C4E9B", lwd = 2.5,add = TRUE)
plot(roc3, col = "#F1b656", lwd = 2.5,add = TRUE)
plot(roc4, col = "#037F77", lwd = 2.5,add = TRUE)
print(roc1)
##
## Call:
## roc.default(response = Newmodelcoxx$Death28d, predictor = Newmodelcoxx$LP)
##
## Data: Newmodelcoxx$LP in 723 controls (Newmodelcoxx$Death28d 0) < 18 cases (Newmodelcoxx$Death28d 1).
## Area under the curve: 0.9368
print(roc2)
##
## Call:
## roc.default(response = Newmodelcoxx$Death28d, predictor = Newmodelcoxx$riskgroupnhl)
##
## Data: Newmodelcoxx$riskgroupnhl in 723 controls (Newmodelcoxx$Death28d 0) < 18 cases (Newmodelcoxx$Death28d 1).
## Area under the curve: 0.645
print(roc3)
##
## Call:
## roc.default(response = Newmodelcoxx$Death28d, predictor = Newmodelcoxx$riskgroupRB)
##
## Data: Newmodelcoxx$riskgroupRB in 723 controls (Newmodelcoxx$Death28d 0) < 18 cases (Newmodelcoxx$Death28d 1).
## Area under the curve: 0.595
print(roc4)
##
## Call:
## roc.default(response = Newmodelcoxx$Death28d, predictor = Newmodelcoxx$meld)
##
## Data: Newmodelcoxx$meld in 723 controls (Newmodelcoxx$Death28d 0) < 18 cases (Newmodelcoxx$Death28d 1).
## Area under the curve: 0.9358
legend(x = 0.75, y = 0.35, #"bottomright"
legend = c("new model AUC: 0.937 ",
"nHy's Law AUC: 0.645 ",
"Robles-Diaz AUC: 0.595 ",
"MELD AUC: 0.936 "),
col = c("#C72228", "#0C4E9B", "#F1b656", "#037F77"),
lty = 1, lwd = 2.5,cex = 0.8,
y.intersp = 1)
delong_test1 <- roc.test(roc1, roc2, method = "delong")
delong_test2 <- roc.test(roc1, roc3, method = "delong")
delong_test3 <- roc.test(roc1, roc4, method = "delong")
delong_results <- data.frame(
Comparison = c("roc1 vs roc2", "roc1 vs roc3", "roc1 vs roc4"),
Statistic = c(delong_test1$statistic, delong_test2$statistic, delong_test3$statistic),
P_Value = c(delong_test1$p.value, delong_test2$p.value, delong_test3$p.value),
Method = c(delong_test1$method, delong_test2$method, delong_test3$method)
)
delong_results$Significance <- ifelse(delong_results$P_Value < 0.001, "***",
ifelse(delong_results$P_Value < 0.01, "**",
ifelse(delong_results$P_Value < 0.05, "*", "ns")))
print(delong_results)
## Comparison Statistic P_Value
## 1 roc1 vs roc2 4.8454318 1.263366e-06
## 2 roc1 vs roc3 6.9616729 3.362556e-12
## 3 roc1 vs roc4 0.0657481 9.475784e-01
## Method Significance
## 1 DeLong's test for two correlated ROC curves ***
## 2 DeLong's test for two correlated ROC curves ***
## 3 DeLong's test for two correlated ROC curves ns
#DCA
library(rmda)
dca_newmodel <- decision_curve(Death28d ~ LP, data = Newmodelcoxx, family = binomial)
dca_nHysLaw <- decision_curve(Death28d ~ riskgroupnhl, data = Newmodelcoxx, family = binomial)
dca_RB <- decision_curve(Death28d ~ riskgroupRB, data = Newmodelcoxx, family = binomial)
dca_meld <- decision_curve(Death28d ~ meld, data = Newmodelcoxx, family = binomial)
plot_decision_curve(
list(dca_newmodel, dca_nHysLaw, dca_RB, dca_meld),
curve.names = c("newmodel","nHy's Law","Robles-Diaz","MELD"),
col = c("#C72228", "#0C4E9B", "#F1b656","#037F77"),
lty = c(1,1,1,1),
lwd = c(3,3,3,3),
confidence.intervals = FALSE,
legend.position = "none",
standardize = FALSE,
cost.benefit.axis = FALSE
)
## Note: When multiple decision curves are plotted, decision curves for 'All' are calculated using the prevalence from the first DecisionCurve object in the list provided.
legend("bottomright",
legend = c("newmodel","nHy's Law","Robles-Diaz","MELD"),
col = c("#C72228", "#0C4E9B", "#F1b656", "#037F77"),
lty = 1,
lwd = 3,
cex = 0.8)
