diff --git a/NAMESPACE b/NAMESPACE
index dae91b42b..cb61e2870 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -56,6 +56,8 @@ S3method(print_md,visualisation_matrix)
S3method(reshape_grouplevel,data.frame)
S3method(reshape_grouplevel,default)
S3method(reshape_grouplevel,estimate_grouplevel)
+S3method(residualize_over_grid,data.frame)
+S3method(residualize_over_grid,estimate_means)
S3method(smoothing,data.frame)
S3method(smoothing,numeric)
S3method(standardize,estimate_contrasts)
@@ -95,6 +97,7 @@ export(pool_slopes)
export(print_html)
export(print_md)
export(reshape_grouplevel)
+export(residualize_over_grid)
export(smoothing)
export(standardize)
export(unstandardize)
diff --git a/R/residualize_over_grid.R b/R/residualize_over_grid.R
new file mode 100644
index 000000000..d42f1a9d2
--- /dev/null
+++ b/R/residualize_over_grid.R
@@ -0,0 +1,171 @@
+#' @title Compute partial residuals from a data grid
+#' @name residualize_over_grid
+#'
+#' @description This function computes partial residuals based on a data grid,
+#' where the data grid is usually a data frame from all combinations of factor
+#' variables or certain values of numeric vectors. This data grid is usually used
+#' as `newdata` argument in `predict()`, and can be created with
+#' [`insight::get_datagrid()`].
+#'
+#' @param grid A data frame representing the data grid, or an object of class
+#' `estimate_means` or `estimate_predicted`, as returned by the different
+#' `estimate_*()` functions.
+#' @param model The model for which to compute partial residuals. The data grid
+#' `grid` should match to predictors in the model.
+#' @param predictor_name The name of the focal predictor, for which partial residuals
+#' are computed.
+#' @param ... Currently not used.
+#'
+#' @section Partial Residuals:
+#' For **generalized linear models** (glms), residualized scores are computed as
+#' `inv.link(link(Y) + r)` where `Y` are the predicted values on the response
+#' scale, and `r` are the *working* residuals.
+#'
+#' For (generalized) linear **mixed models**, the random effect are also
+#' partialled out.
+#'
+#' @references
+#' Fox J, Weisberg S. Visualizing Fit and Lack of Fit in Complex Regression
+#' Models with Predictor Effect Plots and Partial Residuals. Journal of
+#' Statistical Software 2018;87.
+#'
+#' @return A data frame with residuals for the focal predictor.
+#'
+#' @examplesIf requireNamespace("marginaleffects", quietly = TRUE)
+#' set.seed(1234)
+#' x1 <- rnorm(200)
+#' x2 <- rnorm(200)
+#' # quadratic relationship
+#' y <- 2 * x1 + x1^2 + 4 * x2 + rnorm(200)
+#'
+#' d <- data.frame(x1, x2, y)
+#' model <- lm(y ~ x1 + x2, data = d)
+#'
+#' pr <- estimate_means(model, c("x1", "x2"))
+#' head(residualize_over_grid(pr, model))
+#' @export
+residualize_over_grid <- function(grid, model, ...) {
+ UseMethod("residualize_over_grid")
+
+}
+
+
+#' @rdname residualize_over_grid
+#' @export
+residualize_over_grid.data.frame <- function(grid, model, predictor_name, ...) {
+
+ old_d <- insight::get_predictors(model)
+ fun_link <- insight::link_function(model)
+ inv_fun <- insight::link_inverse(model)
+ predicted <- grid[[predictor_name]]
+ grid[[predictor_name]] <- NULL
+
+ is_fixed <- sapply(grid, function(x) length(unique(x))) == 1
+ grid <- grid[, !is_fixed, drop = FALSE]
+ old_d <- old_d[, colnames(grid)[colnames(grid) %in% colnames(old_d)], drop = FALSE]
+
+ if (!.is_grid(grid)) {
+ insight::format_error("Grid for partial residuals must be a fully crossed grid.")
+ }
+
+ # for each var
+ best_match <- NULL
+
+ for (p in colnames(old_d)) {
+ if (is.factor(old_d[[p]]) || is.logical(old_d[[p]]) || is.character(old_d[[p]])) {
+ grid[[p]] <- as.character(grid[[p]])
+ old_d[[p]] <- as.character(old_d[[p]])
+ } else {
+ grid[[p]] <- .validate_num(grid[[p]])
+ }
+
+ # if factor / logical / char in old data, find where it is equal
+ # if numeric in old data, find where it is closest
+ best_match <- .closest(old_d[[p]], grid[[p]], best_match = best_match)
+ }
+
+ idx <- apply(best_match, 2, which)
+ idx <- sapply(idx, "[", 1)
+
+ # extract working residuals
+ res <- .safe(stats::residuals(model, type = "working"))
+
+ # if failed, and model linear, extract response residuals
+ if (is.null(res)) {
+ minfo <- insight::model_info(model)
+ if (minfo$is_linear) {
+ res <- .safe(insight::get_residuals(model, type = "response"))
+ }
+ }
+
+ if (is.null(res)) {
+ insight::format_alert("Could not extract residuals.")
+ return(NULL)
+ }
+
+ my_points <- grid[idx, , drop = FALSE]
+ my_points[[predictor_name]] <- inv_fun(fun_link(predicted[idx]) + res) # add errors
+
+ my_points
+}
+
+
+#' @export
+residualize_over_grid.estimate_means <- function(grid, model, ...) {
+ new_d <- as.data.frame(grid)
+
+ relevant_columns <- unique(c(
+ attributes(grid)$trend,
+ attributes(grid)$contrast,
+ attributes(grid)$focal_terms,
+ attributes(grid)$coef_name
+ ))
+
+ new_d <- new_d[colnames(new_d) %in% relevant_columns]
+
+ residualize_over_grid(new_d, model, predictor_name = attributes(grid)$coef_name, ...)
+}
+
+
+# utilities --------------------------------------------------------------------
+
+
+.is_grid <- function(df) {
+ unq <- lapply(df, unique)
+
+ if (prod(lengths(unq)) != nrow(df)) {
+ return(FALSE)
+ }
+
+ df2 <- do.call(expand.grid, args = unq)
+ df2$..1 <- 1
+
+ res <- merge(df, df2, by = colnames(df), all = TRUE)
+
+ sum(res$..1) == sum(df2$..1)
+}
+
+
+.closest <- function(x, target, best_match) {
+ if (is.numeric(x)) {
+ # AD <- outer(x, target, FUN = function(x, y) abs(x - y))
+ AD <- abs(outer(x, target, FUN = `-`))
+ idx <- apply(AD, 1, function(x) x == min(x))
+ } else {
+ idx <- t(outer(x, target, FUN = `==`))
+ }
+
+ if (is.matrix(best_match)) {
+ idx <- idx & best_match
+ }
+
+ idx
+}
+
+
+.validate_num <- function(x) {
+ if (!is.numeric(x)) {
+ x <- as.numeric(as.character(x))
+ }
+ x
+}
diff --git a/R/visualisation_recipe.R b/R/visualisation_recipe.R
index 98672830e..cd114d66b 100644
--- a/R/visualisation_recipe.R
+++ b/R/visualisation_recipe.R
@@ -147,6 +147,7 @@
#' @export
visualisation_recipe.estimate_predicted <- function(x,
show_data = FALSE,
+ show_residuals = FALSE,
point = NULL,
line = NULL,
pointrange = NULL,
@@ -170,6 +171,7 @@ visualisation_recipe.estimate_predicted <- function(x,
.visualization_recipe(
x,
show_data = show_data,
+ show_residuals = show_residuals,
point = point,
line = line,
pointrange = pointrange,
diff --git a/R/visualisation_recipe_internal.R b/R/visualisation_recipe_internal.R
index 39f6139fe..52f0acf23 100644
--- a/R/visualisation_recipe_internal.R
+++ b/R/visualisation_recipe_internal.R
@@ -246,6 +246,7 @@
#' @keywords internal
.visualization_recipe <- function(x,
show_data = TRUE,
+ show_residuals = FALSE,
point = NULL,
line = NULL,
pointrange = NULL,
@@ -300,6 +301,15 @@
}
+ # add residual data as next lowest layer
+ if (show_residuals) {
+ layers[[paste0("l", l)]] <- .visualization_recipe_residuals(x, aes)
+ # Update with additional args
+ if (!is.null(point)) layers[[paste0("l", l)]] <- utils::modifyList(layers[[paste0("l", l)]], point)
+ l <- l + 1
+ }
+
+
# intercept line for slopes ----------------------------------
if (inherits(x, "estimate_slopes")) {
layers[[paste0("l", l)]] <- insight::compact_list(list(
@@ -524,3 +534,44 @@
out
}
+
+
+# residuals ----------------------------------------------------------------
+
+
+#' @keywords internal
+.visualization_recipe_residuals <- function(x, aes) {
+ model <- attributes(x)$model
+ residual_data <- residualize_over_grid(x, model)
+
+ # Default changes for binomial models
+ shape <- 16
+ stroke <- 0
+ if (insight::model_info(model)$is_binomial) {
+ shape <- "|"
+ stroke <- 1
+ }
+
+ out <- list(
+ geom = "point",
+ data = residual_data,
+ aes = list(
+ y = y,
+ x = aes$x,
+ color = aes$color,
+ alpha = aes$alpha
+ ),
+ height = 0,
+ shape = shape,
+ stroke = stroke
+ )
+
+ # set default alpha, if not mapped by aes
+ if (is.null(aes$alpha)) {
+ out$alpha <- 1 / 3
+ } else {
+ out$alpha <- NULL
+ }
+
+ out
+}
diff --git a/man/residualize_over_grid.Rd b/man/residualize_over_grid.Rd
new file mode 100644
index 000000000..29455dcb1
--- /dev/null
+++ b/man/residualize_over_grid.Rd
@@ -0,0 +1,64 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/residualize_over_grid.R
+\name{residualize_over_grid}
+\alias{residualize_over_grid}
+\alias{residualize_over_grid.data.frame}
+\title{Compute partial residuals from a data grid}
+\usage{
+residualize_over_grid(grid, model, ...)
+
+\method{residualize_over_grid}{data.frame}(grid, model, predictor_name, ...)
+}
+\arguments{
+\item{grid}{A data frame representing the data grid, or an object of class
+\code{estimate_means} or \code{estimate_predicted}, as returned by the different
+\verb{estimate_*()} functions.}
+
+\item{model}{The model for which to compute partial residuals. The data grid
+\code{grid} should match to predictors in the model.}
+
+\item{...}{Currently not used.}
+
+\item{predictor_name}{The name of the focal predictor, for which partial residuals
+are computed.}
+}
+\value{
+A data frame with residuals for the focal predictor.
+}
+\description{
+This function computes partial residuals based on a data grid,
+where the data grid is usually a data frame from all combinations of factor
+variables or certain values of numeric vectors. This data grid is usually used
+as \code{newdata} argument in \code{predict()}, and can be created with
+\code{\link[insight:get_datagrid]{insight::get_datagrid()}}.
+}
+\section{Partial Residuals}{
+
+For \strong{generalized linear models} (glms), residualized scores are computed as
+\code{inv.link(link(Y) + r)} where \code{Y} are the predicted values on the response
+scale, and \code{r} are the \emph{working} residuals.
+
+For (generalized) linear \strong{mixed models}, the random effect are also
+partialled out.
+}
+
+\examples{
+\dontshow{if (requireNamespace("marginaleffects", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
+set.seed(1234)
+x1 <- rnorm(200)
+x2 <- rnorm(200)
+# quadratic relationship
+y <- 2 * x1 + x1^2 + 4 * x2 + rnorm(200)
+
+d <- data.frame(x1, x2, y)
+model <- lm(y ~ x1 + x2, data = d)
+
+pr <- estimate_means(model, c("x1", "x2"))
+head(residualize_over_grid(pr, model))
+\dontshow{\}) # examplesIf}
+}
+\references{
+Fox J, Weisberg S. Visualizing Fit and Lack of Fit in Complex Regression
+Models with Predictor Effect Plots and Partial Residuals. Journal of
+Statistical Software 2018;87.
+}
diff --git a/man/visualisation_recipe.estimate_predicted.Rd b/man/visualisation_recipe.estimate_predicted.Rd
index 3d6a495fd..52ee38ee2 100644
--- a/man/visualisation_recipe.estimate_predicted.Rd
+++ b/man/visualisation_recipe.estimate_predicted.Rd
@@ -15,6 +15,7 @@
\method{visualisation_recipe}{estimate_predicted}(
x,
show_data = FALSE,
+ show_residuals = FALSE,
point = NULL,
line = NULL,
pointrange = NULL,