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

# Load the RHC data set
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.00

Note 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) +

# 2. fit the propensity score model

# 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.rhcoverlap[,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 Treated.means ## 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.00 Step 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 = rhcswang1
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 sum(rhc$OW[Z==0])
## [1] 3551
surv_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:
cat("Untreated Unadjusted survival probability =", survival_probability_0, "\n") 
## Untreated Unadjusted survival probability = 0.53
cat("Treated Unadjusted survival probability =", survival_probability_1, "\n") 
## Treated Unadjusted survival probability = 0.46
cat("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:
cat("Untreated IPW Adjusted survival probability =", survival_probability_0, "\n") 
## Untreated IPW Adjusted survival probability = 0.52
cat("Treated IPW Adjusted survival probability =", survival_probability_1, "\n") 
## Treated IPW Adjusted survival probability = 0.48
cat("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:
cat("Untreated OW Adjusted survival probability =", survival_probability_0, "\n") 
## Untreated OW Adjusted survival probability = 0.52
cat("Treated OW Adjusted survival probability =", survival_probability_1, "\n") 
## Treated OW Adjusted survival probability = 0.47
cat("Survival Probability Causal Effect at time 180 =", survival_probability_1-survival_probability_0)
## Survival Probability Causal Effect at time 180 = -0.054

Step 3: Weighted estimation of the Marginal HR

Notice that Coxph gives only one standard error estimate when there are no weights.

coxph(Surv(survtime,death) ~ swang1,  data = rhc)
## 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= 3722

Notice 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.

coxph(Surv(survtime,death) ~ swang1,  data = rhc, weights = IPW)
## 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= 3722

Look 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.

rhc$IPWraw = Z/PS + (1-Z)/(1-PS) coxph(Surv(survtime,death) ~ swang1, data = rhc, weights = IPWraw) ## 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= 3722 Here 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. coxph(Surv(survtime,death) ~ swang1, data = rhc, weights = OW) ## 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= 3722 Step 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.2

Step 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.

tmax=365
print(surv, print.rmean=TRUE, rmean=tmax)
## 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
print(surv_ipw, print.rmean=TRUE, rmean=tmax)
## 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
print(surv_ow, print.rmean=TRUE, rmean=tmax)
## 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 =  365

Additional 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_limitswang1==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.trtIPW
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.037 The 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.052

Treatment 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