R/psrwe_powerprior.R

In psrwe: PS-Integrated Methods for Incorporating RWE in Clinical Studies

#' Call STAN models
#'
#' Call STAN models. Called by \code{psrwe_powerp}.
#'
#' @param lst_data List of study data to be passed to STAN
#' @param stan_mdl STAN model including
#'    \describe{
#'     \item{powerps}{PS-power prior model for continuous outcomes}
#'     \item{powerpsbinary}{PS-power prior model for binary outcomes}
#'     \item{powerp}{Power prior model}
#' }
#'
#' @param chains STAN parameter. Number of Markov chainsm
#' @param iter STAN parameter. Number of iterations
#' @param warmup STAN parameter. Number of burnin.
#' @param control STAN parameter. See \code{rstan::stan} for details.
#' @param ... other options to call STAN sampling such as \code{thin},
#'     \code{algorithm}. See \code{rstan::sampling} for details.#'
#'
#' @return Result from STAN sampling
#'
#' @export
#'
rwe_stan <- function(lst_data,
                     stan_mdl = c("powerps", "powerpsbinary", "powerp"),
                     chains = 4, iter = 2000, warmup = 1000,
                     control = list(adapt_delta = 0.95), ...) {

    stan_mdl <- match.arg(stan_mdl)
    stan_rst <- rstan::sampling(stanmodels[[stan_mdl]],
                                data    = lst_data,
                                chains  = chains,
                                iter    = iter,
                                warmup  = warmup,
                                control = control,
                                ...)

    stan_rst
}


#' Get posterior samples based on PS-power prior approach
#'
#' Draw posterior samples of the parameters of interest for the PS-power prior
#' approach
#'
#'
#' @param dta_psbor A class \code{PSRWE_BOR} object generated by
#'     \code{\link{psrwe_borrow}}.
#' @param v_outcome Column name corresponding to the outcome.
#' @param outcome_type Type of outcomes: \code{continuous} or \code{binary}.
#' @param prior_type Whether treat power parameter as fixed (\code{fixed}) or
#'     fully Bayesian (\code{random}).
#' @param seed Random seed.
#' @param ... extra parameters for calling function \code{\link{rwe_stan}}.
#'
#' @return A class \code{PSRWE_RST} list with the following objects
#'
#' \describe{
#'
#'     \item{Observed}{Observed mean and SD of the outcome by group, arm
#'     and stratum}
#'
#'     \item{Control}{A list of estimated mean and SD of the outcome by stratum
#'     in the control arm}
#'
#'     \item{Treatment}{A list of estimated mean and SD of the outcome by
#'     stratum in the treatment arm for RCT}
#'
#'     \item{Effect}{A list of estimated mean and SD of the treatment effect by
#'     stratum for RCT}
#'
#'     \item{Borrow}{Borrowing information from \code{dta_psbor}}
#'
#'     \item{stan_rst}{Result from STAN sampling}
#' }
#'
#' @examples
#'
#' \donttest{
#' data(ex_dta)
#' dta_ps <- psrwe_est(ex_dta,
#'        v_covs = paste("V", 1:7, sep = ""),
#'        v_grp = "Group",
#'        cur_grp_level = "current")
#' ps_borrow <- psrwe_borrow(total_borrow = 30, dta_ps)
#' rst <- psrwe_powerp(ps_borrow, v_outcome = "Y_Con", seed = 123)}
#'
#' @export
#'
psrwe_powerp <- function(dta_psbor, v_outcome = "Y",
                          outcome_type = c("continuous", "binary"),
                          prior_type = c("fixed", "random"),
                          ..., seed = NULL) {

    ## check
    stopifnot(inherits(dta_psbor,
                       what = get_rwe_class("PSDIST")))

    type       <- match.arg(outcome_type)
    prior_type <- match.arg(prior_type)
    stopifnot(v_outcome %in% colnames(dta_psbor$data))

    ## save the seed from global if any then set random seed
    old_seed <- NULL
    if (!is.null(seed)) {
        if (exists(".Random.seed", envir = .GlobalEnv)) {
            old_seed <- get(".Random.seed", envir = .GlobalEnv)
        }
        set.seed(seed)
    }

    ## observed
    rst_obs <- get_observed(dta_psbor$data, v_outcome)

    ## prepare stan data
    lst_dta <- get_stan_data(dta_psbor, v_outcome, prior_type)

    ## sampling
    stan_mdl <- if_else("continuous" == type,
                        "powerps",
                        "powerpsbinary")

    ctl_post   <- rwe_stan(lst_data = lst_dta$ctl, stan_mdl = stan_mdl, ...)
    ctl_thetas <- extract(ctl_post, "thetas")$thetas

    is_rct     <- dta_psbor$is_rct
    trt_post   <- NULL
    trt_thetas <- NULL
    if (is_rct) {
        trt_post <- rwe_stan(lst_data = lst_dta$trt,
                             stan_mdl = stan_mdl, ...)
        trt_thetas <- extract(trt_post, "thetas")$thetas
    }

    ## summary
    rst_trt    <- NULL
    rst_effect <- NULL
    if (is_rct) {
        rst_trt    <- get_post_theta(trt_thetas, dta_psbor$Borrow$N_Cur_TRT)
        rst_effect <- get_post_theta(trt_thetas - ctl_thetas,
                                     dta_psbor$Borrow$N_Current)
        n_ctl      <- dta_psbor$Borrow$N_Cur_CTL
    } else {
        n_ctl      <- dta_psbor$Borrow$N_Current
    }
    rst_ctl <- get_post_theta(ctl_thetas, n_ctl)

    ## reset the orignal seed back to the global or
    ## remove the one set within this session earlier.
    if (!is.null(seed)) {
        if (!is.null(old_seed)) {
            invisible(assign(".Random.seed", old_seed, envir = .GlobalEnv))
        } else {
            invisible(rm(list = c(".Random.seed"), envir = .GlobalEnv))
        }
    }

    ## return
    rst <-  list(stan_rst = list(ctl_post = ctl_post,
                                 trt_post = trt_post),
                 Observed  = rst_obs,
                 Control   = rst_ctl,
                 Treatment = rst_trt,
                 Effect    = rst_effect,
                 Borrow    = dta_psbor$Borrow,
                 Total_borrow = dta_psbor$Total_borrow,
                 Method       = "ps_pp",
                 Outcome_type = type,
                 is_rct       = is_rct)

    class(rst) <- get_rwe_class("ANARST")
    rst
}


