Package 'gratia'

Description

Add a confidence interval to an existing object

Usage

add_confint(object, coverage = 0.95, ...)

## S3 method for class 'smooth_estimates'
add_confint(object, coverage = 0.95, ...)

## S3 method for class 'parametric_effects'
add_confint(object, coverage = 0.95, ...)

## Default S3 method:
add_confint(object, coverage = 0.95, ...)

Arguments

Details

For smooth estimates, this function computes Gaussian pointwise intervals from the saved estimates and standard errors. It replaces any existing limits and updates their interval metadata. To retain stored bootstrap, shape-enforced or simultaneous limits, draw the estimates directly with draw() without calling add_confint().

Description

Add a constant to estimated values

Usage

add_constant(object, constant = NULL, ...)

## S3 method for class 'smooth_estimates'
add_constant(object, constant = NULL, ...)

## S3 method for class 'smooth_samples'
add_constant(object, constant = NULL, ...)

## S3 method for class 'mgcv_smooth'
add_constant(object, constant = NULL, ...)

## S3 method for class 'parametric_effects'
add_constant(object, constant = NULL, ...)

## S3 method for class 'tbl_df'
add_constant(object, constant = NULL, column = NULL, ...)

## S3 method for class 'evaluated_parametric_term'
add_constant(object, constant = NULL, ...)

Arguments

Value

Returns object but with the estimate shifted by the addition of the supplied constant.

Author(s)

Gavin L. Simpson

Description

Add fitted values from a model to a data frame

Usage

add_fitted(data, model, value = ".fitted", ...)

## S3 method for class 'gamm4'
add_fitted(data, model, ...)

## S3 method for class 'scam'
add_fitted(data, model, ...)

Arguments

Value

A data frame (tibble) formed from data and fitted values from model, stored in ".fitted" by default.

See Also

gamm4_support

Description

Adds draws from the posterior distribution of model to the data object using one of fitted_samples(), predicted_samples(), or posterior_samples().

Usage

add_fitted_samples(
  data,
  model,
  ...,
  n_draws = NULL,
  seed = NULL,
  chains = NULL
)

add_predicted_samples(data, model, ..., n_draws = NULL, seed = NULL)

add_posterior_samples(
  data,
  model,
  ...,
  n_draws = NULL,
  seed = NULL,
  chains = NULL
)

add_smooth_samples(
  data,
  model,
  ...,
  n_draws = NULL,
  seed = NULL,
  select = NULL,
  chains = NULL
)

Arguments

Examples

load_mgcv()

df <- data_sim("eg1", n = 400, seed = 42)

m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")

# add fitted samples (posterior draws of the expected value of the response)
# note that there are 800 rows in the output: 400 data by `n_draws = 2` samples.
df |>
  add_fitted_samples(m, n_draws = 2, seed = 84)

# add posterior draws from smooth s(x2)
df |>
  add_smooth_samples(m, n_draws = 2, seed = 2, select= "s(x2)")

Description

Add fitted values from a GAM to a data frame

Usage

## S3 method for class 'gam'
add_fitted(
  data,
  model,
  value = ".fitted",
  type = "response",
  ...,
  envir = NULL,
  marginalise = NULL,
  integration = c("quadrature", "monte_carlo"),
  n_quad = 40L,
  n_integrate = 5000L,
  seed = NULL
)

Arguments

Value

A data frame (tibble) formed from data and predictions from model.

Examples

load_mgcv()
df <- data_sim("eg1", seed = 1)
df <- df[, c("y", "x0", "x1", "x2", "x3")]
m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")

# add fitted values to our data
add_fitted(df, m)

# with type = "terms" or "iterms"
add_fitted(df, m, type = "terms")

Description

Add partial residuals

Usage

add_partial_residuals(data, model, ...)

## S3 method for class 'gam'
add_partial_residuals(data, model, select = NULL, partial_match = FALSE, ...)

## S3 method for class 'gamm4'
add_partial_residuals(data, model, ...)

## S3 method for class 'scam'
add_partial_residuals(data, model, ...)

Arguments

Details

Partial residuals require working residuals and working weights from the fitted model. Families that do not provide these, including the location-scale families in mgcv, are not supported, even for terms in the location predictor.

See Also

gamm4_support

Examples

load_mgcv()

df <- data_sim("eg1", seed = 1)
df <- df[, c("y", "x0", "x1", "x2", "x3")]
m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")

## add partial residuals
add_partial_residuals(df, m)

## add partial residuals for selected smooths
add_partial_residuals(df, m, select = "s(x0)")

Description

Add residuals from a model to a data frame

Usage

add_residuals(data, model, value = ".residual", ...)

## S3 method for class 'gamm4'
add_residuals(data, model, ...)

## S3 method for class 'scam'
add_residuals(data, model, ...)

Arguments

Value

A data frame (tibble) formed from data and residuals from model.

See Also

gamm4_support

Description

Add residuals from a GAM to a data frame

Usage

## S3 method for class 'gam'
add_residuals(data, model, value = ".residual", type = "deviance", ...)

Arguments

Value

A data frame (tibble) formed from data and residuals from model.

Examples

load_mgcv()

df <- data_sim("eg1", seed = 1)
df <- df[, c("y", "x0", "x1", "x2", "x3")]
m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")

##
add_residuals(df, m)

Description

Add indicators of significant change after SiZeR

Usage

add_sizer(object, type = c("change", "sizer"), ...)

## S3 method for class 'derivatives'
add_sizer(object, type = c("change", "sizer"), ...)

## S3 method for class 'smooth_estimates'
add_sizer(object, type = c("change", "sizer"), derivatives = NULL, ...)

Arguments

Examples

load_mgcv()

df <- data_sim("eg1", n = 400, dist = "normal", scale = 2, seed = 42)
m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")

## first derivatives of all smooths using central finite differences
d <- derivatives(m, type = "central") |>
  add_sizer()

# default adds a .change column
names(d)

Description

Calculate one-dimensional accumulated local effects (ALE) of a numeric covariate, combining all model terms involving that covariate.

Usage

ale(object, ...)

## Default S3 method:
ale(object, ...)

## S3 method for class 'gam'
ale(
  object,
  focal,
  data = NULL,
  scale = c("link", "response", "linear_predictor"),
  n_bins = 40L,
  breaks = NULL,
  parameter = NULL,
  which_eta = NULL,
  center = TRUE,
  interval = c("confidence", "simultaneous"),
  ci_level = 0.95,
  unconditional = FALSE,
  frequentist = FALSE,
  n_sim = 10000L,
  n_cores = 1L,
  seed = NULL,
  block_size = 10000L,
  exclude = NULL,
  envir = NULL,
  ...
)

## S3 method for class 'bam'
ale(object, ...)

## S3 method for class 'scam'
ale(object, ..., unconditional = FALSE)

Arguments

Details

Within each interval, predict at its two boundaries for every reference observation in that interval, retaining its other covariates. Average the upper-minus-lower prediction differences, then accumulate these averages. Intervals are right closed, with the minimum included in the first interval. The returned curve is linearly interpolated between boundaries. Centring subtracts its mean evaluated at the reference observations, rather than an unweighted mean over boundaries. Fitting weights do not weight these averages.

On the response scale, transform endpoint predictions before differencing. Transforming a link-scale ALE curve does not produce response-scale ALE. Response-scale effects are differences in natural parameter units, not predicted parameter values, and can be negative.

Distributional families gaulss, gammals, gumbls, gevlss, shash, twlss, and ziplss are supported. Natural units mean standard deviation for Gaussian scale, positive scale or dispersion for logged scales, natural Tweedie power, positive sinh-arcsinh tail parameter, and Poisson rate or probability of presence for ziplss. Location parameters need not be response means. Overall means combining several parameters, multinomial probabilities, multivariate covariance parameters and grouped-family outputs are unsupported. For other families, link-scale ALE is available for a single linear predictor; response-scale ALE requires a supported scalar mean family.

Supported response-scale parameters are:

