Nothing
.predictionInterval_ <- function(thetao,
seo,
ser,
tau = 0,
conf.level = 0.95,
designPrior = c("predictive", "conditional", "EB")) {
stopifnot(is.numeric(thetao),
length(thetao) == 1,
is.finite(thetao),
is.numeric(seo),
length(seo) == 1,
is.finite(seo),
0 < seo,
is.numeric(ser),
length(ser) == 1,
is.finite(ser),
0 < ser,
is.numeric(tau),
length(tau) == 1,
is.finite(tau),
0 <= tau,
is.numeric(conf.level),
length(conf.level) == 1,
is.finite(conf.level),
0 <= conf.level, conf.level <= 1,
!is.null(designPrior))
designPrior <- match.arg(designPrior)
## determine parameters of predictive distribution of yr
if (designPrior == "conditional") {
mu <- thetao
sigma <- ser
} else if (designPrior == "predictive") {
mu <- thetao
sigma <- sqrt(seo^2 + ser^2 + 2 * tau^2)
} else { ## designPrior == "EB"
s <- pmax(1 - (seo^2 + tau^2) / thetao^2, 0)
mu <- s * thetao
sigma <- sqrt(s * (seo^2 + tau^2) + ser^2 + tau^2)
}
## compute prediction interval
lower <- stats::qnorm(p = (1 - conf.level) / 2, mean = mu, sd = sigma)
upper <- stats::qnorm(p = (1 + conf.level) / 2, mean = mu, sd = sigma)
result <- cbind(lower = lower, mean = mu, upper = upper)
return(result)
}
.predictionInterval__ <- Vectorize(.predictionInterval_)
#' Prediction interval for effect estimate of replication study
#'
#' Computes a prediction interval for the effect estimate of the replication study.
#' @param thetao Numeric vector of effect estimates from original studies.
#' @param seo Numeric vector of standard errors of the original effect estimates.
#' @param ser Numeric vector of standard errors of the replication effect estimates.
#' @param tau Between-study heterogeneity standard error.
#' Default is \code{0} (no heterogeneity).
#' Is only taken into account when \code{designPrior} is "predictive" or "EB".
#' @param conf.level The confidence level of the prediction intervals. Default is 0.95.
#' @param designPrior Either "predictive" (default), "conditional", or "EB".
#' If "EB", the contribution of the original study to the predictive distribution is
#' shrunken towards zero based on the evidence in the original study (with empirical Bayes).
#' @details This function computes a prediction interval and a mean estimate under a
#' specified predictive distribution of the replication effect estimate. Setting
#' \code{designPrior = "conditional"} is not recommended since this ignores the
#' uncertainty of the original effect estimate. See Patil, Peng, and Leek (2016)
#' and Pawel and Held (2020) for details.
#' @return A data frame with the following columns
#' \item{lower}{Lower limit of prediction interval,}
#' \item{mean}{Mean of predictive distribution,}
#' \item{upper}{Upper limit of prediction interval.}
#' @details \code{predictionInterval} is the vectorized version of \code{.predictionInterval_}.
#' \code{\link[base]{Vectorize}} is used to vectorize the function.
#' @references
#' Patil, P., Peng, R. D., Leek, J. T. (2016).
#' What should researchers expect when they replicate studies? A statistical view of
#' replicability in psychological science. \emph{Perspectives on Psychological Science},
#' \bold{11}, 539-544. \doi{10.1177/1745691616646366}
#'
#' Pawel, S., Held, L. (2020). Probabilistic forecasting of replication studies.
#' \emph{PLoS ONE}. \bold{15}, e0231416. \doi{10.1371/journal.pone.0231416}
#' @author Samuel Pawel
#' @examples
#' predictionInterval(thetao = c(1.5, 2, 5), seo = 1, ser = 0.5, designPrior = "EB")
#'
#' # compute prediction intervals for replication projects
#' data("RProjects", package = "ReplicationSuccess")
#' parOld <- par(mfrow = c(2, 2))
#' for (p in unique(RProjects$project)) {
#' data_project <- subset(RProjects, project == p)
#' PI <- predictionInterval(thetao = data_project$fiso, seo = data_project$se_fiso,
#' ser = data_project$se_fisr)
#' PI <- tanh(PI) # transforming back to correlation scale
#' within <- (data_project$rr < PI$upper) & (data_project$rr > PI$lower)
#' coverage <- mean(within)
#' color <- ifelse(within == TRUE, "#333333B3", "#8B0000B3")
#' study <- seq(1, nrow(data_project))
#' plot(data_project$rr, study, col = color, pch = 20,
#' xlim = c(-0.5, 1), xlab = expression(italic(r)[r]),
#' main = paste0(p, ": ", round(coverage*100, 1), "% coverage"))
#' arrows(PI$lower, study, PI$upper, study, length = 0.02, angle = 90,
#' code = 3, col = color)
#' abline(v = 0, lty = 3)
#' }
#' par(parOld)
#' @export
predictionInterval <- function(thetao,
seo,
ser,
tau = 0,
conf.level = 0.95,
designPrior = "predictive") {
res <- .predictionInterval__(thetao = thetao, seo = seo, ser = ser, tau = tau,
conf.level = conf.level, designPrior = designPrior)
res <- matrix(res, ncol = 3, byrow = TRUE)
res <- data.frame(res)
colnames(res) <- c("lower", "mean", "upper")
res
}
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.