Competing death and Fracture risk assessment
Thach Tran
10/10/2021
Competing Death and Fracture Risk Assessment
using the “standard” longitudinal study design:
library(survival); library(tableone); library(prodlim); library(riskRegression); library(CalibrationCurves)## riskRegression version 2023.09.08## Loading required package: rms## Loading required package: Hmisc##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:base':
##
## format.pval, units## Loading required package: ggplot2Aim 1: Development of prediction models using different methods
wpps<- read.csv("C:\\Thach\\Research projects\\MSM and competing risk\\Data_Analysis\\wpp_s60.csv")
head(wpps)## Fx1st_sen BMI_0 neuro_0 rheum_0 hypert_0 diabetes_0 respi_0 cancer_0 cvd_0
## 1 0 27.34 0 0 0 1 0 0 1
## 2 0 27.34 0 0 1 0 0 0 1
## 3 0 24.06 1 0 1 0 0 0 0
## 4 0 26.13 0 0 1 0 0 0 1
## 5 0 27.34 0 0 1 1 0 0 1
## 6 0 26.13 0 0 0 0 0 0 1
## prior_fx age StudyID sex falls_n death_sen fnbmd kmfu_fx kmfu_death
## 1 0 73.07 3 0 0 1 -0.3250 0.5475702 0.5475702
## 2 0 68.47 9 0 0 1 0.2250 5.0513347 5.0513347
## 3 0 62.10 10 1 0 1 1.3417 26.5735797 26.5735797
## 4 0 75.81 24 1 0 0 2.1417 2.8555784 2.8555784
## 5 0 63.53 26 0 0 0 0.2833 9.2621492 9.2621492
## 6 0 76.24 28 1 1 1 2.8833 20.5338809 20.5338809
## group
## 1 2
## 2 2
## 3 2
## 4 0
## 5 0
## 6 21.1 Baseline characteristics
factorVars <- c("group","sex","falls_n","prior_fx","cvd_0","respi_0","diabetes_0","hypert_0","cancer_0","rheum_0","neuro_0")
Vars <- c("sex","age","fnbmd","falls_n","prior_fx","kmfu_fx","kmfu_death","BMI_0","cvd_0","respi_0","diabetes_0","hypert_0","cancer_0","rheum_0","neuro_0")
table.1d <- CreateTableOne(vars = Vars, strata = "group", data = wpps, factorVars = factorVars)
print(table.1d, nonnormal = c("kmfu_fx","kmfu_death"))## Stratified by group
## 0 1
## n 1045 505
## sex = 1 (%) 646 (61.8) 389 (77.0)
## age (mean (SD)) 68.34 (5.74) 69.43 (6.31)
## fnbmd (mean (SD)) 1.16 (1.16) 1.80 (1.06)
## falls_n (%)
## 0 844 (80.8) 377 (74.7)
## 1 135 (12.9) 91 (18.0)
## 2 66 ( 6.3) 37 ( 7.3)
## prior_fx = 1 (%) 123 (11.8) 88 (17.4)
## kmfu_fx (median [IQR]) 9.65 [6.08, 13.65] 7.20 [3.24, 11.14]
## kmfu_death (median [IQR]) 9.65 [6.08, 13.65] 11.64 [7.93, 16.80]
## BMI_0 (mean (SD)) 27.19 (4.11) 26.56 (3.98)
## cvd_0 = 1 (%) 322 (30.8) 177 (35.0)
## respi_0 = 1 (%) 103 ( 9.9) 73 (14.5)
## diabetes_0 = 1 (%) 122 (11.7) 41 ( 8.1)
## hypert_0 = 1 (%) 539 (51.6) 257 (50.9)
## cancer_0 = 1 (%) 88 ( 8.4) 54 (10.7)
## rheum_0 = 1 (%) 23 ( 2.2) 30 ( 5.9)
## neuro_0 = 1 (%) 60 ( 5.7) 48 ( 9.5)
## Stratified by group
## 2 p test
## n 271
## sex = 1 (%) 120 (44.3) <0.001
## age (mean (SD)) 71.43 (6.93) <0.001
## fnbmd (mean (SD)) 1.31 (1.30) <0.001
## falls_n (%) 0.014
## 0 227 (83.8)
## 1 33 (12.2)
## 2 11 ( 4.1)
## prior_fx = 1 (%) 29 (10.7) 0.004
## kmfu_fx (median [IQR]) 6.50 [3.55, 11.37] <0.001 nonnorm
## kmfu_death (median [IQR]) 6.50 [3.55, 11.37] <0.001 nonnorm
## BMI_0 (mean (SD)) 26.82 (1.45) 0.008
## cvd_0 = 1 (%) 124 (45.8) <0.001
## respi_0 = 1 (%) 38 (14.0) 0.014
## diabetes_0 = 1 (%) 44 (16.2) 0.003
## hypert_0 = 1 (%) 138 (50.9) 0.960
## cancer_0 = 1 (%) 29 (10.7) 0.258
## rheum_0 = 1 (%) 15 ( 5.5) <0.001
## neuro_0 = 1 (%) 17 ( 6.3) 0.0211.2 Descriptive curves from different methods
(1.2.1) Conventional method
(1.2.1.1) Kaplan-Meier curve
table(wpps$Fx1st_sen, wpps$sex)##
## 0 1
## 0 550 766
## 1 116 389mfit.con.f <- survfit(Surv(kmfu_fx, Fx1st_sen) ~ sex, data = wpps)
mfit.con.f## Call: survfit(formula = Surv(kmfu_fx, Fx1st_sen) ~ sex, data = wpps)
##
## n events median 0.95LCL 0.95UCL
## sex=0 666 116 NA 21.4 NA
## sex=1 1155 389 17.2 16.3 18.3table(wpps$death_sen, wpps$sex)##
## 0 1
## 0 490 978
## 1 176 177mfit.con.d <- survfit(Surv(kmfu_death, death_sen) ~ sex, data = wpps)
mfit.con.d## Call: survfit(formula = Surv(kmfu_death, death_sen) ~ sex, data = wpps)
##
## n events median 0.95LCL 0.95UCL
## sex=0 666 176 21.4 19.4 26.6
## sex=1 1155 177 NA 26.6 NA(1.2.1.2) Cox’s PH model
## Conventional model for fracture
cfit.con.f <- coxph(Surv(kmfu_fx, Fx1st_sen) ~ sex + age + fnbmd + falls_n + prior_fx,
data = wpps, id = StudyID)
summary(cfit.con.f)## Call:
## coxph(formula = Surv(kmfu_fx, Fx1st_sen) ~ sex + age + fnbmd +
## falls_n + prior_fx, data = wpps, id = StudyID)
##
## n= 1821, number of events= 505
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex 0.394651 1.483867 0.111300 3.546 0.000391 ***
## age 0.042794 1.043723 0.008007 5.345 9.07e-08 ***
## fnbmd 0.404080 1.497923 0.047510 8.505 < 2e-16 ***
## falls_n 0.009253 1.009296 0.075632 0.122 0.902632
## prior_fx 0.397371 1.487908 0.123750 3.211 0.001322 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 1.484 0.6739 1.1930 1.846
## age 1.044 0.9581 1.0275 1.060
## fnbmd 1.498 0.6676 1.3647 1.644
## falls_n 1.009 0.9908 0.8702 1.171
## prior_fx 1.488 0.6721 1.1675 1.896
##
## Concordance= 0.693 (se = 0.013 )
## Likelihood ratio test= 221.9 on 5 df, p=<2e-16
## Wald test = 228.3 on 5 df, p=<2e-16
## Score (logrank) test = 235.9 on 5 df, p=<2e-16zph.cfit0f <- cox.zph(cfit.con.f)
zph.cfit0f## chisq df p
## sex 0.024096 1 0.88
## age 1.118674 1 0.29
## fnbmd 0.179847 1 0.67
## falls_n 0.933495 1 0.33
## prior_fx 0.000869 1 0.98
## GLOBAL 3.268397 5 0.66### Proportional hazards assumption
windows(width = 14, height = 12)
par(mfrow = c(3,2), oma = c(0,0,2,0))
plot(zph.cfit0f[1])
abline(h=coef(cfit.con.f)[1], lty=2, col=2)
plot(zph.cfit0f[2])
abline(h=coef(cfit.con.f)[2], lty=2, col=2)
plot(zph.cfit0f[3])
abline(h=coef(cfit.con.f)[3], lty=2, col=2)
plot(zph.cfit0f[4])
abline(h=coef(cfit.con.f)[4], lty=2, col=2)
plot(zph.cfit0f[5])
abline(h=coef(cfit.con.f)[5], lty=2, col=2)
mtext("Proportional hazards assumption: Conventional model", line=0, side=3, outer=TRUE, cex=1.8)
## Conventional model for death without a fracture
cfit.con.d <- coxph(Surv(kmfu_death, death_sen) ~ sex + age + fnbmd + falls_n + prior_fx,
data = wpps, id = StudyID)
summary(cfit.con.d)## Call:
## coxph(formula = Surv(kmfu_death, death_sen) ~ sex + age + fnbmd +
## falls_n + prior_fx, data = wpps, id = StudyID)
##
## n= 1821, number of events= 353
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex -0.747923 0.473349 0.116521 -6.419 1.37e-10 ***
## age 0.131206 1.140203 0.009097 14.422 < 2e-16 ***
## fnbmd 0.127560 1.136054 0.054395 2.345 0.019 *
## falls_n -0.262873 0.768839 0.105832 -2.484 0.013 *
## prior_fx 0.097152 1.102028 0.168334 0.577 0.564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 0.4733 2.1126 0.3767 0.5948
## age 1.1402 0.8770 1.1201 1.1607
## fnbmd 1.1361 0.8802 1.0212 1.2639
## falls_n 0.7688 1.3007 0.6248 0.9461
## prior_fx 1.1020 0.9074 0.7923 1.5328
##
## Concordance= 0.749 (se = 0.014 )
## Likelihood ratio test= 273.6 on 5 df, p=<2e-16
## Wald test = 293.5 on 5 df, p=<2e-16
## Score (logrank) test = 319.4 on 5 df, p=<2e-16zph.cfit0d <- cox.zph(cfit.con.d)
zph.cfit0d## chisq df p
## sex 0.018 1 0.893
## age 1.445 1 0.229
## fnbmd 0.485 1 0.486
## falls_n 1.018 1 0.313
## prior_fx 3.454 1 0.063
## GLOBAL 6.833 5 0.233
(1.2.2) Competing risk methods
Data preparation
etime <- with(wpps, ifelse(Fx1st_sen == 0, kmfu_death, kmfu_fx))
event <- with(wpps, ifelse(Fx1st_sen == 0, 2*death_sen, 1))
event <- factor(event, 0:2, labels=c("Event-free", "Fracture", "Death"))
table(event)## event
## Event-free Fracture Death
## 1045 505 271Descriptive Kaplan-Meier curves
mfit.cs <- survfit(Surv(etime, event) ~ sex, data = wpps)
print(mfit.cs)## Call: survfit(formula = Surv(etime, event) ~ sex, data = wpps)
##
## n nevent rmean*
## sex=0, (s0) 666 0 15.916935
## sex=1, (s0) 1155 0 14.892149
## sex=0, Fracture 666 116 5.696660
## sex=1, Fracture 1155 389 10.850160
## sex=0, Death 666 151 7.519875
## sex=1, Death 1155 120 3.391161
## *restricted mean time in state (max time = 29.13347 )mfit.cs$transitions ## to
## from Fracture Death (censored)
## (s0) 505 271 1045
## Fracture 0 0 0
## Death 0 0 0windows(width = 16, height = 12)
par(mfrow = c(1,2), oma = c(0,0,2,0))
plot(mfit.con.f[2], fun = "event", mark.time = FALSE, col = "black", lty = 1, lwd = 2, conf.int = FALSE,
xlab="Time (years)", ylab="Probability", ylim = c(0,1))
lines(mfit.con.d[2], fun = "event", mark.time = FALSE, col = "black", lty = 2, lwd = 2, conf.int = FALSE)
lines(mfit.cs[2,2], mark.time = FALSE, col = "red", lty = 1, lwd = 2, conf.int = FALSE)
lines(mfit.cs[2,3], mark.time = FALSE, col = "red", lty = 2, lwd = 2, conf.int = FALSE)
title(main = "Women", cex.main = 1.2)
plot(mfit.con.f[1], fun = "event", mark.time = FALSE, col = "black", lty = 1, lwd = 2, conf.int = FALSE,
xlab="Time (years)", ylab="Probability", ylim = c(0,1))
lines(mfit.con.d[1], fun = "event", mark.time = FALSE, col = "black", lty = 2, lwd = 2, conf.int = FALSE)
lines(mfit.cs[1,2], mark.time = FALSE, col = "red", lty = 1, lwd = 2, conf.int = FALSE)
lines(mfit.cs[1,3], mark.time = FALSE, col = "red", lty = 2, lwd = 2, conf.int = FALSE)
title(main = "Men", cex.main = 1.2)
legend("topleft", c("Fracture (Conventional)","Death (Conventional)", "Fracture (Competing risk)", "death_sen (Competing risk)"),
col = c("black","black","red","red"), lty = c(1,2,1,2), lwd = 2, bty = 'n', cex = 0.8)
mtext("Risk of fracture and death by different methods", line = 0, side = 3,
outer = TRUE, cex = 1.6)
(1.2.2.1) Cause-specific method
## For descriptive graphs:
csfit <- coxph(Surv(etime, event) ~ sex + age + fnbmd + falls_n + prior_fx,
data = wpps, id = StudyID)
summary(csfit)## Call:
## coxph(formula = Surv(etime, event) ~ sex + age + fnbmd + falls_n +
## prior_fx, data = wpps, id = StudyID)
##
## n= 1821, number of events= 776
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## sex_1:2 0.394651 1.483867 0.111300 0.111066 3.553 0.00038 ***
## age_1:2 0.042794 1.043723 0.008007 0.008289 5.163 2.43e-07 ***
## fnbmd_1:2 0.404080 1.497923 0.047510 0.047542 8.499 < 2e-16 ***
## falls_n_1:2 0.009253 1.009296 0.075632 0.079901 0.116 0.90781
## prior_fx_1:2 0.397371 1.487908 0.123750 0.129847 3.060 0.00221 **
## sex_1:3 -0.733062 0.480436 0.131771 0.139419 -5.258 1.46e-07 ***
## age_1:3 0.127076 1.135503 0.010267 0.010603 11.985 < 2e-16 ***
## fnbmd_1:3 0.031901 1.032415 0.061144 0.070864 0.450 0.65259
## falls_n_1:3 -0.352542 0.702899 0.129767 0.128472 -2.744 0.00607 **
## prior_fx_1:3 0.010446 1.010501 0.205829 0.212640 0.049 0.96082
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex_1:2 1.4839 0.6739 1.1936 1.8447
## age_1:2 1.0437 0.9581 1.0269 1.0608
## fnbmd_1:2 1.4979 0.6676 1.3647 1.6442
## falls_n_1:2 1.0093 0.9908 0.8630 1.1804
## prior_fx_1:2 1.4879 0.6721 1.1536 1.9191
## sex_1:3 0.4804 2.0814 0.3656 0.6314
## age_1:3 1.1355 0.8807 1.1121 1.1593
## fnbmd_1:3 1.0324 0.9686 0.8985 1.1862
## falls_n_1:3 0.7029 1.4227 0.5464 0.9042
## prior_fx_1:3 1.0105 0.9896 0.6661 1.5330
##
## Concordance= 0.711 (se = 0.01 )
## Likelihood ratio test= 417 on 10 df, p=<2e-16
## Wald test = 463.9 on 10 df, p=<2e-16
## Score (logrank) test = 461.1 on 10 df, p=<2e-16, Robust = 291.8 p=<2e-16
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).zph.csfit <- cox.zph(csfit)
zph.csfit## chisq df p
## sex_1:2 0.029118 1 0.865
## age_1:2 1.114354 1 0.291
## fnbmd_1:2 0.182523 1 0.669
## falls_n_1:2 0.887880 1 0.346
## prior_fx_1:2 0.000837 1 0.977
## sex_1:3 0.125561 1 0.723
## age_1:3 1.210546 1 0.271
## fnbmd_1:3 0.120813 1 0.728
## falls_n_1:3 0.861422 1 0.353
## prior_fx_1:3 3.036326 1 0.081
## GLOBAL 8.724266 10 0.558windows(width = 14, height = 12)
par(mfrow=c(2,5), oma=c(0,0,2,0))
plot(zph.csfit[1])
abline(h=coef(csfit)[1], lty=2, col=1)
title(main = "Fracture model", cex.main= 1.5)
plot(zph.csfit[2])
abline(h=coef(csfit)[2], lty=2, col=1)
plot(zph.csfit[3])
abline(h=coef(csfit)[3], lty=2, col=1)
plot(zph.csfit[4])
abline(h=coef(csfit)[4], lty=2, col=1)
plot(zph.csfit[5])
abline(h=coef(csfit)[5], lty=2, col=1)
plot(zph.csfit[6])
abline(h=coef(csfit)[6], lty=2, col=2)
title(main = "Mortality model", cex.main= 1.5)
plot(zph.csfit[7])
abline(h=coef(csfit)[7], lty=2, col=2)
plot(zph.csfit[8])
abline(h=coef(csfit)[8], lty=2, col=2)
plot(zph.csfit[9])
abline(h=coef(csfit)[9], lty=2, col=2)
plot(zph.csfit[10])
abline(h=coef(csfit)[10], lty=2, col=2)
mtext("Sensitivity analysis- PH assumption: Cause-specific model", line=0, side=3, outer=TRUE, cex=1.8)
# For prediction of absolute risk: 'CSC' from riskRegression using the combined coefficients
CSC.fit <- CSC(Hist(etime, event)~ sex + age + fnbmd + falls_n + prior_fx,
data = wpps, method = "breslow")
print(CSC.fit)## CSC(formula = Hist(etime, event) ~ sex + age + fnbmd + falls_n +
## prior_fx, data = wpps, method = "breslow")
##
## Uncensored response of a competing.risks model
##
## No.Observations: 1821
##
## Pattern:
##
## Cause event
## Event-free 1045
## Fracture 505
## Death 271
## unknown 0
##
##
## ----------> Cause: Event-free
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ sex + age + fnbmd +
## falls_n + prior_fx, x = TRUE, y = TRUE, method = "breslow")
##
## n= 1821, number of events= 1045
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex 0.055579 1.057152 0.068079 0.816 0.41428
## age 0.069721 1.072209 0.006156 11.326 < 2e-16 ***
## fnbmd -0.080099 0.923025 0.031174 -2.569 0.01019 *
## falls_n -0.123397 0.883912 0.059071 -2.089 0.03671 *
## prior_fx 0.339550 1.404316 0.099942 3.397 0.00068 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 1.0572 0.9459 0.9251 1.2081
## age 1.0722 0.9327 1.0593 1.0852
## fnbmd 0.9230 1.0834 0.8683 0.9812
## falls_n 0.8839 1.1313 0.7873 0.9924
## prior_fx 1.4043 0.7121 1.1545 1.7082
##
## Concordance= 0.614 (se = 0.01 )
## Likelihood ratio test= 138.7 on 5 df, p=<2e-16
## Wald test = 152.3 on 5 df, p=<2e-16
## Score (logrank) test = 154.3 on 5 df, p=<2e-16
##
##
##
## ----------> Cause: Fracture
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ sex + age + fnbmd +
## falls_n + prior_fx, x = TRUE, y = TRUE, method = "breslow")
##
## n= 1821, number of events= 505
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex 0.394636 1.483844 0.111300 3.546 0.000392 ***
## age 0.042792 1.043720 0.008007 5.344 9.07e-08 ***
## fnbmd 0.404020 1.497834 0.047507 8.504 < 2e-16 ***
## falls_n 0.009284 1.009327 0.075632 0.123 0.902308
## prior_fx 0.397346 1.487871 0.123748 3.211 0.001323 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 1.484 0.6739 1.1930 1.846
## age 1.044 0.9581 1.0275 1.060
## fnbmd 1.498 0.6676 1.3647 1.644
## falls_n 1.009 0.9908 0.8703 1.171
## prior_fx 1.488 0.6721 1.1674 1.896
##
## Concordance= 0.693 (se = 0.013 )
## Likelihood ratio test= 221.9 on 5 df, p=<2e-16
## Wald test = 228.3 on 5 df, p=<2e-16
## Score (logrank) test = 235.8 on 5 df, p=<2e-16
##
##
##
## ----------> Cause: Death
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ sex + age + fnbmd +
## falls_n + prior_fx, x = TRUE, y = TRUE, method = "breslow")
##
## n= 1821, number of events= 271
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex -0.73304 0.48044 0.13177 -5.563 2.65e-08 ***
## age 0.12707 1.13550 0.01027 12.377 < 2e-16 ***
## fnbmd 0.03190 1.03241 0.06114 0.522 0.6019
## falls_n -0.35250 0.70293 0.12977 -2.716 0.0066 **
## prior_fx 0.01044 1.01049 0.20583 0.051 0.9596
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 0.4804 2.0814 0.3711 0.6220
## age 1.1355 0.8807 1.1129 1.1586
## fnbmd 1.0324 0.9686 0.9158 1.1639
## falls_n 0.7029 1.4226 0.5451 0.9065
## prior_fx 1.0105 0.9896 0.6750 1.5126
##
## Concordance= 0.745 (se = 0.016 )
## Likelihood ratio test= 195.1 on 5 df, p=<2e-16
## Wald test = 210 on 5 df, p=<2e-16
## Score (logrank) test = 225.2 on 5 df, p=<2e-16(1.2.2.2) Fine-Gray method
fgdata.fx <- finegray(Surv(etime, event) ~ ., data = wpps, etype = "Fracture")
fgdata.death <- finegray(Surv(etime, event) ~ ., data = wpps, etype = "Death")
mfit.fg <- survfit(Surv(fgstart, fgstop, fgstatus) ~ sex, data = fgdata.fx)
fgfit.f <- coxph(Surv(fgstart, fgstop, fgstatus) ~ sex + age + fnbmd + falls_n + prior_fx,
data = fgdata.fx, weight= fgwt)
summary(fgfit.f) ## Call:
## coxph(formula = Surv(fgstart, fgstop, fgstatus) ~ sex + age +
## fnbmd + falls_n + prior_fx, data = fgdata.fx, weights = fgwt)
##
## n= 49711, number of events= 505
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## sex 0.532712 1.703547 0.110046 0.104024 5.121 3.04e-07 ***
## age 0.011767 1.011836 0.007528 0.006885 1.709 0.08744 .
## fnbmd 0.365037 1.440567 0.044831 0.040265 9.066 < 2e-16 ***
## falls_n 0.103195 1.108708 0.074308 0.071336 1.447 0.14801
## prior_fx 0.338834 1.403311 0.122979 0.117573 2.882 0.00395 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 1.704 0.5870 1.3893 2.089
## age 1.012 0.9883 0.9983 1.026
## fnbmd 1.441 0.6942 1.3313 1.559
## falls_n 1.109 0.9020 0.9640 1.275
## prior_fx 1.403 0.7126 1.1145 1.767
##
## Concordance= 0.683 (se = 0.012 )
## Likelihood ratio test= 176.2 on 5 df, p=<2e-16
## Wald test = 224.6 on 5 df, p=<2e-16
## Score (logrank) test = 177.7 on 5 df, p=<2e-16, Robust = 165.7 p=<2e-16
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).zph.fgfit.f <- cox.zph(fgfit.f)
zph.fgfit.f## chisq df p
## sex 1.081 1 0.29855
## age 11.796 1 0.00059
## fnbmd 4.834 1 0.02790
## falls_n 0.453 1 0.50082
## prior_fx 1.563 1 0.21121
## GLOBAL 14.964 5 0.01052windows(width = 14, height = 12)
par(mfrow=c(3,2), oma=c(0,0,2,0))
plot(zph.fgfit.f[1])
abline(h=coef(fgfit.f)[1], lty=2, col=2)
plot(zph.fgfit.f[2])
abline(h=coef(fgfit.f)[2], lty=2, col=2)
plot(zph.fgfit.f[3])
abline(h=coef(fgfit.f)[3], lty=2, col=2)
plot(zph.fgfit.f[4])
abline(h=coef(fgfit.f)[4], lty=2, col=2)
plot(zph.fgfit.f[5])
abline(h=coef(fgfit.f)[5], lty=2, col=2)
mtext("Proportional hazards assumption: Fine-Gray model", line=0, side=3, outer=TRUE, cex=1.8)
(1.2.2.3) Multi-state model
msm.data <- tmerge(wpps, wpps, id = StudyID, death_sen = event(kmfu_death, death_sen),
Fx1st_sen = event(kmfu_fx, Fx1st_sen))
msm.data <- tmerge(msm.data, msm.data, StudyID, enum = cumtdc(tstart))
head(msm.data)## Fx1st_sen BMI_0 neuro_0 rheum_0 hypert_0 diabetes_0 respi_0 cancer_0 cvd_0
## 1 0 27.34 0 0 0 1 0 0 1
## 2 0 27.34 0 0 1 0 0 0 1
## 3 0 24.06 1 0 1 0 0 0 0
## 4 0 26.13 0 0 1 0 0 0 1
## 5 0 27.34 0 0 1 1 0 0 1
## 6 0 26.13 0 0 0 0 0 0 1
## prior_fx age StudyID sex falls_n death_sen fnbmd kmfu_fx kmfu_death
## 1 0 73.07 3 0 0 1 -0.3250 0.5475702 0.5475702
## 2 0 68.47 9 0 0 1 0.2250 5.0513347 5.0513347
## 3 0 62.10 10 1 0 1 1.3417 26.5735797 26.5735797
## 4 0 75.81 24 1 0 0 2.1417 2.8555784 2.8555784
## 5 0 63.53 26 0 0 0 0.2833 9.2621492 9.2621492
## 6 0 76.24 28 1 1 1 2.8833 20.5338809 20.5338809
## group tstart tstop enum
## 1 2 0 0.5475702 1
## 2 2 0 5.0513347 1
## 3 2 0 26.5735797 1
## 4 0 0 2.8555784 1
## 5 0 0 9.2621492 1
## 6 2 0 20.5338809 1with(msm.data, table(death_sen, Fx1st_sen))## Fx1st_sen
## death_sen 0 1
## 0 1468 436
## 1 353 0temp <- with(msm.data, ifelse(death_sen==1, 2, Fx1st_sen))
msm.data$event <- factor(temp, 0:2, labels=c("Event-free", "Fracture", "Death"))
mfit.msm1 <- survfit(Surv(tstart, tstop, event) ~ sex, data = msm.data, id = StudyID)
mfit.msm1## Call: survfit(formula = Surv(tstart, tstop, event) ~ sex, data = msm.data,
## id = StudyID)
##
## n nevent rmean*
## sex=0, (s0) 760 0 16.584865
## sex=1, (s0) 1497 0 15.877738
## sex=0, Fracture 760 94 2.634750
## sex=1, Fracture 1497 342 7.172635
## sex=0, Death 760 176 9.913855
## sex=1, Death 1497 177 6.083096
## *restricted mean time in state (max time = 29.13347 )mfit.msm1$transitions## to
## from Fracture Death (censored)
## (s0) 436 271 1114
## Fracture 0 82 354
## Death 0 0 0temp2 <- with(msm.data, ifelse(enum==2 & event=='Death', 4, as.numeric(event)))
table(temp2)## temp2
## 1 2 3 4
## 1468 436 271 82temp3 <- factor(temp2, labels=c("Event-free", "Fracture", "Death w/o fracture", "Death post fracture"))
table(temp3)## temp3
## Event-free Fracture Death w/o fracture Death post fracture
## 1468 436 271 82mfit.msm2 <- survfit(Surv(tstart, tstop, temp3) ~ sex, data = msm.data, id = StudyID)
print(mfit.msm2)## Call: survfit(formula = Surv(tstart, tstop, temp3) ~ sex, data = msm.data,
## id = StudyID)
##
## n nevent rmean*
## sex=0, (s0) 760 0 16.584865
## sex=1, (s0) 1497 0 15.877738
## sex=0, Fracture 760 94 2.634750
## sex=1, Fracture 1497 342 7.172635
## sex=0, Death w/o fracture 760 151 7.824149
## sex=1, Death w/o fracture 1497 120 3.587283
## sex=0, Death post fracture 760 25 2.089707
## sex=1, Death post fracture 1497 57 2.495813
## *restricted mean time in state (max time = 29.13347 )msmfit <- coxph(Surv(tstart, tstop, event) ~ sex + age +fnbmd + falls_n + prior_fx,
data = msm.data, id = StudyID)
summary(msmfit)## Call:
## coxph(formula = Surv(tstart, tstop, event) ~ sex + age + fnbmd +
## falls_n + prior_fx, data = msm.data, id = StudyID)
##
## n= 2257, number of events= 789
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## sex_1:2 0.484042 1.622620 0.121971 0.123186 3.929 8.52e-05 ***
## age_1:2 0.030720 1.031197 0.008691 0.009045 3.397 0.000682 ***
## fnbmd_1:2 0.402157 1.495046 0.050721 0.051157 7.861 3.80e-15 ***
## falls_n_1:2 0.025733 1.026067 0.081235 0.086838 0.296 0.766977
## prior_fx_1:2 0.370094 1.447871 0.133482 0.140989 2.625 0.008665 **
## sex_1:3 -0.731544 0.481165 0.131667 0.139269 -5.253 1.50e-07 ***
## age_1:3 0.127881 1.136418 0.010263 0.010593 12.073 < 2e-16 ***
## fnbmd_1:3 0.034614 1.035220 0.061108 0.070760 0.489 0.624722
## falls_n_1:3 -0.348446 0.705784 0.129952 0.128053 -2.721 0.006506 **
## prior_fx_1:3 0.014681 1.014789 0.205793 0.212689 0.069 0.944970
## sex_2:3 -1.170363 0.310254 0.275851 0.267189 -4.380 1.19e-05 ***
## age_2:3 0.157412 1.170478 0.021315 0.020866 7.544 4.56e-14 ***
## fnbmd_2:3 0.343225 1.409485 0.124269 0.126827 2.706 0.006805 **
## falls_n_2:3 -0.051889 0.949435 0.183594 0.188336 -0.276 0.782924
## prior_fx_2:3 0.145217 1.156291 0.300550 0.297771 0.488 0.625775
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex_1:2 1.6226 0.6163 1.2746 2.0657
## age_1:2 1.0312 0.9697 1.0131 1.0496
## fnbmd_1:2 1.4950 0.6689 1.3524 1.6527
## falls_n_1:2 1.0261 0.9746 0.8655 1.2164
## prior_fx_1:2 1.4479 0.6907 1.0983 1.9087
## sex_1:3 0.4812 2.0783 0.3662 0.6322
## age_1:3 1.1364 0.8800 1.1131 1.1603
## fnbmd_1:3 1.0352 0.9660 0.9012 1.1892
## falls_n_1:3 0.7058 1.4169 0.5491 0.9071
## prior_fx_1:3 1.0148 0.9854 0.6689 1.5396
## sex_2:3 0.3103 3.2232 0.1838 0.5238
## age_2:3 1.1705 0.8544 1.1236 1.2193
## fnbmd_2:3 1.4095 0.7095 1.0993 1.8072
## falls_n_2:3 0.9494 1.0533 0.6564 1.3733
## prior_fx_2:3 1.1563 0.8648 0.6451 2.0727
##
## Concordance= 0.711 (se = 0.01 )
## Likelihood ratio test= 459.4 on 15 df, p=<2e-16
## Wald test = 521.4 on 15 df, p=<2e-16
## Score (logrank) test = 501.8 on 15 df, p=<2e-16, Robust = 275.5 p=<2e-16
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).zph.msmfit <- cox.zph(msmfit)
zph.msmfit## chisq df p
## sex_1:2 0.00323 1 0.955
## age_1:2 0.04012 1 0.841
## fnbmd_1:2 0.81804 1 0.366
## falls_n_1:2 1.65837 1 0.198
## prior_fx_1:2 0.09868 1 0.753
## sex_1:3 0.08792 1 0.767
## age_1:3 1.43172 1 0.231
## fnbmd_1:3 0.07019 1 0.791
## falls_n_1:3 1.12477 1 0.289
## prior_fx_1:3 3.14762 1 0.076
## sex_2:3 0.10321 1 0.748
## age_2:3 0.12452 1 0.724
## fnbmd_2:3 3.10354 1 0.078
## falls_n_2:3 0.05893 1 0.808
## prior_fx_2:3 1.43041 1 0.232
## GLOBAL 15.09513 15 0.445windows(width = 14, height = 12)
par(mfrow=c(3,2), oma=c(0,0,2,0))
plot(zph.msmfit[1])
abline(h=coef(msmfit)[1], lty=2, col=2)
plot(zph.msmfit[2])
abline(h=coef(msmfit)[2], lty=2, col=2)
plot(zph.msmfit[3])
abline(h=coef(msmfit)[3], lty=2, col=2)
plot(zph.msmfit[4])
abline(h=coef(msmfit)[4], lty=2, col=2)
plot(zph.msmfit[5])
abline(h=coef(msmfit)[5], lty=2, col=2)
mtext("Proportional hazards assumption: Multistate model", line=0, side=3, outer=TRUE, cex=1.8)
(1.2.2.4) Survival curves for different methods
windows(width = 16, height = 12)
par(mfrow=c(1,2), oma=c(0,0,2,0))
plot(mfit.con.f[2], col = "black", lty = 1, fun = "event", mark.time = FALSE, lwd = 2, conf.int = FALSE,
xlab="Time to event (years)", ylab="Probability of event", ylim = c(0,1))
lines(mfit.con.d[2], fun = "event", mark.time = FALSE, col = "black", lty = 2, lwd = 2, conf.int = FALSE)
lines(mfit.cs[2,2], fun = "event", mark.time = FALSE, col = "red", lty = 1, lwd = 2, conf.int = FALSE)
lines(mfit.fg[2], fun = "event", mark.time = FALSE, col = "blue", lty = 1, lwd = 2, conf.int = FALSE)
lines(mfit.msm2[2,2], fun = "event", mark.time = FALSE, col = "green", lty = 1, lwd = 2, conf.int = FALSE)
lines(mfit.msm2[2,4], fun = "event", mark.time = FALSE, col = "green", lty = 3, lwd = 2, conf.int = FALSE)
title(main= "Women")
plot(mfit.con.f[1], col = "black", lty = 1, fun = "event", mark.time = FALSE, lwd = 2, conf.int = FALSE,
xlab="Time to event (years)", ylab="Probability in State", ylim = c(0,1))
lines(mfit.con.d[1], fun = "event", mark.time = FALSE, col = "black", lty = 2, lwd = 2, conf.int = FALSE)
lines(mfit.cs[1,2], fun = "event", mark.time = FALSE, col = "red", lty = 1, lwd = 2, conf.int = FALSE)
lines(mfit.fg[1], fun = "event", mark.time = FALSE, col = "blue", lty = 1, lwd = 2, conf.int = FALSE)
lines(mfit.msm2[1,2], fun = "event", mark.time = FALSE, col = "green", lty = 1, lwd = 2, conf.int = FALSE)
lines(mfit.msm2[1,4], fun = "event", mark.time = FALSE, col = "green", lty = 3, lwd = 2, conf.int = FALSE)
legend("topleft", c("Fracture (Conventional)", "Death (Conventional)", "Fracture (Cause-specific)",
"Fracture (Fine-Gray)", "Fracture (Multistate)", "Death post fracture (Multistate)"),
col = c("black","black","red","blue","green","green"), lty = c(1,2,1,1,1,3), lwd = 2, bty = 'n', cex = 0.8)
title(main= "Men")
mtext("Difference in fracture incidence by competing risk models", line = 0, side = 3, outer = TRUE, cex = 2)
(1.3.3) Baseline hazards from different methods
mean(wpps$sex); mean(wpps$age); mean(wpps$fnbmd); mean(wpps$falls_n); mean(wpps$prior_fx)## [1] 0.6342669## [1] 69.10359## [1] 1.357313## [1] 0.2674355## [1] 0.1317957subj.base.mean <- expand.grid(sex = 0.63, age = 69.1, fnbmd = 1.36, falls_n = 0.28, prior_fx = 0.13)
subj.base.mean## sex age fnbmd falls_n prior_fx
## 1 0.63 69.1 1.36 0.28 0.13# Conventional #
summary(cfit.con.f)## Call:
## coxph(formula = Surv(kmfu_fx, Fx1st_sen) ~ sex + age + fnbmd +
## falls_n + prior_fx, data = wpps, id = StudyID)
##
## n= 1821, number of events= 505
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex 0.394651 1.483867 0.111300 3.546 0.000391 ***
## age 0.042794 1.043723 0.008007 5.345 9.07e-08 ***
## fnbmd 0.404080 1.497923 0.047510 8.505 < 2e-16 ***
## falls_n 0.009253 1.009296 0.075632 0.122 0.902632
## prior_fx 0.397371 1.487908 0.123750 3.211 0.001322 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 1.484 0.6739 1.1930 1.846
## age 1.044 0.9581 1.0275 1.060
## fnbmd 1.498 0.6676 1.3647 1.644
## falls_n 1.009 0.9908 0.8702 1.171
## prior_fx 1.488 0.6721 1.1675 1.896
##
## Concordance= 0.693 (se = 0.013 )
## Likelihood ratio test= 221.9 on 5 df, p=<2e-16
## Wald test = 228.3 on 5 df, p=<2e-16
## Score (logrank) test = 235.9 on 5 df, p=<2e-16con.fx.basem <- survfit(cfit.con.f, newdata = subj.base.mean)
summary(con.fx.basem, time = c(1,5,10,15,20,25))## Call: survfit(formula = cfit.con.f, newdata = subj.base.mean)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1758 43 0.982 0.00285 0.976 0.987
## 5 1326 145 0.908 0.00697 0.895 0.922
## 10 729 155 0.783 0.01196 0.759 0.806
## 15 313 94 0.642 0.01737 0.609 0.677
## 20 97 50 0.467 0.02654 0.418 0.522
## 25 32 14 0.347 0.03539 0.284 0.423# Cause-specific (not needed as the 'CSC' was used) #
# Fine-Gray #
summary(fgfit.f) ## Call:
## coxph(formula = Surv(fgstart, fgstop, fgstatus) ~ sex + age +
## fnbmd + falls_n + prior_fx, data = fgdata.fx, weights = fgwt)
##
## n= 49711, number of events= 505
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## sex 0.532712 1.703547 0.110046 0.104024 5.121 3.04e-07 ***
## age 0.011767 1.011836 0.007528 0.006885 1.709 0.08744 .
## fnbmd 0.365037 1.440567 0.044831 0.040265 9.066 < 2e-16 ***
## falls_n 0.103195 1.108708 0.074308 0.071336 1.447 0.14801
## prior_fx 0.338834 1.403311 0.122979 0.117573 2.882 0.00395 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 1.704 0.5870 1.3893 2.089
## age 1.012 0.9883 0.9983 1.026
## fnbmd 1.441 0.6942 1.3313 1.559
## falls_n 1.109 0.9020 0.9640 1.275
## prior_fx 1.403 0.7126 1.1145 1.767
##
## Concordance= 0.683 (se = 0.012 )
## Likelihood ratio test= 176.2 on 5 df, p=<2e-16
## Wald test = 224.6 on 5 df, p=<2e-16
## Score (logrank) test = 177.7 on 5 df, p=<2e-16, Robust = 165.7 p=<2e-16
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).fg.fx.basem <- survfit(fgfit.f, newdata = subj.base.mean)
summary(fg.fx.basem, time = c(1,5,10,15,20,25))## Call: survfit(formula = fgfit.f, newdata = subj.base.mean)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1778 43 0.980 0.00301 0.974 0.986
## 5 1418 145 0.907 0.00679 0.894 0.921
## 10 729 155 0.800 0.01067 0.779 0.821
## 15 313 94 0.694 0.01442 0.666 0.723
## 20 162 50 0.579 0.01997 0.541 0.619
## 25 32 14 0.513 0.02462 0.467 0.563# MSM #
summary(msmfit)## Call:
## coxph(formula = Surv(tstart, tstop, event) ~ sex + age + fnbmd +
## falls_n + prior_fx, data = msm.data, id = StudyID)
##
## n= 2257, number of events= 789
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## sex_1:2 0.484042 1.622620 0.121971 0.123186 3.929 8.52e-05 ***
## age_1:2 0.030720 1.031197 0.008691 0.009045 3.397 0.000682 ***
## fnbmd_1:2 0.402157 1.495046 0.050721 0.051157 7.861 3.80e-15 ***
## falls_n_1:2 0.025733 1.026067 0.081235 0.086838 0.296 0.766977
## prior_fx_1:2 0.370094 1.447871 0.133482 0.140989 2.625 0.008665 **
## sex_1:3 -0.731544 0.481165 0.131667 0.139269 -5.253 1.50e-07 ***
## age_1:3 0.127881 1.136418 0.010263 0.010593 12.073 < 2e-16 ***
## fnbmd_1:3 0.034614 1.035220 0.061108 0.070760 0.489 0.624722
## falls_n_1:3 -0.348446 0.705784 0.129952 0.128053 -2.721 0.006506 **
## prior_fx_1:3 0.014681 1.014789 0.205793 0.212689 0.069 0.944970
## sex_2:3 -1.170363 0.310254 0.275851 0.267189 -4.380 1.19e-05 ***
## age_2:3 0.157412 1.170478 0.021315 0.020866 7.544 4.56e-14 ***
## fnbmd_2:3 0.343225 1.409485 0.124269 0.126827 2.706 0.006805 **
## falls_n_2:3 -0.051889 0.949435 0.183594 0.188336 -0.276 0.782924
## prior_fx_2:3 0.145217 1.156291 0.300550 0.297771 0.488 0.625775
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex_1:2 1.6226 0.6163 1.2746 2.0657
## age_1:2 1.0312 0.9697 1.0131 1.0496
## fnbmd_1:2 1.4950 0.6689 1.3524 1.6527
## falls_n_1:2 1.0261 0.9746 0.8655 1.2164
## prior_fx_1:2 1.4479 0.6907 1.0983 1.9087
## sex_1:3 0.4812 2.0783 0.3662 0.6322
## age_1:3 1.1364 0.8800 1.1131 1.1603
## fnbmd_1:3 1.0352 0.9660 0.9012 1.1892
## falls_n_1:3 0.7058 1.4169 0.5491 0.9071
## prior_fx_1:3 1.0148 0.9854 0.6689 1.5396
## sex_2:3 0.3103 3.2232 0.1838 0.5238
## age_2:3 1.1705 0.8544 1.1236 1.2193
## fnbmd_2:3 1.4095 0.7095 1.0993 1.8072
## falls_n_2:3 0.9494 1.0533 0.6564 1.3733
## prior_fx_2:3 1.1563 0.8648 0.6451 2.0727
##
## Concordance= 0.711 (se = 0.01 )
## Likelihood ratio test= 459.4 on 15 df, p=<2e-16
## Wald test = 521.4 on 15 df, p=<2e-16
## Score (logrank) test = 501.8 on 15 df, p=<2e-16, Robust = 275.5 p=<2e-16
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).msm.basem <- survfit(msmfit, newdata = subj.base.mean)
summary(msm.basem, time = c(1,5,10,15,20,25))## Call: survfit(formula = msmfit, newdata = subj.base.mean)
##
## data 1
## time n.risk n.event Pr((s0)) Pr(Fracture) Pr(Death)
## 1 1799 65 0.974 0.0184 0.00735
## 5 1452 236 0.871 0.0839 0.04522
## 10 863 244 0.712 0.1653 0.12300
## 15 414 135 0.547 0.2368 0.21614
## 20 143 74 0.359 0.2787 0.36196
## 25 52 28 0.239 0.1786 0.58211Aim 2 Quantification of predictive performance
wpp_val <- read.csv("C:\\Thach\\Research projects\\MSM and competing risk\\Data_Analysis\\wpp_s40.csv")
head(wpp_val)## Fx1st_sen BMI_0 neuro_0 rheum_0 hypert_0 diabetes_0 respi_0 cancer_0 cvd_0
## 1 1 26.13 0 0 1 0 0 0 0
## 2 0 30.48 0 0 1 1 0 0 0
## 3 0 32.30 0 0 1 0 0 0 0
## 4 0 38.57 0 0 1 0 0 0 0
## 5 0 27.34 0 0 1 0 0 0 1
## 6 1 25.59 1 0 1 0 0 0 0
## prior_fx age StudyID sex falls_n death_sen fnbmd kmfu_fx kmfu_death
## 1 0 67.39 8 1 0 0 0.2833 19.789185 19.876797
## 2 0 64.55 27 1 0 0 1.1500 19.301848 19.301848
## 3 0 64.19 29 1 0 0 0.2083 2.603696 2.603696
## 4 0 60.76 32 1 0 0 0.3333 2.168378 2.168378
## 5 0 62.22 36 0 0 1 -0.9583 17.768652 17.768652
## 6 0 67.67 40 0 0 0 2.4583 12.854209 12.646133
## group fx_1y fx_5y fx_10y fx_15y fx_20y fx_25y
## 1 1 0 0 0 0 1 1
## 2 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0
## 5 2 0 0 0 0 0 0
## 6 1 0 0 0 1 1 1(2.1) Calculating absolute risks from different methods
Conventional, FG and MSM: using baseline hazards and coefficients CS: using a ‘black box’ from CSC command
# Conventional
wpp_val$pifx_con = 0.39465*(wpp_val$sex - 0.63) + 0.04279*(wpp_val$age - 69.1) + 0.40408*(wpp_val$fnbmd - 1.36) + 0.00925*(wpp_val$falls_n - 0.28) + 0.3974*(wpp_val$prior_fx - 0.13);
wpp_val$riskfx_con1y = 1- 0.981673**exp(wpp_val$pifx_con)
wpp_val$riskfx_con5y = 1- 0.908361**exp(wpp_val$pifx_con)
wpp_val$riskfx_con10y = 1- 0.782545**exp(wpp_val$pifx_con)
wpp_val$riskfx_con15y = 1- 0.642094**exp(wpp_val$pifx_con)
wpp_val$riskfx_con20y = 1- 0.466764**exp(wpp_val$pifx_con)
wpp_val$riskfx_con25y = 1- 0.346552**exp(wpp_val$pifx_con)
# Cause-specific
CSC.risk <- predict(CSC.fit, newdata = wpp_val, times = c(1,5,10,15,20,25), cause = "Fracture")
CSC.dat <- data.frame(CSC.risk$absRisk)
# Fine-Gray
wpp_val$pifx_fg = 0.5327*(wpp_val$sex - 0.63) + 0.011767*(wpp_val$age - 69.1) + 0.36504*(wpp_val$fnbmd - 1.36) + 0.1032*(wpp_val$falls_n - 0.28) + 0.3388*(wpp_val$prior_fx - 0.13);
wpp_val$riskfx_fg1y = 1- 0.980345**exp(wpp_val$pifx_fg)
wpp_val$riskfx_fg5y = 1- 0.907297**exp(wpp_val$pifx_fg)
wpp_val$riskfx_fg10y = 1- 0.799797**exp(wpp_val$pifx_fg)
wpp_val$riskfx_fg15y = 1- 0.693884**exp(wpp_val$pifx_fg)
wpp_val$riskfx_fg20y = 1- 0.578737**exp(wpp_val$pifx_fg)
wpp_val$riskfx_fg25y = 1- 0.512648**exp(wpp_val$pifx_fg)
# Multistate
wpp_val$pifx_msm = 0.48404*(wpp_val$sex - 0.63) + 0.03072*(wpp_val$age - 69.1) + 0.40216*(wpp_val$fnbmd - 1.36) + 0.0257*(wpp_val$falls_n - 0.28) + 0.3701*(wpp_val$prior_fx - 0.13);
wpp_val$riskfx_msm1y = 1- (1- 0.0183764)**exp(wpp_val$pifx_msm)
wpp_val$riskfx_msm5y = 1- (1- 0.0839146)**exp(wpp_val$pifx_msm)
wpp_val$riskfx_msm10y = 1- (1- 0.1652715)**exp(wpp_val$pifx_msm)
wpp_val$riskfx_msm15y = 1- (1- 0.2367901)**exp(wpp_val$pifx_msm)
wpp_val$riskfx_msm20y = 1- (1- 0.2786909)**exp(wpp_val$pifx_msm)
wpp_val$riskfx_msm25y = 1- (1- 0.1785696)**exp(wpp_val$pifx_msm) (2.2) Baseline characteristics of the validation cohort
table.1v <- CreateTableOne(vars = Vars, strata = "group", data = wpp_val, factorVars = factorVars)
print(table.1v, nonnormal = c("kmfu_fx","kmfu_death"))## Stratified by group
## 0 1
## n 554 259
## sex = 1 (%) 345 (62.3) 203 (78.4)
## age (mean (SD)) 68.08 (5.36) 69.81 (6.74)
## fnbmd (mean (SD)) 1.14 (1.17) 1.70 (1.06)
## falls_n (%)
## 0 457 (82.5) 187 (72.2)
## 1 66 (11.9) 51 (19.7)
## 2 31 ( 5.6) 21 ( 8.1)
## prior_fx = 1 (%) 59 (10.6) 50 (19.3)
## kmfu_fx (median [IQR]) 9.47 [5.68, 14.15] 7.21 [2.85, 12.49]
## kmfu_death (median [IQR]) 9.47 [5.68, 14.15] 12.68 [7.43, 17.15]
## BMI_0 (mean (SD)) 27.13 (4.43) 26.69 (3.57)
## cvd_0 = 1 (%) 175 (31.6) 99 (38.2)
## respi_0 = 1 (%) 52 ( 9.4) 29 (11.2)
## diabetes_0 = 1 (%) 76 (13.7) 25 ( 9.7)
## hypert_0 = 1 (%) 286 (51.6) 134 (51.7)
## cancer_0 = 1 (%) 44 ( 7.9) 28 (10.8)
## rheum_0 = 1 (%) 12 ( 2.2) 10 ( 3.9)
## neuro_0 = 1 (%) 40 ( 7.2) 21 ( 8.1)
## Stratified by group
## 2 p test
## n 122
## sex = 1 (%) 57 (46.7) <0.001
## age (mean (SD)) 72.64 (7.67) <0.001
## fnbmd (mean (SD)) 1.45 (1.48) <0.001
## falls_n (%) 0.004
## 0 102 (83.6)
## 1 10 ( 8.2)
## 2 10 ( 8.2)
## prior_fx = 1 (%) 19 (15.6) 0.003
## kmfu_fx (median [IQR]) 5.93 [2.34, 9.21] <0.001 nonnorm
## kmfu_death (median [IQR]) 5.93 [2.34, 9.21] <0.001 nonnorm
## BMI_0 (mean (SD)) 26.95 (1.40) 0.332
## cvd_0 = 1 (%) 64 (52.5) <0.001
## respi_0 = 1 (%) 16 (13.1) 0.416
## diabetes_0 = 1 (%) 18 (14.8) 0.207
## hypert_0 = 1 (%) 63 (51.6) 1.000
## cancer_0 = 1 (%) 11 ( 9.0) 0.407
## rheum_0 = 1 (%) 7 ( 5.7) 0.085
## neuro_0 = 1 (%) 6 ( 4.9) 0.529(2.3) Assessment of predictive performance
(2.3.1) Prediting 5-year fracture risk by different methods
windows(width = 16, height = 12)
par(mfrow=c(2,2), oma=c(0,0,2,0))
val.prob.ci.2(wpp_val$riskfx_con5y, wpp_val$fx_5y, pl = T, logistic.cal = F, g = 10,
xlab="Predicted risk", ylab="Observed risk", ylim = c(0, 0.45), xlim = c(0,0.45), cex.axis = .8, cex.lab = 1,
statloc = c(-0.01, 0.35), cex = .8)## Call:
## val.prob.ci.2(p = wpp_val$riskfx_con5y, y = wpp_val$fx_5y, pl = T,
## logistic.cal = F, xlab = "Predicted risk", ylab = "Observed risk",
## xlim = c(0, 0.45), ylim = c(0, 0.45), g = 10, statloc = c(-0.01,
## 0.35), cex = 0.8, cex.axis = 0.8, cex.lab = 1)
##
## A 95% confidence interval is given for the calibration intercept, calibration slope and c-statistic.
##
## Dxy C (ROC) R2 D D:Chi-sq
## 4.430243e-01 7.215121e-01 1.198782e-01 6.049030e-02 5.755843e+01
## D:p U U:Chi-sq U:p Q
## 3.286260e-14 -1.794149e-03 3.224707e-01 8.510917e-01 6.228445e-02
## Brier Intercept Slope Emax Brier scaled
## 9.031487e-02 -6.019106e-02 9.835171e-01 2.353082e-02 7.073985e-02
## Eavg ECI
## 2.004853e-02 1.980804e-01title(main= "A. Conventional Cox's PH model", adj = 0, cex.main = 1.4)
val.prob.ci.2(CSC.dat$X2, wpp_val$fx_5y, pl = T, logistic.cal = F, g = 10,
xlab="Predicted risk", ylab="Observed risk", ylim = c(0, 0.45), xlim = c(0,0.45), cex.axis = .8, cex.lab = 1,
statloc = c(-0.01, 0.35), cex = .8)## Call:
## val.prob.ci.2(p = CSC.dat$X2, y = wpp_val$fx_5y, pl = T, logistic.cal = F,
## xlab = "Predicted risk", ylab = "Observed risk", xlim = c(0,
## 0.45), ylim = c(0, 0.45), g = 10, statloc = c(-0.01,
## 0.35), cex = 0.8, cex.axis = 0.8, cex.lab = 1)
##
## A 95% confidence interval is given for the calibration intercept, calibration slope and c-statistic.
##
## Dxy C (ROC) R2 D D:Chi-sq
## 4.473319e-01 7.236659e-01 1.218747e-01 6.154834e-02 5.854769e+01
## D:p U U:Chi-sq U:p Q
## 1.987299e-14 -1.455629e-03 6.389870e-01 7.265169e-01 6.300396e-02
## Brier Intercept Slope Emax Brier scaled
## 9.024284e-02 5.592022e-02 1.090645e+00 6.288579e-02 7.148100e-02
## Eavg ECI
## 1.621784e-02 9.848144e-02title(main= "B. Cause-specific model", adj = 0, cex.main = 1.4)
val.prob.ci.2(wpp_val$riskfx_fg5y, wpp_val$fx_5y, pl = T, logistic.cal = F, g = 10,
xlab="Predicted risk", ylab="Observed risk", ylim = c(0, 0.45), xlim = c(0,0.45), cex.axis = .8, cex.lab = 1,
statloc = c(-0.01, 0.35), cex = .8)## Call:
## val.prob.ci.2(p = wpp_val$riskfx_fg5y, y = wpp_val$fx_5y, pl = T,
## logistic.cal = F, xlab = "Predicted risk", ylab = "Observed risk",
## xlim = c(0, 0.45), ylim = c(0, 0.45), g = 10, statloc = c(-0.01,
## 0.35), cex = 0.8, cex.axis = 0.8, cex.lab = 1)
##
## A 95% confidence interval is given for the calibration intercept, calibration slope and c-statistic.
##
## Dxy C (ROC) R2 D D:Chi-sq
## 4.352565e-01 7.176282e-01 1.140713e-01 5.741938e-02 5.468712e+01
## D:p U U:Chi-sq U:p Q
## 1.413314e-13 -2.580716e-04 1.758703e+00 4.150520e-01 5.767745e-02
## Brier Intercept Slope Emax Brier scaled
## 9.087069e-02 3.526537e-03 1.229851e+00 1.203697e-01 6.502097e-02
## Eavg ECI
## 1.679280e-02 6.927981e-02title(main= "C. Fine-Gray model", adj = 0, cex.main = 1.4)
val.prob.ci.2(wpp_val$riskfx_msm5y, wpp_val$fx_5y, pl = T, logistic.cal = F, g = 10,
xlab="Predicted risk", ylab="Observed risk", ylim = c(0, 0.45), xlim = c(0,0.45), cex.axis = .8, cex.lab = 1,
statloc = c(-0.01, 0.35), cex = .8, legendloc = c(0.2,0.1), cex.leg = .8)## Call:
## val.prob.ci.2(p = wpp_val$riskfx_msm5y, y = wpp_val$fx_5y, pl = T,
## logistic.cal = F, xlab = "Predicted risk", ylab = "Observed risk",
## xlim = c(0, 0.45), ylim = c(0, 0.45), g = 10, legendloc = c(0.2,
## 0.1), statloc = c(-0.01, 0.35), cex = 0.8, cex.leg = 0.8,
## cex.axis = 0.8, cex.lab = 1)
##
## A 95% confidence interval is given for the calibration intercept, calibration slope and c-statistic.
##
## Dxy C (ROC) R2 D D:Chi-sq
## 4.459078e-01 7.229539e-01 1.203610e-01 6.074603e-02 5.779754e+01
## D:p U U:Chi-sq U:p Q
## 2.908784e-14 -1.451570e-03 6.427818e-01 7.251398e-01 6.219760e-02
## Brier Intercept Slope Emax Brier scaled
## 9.032284e-02 6.503964e-02 1.079130e+00 5.894907e-02 7.065784e-02
## Eavg ECI
## 1.712531e-02 1.148505e-01title(main= "D. Multistate model", adj = 0, cex.main = 1.4)
mtext("Prediction of 5-year fracture risk", line = 0, side = 3, outer = TRUE, cex = 2)
(B2.3.2) Prediting 10-year fracture risk by different methods
windows(width = 16, height = 12)
par(mfrow=c(2,2), oma=c(0,0,2,0))
val.prob.ci.2(wpp_val$riskfx_con10y, wpp_val$fx_10y, pl = T, logistic.cal = F, g = 10,
xlab="Predicted risk", ylab="Observed risk", ylim = c(0, 0.8), xlim = c(0,0.8), cex.axis = .8, cex.lab = 1,
statloc = c(-0.02, 0.62), cex = .8, legendloc = c(0,0))## Call:
## val.prob.ci.2(p = wpp_val$riskfx_con10y, y = wpp_val$fx_10y,
## pl = T, logistic.cal = F, xlab = "Predicted risk", ylab = "Observed risk",
## xlim = c(0, 0.8), ylim = c(0, 0.8), g = 10, legendloc = c(0,
## 0), statloc = c(-0.02, 0.62), cex = 0.8, cex.axis = 0.8,
## cex.lab = 1)
##
## A 95% confidence interval is given for the calibration intercept, calibration slope and c-statistic.
##
## Dxy C (ROC) R2 D D:Chi-sq
## 3.816175e-01 6.908087e-01 1.078817e-01 6.671796e-02 6.338130e+01
## D:p U U:Chi-sq U:p Q
## 1.665335e-15 4.126414e-02 4.058197e+01 1.540767e-09 2.545382e-02
## Brier Intercept Slope Emax Brier scaled
## 1.426678e-01 -5.153542e-01 7.517738e-01 2.278306e-01 2.295429e-02
## Eavg ECI
## 8.080606e-02 1.026543e+00title(main= "A. Conventional Cox's PH model", adj = 0, cex.main = 1.4)
val.prob.ci.2(CSC.dat$X3, wpp_val$fx_10y, pl = T, logistic.cal = F, g = 10,
xlab="Predicted risk", ylab="Observed risk", ylim = c(0, 0.8), xlim = c(0,0.8), cex.axis = .8, cex.lab = 1,
statloc = c(-0.02, 0.62), cex = .8, legendloc = c(0,0))## Call:
## val.prob.ci.2(p = CSC.dat$X3, y = wpp_val$fx_10y, pl = T, logistic.cal = F,
## xlab = "Predicted risk", ylab = "Observed risk", xlim = c(0,
## 0.8), ylim = c(0, 0.8), g = 10, legendloc = c(0, 0),
## statloc = c(-0.02, 0.62), cex = 0.8, cex.axis = 0.8, cex.lab = 1)
##
## A 95% confidence interval is given for the calibration intercept, calibration slope and c-statistic.
##
## Dxy C (ROC) R2 D D:Chi-sq
## 3.845003e-01 6.922501e-01 1.109204e-01 6.869552e-02 6.523031e+01
## D:p U U:Chi-sq U:p Q
## 6.661338e-16 -1.009903e-03 1.055741e+00 5.898598e-01 6.970542e-02
## Brier Intercept Slope Emax Brier scaled
## 1.348481e-01 -8.904794e-02 1.024066e+00 1.630391e-02 7.650677e-02
## Eavg ECI
## 2.138634e-02 9.658921e-02title(main= "B. Cause-specific model", adj = 0, cex.main = 1.4)
val.prob.ci.2(wpp_val$riskfx_fg10y, wpp_val$fx_10y, pl = T, logistic.cal = F, g = 10,
xlab="Predicted risk", ylab="Observed risk", ylim = c(0, 0.8), xlim = c(0,0.8), cex.axis = .8, cex.lab = 1,
statloc = c(-0.02, 0.62), cex = .8, legendloc = c(0,0)) ## Call:
## val.prob.ci.2(p = wpp_val$riskfx_fg10y, y = wpp_val$fx_10y, pl = T,
## logistic.cal = F, xlab = "Predicted risk", ylab = "Observed risk",
## xlim = c(0, 0.8), ylim = c(0, 0.8), g = 10, legendloc = c(0,
## 0), statloc = c(-0.02, 0.62), cex = 0.8, cex.axis = 0.8,
## cex.lab = 1)
##
## A 95% confidence interval is given for the calibration intercept, calibration slope and c-statistic.
##
## Dxy C (ROC) R2 D D:Chi-sq
## 3.816175e-01 6.908087e-01 1.094833e-01 6.775979e-02 6.435541e+01
## D:p U U:Chi-sq U:p Q
## 9.992007e-16 1.345834e-02 1.458355e+01 6.811177e-04 5.430145e-02
## Brier Intercept Slope Emax Brier scaled
## 1.372787e-01 -3.294245e-01 9.993744e-01 8.236295e-02 5.986099e-02
## Eavg ECI
## 5.656463e-02 4.301090e-01title(main= "C. Fine-Gray model", adj = 0, cex.main = 1.4)
val.prob.ci.2(wpp_val$riskfx_msm10y, wpp_val$fx_10y, pl = T, logistic.cal = F, g = 10,
xlab="Predicted risk", ylab="Observed risk", ylim = c(0, 0.8), xlim = c(0,0.8), cex.axis = .8, cex.lab = 1,
statloc = c(-0.02, 0.62), cex = .8, legendloc = c(0.45,0.2), cex.leg = .8)## Call:
## val.prob.ci.2(p = wpp_val$riskfx_msm10y, y = wpp_val$fx_10y,
## pl = T, logistic.cal = F, xlab = "Predicted risk", ylab = "Observed risk",
## xlim = c(0, 0.8), ylim = c(0, 0.8), g = 10, legendloc = c(0.45,
## 0.2), statloc = c(-0.02, 0.62), cex = 0.8, cex.leg = 0.8,
## cex.axis = 0.8, cex.lab = 1)
##
## A 95% confidence interval is given for the calibration intercept, calibration slope and c-statistic.
##
## Dxy C (ROC) R2 D D:Chi-sq
## 3.857458e-01 6.928729e-01 1.114422e-01 6.903546e-02 6.554816e+01
## D:p U U:Chi-sq U:p Q
## 5.551115e-16 1.192882e-03 3.115345e+00 2.106258e-01 6.784258e-02
## Brier Intercept Slope Emax Brier scaled
## 1.350592e-01 -1.286991e-01 8.827999e-01 8.454544e-02 7.506117e-02
## Eavg ECI
## 2.930845e-02 2.239376e-01title(main= "D. Multistate model", adj = 0, cex.main = 1.4)
mtext("Prediction of 10-year fracture risk", line = 0, side = 3, outer = TRUE, cex = 2)