• gaulss: location (mean) and scale (standard deviation);

• gammals: location (mean) and scale (dispersion);

• gumbls and gevlss: location and positive scale, plus GEV shape;

• shash: location, positive scale, skewness and positive tail parameter (named kurtosis to match gratia's parameter names);

• twlss: location (mean), power and scale (dispersion);

• ziplss: location (Poisson rate) and pi (probability of presence).

Raw-data recovery can require envir when a model was fitted locally and its formula no longer retains the environment containing the fitting data. Alternatively, supply reference data explicitly.

Coefficient gradients are averaged, accumulated and centred in the same way as the effects, retaining covariance across intervals and the centring constant. Link-scale gradients are exact; response-scale intervals use the delta method. Both interval types condition on the reference observations and boundaries. Simultaneous coverage concerns the returned boundaries, not all points between them. No uncertainty for the reference distribution is added.

Formula offsets are included and reevaluated at interval boundaries. Offsets supplied only through the fitting function's offset argument are omitted, following explicit-newdata prediction semantics. BAM predictions use discrete = FALSE, including for models fitted with discrete = TRUE.

Prediction blocks use an existing pool of workers configured with mirai::daemons(), when mirai and carrier are installed. Set block_size small enough to give each worker a block; serial and parallel calculations retain the same block boundaries and accumulation order. Worker startup and model transfer can outweigh the benefit for small models or reference data. The caller owns the pool; ale() does not start or stop workers. n_cores controls coefficient simulation for simultaneous intervals, not the number of prediction workers. Workers use their configured BLAS thread limits.

Data recovery and validation, centring, covariance calculations and interval simulation remain in the caller. Simultaneous bands use one shared set of coefficient deviations for all boundaries, preserving results and the caller's random-number state across worker counts. Models requiring custom prediction methods, family callbacks, contrasts or formula functions that workers cannot safely evaluate retain serial prediction. Standard arithmetic and common base and stats transformations can run on workers; local helpers, stochastic expressions, custom S3 transformation methods and other package functions retain serial prediction. Method checks do not invoke active or lazy caller bindings.

This implementation is written independently from the published ALE definition, without using code from mgcViz or ALEPlot.

Value

A tibble of class "ale", with .focal and .parameter identifying the curve, .value giving boundary coordinates, .ale the effects, .se their standard errors, .crit the interval critical value, and .lower_ci and .upper_ci the interval limits. .bin identifies the interval ending at each boundary and .n counts its reference rows. At the initial boundary both are zero. The ale_info attribute records computation settings without retaining the model or reference data.

References

Apley, D. W. and Zhu, J. (2020). Visualizing the effects of predictor variables in black box supervised learning models. Journal of the Royal Statistical Society: Series B, 82, 1059–1086. doi:10.1111/rssb.12377.

See Also

draw.ale(), conditional_values(), response_derivatives()

Examples

load_mgcv()
df <- data_sim("eg1", n = 200, seed = 42)
m <- gam(y ~ s(x0) + te(x1, x2, k = c(5, 5)) + s(x3),
  data = df, method = "REML")
a <- ale(m, focal = "x1")
draw(a)

# Standard deviation in its natural units, rather than reciprocal SD
ml <- gam(list(y ~ s(x1), ~ x2), data = df, family = gaulss())
ale(ml, "x2", parameter = "scale", scale = "response") |> draw()

Description

Model diagnostic plots

Usage

appraise(model, ...)

## S3 method for class 'gam'
appraise(
  model,
  method = NULL,
  use_worm = FALSE,
  n_uniform = 10,
  n_simulate = 50,
  seed = NULL,
  type = NULL,
  n_bins = c("sturges", "scott", "fd"),
  ncol = NULL,
  nrow = NULL,
  guides = "keep",
  level = 0.9,
  ci_col = "black",
  ci_alpha = 0.2,
  point_col = "grey20",
  point_alpha = 1,
  line_col = "red",
  by_family = NULL,
  rank = NULL,
  scalable = NULL,
  scalable_threshold = 10000L,
  qq_bins = 1000L,
  density_bins = c(50L, 50L),
  ppcheck = FALSE,
  n_ppc = 50,
  density_method = c("kde", "histogram"),
  bandwidth = NULL,
  conditioning = NULL,
  ppc_batch_size = 10L,
  ...,
  n_dharma = 250L
)

## S3 method for class 'lm'
appraise(model, ...)

## S3 method for class 'gamlss'
appraise(model, ..., type = "quantile", seed = NULL, conditioning = NULL)

## Default S3 method:
appraise(
  model,
  ...,
  type = NULL,
  seed = NULL,
  conditioning = NULL,
  n_dharma = 250L
)

## S3 method for class 'gamm4'
appraise(model, conditioning = NULL, ...)

## S3 method for class 'gjrm'
appraise(
  model,
  ...,
  type = "quantile",
  seed = NULL,
  conditioning = NULL,
  response = NULL,
  ppcheck = FALSE,
  by_family = NULL,
  n_dharma = 250L
)

## S3 method for class 'gamm'
appraise(model, conditioning = NULL, ...)

## S3 method for class 'gratia_diagnostic_proxy'
appraise(model, type = "quantile", ...)

## S3 method for class 'lme'
appraise(model, conditioning = NULL, ...)

## S3 method for class 'merMod'
appraise(model, conditioning = NULL, ...)

## S3 method for class 'sdmTMB'
appraise(
  model,
  conditioning = NULL,
  component = "response",
  field_method = c("mle-mvn", "mle-eb", "mle-mvn-integrated"),
  ...
)

Arguments

Details

Use type = "dharma" for simulation-based PIT residuals calculated internally using gratia's response generators; no DHARMa or mgcViz dependency is needed. These residuals lie in ⁠[0, 1]⁠, including exact endpoints, and do not require a response CDF. n_dharma controls their simulation count. All residual panels share one calculated vector. QQ and worm plots require method = NULL and use uniform reference quantiles, with approximate pointwise bands for independent responses. Bands are omitted for dependent responses. Censored and multivariate outcomes are unsupported. See dharma_residuals() for the calculation.

Use type = "surrogate" for mgcv::ocat() and mgcv::multinom() GAMs. The former uses a scalar latent logistic residual; the latter retains all category utilities, including category 0, as separate panels with centred Gumbel references. QQ and worm plots require method = NULL. Each appraisal shares one residual draw across panels, with the same draw in ordinary and scalable displays. Multinomial residuals are plotted against their own category probabilities; observed-versus-fitted panels compare category indicators with those probabilities. Pointwise bands are calculated separately for each utility. See surrogate_residuals().

With type = NULL, ordinary lm/glm/GAM fits use deviance residuals. Models fitted with mgcv::gfam() use randomised quantile (Dunn–Smyth) residuals when every active component supplies a supported CDF and probability limits for discrete or censored observations. Otherwise, a message identifies unsupported components and all components use deviance residuals; requesting type = "quantile" explicitly is an error when any component is unsupported. Components never mix residual scales. Invalid encodings or numerical CDF failures are errors.

Quantile residuals use each component's fitted parameters, dispersion and weights, including binomial trial counts and Gaussian precision weights. Gamma (including precision weights), negative binomial, beta, scaled-t, Tweedie, inverse Gaussian, ordinal and zero-inflated Poisson components are supported. Unrecognised custom-family bounds remain unsupported. A custom CDF must also support the probabilities just below discrete values and at censoring limits. Replacements for recognised family CDFs must use the same arguments and probability conventions. Scalar uncensored cnorm(), clog() and cpois() components are supported; genuine censoring requires recoverable, aligned censoring bounds. Matrix censored responses rejected by gfam() and upstream fitting failures are not repaired.

Use response = NULL to check every grouped-family component or a single integer to select a component by its position in the mgcv::gfam() family list. Repeated families remain distinct. Selection filters fitting observations without refitting the model and retains original component labels. by_family independently controls the layout of the selected data. With an explicit seed, selected residuals match the corresponding rows of an all-component calculation whenever all components support that calculation.

by_family = NULL pools quantile residuals with colours identifying indexed family components, including repeated uses of the same family. Traditional residuals use component facets with free axes by default. The pooled residual-versus-predictor plot uses within-component predictor ranks: (rank(eta, ties.method = "average") - 0.5) / n, where n counts the finite observations included for that component; constants map to 0.5. Facets retain original linear predictors. Use rank to override this choice independently of residual type or layout. The links and original predictor values can have different meanings across components; ranks compare positions within components, not effects on a common response scale. Pooled histograms use common bin edges and component-specific densities.

Quantile QQ and worm plots compare to fixed standard-normal quantiles, with an identity line or zero line. Independent responses have exact pointwise reference bands for ordered residuals. For dependent responses, bands use simulations that preserve dependence when a joint generator is available; otherwise they are omitted with an explanation. Bands condition on fitted parameters, exclude estimation uncertainty and do not provide simultaneous coverage. References calculated from simulation ranks also use simulated bands; see diagnostic_support. Worm plots subtract the same theoretical quantiles from residuals and reference bands. appraise() computes one seeded residual vector shared by its residual panels and always separates observed-versus-fitted values by component, because response units differ. Standalone observed_fitted_plot() allows explicit pooling. Ordinal panels use expected category codes.

Scalable plots aggregate all observations before rendering. Residual heatmaps show conditional density; observed-versus-fitted heatmaps show counts. Point styling applies to QQ/worm points, and heatmap fill scales can be changed with ggplot scales. QQ references average draws as they are generated or compute exact pooled quantiles using temporary files above the configured memory threshold; see diagnostic_support. Supported univariate GJRM::gamlss models use quantile residuals by default.

If the response expectation does not exist, the observed-versus-fitted panel explains why it is unavailable. Other supported panels remain usable. Prepared diagnostic objects can be given a new conditioning or field choice; the model is prepared again so all panels use the requested choice.

Note

The wording used in mgcv::qq.gam() uses direct in reference to the simulated residuals method (method = "simulated"). To avoid confusion, method = "direct" is deprecated in favour of method = "uniform".

See Also

The plots are produced by functions qq_plot(), residuals_linpred_plot(), residuals_hist_plot(), observed_fitted_plot(), and optionally posterior_predictive_plot().

gamm4_support

Examples

load_mgcv()
# Simulate an additive Gaussian response and fit a GAM
dat <- data_sim("eg1", n = 200, scale = 2, seed = 2)
mod <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat, method = "REML")

