Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

17_getBetaMix_dbetaMix #23

Merged
merged 29 commits into from
Nov 27, 2023
Merged
Show file tree
Hide file tree
Changes from 19 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@ export(dbetaMix)
export(dbetabinom)
export(dbetabinomMix)
export(dbetadiff)
export(getBetamixPost)
export(logit)
export(myPlot)
export(myPlotDiff)
Expand Down
82 changes: 35 additions & 47 deletions R/dbetabinom.R
Original file line number Diff line number Diff line change
Expand Up @@ -35,22 +35,19 @@ dbetabinom <- function(x, m, a, b, log = FALSE) {
}
}


#' Beta-Mixture-Binomial Density Function
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' Calculates the density function for a mixture of beta-binomial distributions.
#'
#' @typed x : numeric
#' number of successes.
#' @typed m : number
#' number of trials.
#' @inheritParams dbetabinom
#' @typed par : matrix
#' the beta parameters matrix, with K rows and 2 columns,
#' corresponding to the beta parameters of the K components.
#' @typed weights : numeric
#' the mixture weights of the beta mixture prior.
#' the mixture weights of the beta mixture prior of length K.
#' Each element corresponds to the row of beta parameters in `par`.
#' @typed log : flag
#' whether to return the log density value (not default).
#' @return The (log) density values of the mixture of beta-binomial distributions at `x`.
Expand All @@ -72,75 +69,66 @@ dbetabinomMix <- Vectorize(dbetabinomMix, vectorize.args = "x")

#' Compute Beta-Mixture-Binomial Posterior Distribution
#'
#' Computes the posterior parameters of a beta-mixture-binomial distribution.
#' A helper function that computes the posterior parameters of a beta-mixture-binomial distribution.
#'
#' @inheritParams dbetabinom
#' @inheritParams dbetabinomMix
#'
#' @param x number of successes
#' @param n number of patients
#' @param par the beta parameters matrix, with K rows and 2 columns,
#' corresponding to the beta parameters of the K components
#' @param weights the mixture weights
#' @return a list with the updated beta parameters and weights
#' @return A list with the updated beta parameters and weights.
#'
#' @importFrom stats dbeta dbinom
#'
#' @example examples/getBetamixPost.R
#' @export

