#' coefficient of determination
#'
#' @export
#'
R2 <- function(x, ...) {
UseMethod("R2")
}
#' Default S3 method for coefficient of determination
#' @param x observed data
#' @param xpred predictions
#' @export
R2.default <- function(x, xpred) {
stopifnot(length(x) == length(xpred))
1 - sum((x - xpred)^2) / sum((x - mean(x))^2)
}
#' Coefficient of dermination for lm objects
#' Optionally gives adjusted R2.
#' @param x an object of class "lm"
#' @export
R2.lm <- function(x, adjust = FALSE) {
if(adjust)
out <- summary(x)$adj.r.squared
else
out <- summary(x)$r.squared
out
}
#' Coefficient of dermination for gam objects
#' Optionally gives adjusted R2.
#' @param x an object of class "gam"
#' @export
R2.gam <- function(x, adjust = FALSE) {
if (!requireNamespace("mgcv", quietly = TRUE)) {
stop("mgcv needed for this function to work. Please install it.",
call. = FALSE)
}
resid <- mgcv::residuals.gam(x, type = "response")
if(adjust)
out <- summary(x)$r.sq
else
out <- 1 - sum(resid^2) / sum((x$y - mean(x$y))^2)
out
}
#' Leave-one-out coefficient of determination
#' Similar to R2, but uses leave-one-out estimates for response and mean.
#' Uses GCV and OCV to estimate MSE of prediction.
#' @param object an object with a class for `deviance`, `residuals`, and `hatvalues`.
#' @export
#'
Q2 <- function(object, method = c("gcv", "ocv")) {
method = match.arg(method)
mse <- ifelse(method == "gcv", gcv(object), ocv(object))
1 - length(object$y) * mse / tss_loo(object$y)
}
#' mean absolute error
#' @param obs observed data
#' @param pred predictions
#'
#' @export
mae <- function(obs, pred) {
stopifnot(length(obs) == length(pred))
mean(abs(obs - pred))
}
#' sum of squared errors
#' @param obs observed data
#' @param pred predictions
#'
#' @export
sse <- function(obs, pred) {
stopifnot(length(obs) == length(pred))
sum((obs - pred)^2)
}
#' leave-one-out sum of squares
#' useful for calculating Q2
#' @export
tss_loo <- function(x) {
omitted_mean <- function(ind) mean(x[-ind])
omeans <- vapply(1:length(x), omitted_mean, numeric(1))
sum((x - omeans)^2)
}
# gam functions --------------------------------------------------------------
#' Diagonals of GAM influence matrix
#' @importFrom mgcv magic.post.proc
#' @param X A model matrix, potentially with unobserved prediction values
#' @param object a gam model
hatvalues.gam <- function(object, X = NULL) {
if (is.null(X))
return(object$hat)
return(magic.post.proc(X = X, object = object)$hat)
}
#' Studentized deleted residuals, given a model object
#' @param object a model object that undertstands `residuals` and `hatvalues`
#' @importFrom stats rstudent
#' @export
rstudent.gam <- function(object, ...) {
r <- delResid(object)
s <- sqrt(delMSE(object) / (1 - hatvalues(object)))
out <- r / s
out
}
#' Deleted residuals, given a model object
#' @param object a model object that undertstands `residuals` and `hatvalues`
#' @export
delResid <- function(object) {
residuals(object = object, type = "response") / (1 - hatvalues(object))
}
#' Mean-squared error for a model
#' @param model a model object that understands `residuals`
#' @export
MSE <- function(model) {
errs <- residuals(model, type = "response")
out <- sum(errs^2) / model$df.residual
out
}
#' Deleted MSE
#' Useful for calculating deleted statistics like dffits, studentized deleted residuals
#' @param model object, e.g. lm or gam
#' @export
delMSE <- function(model) {
n_p <- model$df.residual
errs <- residuals(model, type = "response")
numer <- MSE(model) * n_p - errs^2 / (1 - hatvalues(model))
out <- numer / (n_p - 1)
out
}
#' DFFITS for a given model
#' @export
dffits.gam <- function(model) {
numer <- delResid(model) - residuals(model, type = "response")
out <- numer / sqrt(delMSE(model) * hatvalues(model))
out
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.