# Inspect residuals, fitted values and the residual distribution
appraise(mod, seed = 42)

Description

Coefficient matrices retain their chain mapping in an attribute. Sampled quantities use their .chain, .iteration, and .draw columns. A plain matrix without a mapping is treated as one chain; chain boundaries are never inferred from its number of rows. Filtering must leave the same consecutive iterations in every chain. Missing observations remain missing variables.

Usage

as_mcmc_draws(
  x,
  format = c("draws_array", "draws_df", "mcmc.list", "array"),
  value = NULL
)

Arguments

Value

An iteration by chain by variable array, a posterior draws object, or a coda list of chains. Variable names identify sampled quantities and their observation rows, parameters, or smooths.

Description

Generic function for assembling plot objects created from R objects, using the ggplot2 package.

Usage

assemble(object, ..., label_control = NULL)

## S3 method for class 'gamm4'
assemble(object, ...)

Arguments

Value

A list of ggplot2::ggplot() objects.

Label layout

Pass label_control = draw_label_control(...) to draw() to wrap discrete axis labels and legend text, stagger discrete tick labels, or adjust legend and colour-bar layout. Settings apply to each panel, including unwrapped plot lists. See draw_label_control() for the available settings.

Author(s)

Gavin L. Simpson

See Also

gamm4_support

Description

Plots smooth and parametric contributions from a fitted GAM model using ggplot2::ggplot(). For models fitted by mgcv::gamm(), assemble() and draw() use the fitted gam component.

Usage

## S3 method for class 'gam'
assemble(
  object,
  data = NULL,
  select = NULL,
  parametric = NULL,
  terms = NULL,
  residuals = FALSE,
  scales = c("free", "fixed"),
  ci_level = 0.95,
  n = 100,
  n_2d = 50,
  n_3d = 16,
  n_4d = 4,
  unconditional = FALSE,
  overall_uncertainty = TRUE,
  constant = NULL,
  fun = NULL,
  dist = 0.1,
  rug = TRUE,
  distinct_rug = TRUE,
  contour = TRUE,
  grouped_by = deprecated(),
  ci_alpha = 0.2,
  ci_col = "black",
  smooth_col = "black",
  resid_col = "steelblue3",
  contour_col = "black",
  n_contour = NULL,
  partial_match = FALSE,
  discrete_colour = NULL,
  discrete_fill = NULL,
  continuous_colour = NULL,
  continuous_fill = NULL,
  position = "identity",
  angle = NULL,
  ncol = NULL,
  nrow = NULL,
  guides = "keep",
  widths = NULL,
  heights = NULL,
  crs = NULL,
  default_crs = NULL,
  lims_method = "geometry_bbox",
  wrap = TRUE,
  caption = TRUE,
  envir = NULL,
  surface = c("raster", "tile", "3d"),
  surface_3d = list(),
  ci_3d = FALSE,
  ci_3d_alpha = 0.15,
  ci_3d_col = "grey60",
  back_transform = FALSE,
  re_levels = NULL,
  mrf_plot = c("effects", "qq"),
  interval = c("confidence", "simultaneous"),
  n_sim = 10000,
  n_cores = 1,
  seed = NULL,
  interval_method = c("auto", "gaussian", "bootstrap", "shape"),
  frequentist = FALSE,
  shape_source = c("gaussian", "bootstrap"),
  ...,
  variables = NULL,
  simultaneous_method = c("simulation", "tube"),
  simultaneous_domain = NULL,
  tube_control = list(),
  cyclic = c("cartesian", "polar"),
  polar_args = list(),
  grouped = c("sz", "fs"),
  show_ci = c("factor-by", "sz"),
  time_varying_style = c("ribbon", "pointrange")
)

## S3 method for class 'gamm'
assemble(object, ...)

## S3 method for class 'scam'
assemble(object, ..., overall_uncertainty = FALSE, centre = c("model", "plot"))

Arguments

Details

Smooth and parametric plot subtitles identify the model part or distributional parameter as described in draw.gam(). Custom labels, including subtitle, can be passed through .... Use cyclic = "polar" to draw one-dimensional cyclic smooths with polar coordinates; polar_args customizes the coordinate display as described in draw.gam().

Value

A list of ggplot2::ggplot() objects.

Label layout

Pass label_control = draw_label_control(...) to draw() to wrap discrete axis labels and legend text, stagger discrete tick labels, or adjust legend and colour-bar layout. Settings apply to each panel, including unwrapped plot lists. See draw_label_control() for the available settings.

Note

Internally, plots of each smooth are created using ggplot2::ggplot() and composed into a single plot using patchwork::wrap_plots(). As a result, it is not possible to use + to add to the plots in the way one might typically work with ggplot() plots. Instead, use the & operator; see the examples.

Author(s)

Gavin L. Simpson

Examples

load_mgcv()

# simulate some data
df1 <- data_sim("eg1", n = 400, dist = "normal", scale = 2, seed = 2)
# fit GAM
m1 <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df1, method = "REML")

# assemble plots of all smooths
assemble(m1)

# can be passed to patchwork::wrap_plots()
suppressPackageStartupMessages(library("patchwork"))
assemble(m1) |> wrap_plots(ncol = 2, nrow = 2)

# Prepare individual plots from a GAMM
m2 <- gamm(y ~ s(x0) + x1, data = df1, method = "REML")
gamm_plots <- assemble(m2)
gamm_plots[[1]]

Description

Basis expansions from a definition of a smoother using the syntax of mgcv's smooths via mgcv::s()., mgcv::te(), mgcv::ti(), and mgcv::t2(), or directly from a fitted GAM(M).

