Nothing
# inspired by `geom_smooth` and `stat_smooth` in ggplot2
#' @title Model Prediction
#' @name l_predict
#' @description It is entirely for the purpose of plotting fits and intervals on a scatterplot (or histogram).
#' It is a generic function to predict models for \code{loon} smooth layer (a wrap of the function \code{predict}).
#' However, the output is unified.
#' @param model a model object for which prediction is desired
#' @param ... some arguments passed into the function \code{predict}
#' @return
#' A data frame is returned with \code{x} (if \code{newdata} is given) and \code{y}.
#' If the \code{interval} is not \code{none},
#' two more columns, \code{lower} (lower interval) and \code{upper} (upper interval) will be returned.
#'
#' @examples
#' y <- rnorm(10)
#' x <- rnorm(10)
#' model1 <- lm(y ~ x)
#' # formal output
#' pre <- l_predict(model1, newdata = data.frame(x = sort(x)),
#' interval = "conf")
#' head(pre)
#'
#' if(interactive()) {
#' p <- with(cars, l_plot(speed, dist))
#'
#' # Example taken from
#' # https://stackoverflow.com/questions/23852505/how-to-get-confidence-interval-for-smooth-spline
#' #
#' l_predict.smooth.spline <- function(model, interval = c("confidence", "none"),
#' level = 0.95, ...) {
#' # confidence interval of `smooth.spline`
#' interval <- match.arg(interval)
#'
#' res <- (model$yin - model$y)/(1 - model$lev) # jackknife residuals
#' sigma <- sqrt(var(res)) # estimate sd
#' std <- stats::qnorm(level / 2 + 0.5)
#' upper <- model$y + std * sigma * sqrt(model$lev) # upper 95% conf. band
#' lower <- model$y - std * sigma * sqrt(model$lev) # lower 95% conf. band
#'
#' data.frame(y = model$yin, lower = lower, upper = upper)
#' }
#' l <- l_layer_smooth(p, method = "smooth.spline", interval = "confidence")
#' }
#'
#' @export
l_predict <- function(model, ...)
UseMethod("l_predict")
#' @export
#' @rdname l_predict
l_predict.default <- function(model, ...) {
stop("Unknown method")
}
#' @export
#' @rdname l_predict
#' @param newdata optionally, a data frame in which to look for variables with which to predict.
#' If omitted, the fitted linear predictors are used.
#' @param interval type of interval, could be "none", "confidence" or "prediction" (not for \code{glm})
#' @param level confidence level
#' @param ... arguments passed in \code{predict}
l_predict.lm <- function(model, newdata = NULL,
interval = c("none", "confidence", "prediction"),
level = 0.95, ...) {
interval <- match.arg(interval)
se.fit <- TRUE
if(interval == "none") se.fit <- FALSE
pred <- stats::predict(model,
newdata,
se.fit = se.fit,
level = level,
interval = interval, ...)
if (se.fit) {
fit <- as.data.frame(pred$fit)
names(fit) <- c("y", "lower", "upper")
if(is.null(newdata)) {
data.frame(fit)
} else {
data.frame(newdata, fit)
}
} else {
if(is.null(newdata)) {
data.frame(y = as.vector(pred))
} else {
data.frame(newdata, y = as.vector(pred))
}
}
}
#' @rdname l_predict
#' @export
l_predict.nls <- function(model, newdata = NULL,
interval = c("none", "confidence", "prediction"),
level = 0.95, ...) {
l_predict.lm(model, newdata,
interval = interval,
level = level, ...)
}
#' @rdname l_predict
#' @export
l_predict.glm <- function(model, newdata = NULL,
interval = c("none", "confidence"),
level = 0.95, ...) {
interval <- match.arg(interval)
se.fit <- TRUE
if(interval == "none") se.fit <- FALSE
pred <- stats::predict(model, newdata = newdata,
se.fit = se.fit, type = "link", ...)
if (se.fit) {
if(interval == "prediction") {
message("Prediction interval is not provided for `glm`")
}
se <- as.vector(pred$se.fit)
std <- stats::qnorm(level / 2 + 0.5)
y <- model$family$linkinv(as.vector(pred$fit))
lower <- model$family$linkinv(as.vector(pred$fit - std * se))
upper <- model$family$linkinv(as.vector(pred$fit + std * se))
if(is.null(newdata)) {
data.frame(
y = y,
lower = lower,
upper = upper
)
} else {
data.frame(
newdata,
y = y,
lower = lower,
upper = upper
)
}
} else {
if(is.null(newdata)) {
data.frame(y = model$family$linkinv(as.vector(pred)))
} else {
data.frame(newdata, y = model$family$linkinv(as.vector(pred)))
}
}
}
#' @rdname l_predict
#' @export
l_predict.loess <- function(model, newdata = NULL,
interval = c("none", "confidence", "prediction"),
level = 0.95, ...) {
interval <- match.arg(interval)
se.fit <- TRUE
if(interval == "none") se.fit <- FALSE
pred <- stats::predict(model, newdata, se = se.fit, ...)
if (se.fit) {
se <- if(interval == "confidence") {
as.vector(pred$se.fit)
} else {
# predict interval
sqrt(as.vector(pred$se.fit)^2 + as.vector(pred$residual.scale)^2)
}
y <- pred$fit
ci <- se * stats::qt(level / 2 + .5, pred$df)
lower <- y - ci
upper <- y + ci
if(is.null(newdata)) {
data.frame(y = y,
lower = lower,
upper = upper)
} else {
data.frame(newdata, y = y,
lower = lower,
upper = upper)
}
} else {
if(is.null(newdata)) {
data.frame(y = as.vector(pred))
} else {
data.frame(newdata, y = as.vector(pred))
}
}
}
# l_predict.locfit <- function(model, newdata, se, level) {
# pred <- stats::predict(model, newdata, se.fit = se)
#
# if (se) {
# y <- pred$fit
# lower <- y - pred$se.fit
# upper <- y + pred$se.fit
#
# if(is.null(newdata)) {
# data.frame(y, lower, upper)
# } else {
# data.frame(newdata, y, lower, upper)
# }
#
# } else {
# if(is.null(newdata)) {data.frame(y = as.vector(pred))}
# else {data.frame(newdata, y = as.vector(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.