Cox regression with Inverse Propensity Treatment Weighting

Question

A normal Cox Regression is as following:

coxph(formula = Surv(time, status) ~ v1 + v2 + v3, data = x)

I've calculated the Inverse Propensity Treatment Weighting (IPTW) scores with the subsequent Propensity Scores.

Propensity scores can be calculated as following:

ps<-glm(treat~v1+v2+v3, family="binomial", data=x)

Weights used for IPTW are calculated as following:

weight <- ifelse (treat==1, 1/(ps), 1/(1-ps))

Every subject in the dataset can be weighted with aforementioned method (every subject does get a specific weight, calculated as above), but I see no place to put the weights in the 'normal' Cox regression formula.

Is there a Cox regression formula wherein we can assess the calculated weights to each subject and what R package or code is being used for these calculations?

Answer 1

You can do like this using the DIVAT dataset from iptwsurvival package:

##Generate ID
DIVAT$ID<- 1:nrow(DIVAT)

We can calculate the IPTW as the average treatment effect instead as the average treatment effect among treated

DIVAT$p.score <- glm(retransplant ~ age + hla, data = DIVAT, 
                 family = "binomial")$fitted.values

DIVAT$ate.weights <- with(DIVAT, retransplant * 1/p.score + (1-retransplant)* 1/(1-p.score))

Than we can perform a cox regression

####COX without weight
coxph(Surv(times, failures)~ retransplant, data=DIVAT)->fit
summary(fit)

Adding weight is quite easy

###COX with weight naive model
coxph(Surv(times, failures)~ retransplant, data=DIVAT, weights = ate.weights)->fit
summary(fit)
###COX with weight and robust estimation
coxph(Surv(times, failures)~ retransplant + cluster(ID), data=DIVAT, weights = ate.weights)->fit
summary(fit)

However, in this way the estimation of standard error is biased (please see Austin, Peter C. "Variance estimation when using inverse probability of treatment weighting (IPTW) with survival analysis." Statistics in medicine 35.30 (2016): 5642-5655.). Austin suggested to rely on bootstrap estimator. However I'm stacked too, since I'm not able to find a way to perform this kind of analyses. If you found any answer please let me know.

Answer 2

Propensity score weighting method (inverse probability weighting method) R was used for the following statistical analysis.

• Load the following R packages:

  library(ipw)
  library(survival)
  

• Estimate propensity score for each ID in your data frame (base_model), based on variables.

  The propensity score is the probability of assignment of treatment in the presence of given covariates (v). As shown in your data,

  PS estimation
  ps_model <- glm(treatment~v1+v2+v3...., family = binomial, data = base_model)
  summary(ps_model)
  
  # view propensity score values
  pscore <- ps_model$fitted.values
  dataframe$propensityScore <- predict(ps_model, type = "response")
  

• Calculate weights

        #estimate weight for each patient
        base_model$weight.ATE <- ifelse((base_model$treatment=="1"),(1/base_model$propensityScore), (1/(1-base_model$propensityScore)))
  
        base_weight <- ipwpoint(exposure = treatment, family = "binomial", link="logit", numerator = ~1,  denominator =~v1+v2+v3....vn, data = base_model, trunc=0.05) #truncation of 5% for few extreme weights if needed
  

• Survival analysis: Cox regression

        #time to event analysis with weights
        HR5 <- coxph(Surv(time, event)~as.factor(treat_group), weights = weights.trunc, data = base_model) 
        summary(HR5)
  

weights argument was added based on the estimated weights earlier.

• cobalt or tableOne packages of R would help you view balance in characteristics before and after propensity score weighting.

Good luck!