Usage

basis(object, ...)

## S3 method for class 'gam'
basis(
  object,
  select = NULL,
  term = deprecated(),
  data = NULL,
  n = 100,
  n_2d = 50,
  n_3d = 16,
  n_4d = 4,
  partial_match = FALSE,
  envir = NULL,
  ...,
  variables = NULL,
  unweighted = FALSE
)

## S3 method for class 'scam'
basis(
  object,
  select = NULL,
  term = deprecated(),
  data = NULL,
  n = 100,
  n_2d = 50,
  n_3d = 16,
  n_4d = 4,
  partial_match = FALSE,
  ...,
  variables = NULL,
  unweighted = FALSE
)

## S3 method for class 'gamm'
basis(
  object,
  select = NULL,
  term = deprecated(),
  data = NULL,
  n = 100,
  n_2d = 50,
  n_3d = 16,
  n_4d = 4,
  partial_match = FALSE,
  ...,
  variables = NULL,
  unweighted = FALSE
)

## S3 method for class 'gamm4'
basis(
  object,
  select = NULL,
  term = deprecated(),
  data = NULL,
  n = 100,
  n_2d = 50,
  n_3d = 16,
  n_4d = 4,
  partial_match = FALSE,
  ...,
  variables = NULL,
  unweighted = FALSE
)

## Default S3 method:
basis(
  object,
  data,
  knots = NULL,
  constraints = FALSE,
  at = NULL,
  diagonalize = FALSE,
  coefficients = NULL,
  envir = NULL,
  ...,
  variables = NULL
)

Arguments

Details

When evaluating a fitted smooth at explicit data, missing covariates required by that smooth produce NA basis values. Missing values in other model variables do not discard otherwise evaluable rows. The covariate order chosen by variables also determines which automatically generated axes use n_2d, n_3d, and n_4d. Supplied data or at values are unchanged. Set a resolution argument to NULL to use n instead.

Value

A tibble with class basis. Its plot_variables attribute stores covariate orders by smooth label for subsequent plotting. When the basis is shared across factor levels, factor_by_groups records the smooth labels that draw() combines. The returned rows retain each level.

Author(s)

Gavin L. Simpson

Examples

load_mgcv()

df <- data_sim("eg4", n = 400, seed = 42)

bf <- basis(s(x0), data = df)
bf <- basis(s(x2, by = fac, bs = "bs"), data = df, constraints = TRUE)

Description

Extract basis dimension of a smooth

Usage

basis_size(object, ...)

## S3 method for class 'mgcv.smooth'
basis_size(object, ...)

## S3 method for class 'gam'
basis_size(object, ...)

## S3 method for class 'gamm'
basis_size(object, ...)

## S3 method for class 'gamm4'
basis_size(object, ...)

## S3 method for class 'scam'
basis_size(object, ...)

Arguments

See Also

gamm4_support

Examples

load_mgcv()

df <- data_sim("eg1", n = 200, seed = 1)
m <- bam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df)

basis_size(m)

Description

Data generated from a hypothetical study of bird movement along a migration corridor, sampled throughout the year. This dataset consists of simulated sample records of numbers of observed locations of 100 tagged individuals each from six species of bird, at ten locations along a latitudinal gradient, with one observation taken every four weeks. Counts were simulated randomly for each species in each location and week by creating a species-specific migration curve that gave the probability of finding an individual of a given species in a given location, then simulated the distribution of individuals across sites using a multinomial distribution, and subsampling that using a binomial distribution to simulation observation error (i.e. not every bird present at a location would be detected). The data set (bird_move) consists of the variables count, latitude, week and species.

Format

A data frame

Source

Pedersen EJ, Miller DL, Simpson GL, Ross N. 2018. Hierarchical generalized additive models: an introduction with mgcv. PeerJ Preprints 6:e27320v1 doi:10.7287/peerj.preprints.27320v1.

Description

Usage

## S3 method for class 'gamm4'
boundary(x, ...)

boundary(x, ...)

## S3 method for class 'soap.film'
boundary(x, ...)

## S3 method for class 'gam'
boundary(x, select, ...)

Arguments

Value

A list of lists or data frames specifying the loops that define the boundary of the soap film smooth.

See Also

gamm4_support

mgcv::soap

Description

Given a vector indexing the smooths of a GAM, returns a logical vector selecting the requested smooths.

Usage

check_user_select_smooths(
  smooths,
  select = NULL,
  partial_match = FALSE,
  model_name = NULL
)

Arguments

Value

A logical vector the same length as length(smooths) indicating which smooths have been selected.

Author(s)

Gavin L. Simpson

Description

Extract coefficients from a fitted scam model.

Usage

## S3 method for class 'scam'
coef(object, parametrized = TRUE, ...)

Arguments

Description

Compare smooths across models

Usage

compare_smooths(
  model,
  ...,
  select = NULL,
  smooths = deprecated(),
  n = 100,
  n_2d = 50,
  data = NULL,
  unconditional = FALSE,
  overall_uncertainty = TRUE,
  partial_match = FALSE,
  envir = NULL,
  dist = NULL
)

Arguments

Details

Two-variable smooths are evaluated at common locations. Automatic grids span the intersection of the models' observed covariate ranges; factor levels must overlap. Supplied data are evaluated as given, without filling gaps. Each selected smooth label must occur exactly once in every model and use the same covariates, including any by variable. Numeric and logical select refer to the primary model; matching labels are then used in every model.

Estimates are contributions of individual smooths to the linear predictor, on its link scale, with the centring constraints used in each fitted model. They exclude the intercept and other model terms. Different fitting data can change the centring, so a surface difference can include a constant shift. It is the difference between these fitted smooth contributions, rather than a difference between fitted responses.

The result retains one row per model and smooth, with estimates in the data list column. For distributional models, .parameter retains the label for the linear predictor, while .model identifies the input model. Use draw.compare_smooths() for faceted surface comparisons or estimated surface differences. Smooths with more than two covariates can be evaluated but cannot be drawn by this method.

Examples

load_mgcv()
dat <- data_sim("eg1", seed = 2)

## models to compare smooths across - artificially create differences
m1 <- gam(y ~ s(x0, k = 5) + s(x1, k = 5) + s(x2, k = 5) + s(x3, k = 5),
  data = dat, method = "REML"
)
m2 <- gam(y ~ s(x0, bs = "ts") + s(x1, bs = "ts") + s(x2, bs = "ts") +
  s(x3, bs = "ts"), data = dat, method = "REML")

## build comparisons
comp <- compare_smooths(m1, m2)
comp
## notice that the result is a nested tibble

draw(comp)

# Simulate a noisy response over two continuous covariates, x and z
surf <- data_sim("eg2", n = 100, seed = 7)
# Fit two small bases to the same surface data
coarse <- gam(y ~ s(x, z, k = 10), data = surf, method = "REML")
fine <- gam(y ~ s(x, z, k = 15), data = surf, method = "REML")
# Compare surfaces on shared axes and a shared colour scale
surfaces <- compare_smooths(coarse, fine, n_2d = 20)
draw(surfaces)
# Show coarse minus fine, with a colour scale centred on zero
draw(surfaces, show = "difference", pair = c("coarse", "fine"))

Description

Compare full fitted values between pairs of factor levels, conditional on covariate values. Parametric effects and smooths are included unless excluded explicitly. Unlike difference_smooths(), these are differences of predictions, rather than differences of selected smooth functions.

Usage

conditional_differences(
  model,
  by,
  condition = NULL,
  data = NULL,
  scale = c("response", "link", "linear_predictor"),
  ...
)

## S3 method for class 'gam'
conditional_differences(
  model,
  by,
  condition = NULL,
  data = NULL,
  scale = c("response", "link", "linear_predictor"),
  n_vals = NULL,
  ci_level = 0.95,
  complete = TRUE,
  unconditional = FALSE,
  uncertainty = c("delta", "simulation"),
  n_sim = 10000,
  seed = NULL,
  n_cores = 1,
  method = c("gaussian", "mh", "inla", "user"),
  draws = NULL,
  mvn_method = c("mvnfast", "mgcv"),
  burnin = 1000,
  thin = 1,
  t_df = 40,
  rw_scale = 0.25,
  envir = NULL,
  ...,
  interval = c("confidence", "simultaneous"),
  simultaneous_scope = c("contrast", "all")
)

