R/pos2S.R

In RBesT: R Bayesian Evidence Synthesis Tools

#' Probability of Success for 2 Sample Design
#'
#' The \code{pos2S} function defines a 2 sample design (priors, sample
#' sizes & decision function) for the calculation of the probability
#' of success. A function is returned which calculates the calculates
#' the frequency at which the decision function is evaluated to 1 when
#' parameters are distributed according to the given distributions.
#'
#' @template args-boundary2S
#'
#' @details The \code{pos2S} function defines a 2 sample design and
#' returns a function which calculates its probability of success.
#' The probability of success is the frequency with which the decision
#' function is evaluated to 1 under the assumption of a given true
#' distribution of the data implied by a distirbution of the
#' parameters \eqn{\theta_1} and \eqn{\theta_2}.
#'
#' The calculation is analogous to the operating characeristics
#' \code{\link{oc2S}} with the difference that instead of assuming
#' known (point-wise) true parameter values a distribution is
#' specified for each parameter.
#'
#' Calling the \code{pos2S} function calculates the decision boundary
#' \eqn{D_1(y_2)} and returns a function which can be used to evaluate the
#' PoS for different predictive distributions. It is evaluated as
#'
#' \deqn{ \int\int\int f_2(y_2|\theta_2) \, p(\theta_2) \, F_1(D_1(y_2)|\theta_1) \, p(\theta_1) \, dy_2 d\theta_2 d\theta_1. }
#'
#' where \eqn{F} is the distribution function of the sampling
#' distribution and \eqn{p(\theta_1)} and \eqn{p(\theta_2)} specifies
#' the assumed true distribution of the parameters \eqn{\theta_1} and
#' \eqn{\theta_2}, respectively. Each distribution \eqn{p(\theta_1)}
#' and \eqn{p(\theta_2)} is a mixture distribution and given as the
#' \code{mix1} and \code{mix2} argument to the function.
#'
#' For example, in the binary case an integration of the predictive
#' distribution, the BetaBinomial, instead of the binomial
#' distribution will be performed over the data space wherever the
#' decision function is evaluated to 1. All other aspects of the
#' calculation are as for the 2-sample operating characteristics, see
#' \code{\link{oc2S}}.
#'
#' @return Returns a function which when called with two arguments
#' \code{mix1} and \code{mix2} will return the frequencies at
#' which the decision function is evaluated to 1. Each argument is
#' expected to be a mixture distribution representing the assumed true
#' distribution of the parameter in each group.
#'
#' @family design2S
#'
#' @examples
#'
#' # see ?decision2S for details of example
#' priorT <- mixnorm(c(1,   0, 0.001), sigma=88, param="mn")
#' priorP <- mixnorm(c(1, -49, 20   ), sigma=88, param="mn")
#' # the success criteria is for delta which are larger than some
#' # threshold value which is why we set lower.tail=FALSE
#' successCrit  <- decision2S(c(0.95, 0.5), c(0, 50), FALSE)
#'
#' # example interim outcome
#' postP_interim <- postmix(priorP, n=10, m=-50)
#' postT_interim <- postmix(priorT, n=20, m=-80)
#'
#' # assume that mean -50 / -80 were observed at the interim for
#' # placebo control(n=10) / active treatment(n=20) which gives
#' # the posteriors
#' postP_interim
#' postT_interim
#'
#' # then the PoS to succeed after another 20/30 patients is
#' pos_final <- pos2S(postP_interim, postT_interim, 20, 30, successCrit)
#'
#' pos_final(postP_interim, postT_interim)
#'
#' @export
pos2S <- function(prior1, prior2, n1, n2, decision, ...) UseMethod("pos2S")
#' @export
pos2S.default <- function(prior1, prior2, n1, n2, decision, ...) "Unknown density"

#' @templateVar fun pos2S
#' @template design2S-binomial
#' @export
pos2S.betaMix <- function(prior1, prior2, n1, n2, decision, eps, ...) {

    if(missing(eps) & (n1*n2 > 1e7)) {
        warning("Large sample space. Consider setting eps=1e-6.")
    }

    crit_y1 <- decision2S_boundary(prior1, prior2, n1, n2, decision, eps)

    lower.tail <- attr(decision, "lower.tail")

    approx_method <- !missing(eps)

    design_fun <- function(mix1, mix2) {

        ## for each 0:n1 of the possible outcomes, calculate the
        ## probability mass past the boundary (log space) weighted with
        ## the density as the value for 1 occures (due to theta1)

        pred_mix1 <- preddist(mix1, n=n1)
        pred_mix2 <- preddist(mix2, n=n2)

        assert_that(inherits(pred_mix1, "betaBinomialMix"))
        assert_that(inherits(pred_mix2, "betaBinomialMix"))

        ## now get the decision boundary in the needed range
        lim1 <- c(0,n1)
        lim2 <- c(0,n2)

        ## in case we use the approximate method, we restrict the
        ## evaluation of the decision function range
        if(approx_method) {
            lim1 <- qmix(pred_mix1, c(eps/2, 1-eps/2))
            lim2 <- qmix(pred_mix2, c(eps/2, 1-eps/2))
        }

        boundary <- crit_y1(lim2[1]:lim2[2], lim1=lim1)
        res <- rep(-Inf, times=length(boundary))

        for(i in lim2[1]:lim2[2]) {
            y2ind <- i - lim2[1] + 1
            if(boundary[y2ind] == -1) {
                ## decision was always 0
                res[y2ind] <- -Inf
            } else if(boundary[y2ind] == n1+1) {
                ## decision was always 1
                res[y2ind] <- 0
            } else {
                ## calculate for the predictive for dtheta1 the
                ## probability mass past (or before) the boundary
                res[y2ind] <- pmix(pred_mix1, boundary[y2ind], lower.tail=lower.tail, log.p=TRUE)
            }
            ## finally weight with the density according to the occurence
            ## of i due to theta2; the pmax avoids -Inf in a case of Prob==0
            res[y2ind] <- res[y2ind] + dmix(pred_mix2, i, log=TRUE)
        }
        exp(matrixStats::logSumExp(res))
    }
    design_fun
}


