partial_residuals: Augment a model fit with partial residuals for all terms

Description

Construct a data frame containing the model data, partial residuals for all quantitative predictors, and predictor effects, for use in residual diagnostic plots and other analyses. The result is in tidy form (one row per predictor per observation), allowing it to be easily manipulated for plots and simulations.

Usage

partial_residuals(fit, predictors = everything())

Value

Data frame (tibble) containing the model data and residuals in tidy form. There is one row per selected predictor per observation. All predictors are included as columns, plus the following additional columns:

.obs: Row number of this observation in the original model data frame.
.predictor_name: Name of the predictor this row gives the partial residual for.
.predictor_value: Value of the predictor this row gives the partial residual for.
.partial_resid: Partial residual for this predictor for this observation.
.predictor_effect: Predictor effect $\hat \mu(X_{if}, 0)$ for this observation.

Arguments

fit: The model to obtain residuals for. This can be a model fit with lm() or glm(), or any model with a predict() method that accepts a newdata argument.
predictors: Predictors to calculate partial residuals for. Defaults to all predictors, skipping factors. Predictors can be specified using tidyselect syntax; see help("language", package = "tidyselect") and the examples below.

Predictors and regressors

To define partial residuals, we must distinguish between the predictors, the measured variables we are using to fit our model, and the regressors, which are calculated from them. In a simple linear model, the regressors are equal to the predictors. But in a model with polynomials, splines, or other nonlinear terms, the regressors may be functions of the predictors.

For example, in a regression with a single predictor $X$, the regression model $Y = \beta_0 + \beta_1 X + e$ has one regressor, $X$. But if we choose a polynomial of degree 3, the model is $Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \beta_3 X^3$, and the regressors are $\{X, X^2, X^3\}$.

Similarly, if we have predictors $X_1$ and $X_2$ and form a model with main effects and an interaction, the regressors are $\{X_1, X_2, X_1 X_2\}$.

Partial residuals are defined in terms of the predictors, not the regressors, and are intended to allow us to see the shape of the relationship between a particular predictor and the response, and to compare it to how we have chosen to model it with regressors. Partial residuals are not useful for categorical (factor) predictors, and so these are omitted.

Linear models

Consider a linear model where $\mathbb{E}[Y \mid X = x] = \mu(x)$. The mean function $\mu(x)$ is a linear combination of regressors. Let $\hat \mu$ be the fitted model and $\hat \beta_0$ be its intercept.

Choose a predictor $X_f$, the focal predictor, to calculate partial residuals for. Write the mean function as $\mu(X_f, X_o)$, where $X_f$ is the value of the focal predictor, and $X_o$ represents all other predictors.

If $e_i$ is the residual for observation $i$, the partial residual is

$$r_{if} = e_i + (\hat \mu(x_{if}, 0) - \hat \beta_0).$$

Setting $X_o = 0$ means setting all other numeric predictors to 0; factor predictors are set to their first (baseline) level.

Generalized linear models

Consider a generalized linear model where $g(\mathbb{E}[Y \mid X = x]) = \mu(x)$, where $g$ is a link function. Let $\hat \mu$ be the fitted model and $\hat \beta_0$ be its intercept.

Let $e_i$ be the working residual for observation $i$, defined to be

$$e_i = (y_i - g^{-1}(x_i)) g'(x_i).$$

Choose a predictor $X_f$, the focal predictor, to calculate partial residuals for. Write $\mu$ as $\mu(X_f, X_o)$, where $X_f$ is the value of the focal predictor, and $X_o$ represents all other predictors. Hence $\mu(X_f, X_o)$ gives the model's prediction on the link scale.

The partial residual is again

$$r_{if} = e_i + (\hat \mu(x_{if}, 0) - \hat \beta_0).$$

Interpretation

In linear regression, because the residuals $e_i$ should have mean zero in a well-specified model, plotting the partial residuals against $x_f$ should produce a shape matching the modeled relationship $\mu$. If the model is wrong, the partial residuals will appear to deviate from the fitted relationship. Provided the regressors are uncorrelated or approximately linearly related to each other, the plotted trend should approximate the true relationship between $x_f$ and the response.

In generalized linear models, this is approximately true if the link function $g$ is approximately linear over the range of observed $x$ values.

Additionally, the function $\mu(X_f, 0)$ can be used to show the relationship between the focal predictor and the response. In a linear model, the function is linear; with polynomial or spline regressors, it is nonlinear. This function is the predictor effect function, and the estimated predictor effects $\hat \mu(X_{if}, 0)$ are included in this function's output.

Limitations

Factor predictors (as factors, logical, or character vectors) are detected automatically and omitted. However, if a numeric variable is converted to factor in the model formula, such as with y ~ factor(x), the function cannot determine the appropriate type and will raise an error. Create factors as needed in the source data frame before fitting the model to avoid this issue.

References

R. Dennis Cook (1993). "Exploring Partial Residual Plots", Technometrics, 35:4, 351-362. tools:::Rd_expr_doi("10.1080/00401706.1993.10485350")

Cook, R. Dennis, and Croos-Dabrera, R. (1998). "Partial Residual Plots in Generalized Linear Models." Journal of the American Statistical Association 93, no. 442: 730–39. tools:::Rd_expr_doi("10.2307/2670123")

Fox, J., & Weisberg, S. (2018). "Visualizing Fit and Lack of Fit in Complex Regression Models with Predictor Effect Plots and Partial Residuals." Journal of Statistical Software, 87(9). tools:::Rd_expr_doi("10.18637/jss.v087.i09")

Examples

Run this code

fit <- lm(mpg ~ cyl + disp + hp, data = mtcars)
partial_residuals(fit)

# You can select predictors with tidyselect syntax:
partial_residuals(fit, c(disp, hp))

# Predictors with multiple regressors are supported:
fit2 <- lm(mpg ~ poly(disp, 2), data = mtcars)
partial_residuals(fit2)

# Allowing an interaction by number of cylinders is fine, but partial
# residuals are not generated for the factor. Notice the factor must be
# created first, not in the model formula:
mtcars$cylinders <- factor(mtcars$cyl)
fit3 <- lm(mpg ~ cylinders * disp + hp, data = mtcars)
partial_residuals(fit3)

Run the code above in your browser using DataLab