## S3 method for class 'gamm4'
conditional_differences(model, ...)

## S3 method for class 'scam'
conditional_differences(model, ..., unconditional = FALSE)

Arguments

Details

Differences are group 1 minus group 2, with groups ordered by model factor levels. Both groups use the same values of all other covariates. On the response scale the inverse link is applied to each prediction before subtraction. Shared terms can therefore affect response-scale differences even when they cancel on the link scale.

Intervals describe uncertainty in conditional mean differences, not observation noise. Simulation uses fitted_samples() once for the entire prediction grid, retaining dependence between predictions by subtracting within draws. The estimate is always the difference of fitted predictions. Pointwise simulation intervals need not be centred on this estimate.

Simultaneous intervals use the ci_level quantile of the maximum absolute standardized deviation over the selected coverage set. They are symmetric about .diff, with half-width .crit * .se. For uncertainty = "delta", one shared batch of zero-mean Gaussian coefficient deviations is generated from the Bayesian covariance matrix. Response-scale differences use a first-order Taylor approximation. The covariance must be positive semidefinite; singular covariance matrices are supported. Sampler arguments method, draws, mvn_method, burnin, thin, t_df, and rw_scale apply only to uncertainty = "simulation".

For uncertainty = "simulation", deviations are posterior difference draws minus .diff, standardized by their sample standard deviations. Draws are transformed to the response scale before differencing; invalid inverse-link domains or response means cause an error rather than discarded draws. Zero-variance draws that disagree with .diff cannot calibrate a band and cause an error. Deterministic rows agreeing with .diff have zero width.

Simultaneous coverage is approximate posterior coverage on the finite evaluation grid, not between its points or outside its range, nor exact frequentist coverage. Increase n_vals to evaluate a denser grid. Bands are not clipped to the response-difference range (e.g. ⁠[-1, 1]⁠ for probabilities). They need not contain the equal-tailed pointwise simulation intervals.

exclude sets the named term contributions to zero. In particular, excluding random effects does not integrate over their distribution. Offsets follow mgcv::predict.gam() semantics. Supported models are gam and bam fits with a single linear predictor and an ordinary scalar inverse-link mean (including negative binomial, Tweedie, beta regression and scaled t families). SCAM fits with supported single-predictor families also support link- and response-scale differences, with delta or simulation uncertainty and pointwise or simultaneous intervals. SCAM delta intervals use the local covariance of transformed coefficients. Gaussian simulation draws unrestricted fitting coefficients and applies the fitted positivity transformation; user draws must already be in transformed coordinates and satisfy the coefficient constraints. MH sampling is unavailable for SCAM fits, and smoothing-parameter uncertainty is unsupported; use unconditional = FALSE.

Value

A tibble of class "conditional_differences", with .contrast, .diff, .se, .lower_ci, .upper_ci, and evaluation covariates. .se is on the requested scale; for simulation it is the standard deviation of difference draws. A single comparison factor has .level_1 and .level_2 columns. Multiple factors have ⁠<factor>_1⁠ and ⁠<factor>_2⁠ columns. Attributes record by, condition, scale, uncertainty, ci_level, exclude, interval, simultaneous_scope, and plotting channels. Simultaneous results also contain .crit, the calibrated critical value, constant within each coverage set.

See Also

gamm4_support

Examples

load_mgcv()
# Simulate group-specific curves with Gaussian noise and fit factor-by smooths
df <- data_sim("eg4", n = 300, seed = 2)
m <- gam(y ~ fac + s(x2, by = fac) + s(x0), data = df, method = "REML")

# Compare fitted responses over x2, including the parametric group effects
cd <- conditional_differences(m, by = "fac", condition = "x2")
cd
draw(cd)

Description

Construct covariate combinations using conditional_values() conventions, then estimate response moments with fitted_moments(). Other covariates are held at typical values; observations are not averaged over those covariates.

Usage

conditional_moments(
  model,
  condition = NULL,
  data = NULL,
  moments = c("mean", "variance", "cv"),
  n_vals = 100,
  complete = TRUE,
  ...
)

## S3 method for class 'gam'
conditional_moments(
  model,
  condition = NULL,
  data = NULL,
  moments = c("mean", "variance", "cv"),
  n_vals = 100,
  complete = TRUE,
  ...
)

Arguments

Value

A tibble of class "conditional_moments", with the columns returned by fitted_moments() and the plotting channels of conditional_values().

Examples

load_mgcv()
m <- gam(mpg ~ s(wt, k = 5), data = mtcars, method = "REML")
conditional_moments(m, condition = "wt", moments = c("mean", "variance")) |>
  draw()

Description

Generate predicted values from a GAM, conditional upon supplied values of covariates. conditional_values() is modelled after marginaleffects::plot_predictions(), but with an intentionally simpler, more restrictive functionality. The intended use case is for quickly visualizing predicted values from a fitted GAM on the response scale. For more complex model predictions, you are strongly encouraged to use marginaleffects::plot_predictions().

Usage

conditional_values(
  model,
  condition = NULL,
  data = NULL,
  scale = c("response", "link", "linear_predictor"),
  ...
)

## S3 method for class 'gam'
conditional_values(
  model,
  condition = NULL,
  data = NULL,
  scale = c("response", "link", "linear_predictor"),
  n_vals = NULL,
  ci_level = 0.95,
  complete = TRUE,
  ...
)

## S3 method for class 'gamm4'
conditional_values(model, ...)

## S3 method for class 'scam'
conditional_values(model, ..., unconditional = FALSE)

Arguments

Value

A data frame (tibble) of class "conditional_values".

Author(s)

Gavin L. Simpson

See Also

gamm4_support

Examples

load_mgcv()
# Simulate an additive Gaussian response and fit a GAM
df <- data_sim("eg1", n = 200, seed = 2)
m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")

# Vary x2 while holding the other covariates at representative values
cv <- conditional_values(m, condition = "x2")
cv
draw(cv)

Description

Calculates point-wise confidence or simultaneous intervals for the first derivatives of smooth terms in a fitted GAM.

Usage

## S3 method for class 'fderiv'
confint(
  object,
  parm,
  level = 0.95,
  type = c("confidence", "simultaneous"),
  nsim = 10000,
  ncores = 1L,
  ...
)

Arguments

Value

a data frame with components:

1. term; factor indicating to which term each row relates,

2. lower; lower limit of the confidence or simultaneous interval,

3. est; estimated derivative

4. upper; upper limit of the confidence or simultaneous interval.

Author(s)

Gavin L. Simpson

Examples

load_mgcv()

dat <- data_sim("eg1", n = 1000, dist = "normal", scale = 2, seed = 2)
mod <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat, method = "REML")

# new data to evaluate the derivatives at, say over the middle 50% of range
# of each covariate
middle <- function(x, n = 25, coverage = 0.5) {
  v <- (1 - coverage) / 2
  q <- quantile(x, prob = c(0 + v, 1 - v), type = 8)
  seq(q[1], q[2], length = n)
}
new_data <- sapply(dat[c("x0", "x1", "x2", "x3")], middle)
new_data <- data.frame(new_data)
## first derivatives of all smooths...
fd <- fderiv(mod, newdata = new_data)

## point-wise interval
ci <- confint(fd, type = "confidence")
ci

## simultaneous interval for smooth term of x2

x2_sint <- confint(fd,
  parm = "x2", type = "simultaneous",
  nsim = 10000, ncores = 2
)

x2_sint

Description

Calculates point-wise confidence or simultaneous intervals for the smooth terms of a fitted GAM.

Usage

