R/bdpnormal.R
In graemeleehickey/bayesDP: Implementation of the Bayesian Discount Prior Approach for Clinical Trials
Defines functions
posterior_normal
#' @title Bayesian Discount Prior: Gaussian mean values
#'
#' @description \code{bdpnormal} is used for estimating posterior samples from a
#' Gaussian outcome where an informative prior is used. The prior weight is
#' determined using a discount function. This code is modeled after the
#' methodologies developed in Haddad et al. (2017).
#'
#' @param mu_t scalar. Mean of the current treatment group.
#' @param sigma_t scalar. Standard deviation of the current treatment group.
#' @param N_t scalar. Number of observations of the current treatment group.
#' @param mu0_t scalar. Mean of the historical treatment group.
#' @param sigma0_t scalar. Standard deviation of the historical treatment group.
#' @param N0_t scalar. Number of observations of the historical treatment group.
#' @param mu_c scalar. Mean of the current control group.
#' @param sigma_c scalar. Standard deviation of the current control group.
#' @param N_c scalar. Number of observations of the current control group.
#' @param mu0_c scalar. Mean of the historical control group.
#' @param sigma0_c scalar. Standard deviation of the historical control group.
#' @param N0_c scalar. Number of observations of the historical control group.
#' @param discount_function character. Specify the discount function to use.
#' Currently supports \code{weibull}, \code{scaledweibull}, and
#' \code{identity}. The discount function \code{scaledweibull} scales the
#' output of the Weibull CDF to have a max value of 1. The \code{identity}
#' discount function uses the posterior probability directly as the discount
#' weight. Default value is "\code{identity}".
#' @param alpha_max scalar. Maximum weight the discount function can apply.
#' Default is 1. For a two-arm trial, users may specify a vector of two values
#' where the first value is used to weight the historical treatment group and
#' the second value is used to weight the historical control group.
#' @param fix_alpha logical. Fix alpha at alpha_max? Default value is FALSE.
#' @param weibull_shape scalar. Shape parameter of the Weibull discount function
#' used to compute alpha, the weight parameter of the historical data. Default
#' value is 3. For a two-arm trial, users may specify a vector of two values
#' where the first value is used to estimate the weight of the historical
#' treatment group and the second value is used to estimate the weight of the
#' historical control group. Not used when \code{discount_function} =
#' "identity".
#' @param weibull_scale scalar. Scale parameter of the Weibull discount function
#' used to compute alpha, the weight parameter of the historical data. Default
#' value is 0.135. For a two-arm trial, users may specify a vector of two
#' values where the first value is used to estimate the weight of the
#' historical treatment group and the second value is used to estimate the
#' weight of the historical control group. Not used when
#' \code{discount_function} = "identity".
#' @param number_mcmc scalar. Number of Monte Carlo simulations. Default is
#' 10000.
#' @param method character. Analysis method with respect to estimation of the
#' weight parameter alpha. Default method "\code{mc}" estimates alpha for each
#' Monte Carlo iteration. Alternate value "\code{fixed}" estimates alpha once
#' and holds it fixed throughout the analysis. See the the \code{bdpnormal}
#' vignette \cr \code{vignette("bdpnormal-vignette", package="bayesDP")} for
#' more details.
#' @param compare logical. Should a comparison object be included in the fit?
#' For a one-arm analysis, the comparison object is simply the posterior chain
#' of the treatment group parameter. For a two-arm analysis, the comparison
#' object is the posterior chain of the treatment effect that compares
#' treatment and control. If \code{compare=TRUE}, the comparison object is
#' accessible in the \code{final} slot, else the \code{final} slot is
#' \code{NULL}. Default is \code{TRUE}.
#' @details \code{bdpnormal} uses a two-stage approach for determining the
#' strength of historical data in estimation of a mean outcome. In the first
#' stage, a \emph{discount function} is used that that defines the maximum
#' strength of the historical data and discounts based on disagreement with
#' the current data. Disagreement between current and historical data is
#' determined by stochastically comparing the respective posterior
#' distributions under noninformative priors. With Gaussian data, the
#' comparison is the proability (\code{p}) that the current mean is less than
#' the historical mean. The comparison metric \code{p} is then input into the
#' discount function and the final strength of the historical data is returned
#' (alpha).
#'
#' In the second stage, posterior estimation is performed where the discount
#' function parameter, \code{alpha}, is used incorporated in all posterior
#' estimation procedures.
#'
#' To carry out a single arm (OPC) analysis, data for the current treatment
#' (\code{mu_t}, \code{sigma_t}, and \code{N_t}) and historical treatment
#' (\code{mu0_t}, \code{sigma0_t}, and \code{N0_t}) must be input. The results
#' are then based on the posterior distribution of the current data augmented
#' by the historical data.
#'
#' To carry our a two-arm (RCT) analysis, data for the current treatment and
#' at least one of current or historical control data must be input. The
#' results are then based on the posterior distribution of the difference
#' between current treatment and control, augmented by available historical
#' data.
#'
#' For more details, see the \code{bdpnormal} vignette: \cr
#' \code{vignette("bdpnormal-vignette", package="bayesDP")}
#'
#' @return \code{bdpnormal} returns an object of class "bdpnormal". The
#' functions \code{\link[=summary,bdpnormal-method]{summary}} and
#' \code{\link[=print,bdpnormal-method]{print}} are used to obtain and print a
#' summary of the results, including user inputs. The
#' \code{\link[=plot,bdpnormal-method]{plot}} function displays visual outputs
#' as well.
#'
#' An object of class \code{bdpnormal} is a list containing at least
#' the following components:
#' \describe{
#' \item{\code{posterior_treatment}}{
#' list. Entries contain values related to the treatment group:
#' \itemize{
#' \item \code{alpha_discount}
#' numeric. Alpha value, the weighting parameter of the historical data.
#' \item \code{p_hat}
#' numeric. The posterior probability of the stochastic comparison
#' between the current and historical data.
#' \item \code{posterior_mu}
#' vector. A vector of length \code{number_mcmc} containing the posterior
#' mean of the treatment group. If historical treatment data is present,
#' the posterior incorporates the weighted historical data.
#' \item \code{posterior_sigma2}
#' vector. A vector of length \code{number_mcmc} containing the posterior
#' variance of the treatment group. If historical treatment data is present,
#' the posterior incorporates the weighted historical data.
#' \item \code{posterior_flat_mu}
#' vector. A vector of length \code{number_mcmc} containing
#' Monte Carlo samples of the mean of the current treatment group
#' under a flat/non-informative prior, i.e., no incorporation of the
#' historical data.
#' \item \code{posterior_flat_sigma2}
#' vector. A vector of length \code{number_mcmc} containing
#' Monte Carlo samples of the standard deviation of the current treatment group
#' under a flat/non-informative prior, i.e., no incorporation of the
#' historical data.
#' \item \code{prior_mu}
#' vector. If historical treatment data is present, a vector of length
#' \code{number_mcmc} containing Monte Carlo samples of the mean
#' of the historical treatment group under a flat/non-informative prior.
#' \item\code{prior_sigma2}
#' vector. If historical treatment data is present, a vector of length
#' \code{number_mcmc} containing Monte Carlo samples of the standard deviation
#' of the historical treatment group under a flat/non-informative prior.
#' }
#' }
#' \item{\code{posterior_control}}{
#' list. Similar entries as \code{posterior_treament}. Only present if a
#' control group is specified.
#' }
#' \item{\code{final}}{
#' list. Contains the final comparison object, dependent on the analysis type:
#' \itemize{
#' \item One-arm analysis:
#' vector. Posterior chain of the mean.
#' \item Two-arm analysis:
#' vector. Posterior chain of the mean difference comparing treatment and
#' control groups.
#' }
#' }
#' \item{\code{args1}}{
#' list. Entries contain user inputs. In addition, the following elements
#' are output:
#' \itemize{
#' \item \code{arm2}
#' binary indicator. Used internally to indicate one-arm or two-arm
#' analysis.
#' \item \code{intent}
#' character. Denotes current/historical status of treatment and
#' control groups.
#' }
#' }
#' }
#'
#'
#' @seealso \code{\link[=summary,bdpnormal-method]{summary}},
#' \code{\link[=print,bdpnormal-method]{print}},
#' and \code{\link[=plot,bdpnormal-method]{plot}} for details of each of the
#' supported methods.
#'
#' @references
#' Haddad, T., Himes, A., Thompson, L., Irony, T., Nair, R. MDIC Computer
#' Modeling and Simulation working group.(2017) Incorporation of stochastic
#' engineering models as prior information in Bayesian medical device trials.
#' \emph{Journal of Biopharmaceutical Statistics}, 1-15.
#'
#' @examples
#' # One-arm trial (OPC) example
#' fit <- bdpnormal(
#' mu_t = 30, sigma_t = 10, N_t = 50,
#' mu0_t = 32, sigma0_t = 10, N0_t = 50,
#' method = "fixed"
#' )
#' summary(fit)
#' \dontrun{
#' plot(fit)
#' }
#'
#' # Two-arm (RCT) example
#' fit2 <- bdpnormal(
#' mu_t = 30, sigma_t = 10, N_t = 50,
#' mu0_t = 32, sigma0_t = 10, N0_t = 50,
#' mu_c = 25, sigma_c = 10, N_c = 50,
#' mu0_c = 25, sigma0_c = 10, N0_c = 50,
#' method = "fixed"
#' )
#' summary(fit2)
#' \dontrun{
#' plot(fit2)
#' }
#'
#' @rdname bdpnormal
#' @import methods
#' @importFrom stats sd density is.empty.model median model.offset
#' model.response pweibull quantile rbeta rgamma rnorm var vcov
#' @importFrom ggplot2 aes facet_wrap geom_hline geom_line geom_vline ggplot
#' theme_bw xlab ylab
#' @aliases bdpnormal-method
#' @aliases bdpnormal,ANY-method
#' @export bdpnormal
bdpnormal <- setClass("bdpnormal", slots = c(
posterior_treatment = "list",
posterior_control = "list",
final = "list",
args1 = "list"
))
setGeneric(
"bdpnormal",
function(mu_t = NULL,
sigma_t = NULL,
N_t = NULL,
mu0_t = NULL,
sigma0_t = NULL,
N0_t = NULL,
mu_c = NULL,
sigma_c = NULL,
N_c = NULL,
mu0_c = NULL,
sigma0_c = NULL,
N0_c = NULL,
discount_function = "identity",
alpha_max = 1,
fix_alpha = FALSE,
weibull_scale = 0.135,
weibull_shape = 3,
number_mcmc = 10000,
method = "mc",
compare = TRUE) {
standardGeneric("bdpnormal")
}
)
setMethod(
"bdpnormal",
signature(),
function(mu_t = NULL,
sigma_t = NULL,
N_t = NULL,
mu0_t = NULL,
sigma0_t = NULL,
N0_t = NULL,
mu_c = NULL,
sigma_c = NULL,
N_c = NULL,
mu0_c = NULL,
sigma0_c = NULL,
N0_c = NULL,
discount_function = "identity",
alpha_max = 1,
fix_alpha = FALSE,
weibull_scale = 0.135,
weibull_shape = 3,
number_mcmc = 10000,
method = "mc",
compare = TRUE) {
################################################################################
# Check Input #
################################################################################
intent <- c()
if (length(mu_t + sigma_t + N_t) != 0) {
intent <- c(intent, "current treatment")
# cat("Current Treatment\n")
} else {
if (is.null(mu_t) == TRUE) {
cat("mu_t missing\n")
}
if (is.null(sigma_t) == TRUE) {
cat("sigma_t missing\n")
}
if (is.null(N_t) == TRUE) {
cat("N_t missing\n")
}
stop("Current treatment not provided/incomplete.")
}
if (length(mu0_t + sigma0_t + N0_t) != 0) {
intent <- c(intent, "historical treatment")
# cat("Historical Treatment\n")
} else {
if (length(c(mu0_t, sigma0_t, N0_t)) > 0) {
if (is.null(mu0_t) == TRUE) {
cat("mu0_t missing\n")
}
if (is.null(sigma0_t) == TRUE) {
cat("sigma0_t missing\n")
}
if (is.null(N0_t) == TRUE) {
cat("N0_t missing\n")
}
stop("Historical treatment incomplete.")
}
}
if (length(mu_c + sigma_c + N_c) != 0) {
intent <- c(intent, "current control")
# cat("Current Control\n")
} else {
if (length(c(mu_c, sigma_c, N_c)) > 0) {
if (is.null(mu_c) == TRUE) {
cat("mu_c missing\n")
}
if (is.null(sigma_c) == TRUE) {
cat("sigma_c missing\n")
}
if (is.null(N_c) == TRUE) {
cat("N_c missing\n")
}
stop("Current control not provided/incomplete.")
}
}
if (length(mu0_c + sigma0_c + N0_c) != 0) {
intent <- c(intent, "historical control")
# cat("Historical Contro\nl")
} else {
if (length(c(mu0_c, sigma0_c, N0_c)) > 0) {
if (is.null(mu0_c) == TRUE) {
cat("mu0_c missing\n")
}
if (is.null(sigma0_c) == TRUE) {
cat("sigma0_c missing\n")
}
if (is.null(N0_c) == TRUE) {
cat("N0_c missing\n")
}
stop("Historical Control not provided/incomplete.")
}
}
if (!is.null(N_c) | !is.null(N0_c)) {
arm2 <- TRUE
} else {
arm2 <- FALSE
}
# Check that discount_function is input correctly
all_functions <- c("weibull", "scaledweibull", "identity")
function_match <- match(discount_function, all_functions)
if (is.na(function_match)) {
stop("discount_function input incorrectly.")
}
##############################################################################
# Quick check, if alpha_max, weibull_scale, or weibull_shape have length 1,
# repeat input twice
##############################################################################
if (length(alpha_max) == 1) {
alpha_max <- rep(alpha_max, 2)
}
if (length(weibull_scale) == 1) {
weibull_scale <- rep(weibull_scale, 2)
}
if (length(weibull_shape) == 1) {
weibull_shape <- rep(weibull_shape, 2)
}
################################################################################
# Results #
################################################################################
posterior_treatment <- posterior_normal(
mu = mu_t,
sigma = sigma_t,
N = N_t,
mu0 = mu0_t,
sigma0 = sigma0_t,
N0 = N0_t,
discount_function = discount_function,
alpha_max = alpha_max[1],
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
weibull_scale = weibull_scale[1],
weibull_shape = weibull_shape[1],
method = method
)
if (arm2) {
posterior_control <- posterior_normal(
mu = mu_c,
sigma = sigma_c,
N = N_c,
mu0 = mu0_c,
sigma0 = sigma0_c,
N0 = N0_c,
discount_function = discount_function,
alpha_max = alpha_max[2],
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
weibull_scale = weibull_scale[2],
weibull_shape = weibull_shape[2],
method = method
)
} else {
posterior_control <- NULL
}
args1 <- list(
mu_t = mu_t,
sigma_t = sigma_t,
N_t = N_t,
mu0_t = mu0_t,
sigma0_t = sigma0_t,
N0_t = N0_t,
mu_c = mu_c,
sigma_c = sigma_c,
N_c = N_c,
mu0_c = mu0_c,
sigma0_c = sigma0_c,
N0_c = N0_c,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
number_mcmc = number_mcmc,
method = method,
arm2 = arm2,
intent = paste(intent,
collapse = ", ",
compare = compare
)
)
##############################################################################
### Create final (comparison) object
##############################################################################
if (!compare) {
final <- NULL
} else {
if (arm2) {
final <- list()
final$posterior <- posterior_treatment$posterior_mu - posterior_control$posterior_mu
} else {
final <- list()
final$posterior <- posterior_treatment$posterior_mu
}
}
me <- list(
posterior_treatment = posterior_treatment,
posterior_control = posterior_control,
final = final,
args1 = args1
)
class(me) <- "bdpnormal"
return(me)
}
)
################################################################################
# Normal posterior estimation
# 1) Estimate the discount function (if current+historical data both present)
# 2) Estimate the posterior of the augmented data
################################################################################
posterior_normal <- function(mu, sigma, N, mu0, sigma0, N0, discount_function,
alpha_max, fix_alpha, number_mcmc, weibull_scale,
weibull_shape, method) {
# Compute posterior(s) of current (flat) and historical (prior) data
# with non-informative priors
# Current data:
if (!is.null(N)) {
posterior_flat_sigma2 <- 1 / rgamma(number_mcmc, (N - 1) / 2, ((N - 1) * sigma^2) / 2)
s <- (posterior_flat_sigma2 / ((N - 1) + 1))^0.5
posterior_flat_mu <- rnorm(number_mcmc, mu, s)
} else {
posterior_flat_mu <- posterior_flat_sigma2 <- NULL
}
# Historical data:
if (!is.null(N0)) {
prior_sigma2 <- 1 / rgamma(number_mcmc, (N0 - 1) / 2, ((N0 - 1) * sigma0^2) / 2)
s0 <- (prior_sigma2 / ((N0 - 1) + 1))^0.5
prior_mu <- rnorm(number_mcmc, mu0, s0)
} else {
prior_mu <- prior_sigma2 <- NULL
}
##############################################################################
# Discount function
##############################################################################
### Compute stochastic comparison and alpha discount only if both
### N and N0 are present (i.e., current & historical data are present)
if (!is.null(N) & !is.null(N0)) {
### Test of model vs real
if (method == "fixed") {
p_hat <- mean(posterior_flat_mu < prior_mu) # larger is higher failure
p_hat <- 2 * ifelse(p_hat > 0.5, 1 - p_hat, p_hat)
} else if (method == "mc") {
Z <- abs(posterior_flat_mu - prior_mu) / sqrt(s^2 + s0^2)
p_hat <- 2 * (1 - pnorm(Z))
}
### Number of effective sample size given shape and scale discount function
if (fix_alpha == TRUE) {
alpha_discount <- alpha_max
} else {
# Compute alpha discount based on distribution
if (discount_function == "weibull") {
alpha_discount <- pweibull(p_hat,
shape = weibull_shape,
scale = weibull_scale
) * alpha_max
} else if (discount_function == "scaledweibull") {
max_p <- pweibull(1, shape = weibull_shape, scale = weibull_scale)
alpha_discount <- pweibull(p_hat,
shape = weibull_shape,
scale = weibull_scale
) * alpha_max / max_p
} else if (discount_function == "identity") {
alpha_discount <- p_hat * alpha_max
}
}
} else {
alpha_discount <- NULL
p_hat <- NULL
}
##############################################################################
# Posterior augmentation
# - If current or historical data are missing, this will not augment but
# will return the posterior of the non-missing data (with flat prior)
##############################################################################
### If only the historical data is present, compute posterior on historical
if (is.null(N0) & !is.null(N)) {
posterior_sigma2 <- posterior_flat_sigma2
posterior_mu <- rnorm(number_mcmc, posterior_flat_mu, sqrt(posterior_sigma2))
} else if (!is.null(N0) & is.null(N)) {
posterior_sigma2 <- prior_sigma2
posterior_mu <- rnorm(number_mcmc, prior_mu, sqrt(posterior_sigma2))
} else if (!is.null(N0) & !is.null(N)) {
effective_N0 <- N0 * alpha_discount
posterior_mu0 <- prior_sigma2 * N * mu + posterior_flat_sigma2 * effective_N0 * mu0
posterior_mu0 <- posterior_mu0 / (N * prior_sigma2 + posterior_flat_sigma2 * effective_N0)
posterior_sigma2 <- posterior_flat_sigma2 * prior_sigma2
posterior_sigma2 <- posterior_sigma2 / (N * prior_sigma2 + posterior_flat_sigma2 * effective_N0)
posterior_mu <- rnorm(number_mcmc, posterior_mu0, sqrt(posterior_sigma2))
}
return(list(
alpha_discount = alpha_discount,
p_hat = p_hat,
posterior_mu = posterior_mu,
posterior_sigma2 = posterior_sigma2,
posterior_flat_mu = posterior_flat_mu,
posterior_flat_sigma2 = posterior_flat_sigma2,
prior_mu = prior_mu,
prior_sigma2 = prior_sigma2
))
}
graemeleehickey/bayesDP documentation built on Dec. 30, 2024, 9:32 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions posterior_normal
#' @title Bayesian Discount Prior: Gaussian mean values
#'
#' @description \code{bdpnormal} is used for estimating posterior samples from a
#' Gaussian outcome where an informative prior is used. The prior weight is
#' determined using a discount function. This code is modeled after the
#' methodologies developed in Haddad et al. (2017).
#'
#' @param mu_t scalar. Mean of the current treatment group.
#' @param sigma_t scalar. Standard deviation of the current treatment group.
#' @param N_t scalar. Number of observations of the current treatment group.
#' @param mu0_t scalar. Mean of the historical treatment group.
#' @param sigma0_t scalar. Standard deviation of the historical treatment group.
#' @param N0_t scalar. Number of observations of the historical treatment group.
#' @param mu_c scalar. Mean of the current control group.
#' @param sigma_c scalar. Standard deviation of the current control group.
#' @param N_c scalar. Number of observations of the current control group.
#' @param mu0_c scalar. Mean of the historical control group.
#' @param sigma0_c scalar. Standard deviation of the historical control group.
#' @param N0_c scalar. Number of observations of the historical control group.
#' @param discount_function character. Specify the discount function to use.
#' Currently supports \code{weibull}, \code{scaledweibull}, and
#' \code{identity}. The discount function \code{scaledweibull} scales the
#' output of the Weibull CDF to have a max value of 1. The \code{identity}
#' discount function uses the posterior probability directly as the discount
#' weight. Default value is "\code{identity}".
#' @param alpha_max scalar. Maximum weight the discount function can apply.
#' Default is 1. For a two-arm trial, users may specify a vector of two values
#' where the first value is used to weight the historical treatment group and
#' the second value is used to weight the historical control group.
#' @param fix_alpha logical. Fix alpha at alpha_max? Default value is FALSE.
#' @param weibull_shape scalar. Shape parameter of the Weibull discount function
#' used to compute alpha, the weight parameter of the historical data. Default
#' value is 3. For a two-arm trial, users may specify a vector of two values
#' where the first value is used to estimate the weight of the historical
#' treatment group and the second value is used to estimate the weight of the
#' historical control group. Not used when \code{discount_function} =
#' "identity".
#' @param weibull_scale scalar. Scale parameter of the Weibull discount function
#' used to compute alpha, the weight parameter of the historical data. Default
#' value is 0.135. For a two-arm trial, users may specify a vector of two
#' values where the first value is used to estimate the weight of the
#' historical treatment group and the second value is used to estimate the
#' weight of the historical control group. Not used when
#' \code{discount_function} = "identity".
#' @param number_mcmc scalar. Number of Monte Carlo simulations. Default is
#' 10000.
#' @param method character. Analysis method with respect to estimation of the
#' weight parameter alpha. Default method "\code{mc}" estimates alpha for each
#' Monte Carlo iteration. Alternate value "\code{fixed}" estimates alpha once
#' and holds it fixed throughout the analysis. See the the \code{bdpnormal}
#' vignette \cr \code{vignette("bdpnormal-vignette", package="bayesDP")} for
#' more details.
#' @param compare logical. Should a comparison object be included in the fit?
#' For a one-arm analysis, the comparison object is simply the posterior chain
#' of the treatment group parameter. For a two-arm analysis, the comparison
#' object is the posterior chain of the treatment effect that compares
#' treatment and control. If \code{compare=TRUE}, the comparison object is
#' accessible in the \code{final} slot, else the \code{final} slot is
#' \code{NULL}. Default is \code{TRUE}.
#' @details \code{bdpnormal} uses a two-stage approach for determining the
#' strength of historical data in estimation of a mean outcome. In the first
#' stage, a \emph{discount function} is used that that defines the maximum
#' strength of the historical data and discounts based on disagreement with
#' the current data. Disagreement between current and historical data is
#' determined by stochastically comparing the respective posterior
#' distributions under noninformative priors. With Gaussian data, the
#' comparison is the proability (\code{p}) that the current mean is less than
#' the historical mean. The comparison metric \code{p} is then input into the
#' discount function and the final strength of the historical data is returned
#' (alpha).
#'
#' In the second stage, posterior estimation is performed where the discount
#' function parameter, \code{alpha}, is used incorporated in all posterior
#' estimation procedures.
#'
#' To carry out a single arm (OPC) analysis, data for the current treatment
#' (\code{mu_t}, \code{sigma_t}, and \code{N_t}) and historical treatment
#' (\code{mu0_t}, \code{sigma0_t}, and \code{N0_t}) must be input. The results
#' are then based on the posterior distribution of the current data augmented
#' by the historical data.
#'
#' To carry our a two-arm (RCT) analysis, data for the current treatment and
#' at least one of current or historical control data must be input. The
#' results are then based on the posterior distribution of the difference
#' between current treatment and control, augmented by available historical
#' data.
#'
#' For more details, see the \code{bdpnormal} vignette: \cr
#' \code{vignette("bdpnormal-vignette", package="bayesDP")}
#'
#' @return \code{bdpnormal} returns an object of class "bdpnormal". The
#' functions \code{\link[=summary,bdpnormal-method]{summary}} and
#' \code{\link[=print,bdpnormal-method]{print}} are used to obtain and print a
#' summary of the results, including user inputs. The
#' \code{\link[=plot,bdpnormal-method]{plot}} function displays visual outputs
#' as well.
#'
#' An object of class \code{bdpnormal} is a list containing at least
#' the following components:
#' \describe{
#' \item{\code{posterior_treatment}}{
#' list. Entries contain values related to the treatment group:
#' \itemize{
#' \item \code{alpha_discount}
#' numeric. Alpha value, the weighting parameter of the historical data.
#' \item \code{p_hat}
#' numeric. The posterior probability of the stochastic comparison
#' between the current and historical data.
#' \item \code{posterior_mu}
#' vector. A vector of length \code{number_mcmc} containing the posterior
#' mean of the treatment group. If historical treatment data is present,
#' the posterior incorporates the weighted historical data.
#' \item \code{posterior_sigma2}
#' vector. A vector of length \code{number_mcmc} containing the posterior
#' variance of the treatment group. If historical treatment data is present,
#' the posterior incorporates the weighted historical data.
#' \item \code{posterior_flat_mu}
#' vector. A vector of length \code{number_mcmc} containing
#' Monte Carlo samples of the mean of the current treatment group
#' under a flat/non-informative prior, i.e., no incorporation of the
#' historical data.
#' \item \code{posterior_flat_sigma2}
#' vector. A vector of length \code{number_mcmc} containing
#' Monte Carlo samples of the standard deviation of the current treatment group
#' under a flat/non-informative prior, i.e., no incorporation of the
#' historical data.
#' \item \code{prior_mu}
#' vector. If historical treatment data is present, a vector of length
#' \code{number_mcmc} containing Monte Carlo samples of the mean
#' of the historical treatment group under a flat/non-informative prior.
#' \item\code{prior_sigma2}
#' vector. If historical treatment data is present, a vector of length
#' \code{number_mcmc} containing Monte Carlo samples of the standard deviation
#' of the historical treatment group under a flat/non-informative prior.
#' }
#' }
#' \item{\code{posterior_control}}{
#' list. Similar entries as \code{posterior_treament}. Only present if a
#' control group is specified.
#' }
#' \item{\code{final}}{
#' list. Contains the final comparison object, dependent on the analysis type:
#' \itemize{
#' \item One-arm analysis:
#' vector. Posterior chain of the mean.
#' \item Two-arm analysis:
#' vector. Posterior chain of the mean difference comparing treatment and
#' control groups.
#' }
#' }
#' \item{\code{args1}}{
#' list. Entries contain user inputs. In addition, the following elements
#' are output:
#' \itemize{
#' \item \code{arm2}
#' binary indicator. Used internally to indicate one-arm or two-arm
#' analysis.
#' \item \code{intent}
#' character. Denotes current/historical status of treatment and
#' control groups.
#' }
#' }
#' }
#'
#'
#' @seealso \code{\link[=summary,bdpnormal-method]{summary}},
#' \code{\link[=print,bdpnormal-method]{print}},
#' and \code{\link[=plot,bdpnormal-method]{plot}} for details of each of the
#' supported methods.
#'
#' @references
#' Haddad, T., Himes, A., Thompson, L., Irony, T., Nair, R. MDIC Computer
#' Modeling and Simulation working group.(2017) Incorporation of stochastic
#' engineering models as prior information in Bayesian medical device trials.
#' \emph{Journal of Biopharmaceutical Statistics}, 1-15.
#'
#' @examples
#' # One-arm trial (OPC) example
#' fit <- bdpnormal(
#' mu_t = 30, sigma_t = 10, N_t = 50,
#' mu0_t = 32, sigma0_t = 10, N0_t = 50,
#' method = "fixed"
#' )
#' summary(fit)
#' \dontrun{
#' plot(fit)
#' }
#'
#' # Two-arm (RCT) example
#' fit2 <- bdpnormal(
#' mu_t = 30, sigma_t = 10, N_t = 50,
#' mu0_t = 32, sigma0_t = 10, N0_t = 50,
#' mu_c = 25, sigma_c = 10, N_c = 50,
#' mu0_c = 25, sigma0_c = 10, N0_c = 50,
#' method = "fixed"
#' )
#' summary(fit2)
#' \dontrun{
#' plot(fit2)
#' }
#'
#' @rdname bdpnormal
#' @import methods
#' @importFrom stats sd density is.empty.model median model.offset
#' model.response pweibull quantile rbeta rgamma rnorm var vcov
#' @importFrom ggplot2 aes facet_wrap geom_hline geom_line geom_vline ggplot
#' theme_bw xlab ylab
#' @aliases bdpnormal-method
#' @aliases bdpnormal,ANY-method
#' @export bdpnormal
bdpnormal <- setClass("bdpnormal", slots = c(
posterior_treatment = "list",
posterior_control = "list",
final = "list",
args1 = "list"
))
setGeneric(
"bdpnormal",
function(mu_t = NULL,
sigma_t = NULL,
N_t = NULL,
mu0_t = NULL,
sigma0_t = NULL,
N0_t = NULL,
mu_c = NULL,
sigma_c = NULL,
N_c = NULL,
mu0_c = NULL,
sigma0_c = NULL,
N0_c = NULL,
discount_function = "identity",
alpha_max = 1,
fix_alpha = FALSE,
weibull_scale = 0.135,
weibull_shape = 3,
number_mcmc = 10000,
method = "mc",
compare = TRUE) {
standardGeneric("bdpnormal")
}
)
setMethod(
"bdpnormal",
signature(),
function(mu_t = NULL,
sigma_t = NULL,
N_t = NULL,
mu0_t = NULL,
sigma0_t = NULL,
N0_t = NULL,
mu_c = NULL,
sigma_c = NULL,
N_c = NULL,
mu0_c = NULL,
sigma0_c = NULL,
N0_c = NULL,
discount_function = "identity",
alpha_max = 1,
fix_alpha = FALSE,
weibull_scale = 0.135,
weibull_shape = 3,
number_mcmc = 10000,
method = "mc",
compare = TRUE) {
################################################################################
# Check Input #
################################################################################
intent <- c()
if (length(mu_t + sigma_t + N_t) != 0) {
intent <- c(intent, "current treatment")
# cat("Current Treatment\n")
} else {
if (is.null(mu_t) == TRUE) {
cat("mu_t missing\n")
}
if (is.null(sigma_t) == TRUE) {
cat("sigma_t missing\n")
}
if (is.null(N_t) == TRUE) {
cat("N_t missing\n")
}
stop("Current treatment not provided/incomplete.")
}
if (length(mu0_t + sigma0_t + N0_t) != 0) {
intent <- c(intent, "historical treatment")
# cat("Historical Treatment\n")
} else {
if (length(c(mu0_t, sigma0_t, N0_t)) > 0) {
if (is.null(mu0_t) == TRUE) {
cat("mu0_t missing\n")
}
if (is.null(sigma0_t) == TRUE) {
cat("sigma0_t missing\n")
}
if (is.null(N0_t) == TRUE) {
cat("N0_t missing\n")
}
stop("Historical treatment incomplete.")
}
}
if (length(mu_c + sigma_c + N_c) != 0) {
intent <- c(intent, "current control")
# cat("Current Control\n")
} else {
if (length(c(mu_c, sigma_c, N_c)) > 0) {
if (is.null(mu_c) == TRUE) {
cat("mu_c missing\n")
}
if (is.null(sigma_c) == TRUE) {
cat("sigma_c missing\n")
}
if (is.null(N_c) == TRUE) {
cat("N_c missing\n")
}
stop("Current control not provided/incomplete.")
}
}
if (length(mu0_c + sigma0_c + N0_c) != 0) {
intent <- c(intent, "historical control")
# cat("Historical Contro\nl")
} else {
if (length(c(mu0_c, sigma0_c, N0_c)) > 0) {
if (is.null(mu0_c) == TRUE) {
cat("mu0_c missing\n")
}
if (is.null(sigma0_c) == TRUE) {
cat("sigma0_c missing\n")
}
if (is.null(N0_c) == TRUE) {
cat("N0_c missing\n")
}
stop("Historical Control not provided/incomplete.")
}
}
if (!is.null(N_c) | !is.null(N0_c)) {
arm2 <- TRUE
} else {
arm2 <- FALSE
}
# Check that discount_function is input correctly
all_functions <- c("weibull", "scaledweibull", "identity")
function_match <- match(discount_function, all_functions)
if (is.na(function_match)) {
stop("discount_function input incorrectly.")
}
##############################################################################
# Quick check, if alpha_max, weibull_scale, or weibull_shape have length 1,
# repeat input twice
##############################################################################
if (length(alpha_max) == 1) {
alpha_max <- rep(alpha_max, 2)
}
if (length(weibull_scale) == 1) {
weibull_scale <- rep(weibull_scale, 2)
}
if (length(weibull_shape) == 1) {
weibull_shape <- rep(weibull_shape, 2)
}
################################################################################
# Results #
################################################################################
posterior_treatment <- posterior_normal(
mu = mu_t,
sigma = sigma_t,
N = N_t,
mu0 = mu0_t,
sigma0 = sigma0_t,
N0 = N0_t,
discount_function = discount_function,
alpha_max = alpha_max[1],
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
weibull_scale = weibull_scale[1],
weibull_shape = weibull_shape[1],
method = method
)
if (arm2) {
posterior_control <- posterior_normal(
mu = mu_c,
sigma = sigma_c,
N = N_c,
mu0 = mu0_c,
sigma0 = sigma0_c,
N0 = N0_c,
discount_function = discount_function,
alpha_max = alpha_max[2],
fix_alpha = fix_alpha,
number_mcmc = number_mcmc,
weibull_scale = weibull_scale[2],
weibull_shape = weibull_shape[2],
method = method
)
} else {
posterior_control <- NULL
}
args1 <- list(
mu_t = mu_t,
sigma_t = sigma_t,
N_t = N_t,
mu0_t = mu0_t,
sigma0_t = sigma0_t,
N0_t = N0_t,
mu_c = mu_c,
sigma_c = sigma_c,
N_c = N_c,
mu0_c = mu0_c,
sigma0_c = sigma0_c,
N0_c = N0_c,
discount_function = discount_function,
alpha_max = alpha_max,
fix_alpha = fix_alpha,
weibull_scale = weibull_scale,
weibull_shape = weibull_shape,
number_mcmc = number_mcmc,
method = method,
arm2 = arm2,
intent = paste(intent,
collapse = ", ",
compare = compare
)
)
##############################################################################
### Create final (comparison) object
##############################################################################
if (!compare) {
final <- NULL
} else {
if (arm2) {
final <- list()
final$posterior <- posterior_treatment$posterior_mu - posterior_control$posterior_mu
} else {
final <- list()
final$posterior <- posterior_treatment$posterior_mu
}
}
me <- list(
posterior_treatment = posterior_treatment,
posterior_control = posterior_control,
final = final,
args1 = args1
)
class(me) <- "bdpnormal"
return(me)
}
)
################################################################################
# Normal posterior estimation
# 1) Estimate the discount function (if current+historical data both present)
# 2) Estimate the posterior of the augmented data
################################################################################
posterior_normal <- function(mu, sigma, N, mu0, sigma0, N0, discount_function,
alpha_max, fix_alpha, number_mcmc, weibull_scale,
weibull_shape, method) {
# Compute posterior(s) of current (flat) and historical (prior) data
# with non-informative priors
# Current data:
if (!is.null(N)) {
posterior_flat_sigma2 <- 1 / rgamma(number_mcmc, (N - 1) / 2, ((N - 1) * sigma^2) / 2)
s <- (posterior_flat_sigma2 / ((N - 1) + 1))^0.5
posterior_flat_mu <- rnorm(number_mcmc, mu, s)
} else {
posterior_flat_mu <- posterior_flat_sigma2 <- NULL
}
# Historical data:
if (!is.null(N0)) {
prior_sigma2 <- 1 / rgamma(number_mcmc, (N0 - 1) / 2, ((N0 - 1) * sigma0^2) / 2)
s0 <- (prior_sigma2 / ((N0 - 1) + 1))^0.5
prior_mu <- rnorm(number_mcmc, mu0, s0)
} else {
prior_mu <- prior_sigma2 <- NULL
}
##############################################################################
# Discount function
##############################################################################
### Compute stochastic comparison and alpha discount only if both
### N and N0 are present (i.e., current & historical data are present)
if (!is.null(N) & !is.null(N0)) {
### Test of model vs real
if (method == "fixed") {
p_hat <- mean(posterior_flat_mu < prior_mu) # larger is higher failure
p_hat <- 2 * ifelse(p_hat > 0.5, 1 - p_hat, p_hat)
} else if (method == "mc") {
Z <- abs(posterior_flat_mu - prior_mu) / sqrt(s^2 + s0^2)
p_hat <- 2 * (1 - pnorm(Z))
}
### Number of effective sample size given shape and scale discount function
if (fix_alpha == TRUE) {
alpha_discount <- alpha_max
} else {
# Compute alpha discount based on distribution
if (discount_function == "weibull") {
alpha_discount <- pweibull(p_hat,
shape = weibull_shape,
scale = weibull_scale
) * alpha_max
} else if (discount_function == "scaledweibull") {
max_p <- pweibull(1, shape = weibull_shape, scale = weibull_scale)
alpha_discount <- pweibull(p_hat,
shape = weibull_shape,
scale = weibull_scale
) * alpha_max / max_p
} else if (discount_function == "identity") {
alpha_discount <- p_hat * alpha_max
}
}
} else {
alpha_discount <- NULL
p_hat <- NULL
}
##############################################################################
# Posterior augmentation
# - If current or historical data are missing, this will not augment but
# will return the posterior of the non-missing data (with flat prior)
##############################################################################
### If only the historical data is present, compute posterior on historical
if (is.null(N0) & !is.null(N)) {
posterior_sigma2 <- posterior_flat_sigma2
posterior_mu <- rnorm(number_mcmc, posterior_flat_mu, sqrt(posterior_sigma2))
} else if (!is.null(N0) & is.null(N)) {
posterior_sigma2 <- prior_sigma2
posterior_mu <- rnorm(number_mcmc, prior_mu, sqrt(posterior_sigma2))
} else if (!is.null(N0) & !is.null(N)) {
effective_N0 <- N0 * alpha_discount
posterior_mu0 <- prior_sigma2 * N * mu + posterior_flat_sigma2 * effective_N0 * mu0
posterior_mu0 <- posterior_mu0 / (N * prior_sigma2 + posterior_flat_sigma2 * effective_N0)
posterior_sigma2 <- posterior_flat_sigma2 * prior_sigma2
posterior_sigma2 <- posterior_sigma2 / (N * prior_sigma2 + posterior_flat_sigma2 * effective_N0)
posterior_mu <- rnorm(number_mcmc, posterior_mu0, sqrt(posterior_sigma2))
}
return(list(
alpha_discount = alpha_discount,
p_hat = p_hat,
posterior_mu = posterior_mu,
posterior_sigma2 = posterior_sigma2,
posterior_flat_mu = posterior_flat_mu,
posterior_flat_sigma2 = posterior_flat_sigma2,
prior_mu = prior_mu,
prior_sigma2 = prior_sigma2
))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.