Analisi dei dati VALIGA per età

Cut-off di età usando gli alberi

Per valutare il cut-off di età direttamente con gli alberi di classificazione ho costruito dei modelli con tutta la coorte VALIGA, usando come variabile solo l’età contro l’evento combinato. Usando diverse configurazioni gli alberi costruiti sono sempre uguali e la prima divisione operata dall’algoritmo è sempre a 27 anni, che penso potrebbe essere un valido cut-off alternativo

require(rpart)
## Loading required package: rpart
require(partykit)
## Loading required package: partykit
## Loading required package: grid
require(survival)
## Loading required package: survival
## 
## Attaching package: 'survival'
## 
## The following object is masked _by_ '.GlobalEnv':
## 
##     kidney
# Check classification trees subset by age on entire VALIGA cohort 
fit_ctrl_15 <- rpart.control(minbucket=.15*NROW(valiga_all),cp=0)
fit_ctrl_10 <- rpart.control(minbucket=.10*NROW(valiga_all),cp=0)
fit_ctrl_5 <- rpart.control(minbucket=.05*NROW(valiga_all),cp=0)

fit_15 <- rpart(Surv(time=time_comb120_15y,event=combine120_15y) ~ age, data=valiga_all,control=fit_ctrl_15)
#fit_15$cptable
plot(as.party(fit_15),main="Combined tree: minbucket=15% sample")

fit_10 <- rpart(Surv(time=time_comb120_15y,event=combine120_15y) ~ age, data=valiga_all,control=fit_ctrl_10)
#fit_10$cptable
plot(as.party(prune(fit_10,cp=.0025)),main="Combined tree: minbucket=10% sample")

fit_5 <- rpart(Surv(time=time_comb120_15y,event=combine120_15y) ~ age, data=valiga_all,control=fit_ctrl_5)
#fit_5$cptable
plot(as.party(prune(fit_5,cp=.0026)),main="Combined tree: minbucket=5% sample")



Di seguito l’albero di sopravvivenza (a 15 anni per l’evento combinato) sulla coorte VALIGA con età < 27. L’albero identifica come soggetti a maggior rischio i pazienti con un eGFR iniziale < 53 (ca. 23% surv), dopo quelli con un eGFR maggiore e M=1 e MAP > 86 (ca. 50% surv) e in seguito quelli con eGFR > 53, M=0 e proteinuria > 0.4 (ca. 75% surv). Tra i pazienti con eGFR > 53, quelli con M=1 e MAP < 86 e quelli con M=0 e proteinuria < 0.4 non presentavano eventi a 15 anni.

# 27 years of age cut of classification trees
valiga_27 <- subset(valiga_all,age<27)
#fit_ctrl_27_15 <- rpart.control(minbucket=.15*NROW(valiga_27),cp=0)
#fit_ctrl_27_10 <- rpart.control(minbucket=.10*NROW(valiga_27),cp=0)
fit_ctrl_27_5 <- rpart.control(minbucket=.05*NROW(valiga_27),cp=0)
fit <- rpart(Surv(time=time_comb120_15y,event=combine120_15y) ~ M + E + S + T01 + MAP_adj + prot_init + eGFR_120_init + age, data=valiga_27,control=fit_ctrl_27_5)
fit <- prune(fit,cp=.0269)
plot(as.party(fit),main="Combined tree: minbucket=5% sample, age=27 yrs, 15 yrs Survival")



Alberi di sopravvivenza considerando il cut-off di età a 21 anni

In quest’albero costruito con cut-off a 21 anni, i pazienti individuati come a maggior rischio sono quelli con proteinuria basale > 4 g/die (25% surv), quelli con proteiniria minore e M=1 avevano una sopravvivenza leggermente più alta (55%), tra quelli con M=0 rischiavano leggermente di più quelli com proteinuria iniziale > 0.4 (70%), mentre tra chi aveva una proteinuria < 0.4 non si sono registrati eventi.

