# NOT RUN {
## Using a simulated cohort
data(hrData, package = "hrIPW")
hrIPW(hrData, time = "time", status = "status", exposure = "Trt",
variables = paste("X", 1:9, sep = ""), wtype = "ATE-stab")
# Standard error could be compared with the (robust) Lin's standard error
# which does not take into account the propensity score estimation step:
modT <- glm(Trt ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9, data = hrData, family = "binomial")
probT <- predict(modT, type = "response")
hrData$w <- 1/ifelse(hrData$Trt == 1, probT, 1 - probT)
library(survival)
coxph(Surv(time, status) ~ Trt + cluster(id), data = hrData, method = "breslow", weights = w)
# or with the bootstrap-based standard-error (see Austin 2016):
# }
# NOT RUN {
f.boot <- function(data, i, wtype) {
df <- data[i, ]
modT <- glm(Trt ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9, data = df, family = "binomial")
probT <- predict(modT, type = "response")
df$w <- 1/ifelse(df$Trt == 1, probT, 1 - probT)
return(coxph(Surv(time, status) ~ Trt, data = df, weights = w)$coef)
}
library(boot); set.seed(1234)
rcoefs <- boot(data = hrData, statistic = f.boot, R = 500)$t
sd(rcoefs)
# }
# NOT RUN {
## Using the DIVAT data base (package IPWsurvival, to be installed)
data(DIVAT, package = "IPWsurvival")
hrIPW(data = DIVAT, time = "times", status = "failures", exposure = "ecd",
variables = c("age", "hla", "retransplant"), wtype = "ATE-unstab")
# Standard error could be compared with the (robust) Lin's standard error
# which does not take into account the propensity score estimation step:
modT <- glm(ecd ~ age + hla + retransplant, data = DIVAT, family = "binomial")
probT <- predict(modT, type = "response")
DIVAT$w <- 1/ifelse(DIVAT$ecd == 1, probT, 1 - probT)
DIVAT$id <- 1:nrow(DIVAT)
coxph(Surv(times, failures) ~ ecd + cluster(id), data = DIVAT, method = "breslow", weights = w)
# or with the bootstrap-based estimation (see Austin 2016):
# }
# NOT RUN {
f.boot2 <- function(data, i, wtype) {
df <- data[i, ]
modT <- glm(ecd ~ age + hla + retransplant, data = df, family = "binomial")
probT <- predict(modT, type = "response")
df$w <- 1/ifelse(df$ecd == 1, probT, 1 - probT)
return(coxph(Surv(times, failures) ~ ecd, data = df, weights = w)$coef)
}
set.seed(1234)
rcoefs2 <- boot(data = DIVAT, statistic = f.boot2, R = 500)$t
sd(rcoefs2)
# }
Run the code above in your browser using DataLab