Nothing
#' Confidence Intervals for Counterfactual Mean Under Stochastic Intervention
#'
#' @details Compute confidence intervals for estimates produced by
#' \code{\link{txshift}}.
#'
#' @param object An object of class \code{txshift}, as produced by invoking
#' the function \code{\link{txshift}}, for which a confidence interval is to
#' be computed.
#' @param parm A \code{numeric} vector indicating indices of \code{object$est}
#' for which to return confidence intervals.
#' @param level A \code{numeric} indicating the level of the confidence
#' interval to be computed.
#' @param ... Other arguments. Not currently used.
#' @param ci_mult Pre-computed multipliers for generating confidence intervals.
#' The default of \code{NULL} should generally NOT be changed and is only used
#' by the internal machinery for creating simultaneous confidence bands.
#'
#' @importFrom stats qnorm plogis qlogis
#'
#' @method confint txshift
#'
#' @return A named \code{numeric} vector containing the parameter estimate from
#' a \code{txshift} object, alongside lower and upper Wald-style confidence
#' intervals at a specified coverage level.
#'
#' @examples
#' set.seed(429153)
#' n_obs <- 100
#' W <- replicate(2, rbinom(n_obs, 1, 0.5))
#' A <- rnorm(n_obs, mean = 2 * W, sd = 1)
#' Y <- rbinom(n_obs, 1, plogis(A + W + rnorm(n_obs, mean = 0, sd = 1)))
#' txout <- txshift(
#' W = W, A = A, Y = Y, delta = 0.5,
#' estimator = "tmle",
#' g_exp_fit_args = list(
#' fit_type = "hal", n_bins = 5,
#' grid_type = "equal_mass",
#' lambda_seq = exp(-1:-9)
#' ),
#' Q_fit_args = list(
#' fit_type = "glm",
#' glm_formula = "Y ~ ."
#' )
#' )
#' confint(txout)
#' @export
confint.txshift <- function(object,
parm = seq_len(object$psi),
level = 0.95,
...,
ci_mult = NULL) {
# first, let's get Z_(1 - alpha)
if (is.null(ci_mult)) {
ci_mult <- c(-1, 1) * abs(stats::qnorm(p = (1 - level) / 2))
}
if (length(unique(object$.outcome)) > 2) { # assume continuous outcome
# compute the EIF variance multiplier for the CI
# NOTE: the variance value is already scaled by length of observations
sd_eif <- sqrt(object$var)
# compute the interval around the point estimate
ci_psi <- ci_mult * sd_eif + object$psi
} else if (length(unique(object$.outcome)) == 2) { # binary outcome
# for binary outcome case, compute on the logit scale and back-transform
psi_ratio <- stats::qlogis(object$psi)
grad_ratio_delta <- (1 / object$psi) + (1 / (1 - object$psi))
se_eif_logit <- sqrt(grad_ratio_delta^2 * object$var)
ci_psi <- stats::plogis(ci_mult * se_eif_logit + psi_ratio)
} else {
stop("The outcome has fewer than 2 levels: this case is not supported.")
}
# set up output CI object
ci_out <- c(ci_psi[1], object$psi, ci_psi[2])
names(ci_out) <- c("lwr_ci", "est", "upr_ci")
return(ci_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.