confint.ate | R Documentation |
Confidence intervals and confidence Bands for the predicted absolute risk (cumulative incidence function).
## S3 method for class 'ate'
confint(
object,
parm = NULL,
level = 0.95,
n.sim = 10000,
estimator = object$estimator,
contrasts = object$contrasts,
allContrasts = object$allContrasts,
meanRisk.transform = "none",
diffRisk.transform = "none",
ratioRisk.transform = "none",
seed = NA,
ci = object$inference$se,
band = object$inference$band,
p.value = TRUE,
method.band = "maxT-simulation",
alternative = "two.sided",
bootci.method = "perc",
...
)
object |
A |
parm |
not used. For compatibility with the generic method. |
level |
[numeric, 0-1] Level of confidence. |
n.sim |
[integer, >0] the number of simulations used to compute the quantiles for the confidence bands and/or perform adjustment for multiple comparisons. |
estimator |
[character] The type of estimator relative to which the estimates should be displayed. |
contrasts |
[character vector] levels of the treatment variable for which the risks should be assessed and compared. Default is to consider all levels. |
allContrasts |
[2-row character matrix] levels of the treatment variable to be compared. Default is to consider all pairwise comparisons. |
meanRisk.transform |
[character] the transformation used to improve coverage
of the confidence intervals for the mean risk in small samples.
Can be |
diffRisk.transform |
[character] the transformation used to improve coverage
of the confidence intervals for the risk difference in small samples.
Can be |
ratioRisk.transform |
[character] the transformation used to improve coverage
of the confidence intervals for the risk ratio in small samples.
Can be |
seed |
[integer, >0] seed number set when performing simulation for the confidence bands. If not given or NA no seed is set. |
ci |
[logical] should the confidence intervals be computed? |
band |
[logical] should the confidence bands be computed? |
p.value |
[logical] should the p-values/adjusted p-values be computed?
Requires argument |
method.band |
[character] method used to adjust for multiple comparisons.
Can be any element of |
alternative |
[character] a character string specifying the alternative hypothesis,
must be one of |
bootci.method |
[character] Method for constructing bootstrap confidence intervals. Either "perc" (the default), "norm", "basic", "stud", or "bca". |
... |
not used. |
Argument ci
, band
, p.value
, method.band
, alternative
, meanRisk.transform
, diffRisk.transform
, ratioRisk.transform
are only active when the ate
object contains the influence function.
Argument bootci.method
is only active when the ate
object contains bootstrap samples.
Influence function: confidence bands and confidence intervals computed via the influence function are automatically restricted to the interval of definition of the parameter (e.g. [0;1] for the average risk).
Single step max adjustment for multiple comparisons, i.e. accounting for the correlation between the test statistics but not for the ordering of the tests, can be performed setting the arguemnt method.band
to "maxT-integration"
or "maxT-simulation"
. The former uses numerical integration (pmvnorm
and qmvnorm
to perform the adjustment while the latter using simulation. Both assume that the test statistics are jointly normally distributed.
Bootstrap: confidence intervals obtained via bootstrap are computed
using the boot.ci
function of the boot
package.
p-value are obtained using test inversion method
(finding the smallest confidence level such that the interval contain the null hypothesis).
Brice Ozenne
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)))]
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.