library(survival)
library(data.table)
## ## generate data ####
set.seed(10)
d <- sampleData(70,outcome="survival")
d[, X1 := paste0("T",rbinom(.N, size = 2, prob = c(0.51)))]
## table(d$X1)
#### stratified Cox model ####
fit <- coxph(Surv(time,event)~X1 + strata(X2) + X6,
data=d, ties="breslow", x = TRUE, y = TRUE)
#### average treatment effect ####
fit.ate <- ate(fit, treatment = "X1", times = 1:3, data = d,
se = TRUE, iid = TRUE, band = TRUE)
summary(fit.ate)
dt.ate <- as.data.table(fit.ate)
## manual calculation of se
dd <- copy(d)
dd$X1 <- rep(factor("T0", levels = paste0("T",0:2)), NROW(dd))
out <- predictCox(fit, newdata = dd, se = TRUE, times = 1:3, average.iid = TRUE)
term1 <- -out$survival.average.iid
term2 <- sweep(1-out$survival, MARGIN = 2, FUN = "-", STATS = colMeans(1-out$survival))
sqrt(colSums((term1 + term2/NROW(d))^2))
## fit.ate$meanRisk[treatment=="T0",se]
## note
out2 <- predictCox(fit, newdata = dd, se = TRUE, times = 1:3, iid = TRUE)
mean(out2$survival.iid[1,1,])
out$survival.average.iid[1,1]
## check confidence intervals (no transformation)
dt.ate[,.(lower = pmax(0,estimate + qnorm(0.025) * se),
lower2 = lower,
upper = estimate + qnorm(0.975) * se,
upper2 = upper)]
## add confidence intervals computed on the log-log scale
## and backtransformed
outCI <- confint(fit.ate,
meanRisk.transform = "loglog", diffRisk.transform = "atanh",
ratioRisk.transform = "log")
summary(outCI, type = "risk", short = TRUE)
dt.ate[type == "meanRisk", newse := se/(estimate*log(estimate))]
dt.ate[type == "meanRisk", .(lower = exp(-exp(log(-log(estimate)) - 1.96 * newse)),
upper = exp(-exp(log(-log(estimate)) + 1.96 * newse)))]
Run the code above in your browser using DataLab