Nothing
#' @import graphics stats utils
to.Fx <- function(x) {
# A faster version of ecdf(x)(x)
rank(x, ties.method="max")/length(x)
}
#' Constructs a simultaneous credible band
#'
#' \code{sim.cred.band} returns a simultaneous band over a finite set of
#' covariate points given either a sample from the posterior of the
#' regression surface or a function \code{FUN(x, params)} and a sample from
#' the posterior of the parameters.
#'
#' If \code{design} is \code{NULL} (default), it is taken to be the identity
#' matrix of dimension \code{ncol(params)}, so that the rows of params
#' are treated as draws from the posterior \code{FUN(x, params)}.
#'
#' The 'asymptotic' method assumes that the marginal posteriors of
#' the \code{FUN(x, params)} are asymptotically normal and is usually
#' significantly faster and less memory-intensive than the 'quantile'
#' method, which makes no such assumption.
#'
#' @param params A numeric matrix whose rows are draws from the posterior
#' distribution of either the regression surface or the parameter vector.
#' @param design (Optional) A numeric matrix whose rows are covariate points
#' over which the band is to be constructed.
#' @param FUN (Optional) a function of the form \code{function(x, params)}
#' that takes a row of \code{design} and the entire \code{params} matrix and
#' returns a vector of the same length of \code{x} representing the regression
#' surface.
#' @param cred.level Numeric; the credible level.
#' @param method Either "asymptotic" (default) or "quantile"; see details.
#' @param sides One of "both" (default), "upper", or "lower".
#' Which bounds should be constructed?
#' @param est.FUN The function used to produce estimates of the regression
#' surface. Default is \code{mean}.
#' @param var.FUN The function used to quantify the variability of the
#' regression surface posterior. Default is \code{sd}.
#' @param point.estimate If not null, replaces the mean and sets the reference
#' around which the standard error is computed.
#' Useful for bootstrapping methods.
#' Treated as a row of the \code{params} matrix.
#' @param track A numeric vector of indices indicating which rows (default none)
#' of the design matrix should have the sample of the corresponding
#' \code{FUN(x, params)} returned.
#' @param verbose Logical (default \code{FALSE}); print progress?
#'
#' @return An object of class \code{sim.cred.band}, which contains:
#' \describe{
#' \item{\code{upper}}{A numeric vector of upper bounds.}
#' \item{\code{lower}}{A numeric vector of lower bounds.}
#' \item{\code{cred.level}}{As provided.}
#' \item{\code{method}}{As provided.}
#' \item{\code{sides}}{As provided.}
#' \item{\code{est}}{Posterior estimate of the regression surface.}
#' \item{\code{est.FUN}}{As provided.}
#' \item{\code{var}}{Summary of posterior variability of the regression
#' surface.}
#' \item{\code{var.FUN}}{As provided.}
#' \item{\code{W}}{An estimate of the extremal errors.}
#' \item{\code{W.crit}}{The critical quantile of W.}
#' \item{\code{trace}}{The posterior samples of the regression surface
#' indicated by the \code{track} argument.}
#' }
#'
#' @examples
#' ### Sample from regression surface posterior
#' reg.surf.sample <- matrix(rnorm(1000, mean=1:10), ncol=2, byrow=TRUE)
#' sim.cred.band(reg.surf.sample, cred.level=0.80)
#'
#' ### Parametric case
#' design <- cbind(1, 1:10)
#' params <- matrix(rnorm(200, mean=1:2), ncol=2, byrow=TRUE)
#' sim.cred.band(params, design)
#'
#' ### With custom function
#' params.sd <- cbind(1 / rgamma(100, 1), params)
#' FUN.sd <- function(x, params) { params[, -1] %*% t(x) / params[, 1] }
#' sim.cred.band(params.sd, design, FUN.sd)
#'
#' @export
sim.cred.band <- function(params,
design=NULL,
FUN=function(x, params) { params %*% t(x) },
cred.level=0.95,
method=c('asymptotic', 'quantile'),
sides=c('both', 'upper', 'lower'),
est.FUN=mean,
var.FUN=sd,
point.estimate=NULL,
track=numeric(0),
verbose=FALSE) {
# Validate and shape input
params <- data.matrix(params)
M <- nrow(params)
if (is.null(design)) {
N <- ncol(params)
nonpar <- TRUE
} else {
design <- data.matrix(design)
N <- nrow(design)
nonpar <- FALSE
}
est <- var <- numeric(N)
method <- method[1]
stopifnot(method %in% c('asymptotic', 'quantile'))
if (method == 'asymptotic') {
m <- numeric(N)
s <- numeric(N)
}
sides <- sides[1]
stopifnot(sides %in% c('both', 'upper', 'lower'))
sim.cred.band <- list(upper=rep(Inf, N),
lower=rep(-Inf, N),
cred.level=cred.level,
method=method,
sides=sides,
W.crit=NA,
W=NA,
est=rep(NA, N),
est.FUN=est.FUN,
var=rep(NA, N),
var.FUN=var.FUN,
trace=matrix(NA, nrow=M, ncol=length(track)))
colnames(sim.cred.band$trace) <- track
class(sim.cred.band) <- 'sim.cred.band'
# Compute W
W <- rep(-Inf, M)
# This is iterative to avoid storing an NxM matrix (or two)
for (i in 1:N) {
if (verbose && (i %% 100 == 0)) {
cat(i, "/", N, "\n")
}
if (nonpar) {
fx <- params[, i]
} else {
fx <- FUN(design[i, , drop=FALSE], params)
}
if (i %in% track) {
sim.cred.band$trace[, which(track == i)] <- fx
}
est[i] <- est.FUN(fx)
var[i] <- var.FUN(fx)
if (method == 'asymptotic') {
if (is.null(point.estimate)) {
m[i] <- mean(fx)
s[i] <- sd(fx)
} else {
if (nonpar) {
point.fx <- point.estimate[i]
} else {
point.fx <- FUN(design[i, , drop=FALSE], point.estimate)
}
m[i] <- point.fx
s[i] <- sqrt(mean((fx - point.fx) ^ 2))
}
if (sides == "both") {
z <- abs(fx - m[i]) / s[i]
} else if (sides == "upper") {
z <- (fx - m[i]) / s[i]
} else if (sides == "lower") {
z <- (m[i] - fx) / s[i]
}
} else {
if (sides == "both") {
Fx <- to.Fx(fx)
Gx <- 1 - to.Fx(-fx)
z <- pmax(1-Fx, Gx)
} else if (sides == "upper") {
z <- 1 - to.Fx(-fx)
} else if (sides == "lower") {
z <- 1 - to.Fx(fx)
}
}
W <- pmax(W, z)
}
sim.cred.band$W <- W
# Estimate W.crit
W.crit <- quantile(W, cred.level, type=1)
sim.cred.band$W.crit <- W.crit
# Compute bounds
if (method == 'asymptotic') {
if (sides == "both" | sides == "upper") {
sim.cred.band$upper <- m + W.crit * s
}
if (sides == "both" | sides == "lower") {
sim.cred.band$lower <- m - W.crit * s
}
} else {
if (nonpar) {
if (sides %in% c('both', 'upper')) {
sim.cred.band$upper <- apply(params, 2, function(fx) {
-quantile(-fx, prob=1-W.crit, type=1)
})
}
if (sides %in% c('both', 'lower')) {
sim.cred.band$lower <- apply(params, 2, function(fx) {
quantile(fx, prob=1-W.crit, type=1)
})
}
} else {
if (sides == 'both') {
bounds <- apply(design, 1, function(x, params) {
fx <- FUN(t(x), params)
upper <- -quantile(-fx, prob=1-W.crit, type=1)
lower <- quantile( fx, prob=1-W.crit, type=1)
c(lower, upper)
}, params=params)
sim.cred.band$lower <- bounds[1, ]
sim.cred.band$upper <- bounds[2, ]
} else if (sides == 'upper') {
sim.cred.band$upper <- apply(design, 1, function(x, params) {
-quantile(-fx, prob=1-W.crit, type=1)
}, params=params)
} else if (sides == 'lower') {
sim.cred.band$lower <- apply(design, 1, function(x, params) {
quantile( fx, prob=1-W.crit, type=1)
}, params=params)
}
}
}
sim.cred.band$est <- est
sim.cred.band$var <- var
sim.cred.band
}
#' Constructs a credible subset pair
#'
#' \code{credsubs} returns a credible subset pair over a finite set of
#' covariate points given either a sample from the posterior of the
#' regression surface or a function \code{FUN(x, params)} and a sample from
#' the posterior of the parameters.
#'
#' If design is NULL (default), it is taken to be the identity
#' matrix of dimension ncol(params), so that the rows of params
#' are treated as draws from the posterior FUN(x, params).
#'
#' The 'asymptotic' method assumes that the marginal posteriors of
#' the FUN(x, params) are asymptotically normal and is usually
#' significantly faster and less memory-intensive than the 'quantile'
#' method, which makes no such assumption.
#'
#' @param params A numeric matrix whose rows are draws from the posterior
#' distribution of either the regression surface or the
#' parameter vector.
#' @param design (Optional) A numeric matrix whose rows are covariate points
#' over which the band is to be constructed.
#' @param FUN (Optional) a function of the form \code{function(x, params)}
#' that takes a row of \code{design} and the entire \code{params}
#' matrix and returns a vector of the same length of \code{x}
#' representing the regression surface.
#' @param cred.level Numeric; the credible level.
#' @param threshold Numeric; the value of \code{FUN} above which
#' a covariate is included in the target subset.
#' @param method Either "asymptotic" (default) or "quantile"; see details.
#' @param step.down Logical (default \code{TRUE}); should the step-down
#' procedure be used?
#' @param sides One of "both" (default), "exclusive", or "inclusive".
#' Which bounds should be constructed?
#' @param est.FUN The function used to produce estimates of the regression
#' surface. Default is \code{mean}.
#' @param var.FUN The function used to quantify the variability of the
#' regression surface posterior. Default is \code{sd}.
#' @param point.estimate If not null, replaces the mean and sets the reference
#' around which the standard error is computed.
#' Useful for bootstrapping methods.
#' Treated as a row of the \code{params} matrix.
#' @param track A numeric vector of indices indicating which rows (default none)
#' of the design matrix should have the sample of the corresponding
#' \code{FUN(x, params)} returned.
#' @param verbose Logical (default \code{FALSE}); print progress?
#'
#' @return An object of class \code{credsubs}, which contains:
#' \describe{
#' \item{\code{exclusive}}{A logical vector indicating membership in
#' the exclusive credible subset.}
#' \item{\code{inclusive}}{A logical vector indicating membership in
#' the inclusive credible subset.}
#' \item{\code{cred.level}}{As provided.}
#' \item{\code{threshold}}{As provided.}
#' \item{\code{method}}{As provided.}
#' \item{\code{step.down}}{As provided.}
#' \item{\code{sides}}{As provided.}
#' \item{\code{est}}{Posterior estimate of the regression surface.}
#' \item{\code{est.FUN}}{As provided.}
#' \item{\code{var}}{Summary of posterior variability of the regression
#' surface.}
#' \item{\code{var.FUN}}{As provided.}
#' \item{\code{W}}{An estimate of the extremal errors.}
#' \item{\code{W.crit}}{The critical quantile of W.}
#' \item{\code{trace}}{The posterior samples of the regression surface
#' indicated by the \code{track} argument.}
#' }
#'
#' @examples
#' ### Sample from regression surface posterior
#' reg.surf.sample <- matrix(rnorm(1000, mean=1:10), ncol=2, byrow=TRUE)
#' credsubs(reg.surf.sample, cred.level=0.80)
#'
#' ### Parametric case
#' design <- cbind(1, 1:10)
#' params <- matrix(rnorm(200, mean=1:2), ncol=2, byrow=TRUE)
#' credsubs(params, design)
#'
#' ### With custom function
#' params.sd <- cbind(1 / rgamma(100, 1), params)
#' FUN.sd <- function(x, params) { params[, -1] %*% t(x) / params[, 1] }
#' credsubs(params.sd, design, FUN.sd, threshold=1)
#'
#' @export
credsubs <- function(params,
design=NULL,
FUN=function(x, params) { params %*% t(x) },
cred.level=0.95,
threshold=0,
method=c('asymptotic', 'quantile'),
step.down=TRUE,
sides=c('both', 'exclusive', 'inclusive'),
est.FUN=mean,
var.FUN=sd,
point.estimate=NULL,
track=numeric(0),
verbose=FALSE) {
# Validate and shape input
params <- data.matrix(params)
M <- nrow(params)
if (is.null(design)) {
N <- ncol(params)
nonpar <- TRUE
} else {
design <- data.matrix(design)
N <- nrow(design)
nonpar <- FALSE
}
if (verbose) {
cat("Computing credible subgroups over", N, "points using",
M, "posterior draws.\n")
}
method <- method[1]
stopifnot(method %in% c('asymptotic', 'quantile'))
if (method == 'asymptotic') {
m <- numeric(N)
s <- numeric(N)
}
sides <- sides[1]
stopifnot(sides %in% c('both', 'exclusive', 'inclusive'))
scb.sides <- ifelse(sides == 'both', 'both',
ifelse(sides == 'exclusive', 'lower',
'upper'))
credsubs <- list(exclusive=rep(NA, N),
inclusive=rep(NA, N),
cred.level=cred.level,
threshold=threshold,
method=method,
step.down=step.down,
sides=sides,
est=rep(NA, N),
est.FUN=est.FUN,
var=rep(NA, N),
var.FUN=var.FUN,
W=rep(NA, M),
W.crit=NA,
trace=matrix(NA, nrow=M, ncol=length(track)))
class(credsubs) <- 'credsubs'
credsubs$exclusive <- rep(FALSE, N)
credsubs$inclusive <- rep(TRUE, N)
test.set <- 1:N
reject.set <- numeric(0)
repeat {
test.set <- setdiff(test.set, reject.set)
if (nonpar) {
test.par <- test.set
} else {
test.par <- 1:ncol(params)
}
# Runs on first iteration only
if (length(test.set) == N) {
scb.track <- track
} else {
scb.track <- numeric(0)
}
sim.cred.band <- sim.cred.band(params=params[, test.par, drop=FALSE],
design=design[test.set, , drop=FALSE],
FUN=FUN,
cred.level=cred.level,
method=method,
sides=scb.sides,
est.FUN=est.FUN,
var.FUN=var.FUN,
point.estimate=point.estimate,
track=scb.track,
verbose=verbose)
over.set <- test.set[sim.cred.band$lower > threshold]
under.set <- test.set[sim.cred.band$upper < threshold]
credsubs$exclusive[over.set] <- TRUE
credsubs$inclusive[under.set] <- FALSE
credsubs$est[test.set] <- sim.cred.band$est
credsubs$var[test.set] <- sim.cred.band$var
credsubs$W.crit <- sim.cred.band$W.crit
credsubs$W <- sim.cred.band$W
reject.set <- union(over.set, under.set)
# Runs on first iteration only
if (length(test.set) == N) {
credsubs$trace <- sim.cred.band$trace
}
if (verbose) {
cat(length(test.set), "hypotheses tested,",
length(reject.set), "rejected.\n")
}
if (!step.down ||
length(reject.set) == 0 ||
length(test.set) == length(reject.set)) {
break
}
}
credsubs
}
#' Compute the maximum credible levels at which conclusions may be drawn
#'
#' For each covariate point, \code{credsubs.level} computes the maximum
#' credible level at which a conclusion may be drawn at each point, and
#' what that conclusion is.
#'
#' If design is NULL (default), it is taken to be the identity
#' matrix of dimension ncol(params), so that the rows of params
#' are treated as draws from the posterior FUN(x, params).
#'
#' The 'asymptotic' method assumes that the marginal posteriors of
#' the FUN(x, params) are asymptotically normal and is usually
#' significantly faster and less memory-intensive than the 'quantile'
#' method, which makes no such assumption.
#'
#' By default (\code{z.store = "ram"}), the maximum credible level computation
#' stores a potentially very large amount of intermediate computation results
#' in memory. If not enough memory is available, \code{z.store = "disk"}
#' uses the \code{ff} package to store the intermediate results on disk,
#' which can still be fairly quick if the storage is fast (e.g. a local SSD).
#' Alternatively, \code{z.store = "recompute"} discards the intermediate
#' results and recomputes whenever needed. This uses minimal memory, but
#' is usually the slowest option.
#'
#' @param params A numeric matrix whose rows are draws from the posterior
#' distribution of either the regression surface or the
#' parameter vector.
#' @param design (Optional) A numeric matrix whose rows are covariate points
#' over which the band is to be constructed.
#' @param FUN (Optional) a function of the form \code{function(x, params)}
#' that takes a row of \code{design} and the entire \code{params}
#' matrix and returns a vector of the same length of \code{x}
#' representing the regression surface.
#' @param threshold Numeric; the value of \code{FUN} above which
#' a covariate is included in the target subset.
#' @param method Either "asymptotic" (default) or "quantile"; see details.
#' @param step.down Logical (default \code{TRUE}); should the step-down
#' procedure be used?
#' @param sides One of "both" (default), "exclusive", or "inclusive".
#' Which bounds should be constructed?
#' @param est.FUN The function used to produce estimates of the regression
#' surface. Default is \code{mean}.
#' @param var.FUN The function used to quantify the variability of the
#' regression surface posterior. Default is \code{sd}.
#' @param point.estimate If not null, replaces the mean and sets the reference
#' around which the standard error is computed.
#' Useful for bootstrapping methods.
#' Treated as a row of the \code{params} matrix.
#' @param track A numeric vector of indices indicating which rows (default none)
#' of the design matrix should have the sample of the corresponding
#' \code{FUN(x, params)} returned.
#' @param verbose Logical (default \code{FALSE}); print progress?
#' @param z.store How should certain intermediate computations be handled?
#' See details.
#'
#' @return An object of class \code{credsubs.level}, which contains:
#' \describe{
#' \item{\code{level}}{A numeric vector indicating the maximum credible
#' level at which a conclusion may be drawn at each
#' covariate point.}
#' \item{\code{sign}}{A numeric vector indicating the which credible subsets
#' of which each covariate point is a member at the
#' credible level indicated by \code{level}. Exclusive
#' and inclusive: 1, inclusive only: 0, neither: -1.}
#' \item{\code{threshold}}{As provided.}
#' \item{\code{method}}{As provided.}
#' \item{\code{step.down}}{As provided.}
#' \item{\code{sides}}{As provided.}
#' \item{\code{est}}{Posterior estimate of the regression surface.}
#' \item{\code{est.FUN}}{As provided.}
#' \item{\code{var}}{Summary of posterior variability of the regression
#' surface.}
#' \item{\code{var.FUN}}{As provided.}
#' \item{\code{trace}}{The posterior samples of the regression surface
#' indicated by the \code{track} argument.}
#' }
#'
#' @examples
#' ### Sample from regression surface posterior
#' reg.surf.sample <- matrix(rnorm(1000, mean=1:10), ncol=2, byrow=TRUE)
#' credsubs.level(reg.surf.sample)
#'
#' ### Parametric case
#' design <- cbind(1, 1:10)
#' params <- matrix(rnorm(200, mean=1:2), ncol=2, byrow=TRUE)
#' credsubs(params, design)
#'
#' ### With custom function
#' params.sd <- cbind(1 / rgamma(100, 1), params)
#' FUN.sd <- function(x, params) { params[, -1] %*% t(x) / params[, 1] }
#' credsubs(params.sd, design, FUN.sd, threshold=1)
#'
#' @export
credsubs.level <- function(params, design=NULL,
FUN=function(x, params) { params %*% t(x) },
threshold=0,
method=c('asymptotic', 'quantile'),
step.down=TRUE,
sides=c('both', 'exclusive', 'inclusive'),
est.FUN=mean,
var.FUN=sd,
point.estimate=NULL,
track=numeric(0),
verbose=FALSE,
z.store=c("ram", "recompute", "disk")) {
# Validate and shape input
params <- data.matrix(params)
M <- nrow(params)
if (is.null(design)) {
N <- ncol(params)
nonpar <- TRUE
} else {
design <- data.matrix(design)
N <- nrow(design)
nonpar <- FALSE
}
if (verbose) {
cat("Finding maximum credible level at", N, "points using",
M, "posterior draws.\n")
}
method <- method[1]
stopifnot(method %in% c('asymptotic', 'quantile'))
if (method == 'asymptotic') {
m <- numeric(N)
s <- numeric(N)
} else if (method == 'quantile') {
Fxt <- Gxt <- numeric(N)
}
sides <- sides[1]
stopifnot(sides %in% c('both', 'exclusive', 'inclusive'))
credsubs.level <- list(level=rep(NA, N),
sign=rep(NA, N),
threshold=threshold,
method=method,
step.down=step.down,
sides=sides)
class(credsubs.level) <- 'credsubs.level'
if (z.store[1] == "ram") {
z.store <- matrix(0, nrow=M, ncol=N)
recompute.z <- FALSE
} else if (z.store[1] == "disk") {
if (requireNamespace("ff", quietly=TRUE)) {
z.store <- ff::ff(0, dim=c(M, N))
recompute.z <- FALSE
} else {
warning("Package `ff`` required to use function 'credsubs.level'",
"with option z.store='disk' set. Proceeding using",
"z.store='recompute' instead.")
recompute.z <- TRUE
}
} else if (z.store[1] == "recompute") {
recompute.z <- TRUE
} else {
warning("option z.store must be one of 'ram', 'recompute', or 'disk'.",
"Given: ", z.store[1])
return(sim.cred.band)
}
W <- rep(-Inf, M)
max.i <- numeric(M)
m <- s <- est <- numeric(N)
q <- sgn <- rep(NA, N)
test.set <- 1:N
for (i in 1:N) {
if (verbose & i %% 100 == 0) {
cat("prep", i, "/", N, "\n")
}
if (nonpar) {
fx <- params[, i]
} else {
fx <- FUN(design[i, , drop=FALSE], params)
}
if (method == 'asymptotic') {
if (is.null(point.estimate)) {
m[i] <- mean(fx)
s[i] <- sd(fx)
} else {
if (nonpar) {
point.fx <- point.estimate[i]
} else {
point.fx <- FUN(design[i, , drop=FALSE], point.estimate)
}
m[i] <- point.fx
s[i] <- sqrt(mean((fx - point.fx) ^ 2))
}
est[i] <- m[i]
if (sides == "both") {
z <- abs(fx - m[i]) / s[i]
} else if (sides == "inclusive") {
z <- (fx - m[i]) / s[i]
} else if (sides == "exclusive") {
z <- (m[i] - fx) / s[i]
}
} else {
est[i] <- median(fx)
Fxt[i] <- mean(fx <= threshold)
Gxt[i] <- mean(fx < threshold)
if (sides == "both") {
Fx <- to.Fx(fx)
Gx <- 1 - to.Fx(-fx)
z <- pmax(1-Fx, Gx)
} else if (sides == "inclusive") {
z <- 1 - to.Fx(-fx)
} else if (sides == "exclusive") {
z <- 1 - to.Fx(fx)
}
}
if (!recompute.z) {
z.store[, i] <- z
}
max.i <- ifelse(z > W, i, max.i)
W <- pmax(W, z)
}
sgn <- ifelse(est > threshold &
(sides == 'both' | sides == 'exclusive'),
1,
ifelse(est < threshold &
(sides == 'both' | sides == 'inclusive'),
-1,
0))
q[sgn == 0] <- 1
test.set <- which(sgn != 0)
first.max.i <- max.i
prev.q <- 0
recompute.m <- numeric(0)
while (length(test.set) > 0) {
if (verbose & (N - length(test.set)) %% 10 == 0) {
cat("compute", N - length(test.set), "/", N,
"recompute", length(recompute.m), "\n")
}
# Update W values
if (length(recompute.m) > 0) {
if (recompute.z) {
if (method == 'asymptotic') {
rcm <- recompute.m
} else {
rcm <- 1:M
}
if (nonpar) {
fx <- params[rcm, test.set, drop=FALSE]
} else {
fx <- FUN(design[test.set, , drop=FALSE],
params[rcm, , drop=FALSE])
}
if (method == 'asymptotic') {
if (sides == "both") {
Z <- t(abs(t(fx) - m[test.set]) / s[test.set])
} else if (sides == "inclusive") {
Z <- t((t(fx) - m[test.set]) / s[test.set])
} else if (sides == "exclusive") {
Z <- t((m[test.set] - t(fx)) / s[test.set])
}
} else { # method == quantile
if (sides == "both") {
Fx <- t(apply(fx[recompute.m, , drop=FALSE], 1,
function(x, m) {
colMeans(t(t(m) <= x))
}, m=fx))
Gx <- 1 - t(apply(-fx[recompute.m, , drop=FALSE], 1,
function(x, m) {
colMeans(t(t(m) <= x))
}, m=-fx))
Z <- pmax(1-Fx, Gx)
} else if (sides == "inclusive") {
Z <- 1 - t(apply(-fx[recompute.m, , drop=FALSE], 1,
function(x, m) {
colMeans(t(t(m) <= x))
}, m=-fx))
} else if (sides == "exclusive") {
Z <- 1 - t(apply(fx[recompute.m, , drop=FALSE], 1,
function(x, m) {
colMeans(t(t(m) <= x))
}, m=fx))
}
}
} else {
Z <- z.store[recompute.m, test.set, drop=FALSE]
}
max.i[recompute.m] <- apply(Z, 1, function(z) { test.set[which.max(z)] })
W[recompute.m] <- apply(Z, 1, max)
}
# Update empirical distribution of W
FW <- ecdf(W)
# Compute adjusted p-values
if (method == 'asymptotic') {
if (sides == 'both' | sides == 'exclusive') {
lower <- 1 - FW((m[test.set] - threshold) / s[test.set])
} else {
lower <- rep(1, length(test.set))
}
if (sides == 'both' | sides == 'inclusive') {
upper <- 1 - FW((threshold - m[test.set]) / s[test.set])
} else {
upper <- rep(1, length(test.set))
}
} else {
if (sides == 'both' | sides == 'exclusive') {
lower <- 1 - FW(1 - Fxt[test.set])
} else {
lower <- rep(1, length(test.set))
}
if (sides == 'both' | sides == 'inclusive') {
upper <- 1 - FW(Gxt[test.set])
} else {
upper <- rep(1, length(test.set))
}
}
p <- pmin(lower, upper)
if (step.down) {
lowest.p.i <- test.set[which.min(p)]
q[lowest.p.i] <- prev.q <- max(prev.q, p[which(test.set == lowest.p.i)])
} else {
q <- p
break
}
# these are the MCMC draws of W that change
# when the selected covariate point is removed
recompute.m <- which(max.i == lowest.p.i)
test.set <- test.set[which(test.set != lowest.p.i)]
i <- N - length(test.set)
if (verbose && i %% 100 == 0) {
cat("step", i, "/", N, "(", length(recompute.m), "/", M, ")\n")
}
}
credsubs.level$level <- 1 - q
credsubs.level$sign <- sgn
credsubs.level$threshold <- threshold
credsubs.level$step.down <- step.down
credsubs.level
}
#' Build a credible subset calculator
#'
#' This function builds a \code{shiny} application in the specified directory
#' that gives the maximum credible level at an entered covariate point.
#'
#' The calculator creates a subdirectory according to \code{name} in
#' the directory specified by \code{dir}, and places in it files
#' \code{server.R}, \code{ui.R}, and \code{config.RData}. This application
#' requires the \code{shiny} package to run, and can be executed by passing
#' the directory path to \code{run.shiny.calc()}. The produced application
#' directory may be moved from its original location.
#'
#' @param credsubs.level An object of class \code{credsubs.level}.
#' @param cov.space A data frame whose rows are human-readable
#' covariate points corresponding to the entries
#' of \code{credsubs.level$level}.
#' @param name A character string indicating the name of the application.
#' @param dir The directory in which to place the application.
#' @param title A character string to be displayed as the application title.
#' @param instructions A character string to be displayed as instructions.
#' HTML allowed.
#'
#' @export
build.shiny.calc <- function(credsubs.level,
cov.space,
name="calc",
dir=".",
title="Credible Subsets Calculator",
instructions="Select a covariate point.") {
# Shape and validate input
cov.space <- as.data.frame(cov.space)
stopifnot(length(credsubs.level$level) == nrow(cov.space))
full.name <- paste0(dir, "/", name)
if (!dir.exists(full.name)) {
dir.create(paste0(full.name))
}
shiny.dir <- system.file("shiny", "calc", package="credsubs")
file.copy(paste0(shiny.dir, "/server.R"),
paste0(full.name, "/server.R"),
overwrite=TRUE)
file.copy(paste0(shiny.dir, "/ui.R"),
paste0(full.name, "/ui.R"),
overwrite=TRUE)
# Save file
cat(paste0("Building calculator in ",
full.name, "/", "\n"))
cat(paste0("Execute run.shiny.calc('", full.name , "') to run.", "\n"))
save(credsubs.level, cov.space, title, instructions,
file=paste0(full.name, "/config.RData"))
}
#' Run a calculator
#'
#' Runs the specified calculator using the \code{shiny} package.
#' Calculators must be built using \code{build.shiny.calc}.
#'
#' The \code{app.dir} argument need not exactly match the value recommended
#' by \code{build.shiny.calc()}, as long as it points to the correct directory.
#' For example, \code{"./calc/"}, \code{"calc/"}, and \code{"calc"} are
#' all equivalent. If no value is supplied, an example is run.
#'
#' @param app.dir A character string pointing to the application directory.
#'
#' @export
run.shiny.calc <- function(app.dir=system.file("shiny",
"alzheimers",
package = "credsubs")) {
if (requireNamespace("shiny", quietly=TRUE)) {
shiny::runApp(app.dir, display.mode = "normal")
} else {
warning("Package `shiny` needed for this function to work.", " ",
"Please install it.")
}
}
#' @export
plot.sim.cred.band <- function(x, ...) {
sim.cred.band <- x
o <- order(sim.cred.band$est)
plot(sim.cred.band$est[o],
ylim=c(min(sim.cred.band$lower), max(sim.cred.band$upper)),
xlab="Rank", ylab="Value",
...)
points(sim.cred.band$lower[o], ...)
points(sim.cred.band$upper[o], ...)
}
#' @export
plot.credsubs <- function(x, ...) {
credsubs <- x
# Validate method, estimates, and variabilities
stopifnot(credsubs$method[1] == "asymptotic")
stopifnot(identical(credsubs$est.FUN, mean))
stopifnot(identical(credsubs$var.FUN, sd))
plot(credsubs$est - credsubs$threshold, credsubs$var,
col=ifelse(
credsubs$est -
credsubs$W.crit * credsubs$var > credsubs$threshold,
"forestgreen", ifelse(
credsubs$est +
credsubs$W.crit * credsubs$var >= credsubs$threshold,
"gold", "firebrick"
)),
xlab="Mean - threshold", ylab="Standard deviation", ...)
abline(v=0, lty=2)
abline(h=0, lty=2)
abline(b=1/credsubs$W.crit, a=0)
abline(b=-1/credsubs$W.crit, a=0)
}
#' @export
plot.credsubs.level <- function(x, ...) {
credsubs.level <- x
signed.level <- credsubs.level$level * credsubs.level$sign
plot(sort(signed.level), xlab="Rank", ylab="Signed credible level", ...)
abline(h=0, lty=2)
}
#' @export
print.sim.cred.band <- function(x, ...) {
sim.cred.band <- x
cat("Simultaneous credible band\n")
cat("\n")
cat("Credible level:", sim.cred.band$cred.level, "\n")
cat("Method:", sim.cred.band$method, "\n")
cat("Critical W:", sim.cred.band$W.crit, "\n")
cat("\n")
cat("Upper bounds:", head(sim.cred.band$upper),
ifelse(length(sim.cred.band$upper) > 6, "...\n", "\n"))
cat("Lower bounds:", head(sim.cred.band$lower),
ifelse(length(sim.cred.band$lower) > 6, "...\n", "\n"))
}
#' @export
print.credsubs <- function(x, ...) {
credsubs <- x
cat("Credible subsets\n")
cat("\n")
cat("Credible level:", credsubs$cred.level, "\n")
cat("Threshold:", credsubs$threshold, "\n")
cat("Method:", credsubs$method, "\n")
cat("Step-down:", credsubs$step.down, "\n")
cat("Critical W:", credsubs$W.crit, "\n")
cat("\n")
cat("Exclusive credible subset:", head(credsubs$exclusive),
ifelse(length(credsubs$exclusive) > 6, "...\n", "\n"))
cat("Inclusive credible subset:", head(credsubs$inclusive),
ifelse(length(credsubs$inclusive) > 6, "...\n", "\n"))
}
#' @export
print.credsubs.level <- function(x, ...) {
credsubs.level <- x
cat("Maximum credible levels\n")
cat("\n")
cat("Threshold:", credsubs.level$threshold, "\n")
cat("Method:", credsubs.level$method, "\n")
cat("Step-down:", credsubs.level$step.down, "\n")
cat("\n")
cat("Maximum levels:", head(credsubs.level$level),
ifelse(length(credsubs.level$level) > 6, "...\n", "\n"))
cat("Signs:", head(credsubs.level$sign),
ifelse(length(credsubs.level$sign) > 6, "...\n", "\n"))
}
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.