# 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(dataDIVAT2, package = "RISCA")
hrIPW(data = dataDIVAT2, 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 = dataDIVAT2, family = "binomial")
probT <- predict(modT, type = "response")
dataDIVAT2$w <- 1/ifelse(dataDIVAT2$ecd == 1, probT, 1 - probT)
dataDIVAT2$id <- 1:nrow(dataDIVAT2)
coxph(Surv(times, failures) ~ ecd + cluster(id), data = dataDIVAT2, method = "breslow", weights = w)
# or with the bootstrap-based estimation (see Austin 2016):
# }
# NOT RUN {
f.boot <- 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(123)
rcoefs <- boot(data = dataDIVAT2, statistic = f.boot, R = 200)$t
sd(rcoefs)
# }
Run the code above in your browser using DataCamp Workspace