Nothing
#' 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
}
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.