Nothing
#' Inverse estimation for linear, nonlinear, and generalized linear models
#'
#' Provides point and interval estimates for the unknown predictor value that
#' corresponds to an observed value of the response (or vector thereof) or
#' specified value of the mean response. See the references listed below for
#' more details.
#'
#' @param object An object that inherits from class \code{\link[stats]{lm}},
#' \code{\link[stats]{glm}}, \code{\link[stats]{nls}}, or
#' \code{\link[nlme]{lme}}.
#'
#' @param y0 The value of the observed response(s) or specified value of the
#' mean response. For \code{\link[stats]{glm}} objects, \code{y0} should be on
#' the scale of the response variable (e.g., a number between 0 and 1 for
#' binomial families).
#'
#' @param interval The type of interval required.
#'
#' @param level A numeric scalar between 0 and 1 giving the confidence level
#' for the interval to be calculated.
#'
#' @param mean.response Logical indicating whether confidence intervals should
#' correspond to an individual response (\code{FALSE}) or a mean response
#' (\code{TRUE}). For \code{\link[stats]{glm}} objects, this is always
#' \code{TRUE}.
#'
#' @param x0.name For multiple linear regression, a character string giving the
#' name of the predictor variable of interest.
#'
#' @param newdata For multiple linear regression, a \code{data.frame} giving the
#' values of interest for all other predictor variables (i.e., those other than
#' \code{x0.name}).
#'
#' @param data An optional data frame. This is required if \code{object$data} is
#' \code{NULL}.
#'
#' @param nsim Positive integer specifying the number of bootstrap simulations;
#' the bootstrap B (or R).
#'
#' @param boot.type Character string specifying the type of bootstrap to use
#' when \code{interval = "percentile"}. Options are \code{"parametric"} and
#' \code{"nonparametric"}.
#'
#' @param seed Optional argument to \code{\link{set.seed}}.
#'
#' @param progress Logical indicating whether to display a text-based progress
#' bar during the bootstrap simulation.
#'
#' @param lower The lower endpoint of the interval to be searched.
#'
#' @param upper The upper endpoint of the interval to be searched.
#'
#' @param extendInt Character string specifying if the interval
#' \code{c(lower, upper)} should be extended or directly produce an error when
#' the inverse of the prediction function does not have differing signs at the
#' endpoints. The default, \code{"no"}, keeps the search interval and hence
#' produces an error. Can be abbreviated. See the documentation for the
#' \code{base} R function \code{uniroot} for details.
#'
#' @param q1 Optional lower cutoff to be used in forming confidence intervals.
#' Only used when \code{object} inherits from class \code{\link[nlme]{lme}}.
#' Defaults to \code{stats::qnorm((1+level)/2)}.
#'
#' @param q2 Optional upper cutoff to be used in forming confidence intervals.
#' Only used when \code{object} inherits from class \code{\link[nlme]{lme}}.
#' Defaults to \code{stats::qnorm((1-level)/2)}.
#'
#' @param tol The desired accuracy passed on to \code{\link[stats]{uniroot}}.
#' Recommend a minimum of \code{1e-10}.
#'
#' @param maxiter The maximum number of iterations passed on to \code{uniroot}.
#'
#' @param adjust A logical value indicating if an adjustment should be made to
#' the critical value used in calculating the confidence interval.This is useful
#' for when the calibration curve is to be used multiple, say \code{k}, times.
#'
#' @param k The number times the calibration curve is to be used for computing
#' a confidence interval. Only needed when \code{adjust = "Bonferroni"}.
#'
#' @param ... Additional optional arguments. At present, no optional arguments
#' are used.
#'
#' @return Returns an object of class \code{"invest"} or, if
#' \code{interval = "percentile"}, of class \code{c("invest", "bootCal")}. The
#' generic function \code{{plot}} can be used to plot the output
#' of the bootstrap simulation when \code{interval = "percentile"}.
#'
#' An object of class \code{"invest"} containing the following components:
#' \itemize{
#' \item \code{estimate} The estimate of x0.
#' \item \code{lwr} The lower confidence limit for x0.
#' \item \code{upr} The upper confidence limit for x0.
#' \item \code{se} An estimate of the standard error (Wald and percentile
#' intervals only).
#' \item \code{bias} The bootstrap estimate of bias (percentile interval
#' only).
#' \item \code{bootreps} Vector of bootstrap replicates (percentile
#' interval only).
#' \item \code{nsim} The number of bootstrap replicates (percentile
#' interval only).
#' \item \code{interval} The method used for calculating \code{lower} and
#' \code{upper} (only used by \code{{print}} method).
#' }
#'
#' @references
#' Greenwell, B. M. (2014). \emph{Topics in Statistical Calibration}.
#' Ph.D. thesis, Air Force Institute of Technology.
#' URL \url{https://apps.dtic.mil/sti/pdfs/ADA598921.pdf}
#'
#'
#' Greenwell, B. M., and Schubert Kabban, C. M. (2014). investr: An R Package
#' for Inverse Estimation. \emph{The R Journal}, \bold{6}(1), 90--100.
#' URL http://journal.r-project.org/archive/2014-1/greenwell-kabban.pdf.
#'
#' Graybill, F. A., and Iyer, H. K. (1994).
#' \emph{Regression analysis: Concepts and Applications}. Duxbury Press.
#'
#' Huet, S., Bouvier, A., Poursat, M-A., and Jolivet, E. (2004)
#' \emph{Statistical Tools for Nonlinear Regression: A Practical Guide with
#' S-PLUS and R Examples}. Springer.
#'
#' Norman, D. R., and Smith H. (2014).
#' \emph{Applied Regression Analysis}. John Wiley & Sons.
#'
#' Oman, Samuel D. (1998).
#' Calibration with Random Slopes.
#' \emph{Biometrics} \bold{85}(2): 439--449.
#' doi:10.1093/biomet/85.2.439.
#'
#' Seber, G. A. F., and Wild, C. J. (1989)
#' \emph{Nonlinear regression}. Wiley.
#'
#' @rdname invest
#'
#' @export
#'
#' @examples
#' #
#' # Dobson's beetle data (generalized linear model)
#' #
#'
#' # Complementary log-log model
#' mod <- glm(cbind(y, n-y) ~ ldose, data = beetle,
#' family = binomial(link = "cloglog"))
#' plotFit(mod, pch = 19, cex = 1.2, lwd = 2,
#' xlab = "Log dose of carbon disulphide",
#' interval = "confidence", shade = TRUE,
#' col.conf = "lightskyblue")
#'
#' # Approximate 95% confidence intervals and standard error for LD50
#' invest(mod, y0 = 0.5)
#' invest(mod, y0 = 0.5, interval = "Wald")
#'
#' #
#' # Nasturtium example (nonlinear least-squares with replication)
#' #
#'
#' # Log-logistic model
#' mod <- nls(weight ~ theta1/(1 + exp(theta2 + theta3 * log(conc))),
#' start = list(theta1 = 1000, theta2 = -1, theta3 = 1),
#' data = nasturtium)
#' plotFit(mod, lwd.fit = 2)
#'
#' # Compute approximate 95% calibration intervals
#' invest(mod, y0 = c(309, 296, 419), interval = "inversion")
#' invest(mod, y0 = c(309, 296, 419), interval = "Wald")
#'
#' # Bootstrap calibration intervals. In general, nsim should be as large as
#' # reasonably possible (say, nsim = 9999).
#' boo <- invest(mod, y0 = c(309, 296, 419), interval = "percentile",
#' nsim = 300, seed = 101)
#' boo # print bootstrap summary
#' plot(boo) # plot results
#'
#' #
#' # Bladder volume example (random coefficient model)
#' #
#'
#' # Load required packages
#' library(nlme)
#'
#' # Plot data
#' plot(HD^(3/2) ~ volume, data = bladder, pch = 19,
#' col = adjustcolor("black", alpha.f = 0.5))
#'
#' # Fit a random intercept and slope model
#' bladder <- na.omit(bladder)
#' ris <- lme(HD^(3/2) ~ volume, data = bladder, random = ~volume|subject)
#' invest(ris, y0 = 500)
#' invest(ris, y0 = 500, interval = "Wald")
invest <- function(object, y0, ...) {
UseMethod("invest")
}
# #' @rdname invest
#
# #' @export
# invest.default <- function(object, y0, x0.name, newdata, data, lower, upper,
# extendInt = "no", tol = .Machine$double.eps^0.25,
# maxiter = 1000, ...) {
# }
#' @rdname invest
#'
#' @export
invest.lm <- function(object, y0,
interval = c("inversion", "Wald", "percentile", "none"),
level = 0.95, mean.response = FALSE,
x0.name, newdata, data,
boot.type = c("parametric", "nonparametric"), nsim = 999,
seed = NULL, progress = FALSE, lower, upper,
extendInt = "no", tol = .Machine$double.eps^0.25,
maxiter = 1000, adjust = c("none", "Bonferroni"),
k, ...) {
# Extract data, variable names, etc.
.data <- if (!missing(data)) {
data
} else {
eval(object$call$data, envir = parent.frame())
}
yname <- all.vars(stats::formula(object)[[2]])
# Predictor variable(s)
multi <- FALSE
xnames <- intersect(all.vars(stats::formula(object)[[3]]), colnames(.data))
if (length(xnames) != 1) {
multi <- TRUE
if (missing(x0.name)) {
stop("'x0.name' is missing, please select a valid predictor variable")
}
xnames <- xnames[xnames != x0.name]
if (missing(newdata)) {
# FIXME: Should the user be warned here about the default behavior?
newdata <- as.data.frame(lapply(.data[, xnames], FUN = function(x) {
if (is.numeric(x)) {
stats::median(x, na.rm = TRUE)
} else {
names(which.max(table(x, useNA = "always")))
}
}))
}
if (!is.data.frame(newdata)) {
stop("'newdata' must be a data frame")
}
if (nrow(newdata) != 1) {
stop("'newdata' must have a single row")
}
if (ncol(newdata) != length(xnames) || !all(xnames %in% names(newdata))) {
stop(paste0("'newdata' must contain a column for each predictor variable",
" used by ", deparse(substitute(object)),
" (except ", x0.name, ")"))
}
} else {
x0.name <- xnames
}
# End-points for 'uniroot'
if (missing(lower)) lower <- min(.data[, x0.name]) # lower limit default
if (missing(upper)) upper <- max(.data[, x0.name]) # upper limit default
# Define constants
m <- length(y0) # number of unknowns
if(mean.response && m > 1) {
stop("Only one mean response value allowed.")
}
eta <- mean(y0) # mean unknown
n <- length(stats::resid(object)) # in case of missing values
p <- length(stats::coef(object)) # number of regression coefficients
df1 <- n - p # stage 1 degrees of freedom
df2 <- m - 1 # stage 2 degrees of freedom
var1 <- stats::sigma(object)^2 # stage 1 variance estimate
var2 <- if (m == 1) {
0
} else {
stats::var(y0) # stage 2 variance estimate
}
var.pooled <- (df1 * var1 + df2 * var2) / (df1 + df2) # pooled estimate
rat <- var.pooled / var1 # right variance?
# Calculate point estimate by inverting fitted model
x0.est <- computeInverseEstimate(object, multi = multi, x0.name = x0.name,
newdata = newdata, eta = eta, lower = lower,
upper = upper, extendInt = extendInt,
tol = tol, maxiter = maxiter)
# Match arguments
interval <- match.arg(interval)
boot.type <- match.arg(boot.type)
adjust <- match.arg(adjust)
# Return point estimate only
if (interval == "none") {
return(stats::setNames(x0.est, x0.name))
}
# Critical value for confidence interval computations
crit <- if (adjust == "Bonferroni" && m == 1) { # Bonferroni adjustment
stats::qt((level + 2 * k - 1) / (2 * k), n + m - p - 1)
} else { # no adjustment
stats::qt((level + 1) / 2, n + m - p - 1)
}
# Calculate confidence limits
res <- if (interval == "inversion") {
# Inversion confidence/prediction interval
computeInversionInterval(object, multi = multi, x0.name = x0.name,
var.pooled = var.pooled, m = m, rat = rat,
eta = eta, crit = crit, x0.est = x0.est,
newdata = newdata, mean.response = mean.response,
lower = lower, upper = upper,
extendInt = extendInt, tol = tol,
maxiter = maxiter)
} else if (interval == "Wald") {
# Wald-based interval and standard error
computeWaldInterval(object, multi = multi, x0.name = x0.name,
var.pooled = var.pooled, m = m, p = p, eta = eta, crit = crit,
x0.est = x0.est, mean.response = mean.response,
newdata = newdata,
lower = lower, upper = upper, extendInt = extendInt,
tol = tol, maxiter = maxiter)
} else {
# Calculate bootstrap replicates
x0.boot <- computeBootReps(object, nsim = nsim, seed = seed,
boot.type = boot.type, yname = yname,
mean.response = mean.response, y0 = y0,
x0.name = x0.name, lower = lower, upper = upper,
extendInt = extendInt, tol = tol,
maxiter = maxiter, progress = progress)
# Percentiles of bootstrap replicates
perc <- unname(stats::quantile(x0.boot,
probs = c((1 - level) / 2, (1 + level) / 2)))
# Create list of results
boo <- list("estimate" = x0.est, # original estimate
"lower" = perc[1L], # lower percentile
"upper" = perc[2L], # upper percentile
"se" = stats::sd(x0.boot), # standard error
"bias" = mean(x0.boot) - x0.est, # estimated bias
"bootreps" = x0.boot, # bootstrap replicates
"nsim" = nsim, # number of simulations
"level" = level, # desired confidence level
"interval" = "percentile") # type of interval requested
# Assign number of failed bootstrap replications as an attribute
attr(boo, "bootFail") <- attr(x0.boot, "bootFail")
boo
}
# Assign class label(s) and return result
if (interval == "percentile") {
class(res) <- c("invest", "bootCal")
} else {
class(res) <- "invest"
}
res
}
#' @rdname invest
#'
#' @export
invest.glm <- function(object, y0,
interval = c("inversion", "Wald", "percentile", "none"),
level = 0.95, lower, upper,
x0.name, newdata, data, extendInt = "no",
tol = .Machine$double.eps^0.25, maxiter = 1000, ...) {
# NOTE: Currently, this function only works for the case
# mean.response = TRUE.
# Extract data, variable names, etc.
.data <- if (!missing(data)) {
data
} else {
eval(object$call$data, envir = parent.frame())
}
# Calculations should be done on the "link" scale!
eta <- stats::family(object)$linkfun(y0)
# Predictor variable(s)
multi <- FALSE
xnames <- intersect(all.vars(stats::formula(object)[[3]]), colnames(.data))
if (length(xnames) != 1) {
multi <- TRUE
if (missing(x0.name)) {
stop("'x0.name' is missing, please select a valid predictor variable")
}
xnames <- xnames[xnames != x0.name]
if (missing(newdata)) {
stop("'newdata' must be supplied when multiple predictor variables exist!")
}
if (!is.data.frame(newdata)) {
stop("'newdata' must be a data frame")
}
if (nrow(newdata) != 1) {
stop("'newdata' must have a single row")
}
if (ncol(newdata) != length(xnames) || !all(xnames %in% names(newdata))) {
stop(paste0("'newdata' must contain a column for each predictor variable",
" used by ", deparse(substitute(object)),
" (except ", x0.name, ")"))
}
} else {
x0.name <- xnames
}
# End-points for 'uniroot'
if (missing(lower)) {
lower <- min(.data[, x0.name, drop = TRUE]) # lower limit default
}
if (missing(upper)) {
upper <- max(.data[, x0.name, drop = TRUE]) # upper limit default
}
# Calculate point estimate by inverting fitted model
x0.est <- computeInverseEstimate(object, multi = multi, x0.name = x0.name,
newdata = newdata, eta = eta, lower = lower,
upper = upper, extendInt = extendInt,
tol = tol, maxiter = maxiter)
# Match arguments
interval <- match.arg(interval)
# Return point estimate only
if (interval == "none") {
return(stats::setNames(x0.est, x0.name))
}
# Critical value for confidence interval computations
crit <- stats::qnorm((level + 1) / 2) # quantile from standard normal
# Calculate confidence limits
res <- if (interval == "inversion") {
# Invert approximate confidence interval for the mean response--based on
# exercise 5.31 on pg. 207 of Categorical Data Analysis (2nd ed.) by Alan
# Agresti.
computeInversionInterval(object, multi = multi, x0.name = x0.name,
eta = eta, crit = crit, x0.est = x0.est,
newdata = newdata, lower = lower, upper = upper,
extendInt = extendInt, tol = tol,
maxiter = maxiter)
} else if (interval == "Wald") {
# Wald-based interval and standard error
computeWaldInterval(object, multi = multi, x0.name = x0.name,
eta = eta, crit = crit, x0.est = x0.est,
newdata = newdata, lower = lower, upper = upper,
extendInt = extendInt, tol = tol, maxiter = maxiter)
} else {
# Bootstrapping is currently not available for GLM objects
stop("Bootstrap intervals not available for 'glm' objects.", call. = FALSE)
}
# Assign class label and return results
class(res) <- "invest"
res
}
#' @rdname invest
#'
#' @export
invest.nls <- function(object, y0,
interval = c("inversion", "Wald", "percentile", "none"),
level = 0.95, mean.response = FALSE, data,
boot.type = c("parametric", "nonparametric"), nsim = 1,
seed = NULL, progress = FALSE, lower, upper,
extendInt = "no", tol = .Machine$double.eps^0.25,
maxiter = 1000, adjust = c("none", "Bonferroni"), k,
...) {
# No support for the Golub-Pereyra algorithm for partially linear
# least-squares models
if (object$call$algorithm == "plinear") {
stop(paste("The Golub-Pereyra algorithm for partially linear least-squares
models is currently not supported."))
}
# Extract data, variable names, etc.
.data <- if (!missing(data)) data else eval(object$call$data,
envir = parent.frame())
x0.name <- intersect(all.vars(stats::formula(object)[[3]]), colnames(.data))
yname <- all.vars(stats::formula(object)[[2]])
if (length(x0.name) != 1) stop("Only one independent variable allowed.")
if (missing(lower)) lower <- min(.data[, x0.name]) # lower limit default
if (missing(upper)) upper <- max(.data[, x0.name]) # upper limit default
# Set up for inverse estimation
m <- length(y0) # number of unknowns
if(mean.response && m > 1) stop("Only one mean response value allowed.")
eta <- mean(y0) # mean response
n <- length(stats::resid(object)) # sample size
p <- length(stats::coef(object)) # number of parameters
var.pooled <- stats::sigma(object)^2 # residual variance
# Calculate point estimate by inverting fitted model
x0.est <- computeInverseEstimate(object, x0.name = x0.name, eta = eta,
lower = lower, upper = upper,
extendInt = extendInt, tol = tol,
maxiter = maxiter)
# Match arguments
interval <- match.arg(interval)
boot.type <- match.arg(boot.type)
adjust <- match.arg(adjust)
# Return point estimate only
if (interval == "none") {
return(stats::setNames(x0.est, x0.name))
}
# Critical value for confidence interval computations
crit <- if (adjust == "Bonferroni" && m == 1) { # Bonferroni adjustment
stats::qt((level + 2 * k - 1) / (2 * k), n + m - p - 1)
} else { # no adjustment
stats::qt((level + 1) / 2, n+m-p-1)
}
# Calculate confidence limits
res <- if (interval == "inversion") {
# Inversion confidence/prediction interval
computeInversionInterval(object, x0.name = x0.name, var.pooled = var.pooled,
m = m, eta = eta, crit = crit, x0.est = x0.est,
mean.response = mean.response, lower = lower,
upper = upper, extendInt = extendInt, tol = tol,
maxiter = maxiter)
} else if (interval == "Wald") {
# Wald-based interval and standard error
computeWaldInterval(object, x0.name = x0.name, var.pooled = var.pooled,
m = m, p = p, eta = eta, crit = crit, x0.est = x0.est,
mean.response = mean.response, lower = lower,
upper = upper, extendInt = extendInt, tol = tol,
maxiter = maxiter)
} else {
# Calculate bootstrap replicates
x0.boot <- computeBootReps(object, nsim = nsim, seed = seed,
boot.type = boot.type, yname = yname,
mean.response = mean.response, y0 = y0,
x0.name = x0.name, lower = lower, upper = upper,
extendInt = extendInt, tol = tol,
maxiter = maxiter, progress = progress)
# Percentiles of bootstrap replicates
perc <- unname(stats::quantile(x0.boot,
probs = c((1 - level) / 2, (1 + level) / 2)))
# Create list of results
boo <- list("estimate" = x0.est, # original estimate
"lower" = perc[1L], # lower percentile
"upper" = perc[2L], # upper percentile
"se" = stats::sd(x0.boot), # standard error
"bias" = mean(x0.boot) - x0.est, # estimated bias
"bootreps" = x0.boot, # bootstrap replicates
"nsim" = nsim, # number of simulations
"level" = level, # desired confidence level
"interval" = "percentile") # type of interval requested
# Assign number of failed bootstrap replications as an attribute
attr(boo, "bootFail") <- attr(x0.boot, "bootFail")
boo
}
# Assign class label(s) and return result
if (interval == "percentile") {
class(res) <- c("invest", "bootCal")
} else {
class(res) <- "invest"
}
res
}
#' @rdname invest
#'
#' @export
invest.lme <- function(object, y0,
interval = c("inversion", "Wald", "percentile", "none"),
level = 0.95, mean.response = FALSE, data, lower, upper,
q1, q2, extendInt = "no", tol = .Machine$double.eps^0.25,
maxiter = 1000, ...) {
# Extract data, variable names, etc.
.data <- if (!missing(data)) data else object$data
x0.name <- intersect(all.vars(stats::formula(object)[[3]]), colnames(.data))
yname <- all.vars(stats::formula(object)[[2]])
if (length(x0.name) != 1) stop("Only one independent variable allowed.")
if (missing(lower)) lower <- min(.data[, x0.name]) # lower limit default
if (missing(upper)) upper <- max(.data[, x0.name]) # upper limit default
# Set up for inverse estimation
m <- length(y0)
if(mean.response && m > 1) {
stop("Only one mean response value allowed.")
}
eta <- mean(y0)
if (m != 1) {
stop('Only a single unknown allowed for objects of class "lme".')
}
N <- length(stats::resid(object))
p <- length(nlme::fixef(object))
# res.var <- stats::sigma(object)^2 # residual variance
# Calculate point estimate by inverting fitted model
x0.est <- computeInverseEstimate(object, x0.name = x0.name, eta = eta,
lower = lower, upper = upper,
extendInt = extendInt, tol = tol,
maxiter = maxiter)
# Match arguments
interval <- match.arg(interval)
# Return point estimate only
if (interval == "none") {
return(stats::setNames(x0.est, x0.name))
}
# Critical value. Oman (1998. pg. 445) suggests a t(1-alpha/2, N-1) dist.
if (missing(q1)) {
q1 <- stats::qnorm((1-level) / 2)
}
if (missing(q2)) {
q2 <- stats::qnorm((1+level) / 2)
}
# Estimate variance of new response
if (!mean.response) {
var.y0 <- varY(object, newdata = makeData(x0.est, x0.name))
}
# Calculate confidence limits
res <- if (interval == "inversion") {
# Invert approximate confidence interval for the mean response--based on
# exercise 5.31 on pg. 207 of Categorical Data Analysis (2nd ed.) by Alan
# Agresti.
computeInversionInterval(object, x0.name = x0.name, m = m, eta = eta,
q1 = q1, q2 = q2, x0.est = x0.est,
mean.response = mean.response, var.y0 = var.y0,
lower = lower, upper = upper,
extendInt = extendInt, tol = tol,
maxiter = maxiter)
} else if (interval == "Wald") {
# Wald-based interval and standard error
computeWaldInterval(object, x0.name = x0.name, p = p, eta = eta, q1 = q1,
q2 = q2, x0.est = x0.est, mean.response = mean.response,
var.y0 = var.y0, lower = lower, upper = upper,
extendInt = extendInt, tol = tol, maxiter = maxiter)
} else {
# Bootstrapping is currently not available for GLM objects
stop("Bootstrap intervals not available for 'glm' objects.", call. = FALSE)
}
# Assign class label and return results
class(res) <- "invest"
res
}
#' @keywords internal
#'
#' @export
print.invest <- function(x, digits = getOption("digits"), ...) {
if (x$interval == "inversion") print(round(unlist(x[1:3]), digits))
if (x$interval == "Wald") print(round(unlist(x[1:4]), digits))
if (x$interval == "percentile") print(round(unlist(x[1:5]), digits))
invisible(x)
}
#' Plots method for bootstrap calibration
#'
#' The \code{\link{plot}} method for \code{"bootCal"} objects. In
#' particular, this method takes a \code{"bootCal"} object and produces plots
#' for the bootstrap replicates of the inverse estimate.
#'
#' @param x An object that inherits from class \code{"bootCal"}.
#'
#' @param ... Additional optional arguments. At present, no optional arguments
#' are used.
#'
#' @returns \code{x} is returned invisibly.
#'
#' @rdname plot.bootCal
#'
#' @export
plot.bootCal <- function(x, ...) {
# Extract bootstrap replicates and original estimate
t <- x$bootreps # bootstrap replicates
t0 <- x$estimate # original estimate
# Calculate number of histogram breaks (better than default)
if (!is.null(t0)) {
nclass <- min(max(ceiling(length(t)/25), 10), 100)
rg <- range(t)
if (t0 < rg[1L])
rg[1L] <- t0
else if (t0 > rg[2L])
rg[2L] <- t0
rg <- rg + 0.05 * c(-1, 1) * diff(rg)
lc <- diff(rg)/(nclass - 2)
n1 <- ceiling((t0 - rg[1L])/lc)
n2 <- ceiling((rg[2L] - t0)/lc)
bks <- t0 + (-n1:n2) * lc
}
# Plots
op <- graphics::par(no.readonly = TRUE) # the whole list of settable par's
graphics::par(mfrow = c(1, 2))
graphics::hist(t, breaks = bks, probability = TRUE,
main = "",
xlab = "Bootstrap replicate")
stats::qqnorm(t,
main= "",
xlab = "Standard normal quantile",
ylab = "Bootstrap quantile")
stats::qqline(t)
# graphics::par(mfrow = c(1, 1))
graphics::par(op)
invisible(x)
}
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.