#' Get data for STAN
#'
#'
#' @noRd
#'
get_stan_data <- function(dta_psbor, v_outcome, prior_type) {
    f_curd <- function(i, d1, d0 = NULL) {
        cur_d <- c(N1    = length(d1),
                   YBAR1 = mean(cur_d1),
                   YSUM1 = sum(cur_d1))

        if (is.null(d0)) {
            cur_d <- c(cur_d,
                       N0    = 0,
                       YBAR0 = 0,
                       SD0   = 0)
        } else {
            cur_d <- c(cur_d,
                       N0    = length(cur_d0),
                       YBAR0 = mean(cur_d0),
                       SD0   = sd(cur_d0))
        }

        list(stan_d = cur_d,
             y1     = d1,
             inx1   = rep(i, length(d1)))
    }

    is_rct  <- dta_psbor$is_rct
    data    <- dta_psbor$data
    data    <- data[!is.na(data[["_strata_"]]), ]

    strata  <- levels(data[["_strata_"]])
    nstrata <- length(strata)

    ctl_stan_d  <- NULL
    ctl_y1      <- NULL
    ctl_inx1    <- NULL

    trt_stan_d  <- NULL
    trt_y1      <- NULL
    trt_inx1    <- NULL

    for (i in seq_len(nstrata)) {
        cur_01  <- get_cur_d(data, strata[i], v_outcome)
        cur_d1  <- cur_01$cur_d1
        cur_d0  <- cur_01$cur_d0
        cur_d1t <- cur_01$cur_d1t

        ctl_cur    <- f_curd(i, cur_d1, cur_d0)
        ctl_stan_d <- rbind(ctl_stan_d, ctl_cur$stan_d)
        ctl_y1     <- c(ctl_y1,   ctl_cur$y1)
        ctl_inx1   <- c(ctl_inx1, ctl_cur$inx1)

        if (is_rct) {
            trt_cur    <- f_curd(i, cur_d1t)
            trt_stan_d <- rbind(trt_stan_d, trt_cur$stan_d)
            trt_y1     <- c(trt_y1,   trt_cur$y1)
            trt_inx1   <- c(trt_inx1, trt_cur$inx1)
        }
    }

    ctl_lst_data  <- list(S     = nstrata,
                          A     = dta_psbor$Total_borrow,
                          RS    = as.array(dta_psbor$Borrow$Proportion),
                          FIXVS = as.numeric(prior_type == "fixed"),
                          N0    = as.array(ctl_stan_d[, "N0"]),
                          N1    = as.array(ctl_stan_d[, "N1"]),
                          YBAR0 = as.array(ctl_stan_d[, "YBAR0"]),
                          SD0   = as.array(ctl_stan_d[, "SD0"]),
                          TN1   = length(ctl_y1),
                          Y1    = ctl_y1,
                          INX1  = ctl_inx1,
                          YBAR1 = as.array(ctl_stan_d[, "YBAR1"]),
                          YSUM1 = as.array(ctl_stan_d[, "YSUM1"]))

    trt_lst_data <- NULL
    if (is_rct) {
        trt_lst_data  <- list(S     = nstrata,
                              A     = 0,
                              RS    = as.array(dta_psbor$Borrow$Proportion),
                              FIXVS = as.numeric("fixed" == "fixed"),
                              N0    = as.array(trt_stan_d[, "N0"]),
                              N1    = as.array(trt_stan_d[, "N1"]),
                              YBAR0 = as.array(trt_stan_d[, "YBAR0"]),
                              SD0   = as.array(trt_stan_d[, "SD0"]),
                              TN1   = length(trt_y1),
                              Y1    = trt_y1,
                              INX1  = trt_inx1,
                              YBAR1 = as.array(trt_stan_d[, "YBAR1"]),
                              YSUM1 = as.array(trt_stan_d[, "YSUM1"]))
    }


    list(ctl = ctl_lst_data,
         trt = trt_lst_data)
}


