CODEWITHSUNDEEP
r
Carolin V
Carolin V

Reputation: 51

Kaplan Meier curve after IPTW

I want to use IPTW to find the effects of a medication on cardiovascular death (1), which competes with non-cardiovascular death (2) and survival (0). After the IPTW, I would like to do a competing risk analysis to find the effect of the medication on cardiovascular death and plot the resulting Kaplan Meier curve.

This is my start

library(tableone)
library(crr)
library(ipw)
library(sandwich)
library(survey)
treatment<-as.numeric(df$treatment==1)
#propensity score model
psmodel<-glm(treatment ~ age + sex, data=df )

ps<-predict(psmodel, type="response")


#weights
weight<-ifelse(treatment==1,1/(ps),1/(1-ps))


age<-as.numeric(df$age)
sex<-as.numeric(df$sex==1)

cov1 <-cbind(age,sex, weight, treatment)
ftime <-df$Survival
fstatus<-df$Outcome

competingrisk<-crr(ftime, fstatus, cov1,  failcode=2)
summary(competingrisk)

My dataset looks like this but with 500 lines.

structure(list(Outcome = c(1, 1, 1, 2, 2, 2, 2, 0, 0, 1, 1, 0, 
0, 1, 1, 0, 0, 2), Survival = c(7, 13, 14, 8, 9, 15, 14, 16, 
14, 3, 7, 13, 14, 8, 9, 15, 16, 4), treatment = c(1, 0, 0, 1, 
0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0), age = c(59, 58, 57, 
56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42), 
    BMI = c(25, 24, 23, 22, 21, 20, 29, 28, 27, 26, 25, 24, 25, 
    24, 23, 22, 21, 20), sex = c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 
    0, 1, 0, 1, 0, 1, 0, 1)), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -18L))

Upvotes: 0

Views: 1023

Answers (1)

Michael Barrowman
Michael Barrowman

Reputation: 1171

A Kaplan-Meier plot is very naive to multiple covariates, and to competing risks. I highly doubt it is the best way to analyse your data in this instance, however it could be used to investigate and so here is the solution.

If you want to find the four Kaplan-Meier curves for the two events, you can use the following code:

library(survival)

#create a KM model
mod1 <- survfit(Surv(Survival,Outcome==1)~treatment,data=df)

#plot it
plot(mod1,
     conf.int=TRUE, #show confidence intervals
     col=c("red","green") #treatment=0 is red
     )

#see ?plot.survfit for more parameters

mod2 <- survfit(Surv(Survival,Outcome==2)~treatment,data=df)

#plot it
plot(mod2,
     conf.int=TRUE, #show confidence intervals
     col=c("red","green") #treatment=0 is red
)

However, as stated above, KM plots are naive to competing risks (the probabilities do not take the other event into consideration and so at certain time points, you can have that the probabilities of the two events sum to more than 100%).You would likely be better plotting the Cumulative Incidence Function (or CIF), which does take into account the competing risks. The cmprsk package has the cuminc() function for this.

library(cmprsk)
cif <- cuminc(df$Survival,df$Outcome,df$treatment)

plot(cif)
# See ?plot.cuminc for more parameters

However, I find these base plots to be quite unappealing, and so would highly recommend using the survminer package which utilises ggplot2 to create better plots:

library(survminer)

ggsurvplot(mod1)
ggsurvplot(mod2)

ggcompetingrisks(cif)

If you wish to analyse the data and find the effects of a treatment, using the IPTW, you can use the crr() function as in your question and this will return the subdistribution proportional hazard regression coefficients (i.e. the results of a Fine & Gray model). These can be interpreted in much the same way that of a Cox proportional hazard analysis can be (whilst accounting for the competing risks). Bear in mind that the IPTW is not a panacea and therefore these results may still be non-causal.

Once you have the Fine & Gray model results from crr(), you can create a plot of the subdistributions across the two treatments by using the following code (assuming you have made the competingrisk object above)

cov2 <- matrix(c(0,1),ncol=1)
colnames(cov2) <- "treatment"

cr_pred <- predict(competingrisk,cov2)

plot(cr_pred)
# see ?plot.predict.crr for parameters

A few good resources are:

Upvotes: 2

Related Questions