Nothing
#' Compute confidence intervals for drtmle and adaptive_iptw@
#' @param ... Arguments to be passed to method
#' @export
ci <- function(...) {
UseMethod("ci")
}
#' Confidence intervals for drtmle objects
#'
#' @param object An object of class \code{"drtmle"}
#' @param est A vector indicating for which estimators to return a
#' confidence interval. Possible estimators include the TMLE with doubly robust
#' inference (\code{"drtmle"}, recommended), the AIPTW with additional
#' correction for misspecification (\code{"aiptw_c"}, not recommended), the
#' standard TMLE (\code{"tmle"}, recommended only for comparison to "drtmle"),
#' the standard AIPTW (\code{"aiptw"}, recommended only for comparison to
#' "drtmle"), and G-computation (\code{"gcomp"}, not recommended).
#' @param level The nominal coverage probability of the desired confidence
#' interval (should be between 0 and 1). Default computes 95\% confidence
#' intervals.
#' @param contrast Specifies the parameter for which to return confidence
#' intervals. If \code{contrast=NULL}, then confidence intervals for the
#' marginal means are computed. If instead, \code{contrast} is a numeric vector
#' of ones, negative ones, and zeros to define linear combinations of the
#' various means (e.g., to estimate an average treatment effect, see example).
#' Finally, \code{contrast} can be a list with named functions \code{f},
#' \code{f_inv}, \code{h}, and \code{fh_grad}. The first two functions should
#' take as input argument \code{eff}. Respectively, these specify which
#' transformation of the effect measure to compute the confidence interval for
#' and the inverse transformation to put the confidence interval back on the
#' original scale. The function \code{h} defines the contrast to be estimated
#' and should take as input \code{est}, a vector of the same length as
#' \code{object$a_0}, and output the desired contrast. The function
#' \code{fh_grad} is the gradient of the function \code{h}. See examples and
#' vignette for more information.
#' @param ... Other options (not currently used).
#' @importFrom stats qnorm
#' @export
#' @method ci drtmle
#' @return An object of class \code{"ci.drtmle"} with point estimates and
#' confidence intervals of the specified level.
#'
#' @examples
#' # load super learner
#' library(SuperLearner)
#' # simulate data
#' set.seed(123456)
#' n <- 100
#' W <- data.frame(W1 = runif(n), W2 = rnorm(n))
#' A <- rbinom(n, 1, plogis(W$W1 - W$W2))
#' Y <- rbinom(n, 1, plogis(W$W1 * W$W2 * A))
#' \donttest{
#' # fit drtmle with maxIter = 1 to run fast
#' fit1 <- drtmle(
#' W = W, A = A, Y = Y, a_0 = c(1, 0),
#' family = binomial(),
#' stratify = FALSE,
#' SL_Q = c("SL.glm", "SL.mean"),
#' SL_g = c("SL.glm", "SL.mean"),
#' SL_Qr = "SL.npreg",
#' SL_gr = "SL.npreg", maxIter = 1
#' )
#'
#' # get confidence intervals for each mean
#' ci_mean <- ci(fit1)
#'
#' # get confidence intervals for ATE
#' ci_ATE <- ci(fit1, contrast = c(1, -1))
#'
#' # get confidence intervals for risk ratio by
#' # computing CI on log scale and back-transforming
#' myContrast <- list(
#' f = function(eff) {
#' log(eff)
#' },
#' f_inv = function(eff) {
#' exp(eff)
#' },
#' h = function(est) {
#' est[1] / est[2]
#' },
#' fh_grad = function(est) {
#' c(1 / est[1], -1 / est[2])
#' }
#' )
#' ci_RR <- ci(fit1, contrast = myContrast)
#' }
ci.drtmle <- function(object, est = c("drtmle"), level = 0.95,
contrast = NULL, ...) {
if (!inherits(object, "drtmle")) {
stop("ci only works with drtmle objects")
}
out <- vector(mode = "list", length = length(est))
names(out) <- est
# if no contrast then return an CI for each
# covariate-adjusted mean
if (is.null(contrast)) {
for (i in seq_along(est)) {
out[[i]] <- matrix(NA, nrow = length(object$a_0), ncol = 3)
for (j in seq_along(object$a_0)) {
out[[i]][j, ] <-
rep(object[[est[i]]]$est[j], 3) +
stats::qnorm(c(0.5, (1 - level) / 2, (1 + level) / 2)) *
rep(sqrt(object[[est[i]]]$cov[j, j]), 3)
}
row.names(out[[i]]) <- object$a_0
colnames(out[[i]]) <- c("est", "cil", "ciu")
}
} else if (is.numeric(contrast)) {
# check that contrast came in correctly
if (length(contrast) != length(object$a_0)) {
stop("length of contrast vector not equal to length of a_0")
}
if (!all(contrast %in% c(-1, 1, 0))) {
stop("contrast should only be -1, 1, or 0")
}
for (i in seq_along(est)) {
out[[i]] <- matrix(NA, nrow = 1, ncol = 3)
g <- matrix(contrast, nrow = length(contrast))
p <- matrix(object[[est[i]]]$est, nrow = length(contrast))
v <- object[[est[i]]]$cov
thisC <- t(g) %*% p
thisSe <- sqrt(t(g) %*% v %*% g)
out[[i]][1, ] <- rep(thisC, 3) +
stats::qnorm(c(0.5, (1 - level) / 2, (1 + level) / 2)) *
rep(thisSe, 3)
# E[Y(a_0[1])] - E[Y(a_0[2])]
indMinus <- which(contrast == -1)
indPlus <- which(contrast == 1)
plusTerms <- minusTerms <- ""
if (length(indMinus) > 0) {
minusTerms <- paste0("E[Y(", object$a_0[indMinus], ")]")
}
if (length(indPlus) > 0) {
plusTerms <- paste0("E[Y(", object$a_0[indPlus], ")]")
}
thisName <- paste0(
paste0(plusTerms, collapse = "+"),
ifelse(length(indMinus) > 0, "-", ""),
paste0(minusTerms, collapse = "-")
)
row.names(out[[i]]) <- thisName
colnames(out[[i]]) <- c("est", "cil", "ciu")
}
} else if (is.list(contrast)) {
if (!all(c("f", "f_inv", "h", "fh_grad") %in% names(contrast))) {
stop("some function missing in contrast. see ?ci.drtmle for help.")
}
for (i in seq_along(est)) {
out[[i]] <- matrix(NA, nrow = 1, ncol = 3)
thisC <- do.call(contrast$h, args = list(est = object[[est[i]]]$est))
f_thisC <- do.call(contrast$f, args = list(eff = thisC))
grad <- matrix(do.call(
contrast$fh_grad,
args = list(est = object[[est[i]]]$est)
), nrow = length(object$a_0))
v <- object[[est[i]]]$cov
thisSe <- sqrt(t(grad) %*% v %*% grad)
transformCI <- rep(f_thisC, 3) +
stats::qnorm(c(0.5, (1 - level) / 2, (1 + level) / 2)) *
rep(thisSe, 3)
out[[i]][1, ] <- do.call(contrast$f_inv, args = list(eff = transformCI))
row.names(out[[i]]) <- c("user contrast")
colnames(out[[i]]) <- c("est", "cil", "ciu")
}
}
class(out) <- "ci.drtmle"
return(out)
}
#' Confidence intervals for adaptive_iptw objects
#'
#' Estimate confidence intervals for objects of class \code{"adaptive_iptw"}
#'
#' @param object An object of class \code{"adaptive_iptw"}
#' @param est A vector indicating for which estimators to return a
#' confidence interval. Possible estimators include the TMLE IPTW
#' (\code{"iptw_tmle"}, recommended), the one-step IPTW
#' (\code{"iptw_os"}, not recommended), the standard IPTW
#' (\code{"iptw"}, recommended only for comparison to the other two estimators).
#' @param level The nominal coverage probability of the desired confidence
#' interval (should be between 0 and 1). Default computes 95\% confidence
#' intervals.
#' @param contrast Specifies the parameter for which to return confidence
#' intervals. If \code{contrast=NULL}, then confidence intervals for the
#' marginal means are computed. If instead, \code{contrast} is a numeric vector
#' of ones, negative ones, and zeros to define linear combinations of the
#' various means (e.g., to estimate an average treatment effect, see example).
#' Finally, \code{contrast} can be a list with named functions \code{f},
#' \code{f_inv}, \code{h}, and \code{fh_grad}. The first two functions should
#' take as input argument \code{eff}. Respectively, these specify which
#' transformation of the effect measure to compute the confidence interval for
#' and the inverse transformation to put the confidence interval back on the
#' original scale. The function \code{h} defines the contrast to be estimated
#' and should take as input \code{est}, a vector of the same length as
#' \code{object$a_0}, and output the desired contrast. The function
#' \code{fh_grad} is the gradient of the function \code{h}. See examples and
#' vignette for more information.
#' @param ... Other options (not currently used).
#' @export
#' @method ci adaptive_iptw
#' @return An object of class \code{"ci.adaptive_iptw"} with point estimates and
#' confidence intervals of the specified level.
#'
#' @examples
#' # load super learner
#' library(SuperLearner)
#' # fit adaptive_iptw
#' set.seed(123456)
#' n <- 200
#' W <- data.frame(W1 = runif(n), W2 = rnorm(n))
#' A <- rbinom(n, 1, plogis(W$W1 - W$W2))
#' Y <- rbinom(n, 1, plogis(W$W1 * W$W2 * A))
#'
#' fit1 <- adaptive_iptw(
#' W = W, A = A, Y = Y, a_0 = c(1, 0),
#' SL_g = c("SL.glm", "SL.mean", "SL.step"),
#' SL_Qr = "SL.glm"
#' )
#'
#' # get confidence intervals for each mean
#' ci_mean <- ci(fit1)
#'
#' # get confidence intervals for ATE
#' ci_ATE <- ci(fit1, contrast = c(1, -1))
#'
#' # get confidence intervals for risk ratio
#' # by inputting own contrast function
#' # this computes CI on log scale and back transforms
#' myContrast <- list(
#' f = function(eff) {
#' log(eff)
#' },
#' f_inv = function(eff) {
#' exp(eff)
#' },
#' h = function(est) {
#' est[1] / est[2]
#' },
#' fh_grad = function(est) {
#' c(1 / est[1], -1 / est[2])
#' }
#' )
#' ci_RR <- ci(fit1, contrast = myContrast)
ci.adaptive_iptw <- function(object, est = c("iptw_tmle"), level = 0.95,
contrast = NULL, ...) {
if (any(est == "iptw")) {
stop("Theory does not support inference for naive IPTW with super learner.")
}
if (!inherits(object, "adaptive_iptw")) {
stop("ci only works with adaptive_iptw objects")
}
out <- vector(mode = "list", length = length(est))
names(out) <- est
# if no contrast then return an CI for each
# covariate-adjusted mean
if (is.null(contrast)) {
for (i in seq_along(est)) {
out[[i]] <- matrix(NA, nrow = length(object$a_0), ncol = 3)
for (j in seq_along(object$a_0)) {
out[[i]][j, ] <-
rep(object[[est[i]]]$est[j], 3) +
stats::qnorm(c(0.5, (1 - level) / 2, (1 + level) / 2)) *
rep(sqrt(object[[est[i]]]$cov[j, j]), 3)
}
row.names(out[[i]]) <- object$a_0
colnames(out[[i]]) <- c("est", "cil", "ciu")
}
} else if (is.numeric(contrast)) {
# check that contrast came in correctly
if (length(contrast) != length(object$a_0)) {
stop("length of contrast vector not equal to length of a_0")
}
if (!all(contrast %in% c(-1, 1, 0))) {
stop("contrast should only be -1, 1, or 0")
}
for (i in seq_along(est)) {
out[[i]] <- matrix(NA, nrow = 1, ncol = 3)
g <- matrix(contrast, nrow = length(contrast))
p <- matrix(object[[est[i]]]$est, nrow = length(contrast))
v <- object[[est[i]]]$cov
thisC <- t(g) %*% p
thisSe <- sqrt(t(g) %*% v %*% g)
out[[i]][1, ] <- rep(thisC, 3) +
stats::qnorm(c(0.5, (1 - level) / 2, (1 + level) / 2)) * rep(thisSe, 3)
indMinus <- which(contrast == -1)
indPlus <- which(contrast == 1)
plusTerms <- minusTerms <- ""
if (length(indMinus) > 0) {
minusTerms <- paste0("E[Y(", object$a_0[indMinus], ")]")
}
if (length(indPlus) > 0) {
plusTerms <- paste0("E[Y(", object$a_0[indPlus], ")]")
}
thisName <- paste0(
paste0(plusTerms, collapse = "+"),
ifelse(length(indMinus) > 0, "-", ""),
paste0(minusTerms, collapse = "-")
)
row.names(out[[i]]) <- thisName
colnames(out[[i]]) <- c("est", "cil", "ciu")
}
} else if (is.list(contrast)) {
if (!all(c("f", "f_inv", "h", "fh_grad") %in% names(contrast))) {
stop("some function missing in contrast.")
}
for (i in seq_along(est)) {
out[[i]] <- matrix(NA, nrow = 1, ncol = 3)
thisC <- do.call(contrast$h, args = list(est = object[[est[i]]]$est))
f_thisC <- do.call(contrast$f, args = list(eff = thisC))
grad <- matrix(do.call(
contrast$fh_grad,
args = list(est = object[[est[i]]]$est)
), nrow = length(object$a_0))
v <- object[[est[i]]]$cov
thisSe <- sqrt(t(grad) %*% v %*% grad)
transformCI <- rep(f_thisC, 3) +
stats::qnorm(c(0.5, (1 - level) / 2, (1 + level) / 2)) *
rep(thisSe, 3)
out[[i]][1, ] <- do.call(contrast$f_inv, args = list(eff = transformCI))
row.names(out[[i]]) <- c("user contrast")
colnames(out[[i]]) <- c("est", "cil", "ciu")
}
}
class(out) <- "ci.adaptive_iptw"
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.