#' Get results for control, treatment and effect
#'
#'
#' @noRd
#'
get_post_theta <- function(thetas, weights) {
    ws      <- weights / sum(weights)
    samples <- t(thetas)
    overall <- apply(samples, 2, function(x) sum(x * ws))

    means <- apply(samples, 1, mean)
    sds   <- apply(samples, 1, sd)

    mean_overall <- mean(overall)
    sd_overall   <- sd(overall)

    list(Stratum_Samples  = samples,
         Overall_Samples  = overall,
         Stratum_Estimate = cbind(Mean   = means,
                                  StdErr = sds),
         Overall_Estimate = cbind(Mean   = mean_overall,
                                  StdErr = sd_overall))
}


#' @title Summarize overall estimation results
#'
#' @description S3 method summarizing overall estimation results
#'
#'
#' @param object A list of class \code{PSRWE_RST} that is generated using the
#'     \code{\link{psrwe_powerp}}, \code{\link{psrwe_compl}}, or
#'     \code{\link{psrwe_survkm}} function.
#' @param ... Additional parameters.
#'
#' @return A list with data frames for the borrowing and estimation results.
#'
#' @method summary PSRWE_RST
#'
#'
#' @export
#'
summary.PSRWE_RST <- function(object, ...) {

    is_rct <- object$is_rct

    ## overall estimate
    if (is_rct) {
        type <- "Effect"
    } else {
        type <- "Control"
    }

    dtype <- object[[type]]$Overall_Estimate
    if ("ps_km" == object$Method) {
        dtype <- data.frame(dtype) %>%
            filter(object$pred_tp == T)
    }

    Overall <- data.frame(Type   = type,
                          Mean   = dtype[1, "Mean"],
                          StdErr = dtype[1, "StdErr"])

    rst <- list(Overall = Overall)

    ## for CI
    if (exists("CI", object)) {
        citype <- cbind(object$CI[[type]]$Overall_Estimate)

        if (is.data.frame(citype)) {
            if ("ps_km" == object$Method) {
                citype <- cbind(citype,
                                T = object[[type]]$Overall_Estimate$T)

                citype <- data.frame(citype) %>%
                    filter(object$pred_tp == T)
            }

            rst$CI <- data.frame(Lower    = citype[1, "Lower"],
                                 Upper    = citype[1, "Upper"],
                                 ConfLvl  = object$CI$Conf_int,
                                 ConfType = object$CI$Conf_type,
                                 MethodCI = object$CI$Method_ci,
                                 Method   = object$Method)
        }
    }

    ## for hypothesis inference
    if (exists("INFER", object)) {
        hitype <- object$INFER[[type]]$Overall_InferProb

        if ("ps_km" == object$Method) {
            hitype <- cbind(hitype,
                            T = object[[type]]$Overall_Estimate$T)

            hitype <- data.frame(hitype) %>%
                filter(object$pred_tp == T)
        }

        rst$INFER <- data.frame(InferP   = hitype[1, "Infer_prob"],
                                MethodHI = object$INFER$Method_infer,
                                Alter    = object$INFER$Alternative,
                                Mu       = object$INFER$Mu)
    }

    return(rst)
}