## S3 method for class 'gam'
confint(
  object,
  parm,
  level = 0.95,
  data = newdata,
  n = 100,
  type = c("confidence", "simultaneous"),
  nsim = 10000,
  shift = FALSE,
  transform = FALSE,
  unconditional = FALSE,
  frequentist = FALSE,
  overall_uncertainty = TRUE,
  ncores = 1,
  partial_match = FALSE,
  ...,
  newdata = NULL
)

## S3 method for class 'gamm'
confint(object, ...)

## S3 method for class 'gamm4'
confint(object, ...)

## S3 method for class 'scam'
confint(
  object,
  parm = smooths(object),
  ...,
  unconditional = FALSE,
  overall_uncertainty = FALSE
)

Arguments

Details

For SCAM fits, confint() returns intervals for smooth functions, with parm selecting smooth labels. It uses transformed-coefficient Gaussian uncertainty conditional on the fitted smoothing parameters. By default, intercept uncertainty is excluded. These intervals do not enforce shape constraints on their limits.

Value

a tibble with components:

1. .smooth; character indicating to which term each row relates,

2. .type; the type of smooth,

3. .by the name of the by variable if a by smooth, NA otherwise,

4. one or more vectors of values at which the smooth was evaluated, named as per the variables in the smooth,

5. zero or more variables containing values of the by variable,

6. .estimate; estimated value of the smooth,

7. .se; standard error of the estimated value of the smooth,

8. .crit; critical value for the 100 * level% confidence interval; NA for bootstrap limits.

9. .lower_ci; lower limit of the confidence or simultaneous interval,

10. .upper_ci; upper limit of the confidence or simultaneous interval,

Note

overall_uncertainty applies to both pointwise and simultaneous intervals. Each simultaneous band covers the evaluation points of one fitted smooth by default, or its continuous domain with tube calibration; separate smooths do not form a joint coverage set. Native scasm fits default to overall_uncertainty = FALSE and use stored bootstrap draws for pointwise intervals when available, as in smooth_estimates(). Shifting and transforming a smooth preserves its interval limits, including asymmetric bootstrap limits. Shape-enforced intervals describe unshifted smooth components on the model scale; shift = TRUE and effect transformations are unsupported.

Author(s)

Gavin L. Simpson

Examples

load_mgcv()

dat <- data_sim("eg1", n = 1000, dist = "normal", scale = 2, seed = 2)
mod <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat, method = "REML")

# new data to evaluate the smooths at, say over the middle 50% of range
# of each covariate
middle <- function(x, n = 50, coverage = 0.5) {
  v <- (1 - coverage) / 2
  q <- quantile(x, prob = c(0 + v, 1 - v), type = 8)
  seq(q[1], q[2], length = n)
}
new_data <- sapply(dat[c("x0", "x1", "x2", "x3")], middle)
new_data <- data.frame(new_data)

## point-wise interval for smooth of x2
ci <- confint(mod, parm = "s(x2)", type = "confidence", data = new_data)
ci

Description

All combinations of factor levels plus typical values of continuous variables

Usage

data_combos(object, ...)

## S3 method for class 'gam'
data_combos(object, vars = everything(), complete = TRUE, data = NULL, ...)

## S3 method for class 'gamm4'
data_combos(object, ...)

Arguments

See Also

gamm4_support

Description

A tidy reimplementation of the functions implemented in mgcv::gamSim() that can be used to fit GAMs. The sampling distribution can be applied to all the example types.

Usage

data_sim(
  model = "eg1",
  n = 400,
  scale = NULL,
  theta = 3,
  power = 1.5,
  dist = c("normal", "poisson", "binary", "negbin", "tweedie", "gamma", "ocat",
    "ordered categorical"),
  n_cat = 4,
  cuts = c(-1, 0, 5),
  seed = NULL,
  gfam_families = c("binary", "tweedie", "normal"),
  family = NULL,
  predictors = NULL,
  parameters = list(),
  observation = NULL
)

Arguments

Details

data_sim() can simulate data from several underlying models of known true functions. Descriptive names and their original short names are interchangeable, giving identical output for the same arguments and seed. The available options currently are:

• "additive" or "eg1": a four term additive true model. This is the classic Gu & Wahba four univariate term test model. See gw_functions for more details of the underlying four functions.

• "bivariate_surface" or "eg2": a bivariate smooth true model.

• "varying_coefficient" or "eg3": an example containing a continuous by smooth (varying coefficient) true model. The model is $\hat{y}_i = f_2(x_{2i})x_{1i}$ where the function $f_2()$ is $f_2(x) = 0.2 * x^{11} * (10 * (1 - x))^6 + 10 * (10 * x)^3 * (1 - x)^{10}$.

• "factor_by" or "eg4": a factor by smooth true model. The true model contains a factor with 3 levels, where the response for the nth level follows the nth Gu & Wahba function (for $n \in {1, 2, 3}$).

• "additive_factor" or "eg5": an additive plus factor true model. The response is a linear combination of the Gu & Wahba functions 2, 3, 4 (the latter is a null function) plus a factor term with four levels.

• "additive_group_effects" or "eg6": an additive model plus four group intercepts with fixed realised values ⁠3, 6, 9, 12⁠.

• "additive_correlated" or "eg7": a version of the model in "eg1", but where the covariates are correlated.

• "gu_wahba_f2" or "gwf2": a model where the response is Gu & Wahba's $f_2(x_i)$ plus noise.

• "luo_wahba_f6" or "lwf6": a model where the response is Luo & Wahba's "example 6" function $sin(2(4x-2)) + 2 exp(-256(x-0.5)^2)$ plus noise.

• "gfam": simulates data for use with GAMs with family = gfam(families). See example in mgcv::gfam(). If this model is specified then dist is ignored and gfam_families is used to specify which distributions are included in the simulated data. Can be a vector of the distributions allowed by dist, including repeated names. Defaults use Gaussian standard deviation 0.5, Gamma dispersion 0.5, Tweedie dispersion 1.5, negative binomial size 3 and Tweedie power 1.5. Historical grouped ordinal cuts are c(1, 2, 2.5). Explicit scale, theta, power, n_cat and cuts override these settings. A supplied scale is the Gaussian standard deviation and the Gamma or Tweedie dispersion. Identity, logit and log links are used for normal, binary and positive-mean responses, respectively. Ordinal responses use a logistic latent variable.

The random component providing noise or sampling variation can follow one of the distributions, specified via argument dist

• "normal": Gaussian,

• "poisson": Poisson,

• "binary": Bernoulli,

• "negbin": Negative binomial,

• "tweedie": Tweedie,

• "gamma": gamma , and

• "ordered categorical": ordered categorical

Other arguments provide the parameters for the distribution.

Value

A tibble containing responses, covariates and true-function columns. Legacy f is the predictor after its distribution-specific transformation; component columns retain their original function values. With family, f is the first linear predictor, y can be a matrix column, and the "truth" attribute contains eta (all linear predictors), parameters (natural distribution parameters), functions (the original template), response (complete responses before censoring), and family. Grouped families also record index and component-specific parameter lists. Binomial family objects return proportions and a trials column for fitting weights. Cox models return observed times and a status column. An explicit seed restores the caller's RNG state, including its absence; seed = NULL advances the caller's RNG normally.

Shape-constrained examples

Models named "sc_increasing", "sc_decreasing", "sc_convex", "sc_concave", "sc_positive", "sc_increasing_convex", "sc_increasing_concave", "sc_decreasing_convex" and "sc_decreasing_concave" provide univariate reference shapes on ⁠[0, 1]⁠. "sc_additive" combines four terms; "sc_tensor" is increasing in two coordinates; "sc_location_scale" has a concave log-mean and increasing log dispersion. The latter defaults to family = mgcv::gammals(); other shape templates default to Gaussian noise with standard deviation 0.2. See data_sim_shapes for functions, derivatives and mgcv::scasm() recipes.

References

Gu, C., Wahba, G., (1993). Smoothing Spline ANOVA with Component-Wise Bayesian "Confidence Intervals." J. Comput. Graph. Stat. 2, 97–117.

