--- title: "Functions of covariates and evaluation environments" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Functions of covariates and evaluation environments} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library("gratia") library("mgcv") ``` A smooth can depend on a function of a covariate rather than on the raw covariate. The same distinction arises with offsets, such as logged effort in a count model. ```{r model} set.seed(42) dat <- data.frame(x = runif(150, 1, 5), effort = runif(150, 1, 3)) dat$y <- rpois(nrow(dat), dat$effort * exp(sin(log(dat$x)))) m <- gam(y ~ s(log(x), k = 6) + offset(log(effort)), family = poisson(), data = dat, method = "REML") ``` ## Raw data and smooth coordinates Supply raw covariates to evaluate the model at new observations. Gratia computes `log(x)` and, for full predictions, `offset(log(effort))`. An isolated smooth does not require unrelated predictors or the model offset. ```{r raw-data} nd <- data.frame(x = seq(1.2, 4.8, length.out = 8), effort = 2) smooth_estimates(m, data = nd["x"]) fitted_values(m, data = nd) ``` Automatic smooth grids are evenly spaced in the **smooth coordinate**, here `log(x)`, and the output and plot axes use that name. They use stored evaluated covariates, so they do not need to reconstruct an inverse transformation. ```{r coordinates} sm <- smooth_estimates(m, n = 30) draw(sm) ``` An explicitly evaluated column is also accepted: ```{r evaluated-data} evaluated <- data.frame(`log(x)` = log(nd$x), check.names = FALSE) smooth_estimates(m, data = evaluated) ``` If raw inputs and evaluated columns are both supplied, expressions are recomputed from the raw inputs. Stored model frames and automatically prepared grids retain their evaluated coordinates. `data_slice()` works with raw variables. Supply reference `data` if the original raw covariates cannot be recovered; transformed values are not renamed as raw variables or automatically inverted. ```{r raw-slice} data_slice(m, x = evenly(x, n = 8), effort = 2, data = dat) ``` ## Choosing the plotted covariate scale Drawing methods use the model's covariate coordinates by default. For example, `log(x)` is plotted against `log(x)`, with a matching axis label, whether it is inside a smooth or a parametric term. Set `back_transform = TRUE` to display `x` instead. This changes the plotting coordinates and labels; the estimated partial effect and its interval keep their meaning. Use `fun` to transform the partial effect itself. ```{r back-transform} draw(m, back_transform = TRUE) sm_raw <- smooth_estimates(m, back_transform = TRUE) draw(sm_raw, back_transform = TRUE) m_parametric <- gam(y ~ log(x), family = poisson(), data = dat) draw(m_parametric) draw(m_parametric, back_transform = TRUE) pe <- parametric_effects(m_parametric) draw(pe, back_transform = TRUE) ``` The same rule applies to each covariate in tensor products and parametric interactions, including covariates used for conditioning facets. Rugs and partial-residual locations follow the displayed scale. Surface coordinates can become irregularly spaced after changing scale; these plots use rectangles bounded by adjacent coordinate midpoints to avoid shifting raster pixels. ```{r back-transform-surface} dat_surface <- transform(dat, z = runif(nrow(dat), 1, 4)) m_surface <- gam(y ~ te(log(x), sqrt(z)), family = poisson(), data = dat_surface, method = "REML") draw(m_surface, back_transform = TRUE, contour = FALSE) ``` Smooth grids are evenly spaced on the model scale by default. For example, three points for `s(log(x))` with `x` ranging from 1 to 100 correspond to raw values 1, 10, and 100. With `smooth_estimates(..., back_transform = TRUE)`, the grid is evenly spaced in raw `x`: 1, 50.5, and 100. The expressions are then evaluated forwards and estimates computed at those points. `draw(m, back_transform = TRUE)` performs this preparation automatically. Precomputed estimates retain their chosen evaluation points. An object created on a transformed-only grid cannot be drawn on the original scale: its raw coordinates are unavailable. Recompute it with `back_transform = TRUE`, or supply raw prediction rows through `data`. Returned estimates retain model coordinates and auxiliary `.gratia_raw_*` columns, so either drawing scale can be selected without changing the estimates. Subsets and confidence intervals retain these coordinates. There is no list of supported inverse functions. Each smooth axis must depend on one distinct original model covariate, and gratia must be able to recover that covariate and evaluate the expression. For example: ```{r back-transform-alias} foo <- sqrt m_foo <- gam(y ~ s(foo(x), k = 6), family = poisson(), data = dat) sm_foo <- smooth_estimates(m_foo, back_transform = TRUE, envir = environment()) draw(sm_foo, back_transform = TRUE) ``` Raw fitting rows are recovered and checked against the stored fitting observations. If the fitting data are unavailable, or changed functions or constants no longer reproduce the fit, recovery fails. Supply raw prediction rows in `data` and functions or constants in `envir` as needed. Already prepared estimates can be drawn without recovering data or evaluating functions again. When prediction `data` are supplied, unavailable raw fitting observations omit rugs and partial residuals. Prediction rows are not substituted for observed data. `envir` is used for both predictions and observation recovery. Typical values added by `smooth_data(..., include_all = TRUE)` do not identify the original coordinates of a model-scale grid. Prepare that grid with `back_transform = TRUE` to retain raw axes. For numeric `by` smooths, automatic grids evaluate the smooth with multiplier one, even for `s(log(x), by = x)`; the original axis is retained separately. Supplied raw prediction rows use their actual `by` values. Soap-film boundary coordinates are stored on the model scale. A transformed boundary has no retained raw coordinates, so raw-scale drawing raises an error. Draw these soap films with `back_transform = FALSE` to retain their boundaries. Parametric effects likewise retain both coordinate scales by row, including for many-to-one functions. Multi-column basis terms such as `poly(x, 2)` and parametric expressions with no single coordinate per input retain raw axes. The old `parametric_effects(transform = ...)` argument is deprecated. Replace `transform = TRUE` with `back_transform = FALSE`, and `transform = FALSE` with `back_transform = TRUE`. ## Local functions A fitted model may retain the formula expression without retaining the local environment in which it was evaluated. Preserve that environment when fitting if you will need a local function at new observations. ```{r local-function} fit_local <- function(dat) { shift <- 0.5 transform_x <- function(x) log(x + shift) f <- y ~ s(transform_x(x), k = 6) + offset(log(effort)) list(model = gam(f, family = poisson(), data = dat, method = "REML"), envir = environment(f)) } bundle <- fit_local(dat) smooth_estimates(bundle$model, data = nd, envir = bundle$envir) fitted_values(bundle$model, data = nd, envir = bundle$envir) ``` Automatic smooth plots can still use stored `transform_x(x)` values without this environment. New raw observations require the function. Evaluation errors identify the expression and explain how to supply the missing context; stored training values are never substituted for arbitrary new observations. ## Which derivative? For a contribution $f(\log x)$, the default derivative is with respect to $\log x$. Set `wrt = "covariate"` to differentiate with respect to $x$ instead. The latter includes the chain-rule factor $1/x$ because the transformation is reevaluated after each perturbation of $x$. ```{r derivatives} derivatives(m, data = nd, type = "central") derivatives(m, data = nd, type = "central", wrt = "covariate") ``` For a smooth of several raw inputs, such as `s(I(x / z))`, name the raw `focal` variable explicitly. Other inputs are held at the supplied values, or at typical values when generating a slice. Use `I()` for arithmetic that R would otherwise interpret as formula syntax. For tensor smooths, use `partial_derivatives()` and specify the smooth or raw focal coordinate consistently with `wrt`: **`wrt` chooses the differentiation scale, and `focal` names the coordinate to differentiate along on that scale**. For example, `te(log(x), sqrt(z))` represents a contribution $f(\log x, \sqrt z)$. The two choices have the following meanings: | Arguments | Derivative | |:----------|:-----------| | `wrt = "smooth", focal = "log(x)"` | With respect to $\log x$, holding $\sqrt z$ fixed | | `wrt = "covariate", focal = "x"` | With respect to $x$, holding $z$ fixed | ```{r partial-derivatives} set.seed(43) tensor_dat <- data.frame(x = runif(150, 1, 5), z = runif(150, 1, 3)) tensor_dat$y <- sin(log(tensor_dat$x)) * sqrt(tensor_dat$z) + rnorm(nrow(tensor_dat), sd = 0.1) tensor_model <- gam(y ~ te(log(x), sqrt(z), k = c(4, 4)), data = tensor_dat, method = "REML") # Vary x along a slice with z held fixed. tensor_nd <- data.frame(x = seq(1.2, 4.8, length.out = 8), z = 2) pd_smooth <- partial_derivatives(tensor_model, data = tensor_nd, wrt = "smooth", focal = "log(x)", type = "central") pd_covariate <- partial_derivatives(tensor_model, data = tensor_nd, wrt = "covariate", focal = "x", type = "central") # Compare the raw-covariate derivative with the chain-rule calculation. data.frame(x = tensor_nd$x, covariate_derivative = pd_covariate$.partial_deriv, smooth_derivative_over_x = pd_smooth$.partial_deriv / tensor_nd$x) ``` The first call perturbs `log(x)` directly. The second perturbs `x` and recomputes the smooth expressions, incorporating the chain rule numerically. Thus, for these first derivatives, $\partial f / \partial x = (\partial f / \partial \log x) / x$, up to finite-difference approximation error. In `partial_derivatives()`, `wrt = "smooth"` is the default, and omitting `focal` selects the first non-factor smooth coordinate. With `wrt = "covariate"`, supply `focal` whenever several raw covariates are eligible; for this tensor smooth, choose `"x"` or `"z"`. If selecting several smooths, supply one `focal` per smooth in their selected order; the same `wrt` applies to all of them. Supplied `data` should describe a slice with the other independent numeric inputs held fixed, as `z` is above. Raw-scale differentiation requires the raw focal column: stored transformed values alone are insufficient. For an untransformed coordinate such as `x` in `te(x, z)`, both choices of `wrt` give the same interpretation. ### Choosing a finite-difference step The derivative functions use `eps = NULL` by default to choose a step from the fitted range of the differentiation coordinate. This is the range of `log(x)` when differentiating on that smooth coordinate, or of `x` when using `wrt = "covariate"`. The step therefore follows changes in the covariate's units and normally stays the same when prediction data are split into batches. For derivative order $d$, the range is multiplied by $\epsilon_{\mathrm{machine}}^{1/(d+p)}$, where $p=2$ for central differences and $p=1$ for forward or backward differences. This balances the orders of rounding and approximation error; for example, second central differences use the fourth root of machine precision. Constant coordinates use their magnitude instead of their range, with a fallback of one at zero. A lower bound and rounding of the step ensure representable perturbations at large coordinates. This rule does not estimate a universally optimal step: that would also need information about higher derivatives of the fitted function. Very small steps can amplify rounding errors, particularly for second derivatives, while large steps increase approximation error. For unusual function scales or domain boundaries, supply a positive absolute `eps` and check stability at nearby step sizes. See the [FiniteDiff step-size discussion](https://docs.sciml.ai/FiniteDiff/dev/epsilons/) for the error-balance principle.