fit_ctrl_21_10 <- rpart.control(minbucket=.10*NROW(valiga_21),cp=0)
fit <- rpart(Surv(time=time_comb120_15y,event=combine120_15y) ~ M + E + S + T01 + prot_init + eGFR_120_init + age + MAP_adj, data=valiga_21,control=fit_ctrl_21_10)
#fit <- prune(fit,cp=.0158)
plot(as.party(prune(fit,cp=.0158)),main="Combined tree: minbucket=10% sample, pediatric age=21 yrs, 15 yrs Survival")



Alberi di sopravvivenza considerando il cut-off di età a 30 anni

In qust’albero costruito con cut-off a 30 anni, i pazienti individuati come a maggior rischio sono quelli con eGFR basale > 58 ml/min (25% surv), quelli con eGFR maggiore, S=1 e proteinuria basale > 3.9 g/die avevano una sopravvivenza praticamente identica, mentre qell icon proteinuria inferiore o S=0 presentano sopravvivenze nettamente migliori (rispettivamente 70% e 90%).

valiga_30 <- subset(valiga_all,age<30)
fit_ctrl_30_5 <- rpart.control(minbucket=.05*NROW(valiga_30),cp=0)
fit <- rpart(Surv(time=time_comb120_15y,event=combine120_15y) ~ M + E + S + T01 + MAP_adj + prot_init + eGFR_120_init + age, data=valiga_30,control=fit_ctrl_30_5)
plot(as.party(prune(fit,cp=.015)),main="Combined tree: minbucket=5% sample, age=30 yrs, 15 yrs Survival")



Analisi del rischio per età

In questa prima analisi, aggiustata solo per età e trattamento, possiamo notare che il limite superiore dell’intervallo di confidenza delle cure di rischio incrocia lo 0 intorno ai 23 anni.

require(rms)
require(ggplot2)

ddist <- datadist(valiga_all)
options(datadist="ddist")
cox_all <- cph(Surv(time = time_comb120_10y, event = combine120_10y) ~ rcs(age,6) + any_Immunos, data = valiga_all,x=T,y=T)
p <- Predict(cox_all, age, any_Immunos,np=400)
p$Immunos <- factor(p$any_Immunos,labels = c("not Treated","Treated"))
g <- ggplot(as.data.frame(p),aes(x=age))
g + geom_line(aes(y=yhat,colour=Immunos),size=.5) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=Immunos),alpha=.175) + geom_vline(xintercept = c(23,30,18,21,25)) + annotate("text", x = c(17.25,20.25,22.25,24.25,29.25,rep(50,3)), y = c(rep(1,5),.1,-1.25,-.9), label = c("18","21","23","25","30","Hazard = 1","95%CI crossing 0 at 23 (nT) and 24.5 (T)","Line crossing 0 at 30.6 (unTreat) and 31.2 (Treat)"),size=4) + theme_bw() + ylab("log Relative Hazard") + geom_hline(yintercept = 0) + ggtitle("log Relative Hazard of 15 yrs combined event by patient age and treatment")



Nella seconda analisi, aggiustata per età, trattamento, sesso, eGFR e proteinuria, possiamo notare che in questo caso a 23 anni c’è una chiara inversione nel verso della curva di rischio che raggiunge il primo di 2 “massimi relativi” (l’altro intorno ai 40 anni).

cox_all <- cph(Surv(time = time_comb120_10y, event = combine120_10y) ~ rcs(age,6) + gender + any_Immunos + prot_init + eGFR_120_init, data = valiga_all,x=T,y=T)
p <- Predict(cox_all, age, any_Immunos,np=400)
p$Immunos <- factor(p$any_Immunos,labels = c("not Treated","Treated"))
g <- ggplot(as.data.frame(p),aes(x=age))
g + geom_line(aes(y=yhat,colour=Immunos),size=.5) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=Immunos),alpha=.175) + geom_vline(xintercept = 23) + annotate("text", x = 22.25, y = 1, label = c("23"),size=3.5) + theme_bw() + ylab("log Relative Hazard") + ggtitle("log Relative Hazard (of (15 yrs combined event) by Age and Treatment (adjusted for sex, eGFR and Proteinuria")