Luo, Z., Wahba, G., (1997). Hybrid adaptive splines. J. Am. Stat. Assoc. 92, 107–116.

See Also

data_sim_families, data_sim_shapes

Examples

data_sim("additive", n = 100, seed = 1)

# an ordered categorical response
data_sim("additive", n = 100, dist = "ocat", n_cat = 4, cuts = c(-1, 0, 5))

Description

data_sim() accepts instantiated standard and mgcv family objects. The example template supplies covariates and true functions; predictors supplies the linear predictors; parameters supplies generating values that are not themselves modelled by a predictor. The "truth" attribute records all predictors, distribution parameters and complete responses.

Scalar families

The parameter names below belong inside parameters = list(...). Parameters are scalar unless explicitly described as observation-level. Beta draws use the constructor's eps endpoint adjustment, matching mgcv's generator. This adjustment can matter when the beta law places substantial probability very close to zero or one; eps is recorded in the truth.

Dispersion and trials can have length one or n. All other scalar-family parameters are constant across observations. Fixed constructor parameters are used when available; conflicting explicit parameters are errors. Otherwise the corresponding generating values are required, including cut points for ocat(R = ...). For example, both family = mgcv::nb(theta = 3) and ⁠family = mgcv::nb(), parameters = list(theta = 3)⁠ specify size 3. Constructor starting values for an estimated parameter are not used as generating truth. A size range in negbin(c(1, 10)) also requires an explicit generating theta.

For ziP(), lambda = exp(eta) and P(Y > 0) = 1 - exp(-exp(intercept + slope * eta)). Positive observations follow a Poisson distribution conditioned on being positive. The natural slope includes the constructor's lower bound b; the second entry of ziP(theta = ...) is instead a log excess above b. For fixed ordinal spacings, mgcv fixes the first cut at -1; for example ocat(theta = c(1, 5)) supplies cuts c(-1, 0, 5). For free cuts whose first value is a, mgcv's fixed first cut of -1 represents the same probabilities after shifting the predictor by -(a + 1). Include the corresponding component intercept in grouped fitting recipes.

Quasi families do not define a complete sampling distribution and are rejected. Nonzero Box-Cox powers are currently rejected: some Gaussian transformed responses would lie outside the inverse transformation's domain.

Families with several linear predictors

The columns returned by predictors follow mgcv's order and links, including their constructor-specific bounds. These are link-scale predictors, not the mu arguments of mgcv's random-number generators.

These families need no scalar generating parameters except mvn()'s optional positive definite covariance matrix. In shash(), the transformed normal is multiplied by tail * sigma, following mgcv's parameterisation. In gumbls(), the mean is location plus Euler's constant times sigma; the first predictor specifies the location. For gevlss(), the default shape link restricts shape to (-1, 0.5); the identity link allows any finite shape. In ziplss(), positive counts follow a zero-truncated Poisson distribution.

Default predictors

Let g be the template's f before distribution transformations. Scalar Gaussian, t, Gaussian-censoring, logistic-censoring and Box-Cox families use natural location g. Binomial and beta means are plogis((g - 5) / 5); ordinal location is g - mean(g); the ziP() predictor is g / 5. Other scalar means are exp(g / 5). Each mean is transformed with the supplied family's link. If it is outside that link's domain, supply custom predictors rather than relying on a default.

For general families, default predictors are: ⁠(g, 0)⁠ for gaulss() and gumbls(); ⁠(g / 5, 0)⁠ for gammals() and ziplss(); ⁠(g, 0, linkfun(0))⁠ for gevlss() (zero shape under the supplied shape link); ⁠(g, 0, 0, 0)⁠ for shash(); and ⁠(g / 5, 0, 0)⁠ for twlss(). Multivariate normal means repeat g. Multinomial logit j is (g - mean(g)) / 5 * j / K. The sc_location_scale template instead supplies its own log mean and log dispersion for gammals().

Censoring and survival

For cnorm(), clog(), cpois() and bcg(), observation can contain finite left and right detection limits and a positive interval_width. Responses outside detection limits are left or right censored. Other responses are placed in rounding intervals wholly inside those limits; values in intervals crossing a limit are retained exactly. Omitting observation returns complete responses.

Matrix responses use equal entries for exact observations, ⁠(limit, -Inf)⁠ for left censoring, ⁠(limit, Inf)⁠ for right censoring, and ⁠(lower, upper)⁠ for intervals. bcg() instead uses zero in the second entry for left censoring; its limits must be positive. Poisson detection limits must be non-integer, and interval widths must be integer numbers of counts. Its intervals use half-integer endpoints to avoid ambiguity at count values.

cox.ph() uses Weibull cumulative baseline hazard H0(t) = (t / baseline_scale)^baseline_shape, with both baseline parameters defaulting to 1. Event times have survival exp(-H0(t) * exp(eta)). Independent exponential censoring has censor_rate = 0.2; zero disables censoring. Default eta is (g - mean(g)) / 5. Fit using weights = status. True event and censoring times are retained in the truth parameter list.

Grouped families

Supply family = mgcv::gfam(list(...)), optionally with parameters = list(components = list(...), index = ...). components contains one parameter list per family-list position. index is one integer position per row (a scalar selects the same family for every row). Without index, every component is represented at least once and remaining positions are sampled uniformly; n must be at least the number of families. Repeated family types retain separate indices and parameters.

The default shared predictor is g / 5; custom predictors return one column. Scalar response families are supported; general families are rejected. Censored-family components return complete, uncensored responses. Fit using cbind(y, index) as the response and weights = trials if a component is binomial with more than one trial. Row-level dispersion and trial vectors in component parameter lists align with the full generated data, or can have the component's row count. Truth parameters record each component's rows.

See Also

data_sim(), data_sim_shapes, mgcv::family.mgcv()

Examples

# A beta response: precision is explicitly chosen, rather than estimated.
d <- data_sim(family = mgcv::betar(theta = 12), seed = 2)
attr(d, "truth")$parameters$precision

# Binomial proportions: use trials as the fitting weights.
d <- data_sim(family = binomial(), parameters = list(trials = 10), seed = 2)

# Mean and standard deviation have separate predictors.
d <- data_sim("gwf2", family = mgcv::gaulss(), seed = 2,
  predictors = function(d) cbind(d$f, log(0.2 + d$x)))

# Fit a censored normal model to the matrix response.
d <- data_sim("sc_increasing", family = mgcv::cnorm(theta = 0.2), seed = 2,
  observation = list(left = 0.1, right = 0.9))

m <- mgcv::gam(y ~ s(x), family = mgcv::cnorm(), data = d, method = "REML")

Description

data_sim() supplies analytic reference functions for mgcv::scasm(). All covariates are uniform on ⁠[0, 1]⁠. Univariate templates return f1, the true function, and its analytic derivatives d1 and d2. The response contains sampling variation; the derivative columns describe the truth.

Univariate templates

Remove the intercept for the positivity recipe so the constrained smooth represents the whole function. An unconstrained intercept combined with a centred smooth does not express the same positivity requirement.

Several terms and distribution parameters

sc_additive adapts mgcv's example. Its raw terms are gw_f0(x0), 2 * plogis(20 * (x1 - 0.4)), gw_f2(x2) and 10 * (x3 - 0.55)^4. The true predictor is their sum minus 1, divided by 4. The second term is increasing and the fourth is convex. Use predictors = function(d) d$f with family = poisson() to sample mgcv's Poisson example.

sc_tensor uses f = x + z + x * z, increasing in both covariates. Its dx and dz columns give the true partial derivatives.

sc_location_scale uses log mean sin(pi * x) / 2 and log dispersion -2 + z, returned in f and f_scale. It defaults to gammals() and records derivatives in d1, d2 and d_scale. Constraints act on the specified predictors; in particular the concavity concerns the log mean.

The recipes use the default bs = 0 for fitting without bootstrap replicates. To obtain native constraint-respecting bootstrap coefficient draws, fit with a positive bs, such as bs = 100.