getBetamixPost <- function(x, n, par, weights) {
## check the format
stopifnot(
is.matrix(par),
is.numeric(par),
identical(ncol(par), 2L),
all(par > 0),
identical(nrow(par), length(weights)),
all(weights > 0)
)

## renormalize weights
#' @keywords internal
h_getBetamixPost <- function(x, n, par, weights) {
assert_numeric(x, lower = 0, upper = n, finite = TRUE)
assert_numeric(n, lower = 0, finite = TRUE)
assert_matrix(par, min.rows = 1, max.cols = 2, mode = "numeric")
assert_numeric(weights, min = 0, len = nrow(par), finite = TRUE)
# We renormalize weights.
weights <- weights / sum(weights)

## now compute updated parameters
# We now compute updated parameters.
postPar <- par
postPar[, 1] <- postPar[, 1] + x
postPar[, 2] <- postPar[, 2] + n - x
postParProb <- postPar[, 1] / (postPar[, 1] + postPar[, 2])

## compute updated mixture probabilities
# We compute updated mixture probabilities.
tmp <- exp(
stats::dbinom(x, size = n, prob = postParProb, log = TRUE) +
stats::dbeta(postParProb, par[, 1], par[, 2], log = TRUE) -
stats::dbeta(postParProb, postPar[, 1], postPar[, 2], log = TRUE)
)

# We compute the updated weights of the posterior
postWeights <- weights * tmp / sum(weights * tmp)

return(list(
assert_numeric(postWeights)
list(
par = postPar,
weights = postWeights
))
)
}



#' Beta-mixture density function
#' Beta-Mixture density function
#'
#' Calculating beta-mixture density values of support `x`.
audreyyeoCH marked this conversation as resolved.
Show resolved Hide resolved
#'
#' Note that `x` can be a vector.
#' @description `r lifecycle::badge("experimental")`
#'
#' The beta-mixture density can be calculated with the combination of K set of beta parameters and K
#' length of weights.
audreyyeoCH marked this conversation as resolved.
Show resolved Hide resolved
#'
#' @inheritParams dbetabinom
#' @inheritParams dbetabinomMix
#' @typed log : flag
#' whether log values of the beta-mixture density function are returned.
#'
#' @param x the abscissa
#' @param par the beta parameters matrix, with K rows and 2 columns,
#' corresponding to the beta parameters of the K components
#' @param weights the mixture weights of the beta mixture prior
#' @param log return the log value? (not default)
#' @return the (log) density values
#'
#' @export
dbetaMix <- function(x, par, weights, log = FALSE) {
audreyyeoCH marked this conversation as resolved.
Show resolved Hide resolved
ret <- sum(weights * dbeta(x, par[, 1], par[, 2]))
if (log) {
return(log(ret))
log(ret)
} else {
return(ret)
ret
}
}
dbetaMix <- Vectorize(dbetaMix, vectorize.args = "x")
Expand Down
2 changes: 1 addition & 1 deletion R/oc2.R
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ oc2 <- function(method =
## if prior weights of the beta mixture are not supplied
if (missing(weights)) {
weights <- rep(1, nrow(parE))
## (don't need to be normalized, this is done in getBetamixPost)
## (don't need to be normalized, this is done in h_getBetamixPost)
}

## allocation to active and control arms:
Expand Down
2 changes: 1 addition & 1 deletion R/oc3.R
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ oc3 <- function(method =
## if prior weights of the beta mixture are not supplied
if (missing(weights)) {
weights <- rep(1, nrow(parE))
## (don't need to be normalized, this is done in getBetamixPost)
## (don't need to be normalized, this is done in h_getBetamixPost)
}

## allocation to active and control arms:
Expand Down
4 changes: 2 additions & 2 deletions R/postprob.R
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ postprobBeta <- function(x, n, p, a = 1, b = 1) {
#' The mixture weights of the beta mixture prior. Default are
#' uniform weights across mixture components.
#' @typed betamixPost : matrix
#' optional result of `[getBetamixPost()]` in order
#' optional result of `[h_getBetamixPost()]` in order
#' to speed up the computations. If supplied, this is directly used, bypassing
#' the other arguments (except `p` and `log.p` of course).
#' @typed log.p : number
Expand All @@ -85,7 +85,7 @@ postprob <- function(x, n, p, parE = c(1, 1), weights, betamixPost, log.p = FALS
if (missing(weights)) {
weights <- rep(1, nrow(parE))
}
betamixPost <- getBetamixPost(
betamixPost <- h_getBetamixPost(
x = x,
n = n,
par = parE,
Expand Down
4 changes: 2 additions & 2 deletions R/postprobDist.R
Original file line number Diff line number Diff line change
Expand Up @@ -178,8 +178,8 @@ postprobDist <- function(x,
assert_numeric(weightsS, lower = 0, finite = TRUE)
assert_numeric(parE, lower = 0, finite = TRUE)
assert_numeric(parS, lower = 0, finite = TRUE)
activeBetamixPost <- getBetamixPost(x = x, n = n, par = parE, weights = weights)
controlBetamixPost <- getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS)
activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights)
controlBetamixPost <- h_getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS)
assert_names(names(activeBetamixPost), identical.to = c("par", "weights"))
assert_names(names(controlBetamixPost), identical.to = c("par", "weights"))
if (relativeDelta) {
Expand Down
4 changes: 2 additions & 2 deletions R/predprob.R
Original file line number Diff line number Diff line change
Expand Up @@ -59,11 +59,11 @@ predprob <- function(x, n, Nmax, p, thetaT, parE = c(1, 1),
## if prior weights of the beta mixture are not supplied
if (missing(weights)) {
weights <- rep(1, nrow(parE))
## (don't need to be normalized, this is done in getBetamixPost)
## (don't need to be normalized, this is done in h_getBetamixPost)
}

## now compute updated parameters
betamixPost <- getBetamixPost(x = x, n = n, par = parE, weights = weights)
betamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights)

py <- with(
betamixPost,
Expand Down
6 changes: 3 additions & 3 deletions R/predprobDist.R
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ predprobDist <- function(x, n,
## if prior weights of the beta mixture are not supplied
if (missing(weights)) {
weights <- rep(1, nrow(parE))
## (don't need to be normalized, this is done in getBetamixPost)
## (don't need to be normalized, this is done in h_getBetamixPost)
}

## if parS is a vector => situation where there is only one component
Expand All @@ -119,7 +119,7 @@ predprobDist <- function(x, n,

## now compute updated parameters for beta mixture distribution on the
## treatment proportion
activeBetamixPost <- getBetamixPost(x = x, n = n, par = parE, weights = weights)
activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights)

## now with the beta binomial mixture:
py <- with(
Expand Down Expand Up @@ -160,7 +160,7 @@ predprobDist <- function(x, n,
## counts in future SOC patients:
mS <- NmaxControl - nS

controlBetamixPost <- getBetamixPost(
controlBetamixPost <- h_getBetamixPost(
x = xS, n = nS, par = parS,
weights = weightsS
)
Expand Down
2 changes: 1 addition & 1 deletion examples/dbetabinom.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# dbetabinom --
dbetabinom(x = 2, m = 29, a = 0.2, b = 0.4, log = FALSE)

# Can also specify x as a vector.

dbetabinom(x = 1:28, m = 29, a = 0.2, b = 0.4, log = FALSE)
4 changes: 2 additions & 2 deletions examples/getBetamixPost.R
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
## example from Lee and Liu:
getBetamixPost(x = 16, n = 23, par = t(c(0.6, 0.4)), weights = 1)
h_getBetamixPost(x = 16, n = 23, par = t(c(0.6, 0.4)), weights = 1)

getBetamixPost(
h_getBetamixPost(
x = 16, n = 23,
par =
rbind(
Expand Down
4 changes: 1 addition & 3 deletions examples/postprob.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
# Example taken from Lee and Liu (2006)
#
# Example taken from Lee and Liu (2008) :
audreyyeoCH marked this conversation as resolved.
Show resolved Hide resolved
# We observed 16 successes out of 23 patients
# We set a threshold of 0.60
# Assume a beta(0.6,0.4) prior for P_E
Expand All @@ -9,7 +8,6 @@ postprob(x = 16, n = 23, p = 0.60, par = c(0.6, 0.4))

# We could instead specify a mixture prior
# 2 component beta mixture prior, i.e., P_E ~ 0.6*beta(0.6,0.4) + 0.4*beta(1,1) and Pr(P_E > p | data) = 0.823

postprob(
x = 16, n = 23, p = 0.60,
par =
Expand Down
21 changes: 14 additions & 7 deletions man/dbetaMix.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion man/dbetabinom.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

3 changes: 2 additions & 1 deletion man/dbetabinomMix.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

38 changes: 0 additions & 38 deletions man/getBetamixPost.Rd

This file was deleted.

24 changes: 24 additions & 0 deletions man/h_getBetamixPost.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading