Merge branch 'issue-#76-partial-residuals': closes #76

NAMESPACE

@@ -70,6 +70,7 @@ S3method(n_smooths,gam)
 S3method(n_smooths,gamm)
 S3method(parametric_terms,default)
 S3method(parametric_terms,gam)
+S3method(partial_residuals,gam)
 S3method(posterior_samples,default)
 S3method(posterior_samples,gam)
 S3method(predicted_samples,default)
@@ -126,6 +127,7 @@ export(load_mgcv)
 export(n_smooths)
 export(observed_fitted_plot)
 export(parametric_terms)
+export(partial_residuals)
 export(posterior_samples)
 export(predicted_samples)
 export(qq_plot)
@@ -139,6 +141,7 @@ export(smooths)
 export(vars_from_label)
 export(which_smooths)
 importFrom(cowplot,plot_grid)
+importFrom(dplyr,arrange)
 importFrom(dplyr,bind_cols)
 importFrom(dplyr,bind_rows)
 importFrom(dplyr,everything)

NEWS.md

@@ -2,6 +2,12 @@
 
 ## New features
 
+* Partial residuals for models can be computed with `partial_residuals()`. The
+  partial residuals are the weighted residuals of the model added to the
+  contribution of each smooth term (as returned by `predict(model, type = "terms")`.
+
+  Wish of #76 (@noamross)
+
 * Users can now control to some extent what colour or fill scales are used when
   plotting smooths in those `draw()` methods that use them. This is most useful
   to change the fill scale when plotting 2D smooths, or to change the discrete

R/residuals.R

@@ -0,0 +1,55 @@
+##' Partial residuals
+##'
+##' @param object an R object, typically a model. Currently only objects of
+##'   class `"gam"` (or that inherit from that class) are supported.
+##' @param ... arguments passed to other methods.
+##'
+##' @export
+`partial_residuals` <- function(object, ...) {
+    UseMethod("partial_residuals")
+}
+
+##' @param select character, logical, or numeric; which smooths to plot. If
+##'   `NULL`, the default, then all model smooths are drawn. Numeric `select`
+##'   indexes the smooths in the order they are specified in the formula and
+##'   stored in `object`. Character `select` matches the labels for smooths
+##'   as shown for example in the output from `summary(object)`. Logical
+##'   `select` operates as per numeric `select` in the order that smooths are
+##'   stored.
+##' @param partial_match logical; should smooths be selected by partial matches
+##'   with `select`? If `TRUE`, `select` can only be a single string to match
+##'   against.
+##'
+##' @rdname partial_residuals
+##'
+##' @export
+##' 
+##' @importFrom tibble as_tibble
+##' @importFrom dplyr bind_cols arrange
+`partial_residuals.gam` <- function(object, select = NULL, partial_match = FALSE,
+                                    ...) {
+    ## get a vector of labels for smooths
+    sms <- smooths(object)
+    ## which were selected; select = NULL -> all selected
+    take <- check_user_select_smooths(sms, select = select,
+                                      partial_match = partial_match)
+    if (!any(take)) {
+        stop("No smooth label matched 'select'. Try with 'partial_match = TRUE'?",
+             call. = FALSE)
+    }
+    sms <- sms[take] # subset to selected smooths
+
+    ## get the contributions for each selected smooth
+    p_terms <- predict(object, type = "terms", terms = sms)
+    attr(p_terms, "constant") <- NULL # remove intercept attribute
+    ## weight residuals...
+    w_resid <- object$residuals * sqrt(object$weights)
+    ## and compute partial residuals
+    p_resids <- p_terms + w_resid
+
+    ## cast as a tibble --- do something with the column names?
+    ##  - they are non-standard: `s(x)` for example
+    p_resids <- tibble::as_tibble(p_resids)
+
+    p_resids
+}

tests/testthat/test-residuals.R

@@ -0,0 +1,42 @@
+## Test partial_residuals() and related residuals functions
+
+## load packages
+library("testthat")
+library("gratia")
+library("mgcv")
+
+context("Test partial_residuals")
+
+N <- 100L
+data <- data_sim("eg1", n = N, seed = 42)
+## fit the model
+m <-  gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = data, method = 'REML')
+
+test_that("partial_residuals returns a tibble", {
+    expect_silent(p_res <- partial_residuals(m))
+    expect_s3_class(p_res, class = c("tbl_df", "tbl", "data.frame"), exact = TRUE)
+    expect_named(p_res, c("s(x0)", "s(x1)", "s(x2)", "s(x3)"))
+    expect_identical(nrow(p_res), N)
+})
+
+test_that("select works with partial_residuals", {
+    expect_silent(p_res <- partial_residuals(m, select = "s(x1)"))
+    expect_s3_class(p_res, class = c("tbl_df", "tbl", "data.frame"), exact = TRUE)
+    expect_named(p_res, "s(x1)")
+    expect_identical(nrow(p_res), N)
+})
+
+test_that("partial_match selecting works with partial_residuals", {
+    expect_silent(p_res <- partial_residuals(m, select = "x1", partial_match = TRUE))
+    expect_s3_class(p_res, class = c("tbl_df", "tbl", "data.frame"), exact = TRUE)
+    expect_named(p_res, "s(x1)")
+    expect_identical(nrow(p_res), N)
+})
+
+test_that("selecting throws an error if no match", {
+    err_msg <- "No smooth label matched 'select'. Try with 'partial_match = TRUE'?"
+    expect_error(partial_residuals(m, select = "foo", partial_match = TRUE),
+                 err_msg)
+    expect_error(partial_residuals(m, select = "foo", partial_match = FALSE),
+                 err_msg)
+})