#' @templateVar fun pos2S
#' @template design2S-normal
#' @export
pos2S.normMix <- function(prior1, prior2, n1, n2, decision, sigma1, sigma2, eps=1e-6, Ngrid=10, ...) {
    ## distributions of the means of the data generating distributions
    ## for now we assume that the underlying standard deviation
    ## matches the respective reference scales

    if(missing(sigma1)) {
        sigma1 <- RBesT::sigma(prior1)
        message("Using default prior 1 reference scale ", sigma1)
    }
    assert_number(sigma1, lower=0)
    sigma(prior1) <- sigma1

    if(missing(sigma2)) {
        sigma2 <- RBesT::sigma(prior2)
        message("Using default prior 2 reference scale ", sigma2)
    }
    assert_number(sigma2, lower=0)
    sigma(prior2) <- sigma2

    crit_y1 <- decision2S_boundary(prior1, prior2, n1, n2, decision, sigma1, sigma2, eps, Ngrid)

    lower.tail <- attr(decision, "lower.tail")

    design_fun <- function(mix1, mix2) {

        ## get the predictive of the mean
        pred_mix1_mean <- preddist(mix1, n=n1, sigma=sigma1)
        if(n2 == 0) {
            ## gets ignored anyway
            pred_mix2_mean <- preddist(mix2, n=1, sigma=sigma2)
        } else {
            pred_mix2_mean <- preddist(mix2, n=n2, sigma=sigma2)
        }

        assert_that(inherits(pred_mix1_mean, "normMix"))
        assert_that(inherits(pred_mix2_mean, "normMix"))

        lim1 <- qmix(pred_mix1_mean, c(eps/2, 1-eps/2))
        lim2 <- qmix(pred_mix2_mean, c(eps/2, 1-eps/2))
        crit_y1(lim2, lim1)

        ##return(list(crit=crit_y1, m1=pred_dtheta1_mean, m2=pred_dtheta2_mean))

        if(n2 == 0) {
            mean_prior2 <- summary(prior2, probs=c())["mean"]
            return(pmix(pred_mix1_mean, crit_y1(mean_prior2), lower.tail=lower.tail))
        } else {
            return(integrate_density_log(function(x) pmix(pred_mix1_mean, crit_y1(x, lim1=lim1), lower.tail=lower.tail, log.p=TRUE), pred_mix2_mean, logit(eps/2), logit(1-eps/2)))
        }
    }

    design_fun
}

#' @templateVar fun pos2S
#' @template design2S-poisson
#' @export
pos2S.gammaMix <- function(prior1, prior2, n1, n2, decision, eps=1e-6, ...) {
    assert_that(likelihood(prior1) == "poisson")
    assert_that(likelihood(prior2) == "poisson")

    crit_y1 <- decision2S_boundary(prior1, prior2, n1, n2, decision, eps)

    lower.tail <- attr(decision, "lower.tail")

    design_fun <- function(mix1, mix2) {
        assert_that(likelihood(mix1) == "poisson")
        assert_that(likelihood(mix2) == "poisson")

        ## get the predictive of the sum
        pred_mix1_sum <- preddist(mix1, n=n1)
        pred_mix2_sum <- preddist(mix2, n=n2)

        assert_that(inherits(pred_mix1_sum, "gammaPoissonMix"))
        assert_that(inherits(pred_mix2_sum, "gammaPoissonMix"))

        lim1 <- qmix(pred_mix1_sum, c(eps/2, 1-eps/2))
        lim2 <- qmix(pred_mix2_sum, c(eps/2, 1-eps/2))

        ## force lower limit of lim1 to be 0 such that we will get and
        ## answer in most cases; performance wise it should be ok as
        ## we run a O(log(N)) search
        lim1[1] <- 0

        ## ensure that the boundaries are cached
        crit_y1(lim2, lim1=lim1)
        grid <- seq(lim2[1], lim2[2])
        exp(matrixStats::logSumExp(dmix(pred_mix2_sum, grid, log=TRUE)
                        + pmix(pred_mix1_sum, crit_y1(grid, lim1=lim1), lower.tail=lower.tail, log.p=TRUE)))
    }
    design_fun
}