Nothing
#' Predictions from a Fitted Model
#'
#' Generic prediction method for various types of fitted models. \code{predFit}
#' can be used to obtain standard errors of fitted values and
#' adjusted/unadjusted confidence/prediction intervals for objects of class
#' \code{"lm"}, \code{"nls"}, and \code{"glm"}.
#'
#' @param object An object that inherits from class \code{"lm"}, \code{"glm"},
#' \code{"nls"}, or \code{"lme"}.
#' @param newdata An optional data frame in which to look for variables with
#' which to predict. If omitted, the fitted values are used.
#' @param type Character string specifying the type of prediction. Current
#' options are \code{type = "link"} (the default) and
#' \code{type = "response"}.
#' @param se.fit A logical vaue indicating if standard errors are required.
#' Default is \code{FALSE}.
#' @param interval Type of interval to be calculated. Can be one of "none"
#' (default), "confidence", or "prediction". Default is \code{"none"}.
#' @param level A numeric scalar between 0 and 1 giving the confidence level for
#' the intervals (if any) to be calculated. Default is \code{0.95}.
#' @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 k, times.
#' Default is \code{FALSE}.
#' @param k The number times the calibration curve is to be used for computing
#' a confidence/prediction interval. Only needed when
#' \code{adjust = "Bonferroni"}.
#' @param ... Additional optional arguments. At present, no optional arguments
#' are used.
#'
#' @returns If \code{se.fit = FALSE}, then \code{predFit()} returns a vector of
#' predictions or a matrix of predictions and bounds with column names
#' \code{fit}, \code{lwr}, and \code{upr} if \code{interval} is not
#' \code{"none"}. (This function is more so meant for internal use.)
#'
#' If \code{se.fit = TRUE}, then a list with the following components is
#' returned:
#'
#' \itemize{
#'
#' \item \code{fit} a vector or matrix as described above;
#'
#' \item \code{se.fit} a vector containing the standard errors of the
#' predicted means;
#'
#' \item \code{residual.scale} the residual standard deviations;
#'
#' \item \code{df} the residual degrees of freedom.
#'
#' }
#'
#' @details
#' Confidence and prediction intervals for linear models (i.e., \code{"lm"}
#' objects) are obtained according to the usual formulas. Nonlinear and
#' generalized linear models (i.e., \code{"nls"} and \code{"glm"} objects), on
#' the other hand, rely on Taylor-series approximations for the standard errors
#' used in forming the intervals. Approximate standard errors for the fitted
#' values in linear mixed-effects models (i.e., \code{"lme"} objects) can also
#' be computed; however, these rely on the approximate variance-covariance
#' matrix of the fixed-effects estimates and often under estimate the true
#' standard error. More accurate standard errors can be obtained using the
#' parametric bootstrap; see \code{\link[lme4]{bootMer}} for details.
#'
#' For linear and nonlinear models, it is possible to request \emph{adjusted}
#' confidence or prediction intervals using the Bonferroni method
#' (\code{adjust = "Bonferroni"}) or Scheffe's method
#' (\code{adjust = "Scheffe"}). For the Bonferroni adjustment, you must specify
#' a value for \code{k}, the number of intervals for which the coverage is to
#' hold simultaneously. For the Scheffe adjustment, specifying a value for
#' \code{k} is only required when \code{interval = "prediction"}; if
#' \code{interval = "confidence"}, \code{k} is set equal to \eqn{p}, the number
#' of regression parameters. For example, calling \code{plotFit} on \code{"lm"}
#' objects with \code{interval = "confidence"} and \code{adjust = "Scheffe"}
#' will plot the Working-Hotelling band.
#'
#' @export
#' @examples
#' # A linear regression example (see ?datasets::cars)
#' cars.lm <- lm(dist ~ speed + I(speed^2), data = cars)
#' predFit(cars.lm, interval = "confidence")
#'
#' # A nonlinear least squares example (see ?datasets::Puromycin)
#' data(Puromycin, package = "datasets")
#' Puromycin2 <- Puromycin[Puromycin$state == "treated", ][, 1:2]
#' Puro.nls <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin2,
#' start = c(Vm = 200, K = 0.05))
#' conc <- seq(from = 0.02, to = 1.10, length = 101)
#' pred <- predFit(Puro.nls, newdata = data.frame(conc), interval = "prediction")
#' plot(Puromycin2, ylim = c(min(pred[, "lwr"]), max(pred[, "upr"])))
#' lines(conc, pred[, "fit"], lwd = 2)
#' lines(conc, pred[, "lwr"], lty = 2)
#' lines(conc, pred[, "upr"], lty = 2)
predFit <- function(object, ...) {
UseMethod("predFit")
}
#' @rdname predFit
#' @export
predFit.default <- function(object, ...) {
stats::predict(object, ...)
}
#' @rdname predFit
#' @export
predFit.lm <- function(object, newdata, se.fit = FALSE,
interval = c("none", "confidence", "prediction"),
level = 0.95,
adjust = c("none", "Bonferroni", "Scheffe"), k,
...) {
# Make sure se.fit is set to TRUE if intervals are requested
interval <- match.arg(interval)
compute.se.fit <- if (se.fit || (interval != "none")) TRUE else FALSE
# Predicted values and, if requested, standard errors too
if (missing(newdata)) {
pred <- stats::predict(object, se.fit = compute.se.fit)
} else {
# Suppress warning message printed when calling predict.lm with new data
suppressWarnings(
pred <- stats::predict(object, newdata = as.data.frame(newdata),
se.fit = compute.se.fit)
)
}
# Compute results
if (interval == "none") {
# Vector of fitted/predicted values
res <- pred
} else {
# Critical value for interval computations
adjust <- match.arg(adjust)
crit <- if (adjust == "Bonferroni") { # Bonferroni adjustment
stats::qt((level + 2*k - 1) / (2*k), pred$df)
} else if (adjust == "Scheffe") { # Scheffe adjustment
# Working-Hotelling band or adjusted prediction band for k predictions
if (interval == "confidence") {
p <- length(stats::coef(object))
sqrt(p * stats::qf(level, p, pred$df)) # Working-Hotelling band
} else {
sqrt(k * stats::qf(level, k, pred$df)) # need k for prediction
}
} else { # no adjustment
stats::qt((level + 1) / 2, pred$df)
}
# Interval calculations
if (interval == "confidence") { # confidence interval for mean response
lwr <- pred$fit - crit * pred$se.fit
upr <- pred$fit + crit * pred$se.fit
} else { # prediction interval for individual response
lwr <- pred$fit - crit * sqrt(stats::sigma(object)^2 + pred$se.fit^2)
upr <- pred$fit + crit * sqrt(stats::sigma(object)^2 + pred$se.fit^2)
# warning("predictions on current data refer to _future_ responses")
}
# Store results in a matrix
res <- cbind("fit" = pred$fit, "lwr" = lwr, "upr" = upr)
}
# If standard errors of fitted values are requested, convert results to a list
# and store additional information
if (se.fit) {
res <- list("fit" = if (interval != "none") res else pred$fit,
"se.fit" = pred$se.fit,
"df" = pred$df,
"residual.scale" = pred$residual.scale)
}
# Return results
return(res)
}
#' @rdname predFit
#' @export
predFit.glm <- function(object, newdata, type = c("link", "response"),
se.fit = FALSE, interval = c("none", "confidence"),
level = 0.95, ...) {
# Match arguments
type <- match.arg(type)
interval <- match.arg(interval)
# Make sure se.fit is set to TRUE if intervals are requested
if (se.fit || (interval != "none")) {
compute.se.fit <- TRUE
} else {
compute.se.fit <- FALSE
}
# Predicted values and, if requested, standard errors too
if (missing(newdata)) {
pred <- stats::predict(object, se.fit = compute.se.fit)
} else {
# Suppress warning message printed when calling predict.lm with new data
suppressWarnings(
pred <- stats::predict(object, newdata = as.data.frame(newdata),
type = type, se.fit = compute.se.fit)
)
}
# Compute results
if (interval == "none") {
res <- pred
} else {
res <- cbind("fit" = pred$fit,
"lwr" = pred$fit - pred$se.fit * stats::qnorm((level + 1) / 2),
"upr" = pred$fit + pred$se.fit * stats::qnorm((level + 1) / 2))
if (type == "response") {
res <- apply(res, MARGIN = 2, FUN = function(x) {
stats::family(object)$linkinv(x)
})
}
}
# If standard errors of fitted values are requested, convert results to a list
# and store addional information
if (se.fit) {
res <- list("fit" = if (interval != "none") res else pred$fit,
"se.fit" = pred$se.fit,
"df" = pred$df,
"residual.scale" = pred$residual.scale)
}
# Return results
return(res)
}
#' @rdname predFit
#' @export
predFit.nls <- function(object, newdata, se.fit = FALSE,
interval = c("none", "confidence", "prediction"),
level = 0.95,
adjust = c("none", "Bonferroni", "Scheffe"), k,
...) {
# Match arguments
interval <- match.arg(interval)
adjust <- match.arg(adjust)
# Make sure se.fit is set to TRUE if intervals are requested
compute.se.fit <- if (se.fit || (interval != "none")) {
TRUE
} else {
FALSE
}
# No support for the Golub-Pereyra algorithm for partially linear
# least-squares models
if (interval != "none") {
if (!is.null(object$call$algorithm) && object$call$algorithm == "plinear") {
stop(paste0("The Golub-Pereyra algorithm for partially linear least-",
"squares models is currently not supported."), call. = FALSE)
}
}
# Prediction data
newdata <- if (missing(newdata)) {
eval(stats::getCall(object)$data, envir = parent.frame())
} else {
as.data.frame(newdata)
}
if (is.null(newdata)) {
stop("No data available for predictions.", call. = FALSE)
}
# Name of independent variable
xname <- intersect(all.vars(stats::formula(object)[[3]]), colnames(newdata))
# Predicted values
pred <- object$m$predict(newdata)
# Compute standard error
if (compute.se.fit) {
# Assign values to parameter names in current environment
param.names <- names(stats::coef(object))
for (i in 1:length(param.names)) {
assign(param.names[i], stats::coef(object)[i])
}
# Assign values to independent variable name
assign(xname, newdata[, xname])
# Calculate gradient (numerically)
form <- object$m$formula()
rhs <- eval(form[[3]])
if (is.null(attr(rhs, "gradient"))) {
f0 <- attr(stats::numericDeriv(form[[3]], param.names), "gradient")
} else { # self start models should have gradient attribute
f0 <- attr(rhs, "gradient")
}
# Calculate standard error
R1 <- object$m$Rmat()
# v0 <- diag(f0 %*% solve(t(R1) %*% R1) %*% t(f0))
v0 <- diag(f0 %*% tcrossprod(solve(crossprod(R1)), f0)) # slightly faster
se_fit <- sqrt(stats::sigma(object)^2 * v0)
}
# Compute results
if (interval == "none") {
# Vector of fitted/predicted values
res <- pred
} else {
# Adjustment for simultaneous inference
crit <- if (adjust == "Bonferroni") { # Bonferroni adjustment
stats::qt((level + 2*k - 1) / (2*k), stats::df.residual(object))
} else if (adjust == "Scheffe") { # Scheffe adjustment
if (interval == "confidence") {
p <- length(stats::coef(object)) # number of regression parameters
# sqrt(p * stats::qf((level + 1) / 2, p, stats::df.residual(object)))
sqrt(p * stats::qf(level, p, stats::df.residual(object)))
} else {
# sqrt(k * stats::qf((level + 1) / 2, k, stats::df.residual(object)))
sqrt(k * stats::qf(level, k, stats::df.residual(object)))
}
} else { # no adjustment
stats::qt((level + 1) / 2, stats::df.residual(object))
}
# Interval calculations
if (interval == "confidence") { # confidence limits for mean response
lwr <- pred - crit * se_fit # lower limits
upr <- pred + crit * se_fit # upper limits
} else { # prediction limits for individual response
lwr <- pred - crit * sqrt(stats::sigma(object)^2 + se_fit^2) # lower limits
upr <- pred + crit * sqrt(stats::sigma(object)^2 + se_fit^2) # upper limits
}
# Store results in a matrix
res <- cbind("fit" = pred, "lwr" = lwr, "upr" = upr)
}
# If standard errors of fitted values are requested, convert results to a list
# and store additional information
if (se.fit) {
res <- list("fit" = if (interval != "none") res else pred, #res,
"se.fit" = se_fit,
"df" = stats::df.residual(object),
"residual.scale" = stats::sigma(object))
}
# Return results
return(res)
}
#' @rdname predFit
#' @export
predFit.lme <- function(object, newdata, se.fit = FALSE, ...) {
# Prediction data
newdata <- if (missing(newdata)) {
object$data
} else {
as.data.frame(newdata)
}
# Names of independent variables
xname <- intersect(all.vars(stats::formula(object)[[3L]]), colnames(newdata))
# Population predicted values
pred <- stats::predict(object, newdata = newdata, level = 0)
# Approximate standard errors of fitted values
if (se.fit) {
Xmat <- makeX(object, newdata) # fixed-effects design matrix
se <- sqrt(diag(Xmat %*% stats::vcov(object) %*% t(Xmat)))
cbind("fit" = pred, "se.fit" = se)
} else {
pred
}
}
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.