#' @title Print estimation results
#'
#' @description Print summary information of outcome mean estimation results
#'
#' @param x A list of class \code{PSRWE_RST} that is generated using the
#'     \code{\link{psrwe_powerp}}, \code{\link{psrwe_compl}}, or
#'     \code{\link{psrwe_survkm}} function.
#' @param ... Additional parameters
#'
#' @seealso  \code{\link{summary.PSRWE_RST}}
#'
#'
#' @method print PSRWE_RST
#'
#'
#' @export
#'
print.PSRWE_RST <- function(x, ...) {

    rst <- summary(x, ...)
    rst_sum <- rst$Overall
    rst_ci <- rst$CI
    rst_infer <- rst$INFER

    extra_1 <- ""
    if ("ps_km" == x$Method) {
        extra_1 <- paste(" based on the survival probability at time ",
                         x$pred_tp, ",", sep = "")
    }

    if ("Effect" == rst_sum[1, "Type"]) {
        extra_2 <- "treatment effect"
        extra_2_param <- "theta_trt-theta_ctl"
    } else {
        extra_2 <- "point estimate"
        extra_2_param <- "theta"
    }

    ## for CI
    extra_ci <- ""
    if (exists("CI", x)) {
        extra_ci <- paste(" Confidence interval: ",
                          "(",
                          sprintf("%5.3f", rst_ci[1, "Lower"]),
                          ", ",
                          sprintf("%5.3f", rst_ci[1, "Upper"]),
                          ")",
                          " with two-sided level of ",
                          sprintf("%5.3f", rst_ci[1, "ConfLvl"]),
                          " by ",
                          rst_ci[1, "MethodCI"],
                          " method",
                          sep = "")

        if (!is.na(rst_ci[1, "ConfType"])) {
            extra_ci <- paste(extra_ci,
                              " and ",
                              rst_ci[1, "ConfType"],
                              " transformation",
                              sep = "")                              
        }

        extra_ci <- paste(extra_ci, ".", sep = "")
    }

    ## for hypothesis inference
    extra_hi <- ""
    if (exists("INFER", x)) {
        extra_hi <- paste(" Hypothesis inference: ",
                          "H0: ", extra_2_param, " ",
                          ifelse(rst_infer[1, "Alter"] == "less", ">=", "<="),
                          " ", sprintf("%5.3f", rst_infer[1, "Mu"]),
                          " vs. ",
                          "Ha: ", extra_2_param, " ",
                          ifelse(rst_infer[1, "Alter"] == "less", "<", ">"),
                          " ", sprintf("%5.3f", rst_infer[1, "Mu"]),
                          " with the ", rst_infer[1, "MethodHI"],
                          " ",
                          sprintf("%5.5f", rst_infer[1, "InferP"]),
                          ".",
                          sep = "")
    }

    ## combine all above
    ss <- paste("With a total of ", x$Total_borrow,
                " subject borrowed from the RWD,",
                extra_1, " the ", extra_2,
                " is ",
                sprintf("%5.3f", rst_sum[1, "Mean"]),
                " with standard error ",
                sprintf("%5.3f", rst_sum[1, "StdErr"]),
                ".",
                extra_ci,
                extra_hi,
                sep = "")

    cat_paste(ss)
}

#' @title Plot estimation results for power prior approach
#'
#' @description S3 method plotting estimation results
#'
#'
#' @param x A list of class \code{PSRWE_RST} that is generated using the
#'     \code{\link{psrwe_powerp}}, \code{\link{psrwe_compl}}, or
#'     \code{\link{psrwe_survkm}} function.
#' @param ... Additional parameters.
#'
#'
#' @method plot PSRWE_RST
#'
#'
#' @export
#'
plot.PSRWE_RST <- function(x, ...) {
    rst <- switch(x$Method,
                  ps_pp = plot_pp_rst(x, ...),
                  ps_km = plot_km_rst(x, ...),
                  ps_cl = {
                      stop("This method is currently unavailable
                            for composite likelihood analysis.")
                  })

    rst
}