Survival outcome example: RHC data
The right heart catheterization (RHC) dataset is introduced previously in Lab 3. Here, we follow the setting in Lab 3 but treat the outcome as a time-to-event variable subject to right-censoring. This document uses RHC data as an example to demonstrate usage of balancing weights on analyzing causal effect with a time-to-event outcome. Our goal is to estimate \[ \Delta_{w}(t) = S_w^{(1)}(t) - S_w^{(0)}(t), \] where \(S_w^{(z)}\), \(z\in\{0,1\}\), is the counterfactual survival distribution for target population \(w\). Here, \(w\) can be IPW, ATT, and OW, corresponding to the combined population, treated population, and overlap population.
The RHC Data Set
First, let’s load the dataset
# Load the RHC data set
rhc <- read.csv("rhc_demo.csv")[,-c(1)]
rhc$swang1=as.numeric(rhc$swang1=="RHC")
# summary of the treatment and outcome variables
rhc$censor = 1-rhc$death
summary(rhc[,c("swang1","survtime","death","censor")])## swang1 survtime death censor
## Min. :0.00 Min. : 2 Min. :0.00 Min. :0.00
## 1st Qu.:0.00 1st Qu.: 16 1st Qu.:0.00 1st Qu.:0.00
## Median :0.00 Median : 166 Median :1.00 Median :0.00
## Mean :0.38 Mean : 186 Mean :0.65 Mean :0.35
## 3rd Qu.:1.00 3rd Qu.: 232 3rd Qu.:1.00 3rd Qu.:1.00
## Max. :1.00 Max. :1943 Max. :1.00 Max. :1.00Note that the treated group includes patients receiving RHC and the control group includes patients not receiving RHC.
Step 1: Propensity Score Estimation and Balance Check
Firstly, we fit the propensity score model by logistic regression and check the overlap of the data.
# 1. Specify a logistic regression model
PS.formula <- swang1 ~ as.factor(cat1) + as.factor(cat2) + as.factor(ca) + pafi1 +
wtkilo1 + surv2md1 + as.factor(dementhx) + as.factor(gastr) + wblc1 + temp1 +
das2d3pc + age + as.factor(chfhx) + as.factor(sex) + urin1 +
as.factor(hema) + as.factor(chrpulhx) + as.factor(cardiohx) +
as.factor(meta) + as.factor(adld3p)
# 2. fit the propensity score model
PS.model=glm(PS.formula,data=rhc,family=binomial(link="logit"))
# 3. obtain the propensity score
PS = PS.model$fitted.values
# 4. Density plots of estimated propensity scores
bal.rhc <- SumStat(ps.formula = PS.formula, weight = c('IPW','overlap','treated'),
data = rhc)
plot(bal.rhc, type = 'density')## Propensity score for group 0
## Press [enter] to continue## Propensity score for group 1
# 5. Summary or Comparison of the target populations
Untreated.means=cbind(unweighted=round(bal.rhc$unweighted[,1],2), ipw=round(bal.rhc$IPW[,1],2), ow=round(bal.rhc$overlap[,1],2))
Treated.means=cbind(unweighted=round(bal.rhc$unweighted[,2],2), ipw=round(bal.rhc$IPW[,2],2), ow=round(bal.rhc$overlap[,2],2))
Untreated.means## unweighted ipw ow
## as.factor(cat1)CHF 0.07 0.08 0.10
## as.factor(cat1)Coma 0.10 0.07 0.06
## as.factor(cat1)COPD 0.11 0.08 0.04
## as.factor(cat1)MOSF 0.22 0.29 0.32
## as.factor(cat1)Other 0.06 0.05 0.03
## as.factor(cat2)Cirrhosis 0.01 0.01 0.01
## as.factor(cat2)Colon Cancer 0.00 0.00 0.00
## as.factor(cat2)Coma 0.02 0.02 0.01
## as.factor(cat2)Lung Cancer 0.00 0.00 0.00
## as.factor(cat2)MOSF w/Malignancy 0.05 0.04 0.03
## as.factor(cat2)MOSF w/Sepsis 0.11 0.15 0.16
## as.factor(ca)No 0.75 0.76 0.77
## as.factor(ca)Yes 0.18 0.17 0.17
## pafi1 240.63 220.56 209.57
## wtkilo1 65.04 68.01 69.75
## surv2md1 0.61 0.59 0.59
## as.factor(dementhx)1 0.12 0.10 0.08
## as.factor(gastr)Yes 0.15 0.17 0.17
## wblc1 15.26 15.77 15.84
## temp1 37.63 37.60 37.62
## das2d3pc 20.37 20.44 20.64
## age 61.76 61.28 60.93
## as.factor(chfhx)1 0.17 0.18 0.19
## as.factor(sex)Male 0.54 0.55 0.57
## urin1 2049.95 2061.02 2077.39
## as.factor(hema)Yes 0.07 0.06 0.06
## as.factor(chrpulhx)1 0.22 0.19 0.16
## as.factor(cardiohx)1 0.16 0.19 0.20
## as.factor(meta)Yes 0.05 0.05 0.04
## as.factor(adld3p)1 0.06 0.05 0.04
## as.factor(adld3p)2 0.03 0.02 0.02
## as.factor(adld3p)3 0.01 0.01 0.01
## as.factor(adld3p)4 0.01 0.01 0.01
## as.factor(adld3p)5 0.01 0.01 0.01
## as.factor(adld3p)6 0.01 0.01 0.01
## as.factor(adld3p)7 0.01 0.00 0.00## unweighted ipw ow
## as.factor(cat1)CHF 0.10 0.09 0.10
## as.factor(cat1)Coma 0.04 0.07 0.06
## as.factor(cat1)COPD 0.03 0.08 0.04
## as.factor(cat1)MOSF 0.39 0.29 0.32
## as.factor(cat1)Other 0.03 0.05 0.03
## as.factor(cat2)Cirrhosis 0.01 0.01 0.01
## as.factor(cat2)Colon Cancer 0.00 0.00 0.00
## as.factor(cat2)Coma 0.01 0.02 0.01
## as.factor(cat2)Lung Cancer 0.00 0.00 0.00
## as.factor(cat2)MOSF w/Malignancy 0.03 0.04 0.03
## as.factor(cat2)MOSF w/Sepsis 0.19 0.15 0.16
## as.factor(ca)No 0.79 0.75 0.77
## as.factor(ca)Yes 0.15 0.18 0.17
## pafi1 192.43 220.70 209.57
## wtkilo1 72.36 68.07 69.75
## surv2md1 0.57 0.58 0.59
## as.factor(dementhx)1 0.07 0.09 0.08
## as.factor(gastr)Yes 0.19 0.16 0.17
## wblc1 16.27 15.87 15.84
## temp1 37.59 37.56 37.62
## das2d3pc 20.70 20.47 20.64
## age 60.75 61.44 60.93
## as.factor(chfhx)1 0.19 0.18 0.19
## as.factor(sex)Male 0.59 0.56 0.57
## urin1 2056.12 2078.46 2077.39
## as.factor(hema)Yes 0.05 0.06 0.06
## as.factor(chrpulhx)1 0.14 0.19 0.16
## as.factor(cardiohx)1 0.20 0.19 0.20
## as.factor(meta)Yes 0.04 0.04 0.04
## as.factor(adld3p)1 0.03 0.05 0.04
## as.factor(adld3p)2 0.02 0.02 0.02
## as.factor(adld3p)3 0.01 0.01 0.01
## as.factor(adld3p)4 0.00 0.01 0.01
## as.factor(adld3p)5 0.01 0.01 0.01
## as.factor(adld3p)6 0.01 0.01 0.01
## as.factor(adld3p)7 0.00 0.00 0.00Step 2: Weighted estimation of SPCE and Survival Curves
Here are the unadjusted curves without any weighting. For non-parametric curves define the model by a Surv() function but do not use coxph() to fit the model. Just use survfit. Then ggsurvplot makes nice curves.
surv <- survfit(Surv(survtime,death) ~ swang1, data = rhc)
ggsurvplot(surv, data = rhc, conf.int = TRUE,xlim = c(0, 365),
risk.table = TRUE, risk.table.col = "strata", break.time.by = 31,
xlab = "Time in days", ylab = "Survival probability",
title = "Unadjusted survival curves by treatment")
Now we obtained weights. The weights below are normalized and scaled to the original sample size in each group. Scaling to the original sample size allows the number at risk to be on the original scale per the relative weighting. Alternatively you can get the raw number at risk and combine them.
Z = rhc$swang1
rhc$IPW = Z*sum(Z)/PS/sum(Z/PS) + (1-Z)*sum(1-Z)/(1-PS)/sum((1-Z)/(1-PS))
rhc$ATT = Z + (1-Z)*sum(1-Z)*PS/(1-PS)/sum((1-Z)*PS/(1-PS))
rhc$OW = Z*sum(Z)*(1-PS)/sum(Z*(1-PS)) + (1-Z)*sum(1-Z)*(PS)/sum((1-Z)*(PS))
rhc$PS = PS
#shows the scaling
sum(rhc$OW[Z==1])## [1] 2184## [1] 3551surv_ipw <- survfit(Surv(survtime,death) ~ swang1, data = rhc, weights = IPW)
ggsurvplot(surv_ipw, data = rhc, conf.int = TRUE,xlim = c(0, 365),
risk.table = TRUE, risk.table.col = "strata", break.time.by = 31,
xlab = "Time in days", ylab = "Survival probability",
title = "IPW survival curves by treatment")
surv_ow <- survfit(Surv(survtime,death) ~ swang1, data = rhc, weights = OW)
ggsurvplot(surv_ow, data = rhc, conf.int = TRUE,xlim = c(0, 365),
risk.table = TRUE, risk.table.col = "strata", break.time.by = 31,
xlab = "Time in days", ylab = "Survival probability",
title = "OW survival curves by treatment")
Now if we want to pull off the survival probabilities at a particular point in time.
#Extract a particular row from the survfit object
surv_0 = surv$surv[1:658]
surv_1 = surv$surv[659:(659+468-1)]
time_0 = surv$time[1:658]
time_1 = surv$time[659:(659+468-1)]
#in case 30 is not an observed event time we take the immediately prior time
#this is the position of that time in the data
time_position_0 <- which.min(abs(time_0 - 180))
time_position_1 <- which.min(abs(time_1 - 180))
survival_probability_0 <- surv_0[time_position_0]
survival_probability_1 <- surv_1[time_position_1]
# Print results
cat("At time 180:","\n")## At time 180:## Untreated Unadjusted survival probability = 0.53## Treated Unadjusted survival probability = 0.46cat("Survival Probability Causal Effect at time 180 =", survival_probability_1-survival_probability_0)## Survival Probability Causal Effect at time 180 = -0.075#Extract a particular row from the survfit object
surv_0 = surv_ipw$surv[1:658]
surv_1 = surv_ipw$surv[659:(659+468-1)]
time_0 = surv_ipw$time[1:658]
time_1 = surv_ipw$time[659:(659+468-1)]
time_position_0 <- which.min(abs(time_0 - 180))
time_position_1 <- which.min(abs(time_1 - 180))
survival_probability_0 <- surv_0[time_position_0]
survival_probability_1 <- surv_1[time_position_1]
# Print results
cat("At time 180:","\n")## At time 180:## Untreated IPW Adjusted survival probability = 0.52## Treated IPW Adjusted survival probability = 0.48cat("Survival Probability Causal Effect at time 180 =", survival_probability_1-survival_probability_0)## Survival Probability Causal Effect at time 180 = -0.039#Extract a particular row from the survfit object
surv_0 = surv_ow$surv[1:658]
surv_1 = surv_ow$surv[659:(659+468-1)]
time_0 = surv_ow$time[1:658]
time_1 = surv_ow$time[659:(659+468-1)]
time_position_0 <- which.min(abs(time_0 - 180))
time_position_1 <- which.min(abs(time_1 - 180))
survival_probability_0 <- surv_0[time_position_0]
survival_probability_1 <- surv_1[time_position_1]
# Print results
cat("At time 180:","\n")## At time 180:## Untreated OW Adjusted survival probability = 0.52## Treated OW Adjusted survival probability = 0.47cat("Survival Probability Causal Effect at time 180 =", survival_probability_1-survival_probability_0)## Survival Probability Causal Effect at time 180 = -0.054Step 3: Weighted estimation of the Marginal HR
Notice that Coxph gives only one standard error estimate when there are no weights.
## Call:
## coxph(formula = Surv(survtime, death) ~ swang1, data = rhc)
##
## coef exp(coef) se(coef) z p
## swang1 0.15 1.16 0.03 4 6e-06
##
## Likelihood ratio test=20 on 1 df, p=6e-06
## n= 5735, number of events= 3722Notice that Coxph also gives a robust standard error when there are weights. Here, the two variances look pretty similar. We might conclude it doesn’t matter. We would be wrong.
## Call:
## coxph(formula = Surv(survtime, death) ~ swang1, data = rhc, weights = IPW)
##
## coef exp(coef) se(coef) robust se z p
## swang1 0.09 1.09 0.03 0.04 2 0.03
##
## Likelihood ratio test=7 on 1 df, p=0.01
## n= 5735, number of events= 3722Look what happens if my weights have a different scale. These are the raw IPW weights with no normalization and no scaling. The robust standard error is identical (invariant to scale) but the model-based variance is senstive to scale and is now under-estimated. It is too small because the raw IPW weights inflate the sample size by approximately 2*N.
## Call:
## coxph(formula = Surv(survtime, death) ~ swang1, data = rhc, weights = IPWraw)
##
## coef exp(coef) se(coef) robust se z p
## swang1 0.09 1.09 0.02 0.04 2 0.03
##
## Likelihood ratio test=14 on 1 df, p=2e-04
## n= 5735, number of events= 3722Here the results of OW are very similar to IPW, including the standard error. That reflects that the distribution of propensity scores rarely near 0 or 1.
## Call:
## coxph(formula = Surv(survtime, death) ~ swang1, data = rhc, weights = OW)
##
## coef exp(coef) se(coef) robust se z p
## swang1 0.10 1.11 0.03 0.04 3 0.005
##
## Likelihood ratio test=10 on 1 df, p=0.002
## n= 5735, number of events= 3722Step 4: Getting confidence intervals for the Marginal HR
We exemplify this only for IPWraw because it exhibits a major difference between model-based and robust standard errors. We further add the summary() command and robust=TRUE to emphasize that we want robust standard errors to be used in the confidence interval calculation. Notice that if you remove robust=TRUE the confidence intervals don’t change because the function knows that you are using weights. If you set robust=FALSE you will get the anti-conservative confidence interval that is generally wrong.
temp=summary(coxph(Surv(survtime,death) ~ swang1, data = rhc, weights = IPWraw, robust = TRUE))
temp$conf.int## exp(coef) exp(-coef) lower .95 upper .95
## swang1 1.1 0.92 1 1.2Step 5: Weighted estimation of Restricted Mean Survival Time
Set the truncation time tmax (e.g., 365 days) and calculate RMST from the survival function that was already estimated with weights.
## Call: survfit(formula = Surv(survtime, death) ~ swang1, data = rhc)
##
## n events rmean* se(rmean) median 0.95LCL 0.95UCL
## swang1=0 3551 2236 197 2.73 233 209 247
## swang1=1 2184 1486 172 3.50 106 80 136
## * restricted mean with upper limit = 365## Call: survfit(formula = Surv(survtime, death) ~ swang1, data = rhc,
## weights = IPW)
##
## records n events rmean* se(rmean) median 0.95LCL 0.95UCL
## swang1=0 3551 3551 2267 192 2.74 218 184 239
## swang1=1 2184 2184 1488 178 3.46 139 104 194
## * restricted mean with upper limit = 365## Call: survfit(formula = Surv(survtime, death) ~ swang1, data = rhc,
## weights = OW)
##
## records n events rmean* se(rmean) median 0.95LCL 0.95UCL
## swang1=0 3551 3551 2244 194 2.75 225 194 246
## swang1=1 2184 2184 1477 176 3.50 127 98 166
## * restricted mean with upper limit = 365Additional Step of Modeling the Censoring Process (IPCW)
Until now, we have been concerned with confounding, bias that arises due to common causes of treatment and outcome that are almost certainly present in the absence of randomization (we assume so). The preceding code shows how to apply weights to address confounding and estimate various parameters of the survival distribution including hazard ratios, survival curves, and restricted mean survival time. All of the preceding methods assume that censoring C, occurs completely at random (non-informative censoring). This is plausible in many studies where censoring is due to to the end of study follow-up, or differential enrollment times where some individuals enter the study later than others. When censoring is also attributable to patient choices and decision it may not be random. If censoring is thought to related to measured covariates, but not related to outcome conditional on those measured covariates, we can incorporate this type of censoring at random through inverse probability of censoring weights. This section looks at methods that involve both a propensity score for treatment (handling lack of randomization) and a censoring score that handles censoring that is dependent on measured covariates.
We can use the following Cox model to describe the censoring process. \[ P(C^{(z)} \geq t | \mathbf{X}) = \exp\left\{\Lambda_{z} (t) e^{\mathbf{\theta}_z^T\mathbf{X}}\right\}, \] where \(C^{(z)}\) is the potential censoring time under treatment \(z\) (\(z=1,0\) for treated and untreated groups, respectively), \(\Lambda_{z}\) is the treatment-specific baseline cumulative hazard, and \(\mathbf{\theta}_z\) is treatment-specific coefficients corresponding to pre-treatment covariates \(\mathbf{X}\). The partial likelihood approach will be used to estimate coefficients \(\mathbf{\theta}_z\) and the baseline cumulative hazard \(\Lambda_{z}(t)\) will be calculated through the Breslow approach. Because all the parameters are treatment-specific, we will implement two Cox models using the treated and untreated group samples, separately. See the code below.
Step 1: Model the censoring process for IPCW
#0. Excluding rare conditions that cause modeling problems
keep = (rhc$cat2 != "Cirrhosis" & rhc$cat2 != "Colon Cancer" & rhc$cat2 != "Lung Cancer")*1
rhc_limit = rhc[keep==1,]
# 1. Construct the Cox model formula
Censor.formula=Surv(survtime,censor) ~ as.factor(cat1) + as.factor(cat2) +
as.factor(ca) + pafi1 +wtkilo1 + surv2md1 + as.factor(dementhx) +
as.factor(gastr) + wblc1 + temp1 + das2d3pc + age + as.factor(chfhx) +
as.factor(sex) + urin1 + as.factor(hema) + as.factor(chrpulhx) +
as.factor(cardiohx) + as.factor(meta) + as.factor(adld3p)
# 2. Fit censoring model for treated group
rhc.trt=subset(rhc_limit,rhc_limit$swang1==1)
Censor.trt.model = coxph(Censor.formula,data=rhc.trt)
summary(Censor.trt.model)$conf.int## exp(coef) exp(-coef) lower .95 upper .95
## as.factor(cat1)CHF 0.96 1.04 0.647 1.43
## as.factor(cat1)Coma 1.69 0.59 0.874 3.25
## as.factor(cat1)COPD 0.82 1.23 0.459 1.45
## as.factor(cat1)MOSF 1.11 0.90 0.909 1.35
## as.factor(cat1)Other 1.06 0.94 0.543 2.08
## as.factor(cat2)Coma 1.10 0.91 0.433 2.82
## as.factor(cat2)MOSF w/Malignancy 0.30 3.28 0.093 0.99
## as.factor(cat2)MOSF w/Sepsis 0.98 1.03 0.780 1.22
## as.factor(ca)No 1.32 0.76 0.729 2.38
## as.factor(ca)Yes 1.26 0.79 0.681 2.35
## pafi1 1.00 1.00 0.998 1.00
## wtkilo1 1.00 1.00 0.999 1.00
## surv2md1 3.03 0.33 1.429 6.44
## as.factor(dementhx)1 0.85 1.18 0.579 1.24
## as.factor(gastr)Yes 0.94 1.06 0.762 1.17
## wblc1 1.00 1.00 0.989 1.00
## temp1 1.01 0.99 0.961 1.06
## das2d3pc 1.03 0.97 1.016 1.05
## age 0.99 1.01 0.989 1.00
## as.factor(chfhx)1 0.66 1.51 0.495 0.89
## as.factor(sex)Male 1.03 0.97 0.876 1.20
## urin1 1.00 1.00 1.000 1.00
## as.factor(hema)Yes 0.73 1.36 0.409 1.32
## as.factor(chrpulhx)1 0.85 1.18 0.648 1.11
## as.factor(cardiohx)1 0.91 1.10 0.725 1.14
## as.factor(meta)Yes 0.69 1.45 0.449 1.06
## as.factor(adld3p)1 0.99 1.01 0.703 1.41
## as.factor(adld3p)2 1.04 0.96 0.614 1.75
## as.factor(adld3p)3 0.60 1.67 0.221 1.62
## as.factor(adld3p)4 1.23 0.82 0.386 3.89
## as.factor(adld3p)5 0.72 1.40 0.291 1.76
## as.factor(adld3p)6 0.53 1.90 0.129 2.15
## as.factor(adld3p)7 3.02 0.33 0.933 9.77# Fit censoring model for treated group
rhc.con=subset(rhc_limit,rhc_limit$swang1==0)
Censor.con.model = coxph(Censor.formula,data=rhc.con)
summary(Censor.con.model)$conf.int## exp(coef) exp(-coef) lower .95 upper .95
## as.factor(cat1)CHF 0.79 1.26 0.60 1.05
## as.factor(cat1)Coma 1.18 0.85 0.90 1.56
## as.factor(cat1)COPD 0.89 1.12 0.70 1.14
## as.factor(cat1)MOSF 0.85 1.17 0.72 1.01
## as.factor(cat1)Other 0.86 1.16 0.64 1.17
## as.factor(cat2)Coma 1.25 0.80 0.70 2.22
## as.factor(cat2)MOSF w/Malignancy 0.76 1.32 0.45 1.28
## as.factor(cat2)MOSF w/Sepsis 0.81 1.24 0.66 0.99
## as.factor(ca)No 2.11 0.47 1.34 3.31
## as.factor(ca)Yes 1.94 0.52 1.23 3.05
## pafi1 1.00 1.00 1.00 1.00
## wtkilo1 1.00 1.00 1.00 1.00
## surv2md1 1.76 0.57 0.95 3.26
## as.factor(dementhx)1 0.86 1.16 0.70 1.06
## as.factor(gastr)Yes 0.88 1.14 0.73 1.06
## wblc1 1.00 1.00 0.99 1.00
## temp1 1.03 0.97 0.99 1.07
## das2d3pc 1.02 0.98 1.01 1.03
## age 0.99 1.01 0.99 1.00
## as.factor(chfhx)1 0.70 1.42 0.57 0.86
## as.factor(sex)Male 0.82 1.21 0.74 0.92
## urin1 1.00 1.00 1.00 1.00
## as.factor(hema)Yes 0.74 1.34 0.55 1.01
## as.factor(chrpulhx)1 0.90 1.11 0.75 1.08
## as.factor(cardiohx)1 0.93 1.07 0.77 1.13
## as.factor(meta)Yes 0.76 1.32 0.58 0.99
## as.factor(adld3p)1 1.07 0.94 0.87 1.30
## as.factor(adld3p)2 1.13 0.89 0.81 1.57
## as.factor(adld3p)3 1.13 0.89 0.63 2.01
## as.factor(adld3p)4 0.81 1.24 0.47 1.41
## as.factor(adld3p)5 0.83 1.20 0.52 1.34
## as.factor(adld3p)6 1.03 0.97 0.65 1.62
## as.factor(adld3p)7 1.51 0.66 0.77 2.96# 3. Calculate the censoring scores
# censoring scores in the treated group
Kc.trt=CensorScoreWithCox(CoxModel=Censor.trt.model,TimeVec=rhc.trt[,"survtime"])
# censoring scores in the untreated group
Kc.con=CensorScoreWithCox(CoxModel=Censor.con.model,TimeVec=rhc.con[,"survtime"])
# density plot for the censoring scores
df_plot <- data.frame(Group=c(rep("Treated",length(Kc.trt)),rep("Control",length(Kc.con))), Kc=c(Kc.trt,Kc.con))
ggplot(df_plot, aes(x=Kc, fill=Group)) +geom_density(alpha=0.4)+xlab("Censoring scores")
Step 2: IPTW + IPCW estimation of SPCE
We demonstrate how to estimate the treatment effect at 180 days after enrollment to the RHC study. The following code demonstrates how to obtain IPW estimator (\(\hat\Delta_{IPW}(180)\)):
# 1. Calculate the counterfactual survival probability in the treated group
w.trt =rhc.trt$IPW
S1=1-sum(w.trt*rhc.trt[,"death"]*as.numeric(rhc.trt[,"survtime"]<=180)/Kc.trt)/sum(w.trt)
# 2. Calculate the counterfactual survival probability in the untreated group
w.con =rhc.con$IPW
S0=1-sum(w.con*rhc.con[,"death"]*(rhc.con[,"survtime"]<=180)/Kc.con)/sum(w.con)
# 3. Take the difference
Delta.IPW=S1-S0
names(Delta.IPW)="Delta.IPW"
cat(paste("The IPW estimate is",round(Delta.IPW,3)))## The IPW estimate is -0.037The following code demonstrates how to obtain OW estimator (\(\hat\Delta_{OW}(180)\)):
# 1. Calculate the counterfactual survival probability in the treated group
w.trt = rhc.trt$OW
S1=1-sum(w.trt*rhc.trt[,"death"]*as.numeric(rhc.trt[,"survtime"]<=180)/Kc.trt)/sum(w.trt)
# 2. Calculate the counterfactual survival probability in the untreated group
w.con = rhc.con$OW
S0=1-sum(w.con*rhc.con[,"death"]*(rhc.con[,"survtime"]<=180)/Kc.con)/sum(w.con)
# 3. Finally, obtain the treatment effect
Delta.OW=S1-S0
names(Delta.OW)="Delta.OW"
cat(paste("The OW estimate is",round(Delta.OW,3)))## The OW estimate is -0.052Treatment effect curve with IPTW and IPCW
The R code for implementation of all aforementioned balancing weights are summarized in a unified function SurvEffectWithCox in PSweight_survival.R A full tutorial of using SurvEffectWithCox can be found in the Supplementary Material of Cheng et al. (2022) Addressing Extreme Propensity Scores in Estimating Counterfactual Survival Functions via the Overlap Weights. American Journal of Epidemiology, https://doi.org/10.1093/aje/kwac043