R/im_sace.R

In idem: Inference in Randomized Controlled Trials with Death and Missingness

#' Summary of the inference results
#'
#' Summarize survivors only or Survivor Averaged Causal Effect (SACE) based on
#' the imputation and bootstrap analysis
#'
#' @param object A class \code{IDEMINFER} list generated by \code{\link{imInfer}}
#' @param opt Types of the summary
#' \itemize{
#' \item{\code{survivor}: }{Survivors only analysis}
#' \item{\code{SACE}: }{Survivor Averaged Causal Effect}
#' }
#'
#' @param sace.deltas Vector of sensitivity parameters for SACE estimation. If
#'     \code{NULL}, the values will be generated based on the standard
#'     deviations of the estimated differences in the functional outcomes
#'     between the treatment and control arms
#' @param ... Optional arguments for summary
#' @details
#'
#' For SACE, the default sensitivity parameters will be determined by the
#' standard deviation of the treatment effect size on the functional outcomes.
#'
#' @return
#'
#' A class \code{summary.IDEMINFER} list containing
#' \describe{
#' \item{deltas}{imputation sensitivity parameters}
#' \item{n.boot}{number of bootstrap samples in bootstrap analysis}
#' \item{sace.deltas}{SACE sensitivity parameters when \code{opt = SACE}}
#' \item{rst}{A data frame with columns
#' \itemize{
#' \item \code{Delta0}: Imputation sensitivity parameter for control arm,
#' \item \code{Delta1}: Imputation sensitivity parameter for intervention arm
#' \item \code{SACE_Delta}: SACE sensitivity parameter when \code{opt = SACE}
#' \item \code{Effect}: SACE estimate
#' \item \code{LB}: Lower bound of the 95% CI when \code{n.boot > 0} in the \code{IDEMINFER} object
#' \item \code{UB}: Upper bound of the 95% CI when \code{n.boot > 0} in the \code{IDEMINFER} object
#' \item \code{PValue}: p-value when when \code{n.boot > 0} in the \code{IDEMINFER} object
#' }}}
#'
#' @examples
#' \dontrun{
#' rst.abc <- imData(abc, trt="TRT", surv="SURV", outcome=c("Y1","Y2"),
#'                  y0=NULL, endfml="Y2",
#'                  trt.label = c("UC+SBT", "SAT+SBT"),
#'                  cov=c("AGE"), duration=365, bounds=c(0,100));
#' rst.fit  <- imFitModel(rst.abc);
#' rst.imp <- imImpAll(rst.fit, deltas=c(-0.25,0,0.25),
#'                     normal=TRUE, chains = 2, iter = 2000, warmup = 1000);
#' rst.infer <- imInfer(rst.imp, n.boot = 100, effect.quantiles = c(0.25,0.5,0.75));
#' rst.sace  <- summary(rst.infer, opt = "SACE")}
#'
#' @references
#'
#' Chiba Y, VanderWeele TJ (2011). A simple method for principal strata effects when
#' the outcome has been truncated due to death. American Journal of Epidemiology 173(7):745-751.
#'
#' @method summary IDEMINFER
#'
#' @export
#'
summary.IDEMINFER <- function(object, opt = c("survivor", "SACE"), sace.deltas=NULL, ...) {
    opt       <- match.arg(opt);
    surv.data <- object$survivor;
    n.boot    <- length(object$bootstrap);

    if ("survivor" == opt) {
        rst        <- surv.data[, c("Delta0", "Delta1")];
        rst$Effect <- surv.data[, "Diff"];

        if (n.boot > 0) {
            rst$LB     <- surv.data[, "Diff"] - 1.96*surv.data[, "SD"];
            rst$UB     <- surv.data[, "Diff"] + 1.96*surv.data[, "SD"];
            rst$PValue <- get.pval(surv.data[, "Diff"], surv.data[, "SD"]);
        }
        sace.deltas <- NULL;
    } else {
        ##set sace sensitivity parameters
        ##stopifnot(is.null(sace.deltas));
        if(is.null(sace.deltas)) {
            if (n.boot > 0) {
                cur.sd <- mean(surv.data$SD);
            } else {
                cur.sd <- sd(surv.data$Diff);
            }
            sace.deltas <- -seq(0, cur.sd, len = 5);
        } else {
            sace.deltas <- sort(unique(c(0, sace.deltas)), decreasing = TRUE);
        }

        rst  <- NULL;
        for(i in sace.deltas) {
            cur.rst            <- surv.data[, c("Delta0", "Delta1")];
            cur.rst$Effect     <- surv.data$Diff - i;
            cur.rst$SACE_Delta <- i;

            if (n.boot > 0) {
                cur.rst$LB     <- cur.rst$Effect - 1.96*surv.data$SD;
                cur.rst$UB     <- cur.rst$Effect + 1.96*surv.data$SD;
                cur.rst$PValue <- get.pval(cur.rst$Effect, surv.data$SD);
            }

            rst <- rbind(rst, cur.rst);
        }
    }

    ##return
    rownames(rst) <- NULL;
    rtn.rst       <- list(deltas      = object$deltas,
                          sace.deltas = sace.deltas,
                          n.boot      = n.boot,
                          lst.var     = object$lst.var,
                          rst         = data.frame(rst));
    class(rtn.rst) <- c(class(rtn.rst),
                        paste("summary", get.const("TEST.CLASS"), sep="."));
    invisible(rtn.rst);
}

