R/pph.R In BayesCombo: Bayesian Evidence Combination

#' @title Calculates the posterior probability of hypotheses for one study
#'
#' @description The function takes a single effect size and its standard error
#' and calculates the posterior probability of each hypothesis (H<: the effect
#' size is less than 0, H0: the effect size is zero, or H>: the effect size is
#' greater than zero).
#'
#' @details Effect sizes could be a mean difference between groups, regression
#' slope, odds ratio, or other values provided by statistical models. The
#' standard errors represent the uncertainty in the effect size and are provided
#' by most statistical software.
#'
#' Two types of priors need to be specified. The first is the prior for the
#' effect size, which is given by a mean (usually zero) and standard error,
#' which can be specified by the user or calculated automatically. Second, prior
#' probabilities for each hypothesis need to be provided, and an equal
#' probability of 1/3 is used by default.
#'
#' The null hypothesis (usually effect size = 0) can be specified as either a
#' point null or as a range. A point null tests if the effect size is
#' \emph{exactly} zero while the range null tests if the effect size close
#' enough to zero for practical purposes. 'Close enough' is defined by the user
#' as a range on either size of zero.
#'
#' To combine multiple effect sizes use the \code{ev.combo} function.
#'
#' @param beta Effect size.
#'
#' @param se.beta Standard error for the effect.
#'
#' @param beta0 A prior value for the effect size. Default is zero.
#'
#' @param se0 A prior standard error for the effect size. Default is \code{NULL}
#' and is calculated automatically.
#'
#' @param ci Is used to calculate the prior standard error if \code{se0 =
#' NULL}. The default value of 99 calculates the prior standard error so that
#' the 99% confidence intervals of the prior distribution are aligned with the
#' largest (furthest from zero) confidence interval of the data distribution.
#'
#' @param se.mult Standard error multiplier used to increase or decrease the
#' prior SE and used in conjunction with \code{ci} when \code{se0 = NULL}.
#'
#'
#' @param H0 A vector of length two that defines the null hypothesis. If the
#' values are identical (e.g. \code{H0 = c(0,0)}) a point null is used,
#' otherwise the null is defined as the range between the lower and upper value.
#'
#' @param H.priors Prior hypothesis probabilities; default is an equal
#' probability of 1/3, and they are specified in the following order:
#' \code{H<0}, \code{H=0}, \code{H>0}.
#'
#' @param scale Logical. Whether to scale the effect size by its standard
#' error. Standardising has no effect on the calculations but standardised
#' effect sizes may be easier to compare in a forest plot.
#'
#' probabilities. Adjusting prevents a single study from having too much
#' influence on the results when combining multiple studies. For example, if the
#' probability for a hypothesis from one study is zero, then additional studies
#' cannot alter this probability (multiplying anything by zero is still zero).
#'
#' @param epsilon A small value that a posterior hypothesis probability must
#'
#' probability if \code{adjust = TRUE} and one of the posterior hypothesis
#' probabilities is less than \code{epsilon}. The PPHs are then re-scaled to sum
#' to one.
#'
#' @return Object of class \code{pph} which contains the posterior hypothesis
#' probabilities and other calculated values.
#'
#' @export
#' @examples
#' # library(labstats) # need to install from CRAN
#' # plot(time.immob ~ dose, data=fluoxetine)
#' # summary(lm(time.immob ~ dose, data=fluoxetine))
#' x <- pph(beta=-0.25200, se.beta=0.09913) # dose effect from above output
#' x

pph <- function(beta, se.beta, beta0 = 0, se0 = NULL, ci = 99,
se.mult = 1, H0 = c(0,0), H.priors = rep(1/3, 3),
scale = FALSE, adjust = FALSE, epsilon = 1e-6,

if (beta0 < H0[1] | beta0 > H0[2] ) {
stop("beta0 must lie between the lower and upper bounds of H0.")
}

if (length(H0) != 2) {
stop("H0 must have two values. Use c(0,0) for a point null at zero.")
}

if (any(se.beta <= 0)){
stop("se.beta must be greater than zero.")
}

if (scale) { # standardise
beta <- beta/se.beta
se.beta <- se.beta/se.beta
}

if (ci <= 0 | ci >= 100 ) {
stop("ci must lie between zero and one.")
}

if (is.null(se0)){ # calculate if not given
se0 <- prior.se(beta, se.beta, ci) * se.mult
}

post.b <- calc.post.beta(beta, se.beta, beta0, se0)
post.se <- calc.post.se(se.beta, se0)

# BFs
lt <- pnorm(H0[1], post.b, post.se) / pnorm(H0[1], beta0, se.beta)

gt <- pnorm(H0[2], post.b, post.se, lower.tail = FALSE) /
pnorm(H0[2], beta0, se.beta, lower.tail = FALSE)

if (H0[1] == H0[2]){ # point null

null <- dnorm(beta0, post.b, post.se) / dnorm(H0[1], beta0, se0)

} else { # range null

null <- (pnorm(H0[2], post.b, post.se) - pnorm(H0[1], post.b, post.se)) /
(pnorm(H0[2], beta0, se.beta) - pnorm(H0[1], beta0, se.beta))
}

bfs <-c("H<" = lt, "H0" = null, "H>" = gt)

# PPHs
pphs <- (bfs * H.priors) / sum(bfs * H.priors)
names(pphs) <- c("H<", "H0", "H>")

if (adjust & any(pphs < epsilon)){
}

return(
structure(list(beta = beta, se.beta = se.beta, beta0 = beta0,
ci = ci, se0 = se0, post.b = post.b,
post.se = post.se, H.priors = H.priors, BFs = bfs,
pphs = pphs),
class = "PPH")
)
}


Try the BayesCombo package in your browser

Any scripts or data that you put into this service are public.

BayesCombo documentation built on May 1, 2019, 7:50 p.m.