Working with time to event data and am weighting cases by the propensity score. I want to examine how methods of variance estimation affect the HR in a cox proportional hazards model. This is the code so far
library("survival")
library("survminer")
library("WeightIt")
library("sandwich")
library("boot")
library("survey")
data("cancer", package = "survival")
weighting <- weightit(trt ~ risk+ laser + eye, data=diabetic, estimand = "ATE", method = "ps")
#using robust =TRUE
fit_normal <- coxph( Surv(time, status) ~ trt,
data = diabetic, weights = weighting$weights, robust=TRUE )
fit_normal
des <- svydesign(ids = ~1, weights = get.w(weighting),
data = diabetic)
fit.survey <- svycoxph(Surv(time,status)~trt,design=des)
coef(fit.survey)
exp(confint(fit.survey))
exp(confint(fit_normal))
How would I go about bootstrapping to estimate standard errors for the hazards ratio?