#' Print survivors only or SACE analysis results
#'
#' Print survivors only or SACE analysis results
#'
#' @inheritParams print.IDEMINFER
#'
#' @param ... Extra arguments
#'
#' @method print summary.IDEMINFER
#'
#' @export
#'
print.summary.IDEMINFER <- function(x, delta0 = NULL, delta1 = NULL, ...) {

    cat("\nThe imputation sensitivity parameters considered were\n");
    print(x$deltas);
    sace <- x$rst;

    if (is.null(x$sace.deltas)) {
        cat("\n\nThe estimated survivors only treatment effects are\n\n");
        if (!is.null(delta0)) {
            sace <- subset(sace, sace$Delta0 == delta0);
        }
        if (!is.null(delta1) & nrow(sace) > 0) {
            sace <- subset(sace, sace$Delta1 == delta1);
        }
        print(sace);

        cat("\n\nPLEASE BE CAUTIOUS that survivors only analysis is only valid \n");
        cat("when the treatment has no impact on survival.\n")
    } else {
        cat("\n\nThe SACE sensitivity parameters considered were\n");
        print(x$sace.deltas);

        cat("\n\nThe estimated SACE are\n\n");
        if (!is.null(delta0)) {
            sace <- subset(sace, sace$Delta0 == delta0);
        }

        if (!is.null(delta1) & nrow(sace) > 0) {
            sace <- subset(sace, sace$Delta1 == delta1);
        }
        print(sace);
    }
}


#' Plot survivors only and SACE analysis results
#'
#' Generate a plot of survivor only and survivor average causal effect values
#'
#' @inheritParams plot.IDEMINFER
#' @inheritParams print.IDEMINFER
#'
#' @param x A class \code{summary.IDEMSACE} object generated by summary of
#'     \code{IDEMINFER}
#' @param by.sace Logical value. If True, create a contour plot for given SACE
#'     sensitivity parameter. Otherwise, create a plot for treatment effect for
#'     given imputation sensitivity parameters
#' @param sace.delta Single SACE sensitivity parameter
#' @param ... Options for \code{plot}
#'
#' @details
#'
#' The plot function will only generate the contour plot of p-values or
#' treatment effects on functional outcomes for survivors only analyses.
#'
#' For SACE analysis, the plot function generates contour plot of line plot
#' based on the value of \code{by.sace}.
#'
#' @examples
#' \dontrun{
#' rst.abc <- imData(abc, trt="TRT", surv="SURV", outcome=c("Y1","Y2"),
#'                  y0=NULL, endfml="Y2",
#'                  trt.label = c("UC+SBT", "SAT+SBT"),
#'                  cov=c("AGE"), duration=365, bounds=c(0,100));
#' rst.fit  <- imFitModel(rst.abc);
#' rst.imp <- imImpAll(rst.fit, deltas=c(-0.25,0,0.25),
#'                     normal=TRUE, chains = 2, iter = 2000, warmup = 1000);
#' rst.infer <- imInfer(rst.imp, n.boot = 100, effect.quantiles = c(0.25,0.5,0.75));
#' rst.survivors <- summary(rst.infer, opt="survivor");
#' plot(rst.survivors);}
#'
#' @method plot summary.IDEMINFER
#'
#' @export
#'
plot.summary.IDEMINFER <- function(x,
                                   opt = c("pvalue", "effect"),
                                   by.sace = TRUE,
                                   delta0 = 0,
                                   delta1 = 0,
                                   sace.delta = NULL,
                                   ...) {


    opt <- match.arg(opt);
    col.var <- switch(opt,
                      pvalue = "PValue",
                      effect = "Effect")

    if ("PValue" == col.var & 0 == x$n.boot) {
        cat("\nBootstrap analysis was not conducted.");
        cat("\nOnly treatment effect estimation is available.\n");
        col.var <- "Effect";
    }

    sace    <- x$rst;
    lst.var <- x$lst.var;
    trt.len <- NULL;
    get.para(lst.var, environment());

    if (is.null(x$sace.deltas)) {
        ##survivors only analysis
        cat("\nPLEASE BE CAUTIOUS that survivors only analysis is only valid \n");
        cat("when the treatment has no impact on survival.\n")
        plot.contour(sace, trt.len, col.var, ...);
    } else {
        if (by.sace) {
            stopifnot(1 == length(delta0)
                      & 1 == length(delta1)
                      & delta0 %in% x$deltas
                      & delta1 %in% x$deltas);

            sace.data   <- subset(sace, sace$Delta0 == delta0 & sace$Delta1 == delta1);
            sace.deltas <- sace.data$SACE_Delta;

            ylim <- range(sace.data[,c("Effect", "LB", "UB")]);
            ylim <- ylim + c(-0.05*(ylim[2] - ylim[1]),
                             0.05*(ylim[2] - ylim[1]));

            plot(sace.deltas,
                 sace.data$Effect,
                 type = "b",
                 ylim = ylim,
                 xlab='SACE Sensitivity Parameter',
                 ylab='SACE', ...);

            abline(h=0,lty=2,lwd=2, col="gray");
            if (x$n.boot > 0) {
                lines(sace.deltas, sace.data$LB, lty = 2);
                lines(sace.deltas, sace.data$UB, lty = 2);
            }
        } else {
            ##sace contour plot
            stopifnot(1 == length(sace.delta) &
                      sace.delta %in% x$sace.deltas);
            cur.data <- subset(sace, sace$SACE_Delta == sace.delta);
            plot.contour(cur.data, trt.len, col.var, ...);
        }
    }

}