See Also

data_sim(), data_sim_families, mgcv::scasm()

Examples

load_mgcv()
d <- data_sim("sc_increasing", n = 200, seed = 2)
m <- scasm(y ~ s(x, bs = "sc", xt = "m+", k = 8), data = d)

d <- data_sim("sc_positive", n = 200, seed = 2)
m <- scasm(y ~ s(x, bs = "sc", xt = "+", k = 8) - 1, data = d)


d <- data_sim("sc_additive", n = 240, family = poisson(), seed = 2,
  predictors = function(d) d$f)
m <- scasm(y ~ s(x0, k = 6) + s(x1, bs = "sc", xt = "m+", k = 6) +
  s(x2, k = 6) + s(x3, bs = "sc", xt = "c+", k = 6),
  family = poisson(), data = d)

d <- data_sim("sc_tensor", seed = 2)
m <- scasm(y ~ te(x, z, bs = c("sc", "sc"), xt = list("m+", "m+"),
  k = c(6, 6)), data = d)

d <- data_sim("sc_location_scale", seed = 2)
m <- scasm(list(y ~ s(x, bs = "sc", xt = "c-"),
  ~s(z, bs = "sc", xt = "m+")), family = gammals(), data = d)

Description

Prepare a data slice through model covariates

Usage

data_slice(object, ...)

## Default S3 method:
data_slice(object, ...)

## S3 method for class 'data.frame'
data_slice(object, ..., .observed_only = FALSE, .by = NULL)

## S3 method for class 'gam'
data_slice(
  object,
  ...,
  data = NULL,
  envir = NULL,
  .observed_only = FALSE,
  .by = NULL
)

## S3 method for class 'gamm'
data_slice(object, ...)

## S3 method for class 'gamm4'
data_slice(object, ...)

## S3 method for class 'scam'
data_slice(object, ...)

Arguments

Details

A data slice is the data set that results where one (or more covariates) is varied systematically over some or all of its (their) range or at a specified subset of values of interest, while any remaining covariates in the model are held at fixed, representative values. This is known as a reference grid in package emmeans and a data grid in the marginaleffects package.

For GAMs, any covariates not specified via ... will take representative values determined from the data used to fit the model as follows:

• for numeric covariates, the value in the fitting data that is closest to the median value is used,

• for factor covariates, the modal (most frequently observed) level is used, or the first level (sorted as per the vector returned by base::levels() if several levels are observed the same number of times.

These values are already computed when calling gam() or bam() for example and can be found in the var.summary component of the fitted model. Function typical_values() will extract these values for you if you are interested.

Convenience functions evenly(), ref_level(), and level() are provided to help users specify data slices. ref_level(), and level() also ensure that factor covariates have the correct levels, as needed by mgcv::predict.gam() for example.

For an extended discussion of data_slice() and further examples, see vignette("data-slices", package = "gratia").

With .by, unspecified covariates retain their overall representative values. Specify an expression such as z = mean(z) to use a group-specific value instead. Factor levels and ordered status are preserved, and the result is an ungrouped tibble. .observed_only filters exact matches within each group; interpolated values need not match observed values.

Grouping limits the default range of evenly() to each group's marginal range. It does not ensure joint support for multiple continuous covariates. Explicit values or lower and upper bounds can still request extrapolation.

See Also

The convenience functions evenly(), ref_level(), and level(). typical_values() for extracting the representative values used for covariates in the model but not named in the slice.

Examples

load_mgcv()

# simulate some Gaussian data
df <- data_sim("eg1", n = 50, seed = 2)

# fit a GAM with 1 smooth and 1 linear term
m <- gam(y ~ s(x2, k = 7) + x1, data = df, method = "REML")

# Want to predict over f(x2) while holding `x1` at some value.
# Default will use the observation closest to the median for unspecified
# variables.
ds <- data_slice(m, x2 = evenly(x2, n = 50))
ds

# for full control, specify the values you want
ds <- data_slice(m, x2 = evenly(x2, n = 50), x1 = 0.3)

# or provide an expression (function call) which will be evaluated in the
# data frame passed to `data` or `model.frame(object)`
ds <- data_slice(m, x2 = evenly(x2, n = 50), x1 = mean(x1))

Description

Posterior expectations of derivatives from an estimated model

Usage

## S3 method for class 'gamm4'
derivative_samples(model, ...)

derivative_samples(model, ...)

## Default S3 method:
derivative_samples(model, ...)

## S3 method for class 'gamm'
derivative_samples(model, ...)

## S3 method for class 'gam'
derivative_samples(
  model,
  ...,
  data = NULL,
  n_draws = NULL,
  chains = NULL,
  seed = NULL,
  method = c("gaussian", "mh", "inla", "user"),
  n_cores = 1L,
  n_threads = 1L,
  burnin = 1000,
  thin = 1L,
  t_df = 40,
  rw_scale = 0.25,
  frequentist = FALSE,
  unconditional = FALSE,
  draws = NULL,
  diagnostics = FALSE,
  mh_backend = c("mgcv", "gratia", "family"),
  mvn_method = c("mvnfast", "mgcv"),
  block_size = 10000L,
  dense_limit = 5e+06,
  focal = NULL,
  order = 1L,
  type = c("forward", "backward", "central"),
  scale = c("response", "linear_predictor"),
  n_grid = 200L,
  eps = NULL,
  envir = environment(formula(model))
)

## S3 method for class 'scam'
derivative_samples(
  model,
  ...,
  data = NULL,
  n_draws = NULL,
  chains = NULL,
  seed = NULL,
  method = c("gaussian", "mh", "inla", "user"),
  n_cores = 1L,
  n_threads = 1L,
  burnin = 1000,
  thin = 1L,
  t_df = 40,
  rw_scale = 0.25,
  frequentist = FALSE,
  unconditional = FALSE,
  draws = NULL,
  diagnostics = FALSE,
  mh_backend = c("mgcv", "gratia", "family"),
  mvn_method = c("mvnfast", "mgcv"),
  block_size = 10000L,
  dense_limit = 5e+06,
  focal = NULL,
  order = 1L,
  type = c("forward", "backward", "central"),
  scale = c("response", "linear_predictor"),
  n_grid = 200L,
  eps = NULL,
  envir = environment(formula(model))
)

Arguments

Details

For LSS models, response-scale derivatives differentiate the expected response returned by fitted_samples(), retaining the joint uncertainty in every parameter that affects that expectation. Predictor-scale derivatives refer to the first linear predictor. Smooth-component derivatives and derivatives of the response expectation are different quantities.

Rows with unavailable required covariates retain NA derivatives by default; explicit na.action requests follow response_derivatives(). Non-finite derivative calculations at available rows raise an error, including when only some coefficient draws are affected. Response-scale derivatives require an existing response mean; scaled t families need fitted degrees of freedom greater than one.

Value

A tibble, currently with the following variables:

• .derivative: the estimated partial derivative,

• additional columns containing the covariate values at which the derivative was evaluated.

Author(s)

Gavin L. Simpson

See Also

gamm4_support

Examples

load_mgcv()
df <- data_sim("eg1", dist = "negbin", scale = 0.25, seed = 42)

# fit the GAM (note: for execution time reasons using bam())
m <- bam(y ~ s(x0) + s(x1) + s(x2) + s(x3),
  data = df, family = nb(), method = "fREML")

# data slice through data along x2 - all other covariates will be set to
# typical values (value closest to median)
ds <- data_slice(m, x2 = evenly(x2, n = 200))

# samples from posterior of derivatives
fd_samp <- derivative_samples(m,
  data = ds, type = "central",
  focal = "x2", eps = 0.01, seed = 21, n_draws = 100
)

# plot the first 20 posterior draws
if (requireNamespace("ggplot2") && requireNamespace("dplyr")) {
  library("ggplot2")
  fd_samp |>
    dplyr::filter(.draw <= 20) |>
    ggplot(aes(x = x2, y = .derivative, group = .draw)) +
    geom_line(alpha = 0.5)
}