Package {gdam}

Title:	Fast Robust Additive Models using Gamma Divergence
Version:	0.0.1
Date:	2026-09-10
Description:	Provides a computationally efficient method (at least at the time of development!) for robust additive modeling based on the gamma divergence, implementing the methodology of Hui, Ding, and Sugasawa (2026) <doi:10.1007/s11222-026-10960-3>. The method assumes an identity link and normally distributed errors, and applies the gamma divergence to the resulting restricted maximum likelihood function so as to obtain a loss function that is less sensitive to outlying responses. Because 'gdam' uses a model fitted via the 'mgcv' package as the starting point for optimization, and leverages many of its existing techniques for estimation and inference, users can take advantage of many of the smoothing options available in the 'mgcv' package for constructing additive models.
Maintainer:	Francis KC Hui <fhui28@gmail.com>
License:	GPL (≥ 3)
Encoding:	UTF-8
Imports:	colorspace, ggplot2, gratia, Matrix, mgcv, stats
Suggests:	spelling, tidyverse
Language:	en-US
URL:	https://github.com/fhui28/gdam
BugReports:	https://github.com/fhui28/gdam/issues
Config/roxygen2/version:	8.0.0
RoxygenNote:	7.3.3
NeedsCompilation:	no
Packaged:	2026-09-14 09:46:04 UTC; fkch
Author:	Francis KC Hui [aut, cre], Shonosuke Sugasawa [ctb]
Repository:	CRAN
Date/Publication:	2026-09-24 13:30:20 UTC

Add partial residuals from a fitted gdam model

Description

This function is basically a wrapper around the gratia::add_partial_residuals() function for adding partial residuals, but designed to work with gdam objects. Please see the associated help file in the gratia package for more details.

Usage

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

Arguments

Value

A data frame of the same dimension as data with partial residuals added on, while the final column also contains the final weights from the MM algorithm used to fit the gdam model.

Note

Acknowledgments to Gavin Simpson and the gratia package for the original gratia::add_partial_residuals() function.

Author(s)

Francis K.C. Hui fhui28@gmail.com

Model diagnostics plots for gdam objects

Description

This function is basically a wrapper around the gratia::appraise() function, but designed to work with gdam objects. The main difference is that observations (residuals, quantiles etc...) in the diagnostic plots are colored by the final weights from the gdam object.

Usage

## S3 method for class 'gdam'
appraise(
  model,
  use_worm = FALSE,
  n_uniform = 10,
  n_simulate = 50,
  method = "uniform",
  type = "deviance",
  n_bins = "sturges",
  ncol = NULL,
  nrow = NULL,
  guides = "keep",
  level = 0.9,
  ci_col = "black",
  ci_alpha = 0.2,
  line_col = "red",
  ...
)

Arguments

Value

A set of model diagnostic plots as per the original gratia::appraise() function.

Note

Acknowledgments to Gavin Simpson and the gratia package for the original gratia::appraise() function.

Author(s)

Francis K.C. Hui fhui28@gmail.com

Examples

# See the main `gdam` help file for examples.

Fast robust additive models using gamma divergence.

Description

This function fits a robust additive model, assuming an identity link function and a normal distribution for the errors, using a gamma divergence adaption of the restricted maximum likelihood function (REML) as the loss function (hence the name gdam, which stands for gamma divergence additive models).

The functions works by taking an object of class mgcv::gam() and uses at as a starting point for optimization of the gamma divergence, which itself is performed via a minorization-maximization (MM) algorithm.

Usage

gdam(
  gamObject,
  gamma_tuning = 0.3,
  control = list(tol = 1e-06, max_iteration = 1000, rinse_and_repeat = 0)
)

Arguments

Details

See manuscript for more details (currently in preparation).

Value

An object of class gdam with the following components:

• call The matched call.

• coefficients The estimated (smoothing) coefficients of the additive model.

• sp The estimated smoothing parameters of the additive model.

• sigma2 The estimated residual variance of the additive model.

• linear_predictors The linear predictors of the additive model.

• fitted_values The fitted values of the additive model.

• residuals The raw residuals of the additive model.

• final_weights The final weights from the MM algorithm. The weights sum to the total number of observations in the data. Weights comparably less than 1 indicate that the observation is down-weighted in the estimation process i.e., is considered more likely to be an outlying observation

• gamma_tuning The gamma tuning parameter used in the estimation process.

• gamma_divergence The value of the gamma divergence at convergence.

• gamObject The original gamObject used to fit the additive model.

• covariance_matrix The estimated covariance matrix of the parameters, including of the residual variance (but not the smoothing parameters).

• Hscore The conditional Hyvarinen score (H-score) of the additive model. This is a measure that can be used to select the tuning parameter gamma_tuning, when comparing across difference gdam fits. The H-score calculation is based on the proposal of Sugusawa and Yonekura (2021) doi:10.3390/e23091147.

• aic The conditional Akaike Information Criterion (AIC) of the additive model. This is calculated as -2 \times \text{gamma\_divergence} + 2 \times \text{edf}, where edf is the effective degrees of freedom of the additive model. Effectively adapts the idea of Kurara (2024) doi:10.1080/03610926.2022.2155788; see also mgcv::logLik.gam().

Note

Warning

Not all smooth specifications available in mgcv are supported let alone being properly tested e.g., we don't think smooths linked using the "id" argument, linear functional terms, or varying coefficient models as available using the "by" argument are supported; see mgcv::gam.models() for specification of GAMs. Apologies in advance if you encounter any issues with these, but do email the author if encounter any issues or have a feature request.

Author(s)

Francis K.C. Hui fhui28@gmail.com

Examples

library(mgcv)
library(tidyverse)
library(gratia)

# Simulate a additive model and add some outliers.
set.seed(062025)
simdat <- gamSim(eg = 7,
n = 200,
dist = "normal",
scale = 0.5,
verbose = TRUE)

create_outliers <- rbinom(nrow(simdat), size = 1, prob = 0.05)
# Could also the probability homogeneous e.g., vary with one or more covariate values.
table(create_outliers)
simdat$contaminated_y <- simdat$y
simdat$contaminated_y[create_outliers == 1] <- rnorm(sum(create_outliers == 1),
mean = simdat$f[create_outliers==1],
sd = 9)


# Fit the model using mgcv and then fit the robust additive model using gdam
fit_mgcv <- gam(contaminated_y ~ s(x0) + s(x1) + s(x2) + s(x3),
data = simdat,
family = gaussian(),
method = "REML")

fit_gdam <- gdam(gamObject = fit_mgcv, gamma_tuning = 0.1)
fit_gdam

summary(fit_gdam)

# Basic illustration of functions
residuals(object = fit_gdam)

appraise(fit_gdam)

## Automated partial smooth plots
gdam::smooth_estimates(fit_gdam) %>%
gratia::draw() & theme_bw()

## A bit more manual but customized partial smooth plots -- all covariates
gdam::smooth_estimates(fit_gdam) %>%
gratia::add_confint() %>%
pivot_longer(x0:x3, values_to = "x", names_to = "covariate") %>%
filter(!is.na(x)) %>%
ggplot() +
geom_line(aes(x = x, y = .estimate), lwd = 1.2) +
geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci, x = x), alpha = 0.2) +
geom_rug(aes(x = x),
data = simdat %>% pivot_longer(x0:x3, values_to = "x", names_to = "covariate"),
sides = "b") +
facet_wrap(. ~ covariate, scale = "free", nrow = 2) +
theme_bw()

## A bit more manual but customized partial smooth plots -- one covariate
## but adding partial residuals
simdat_expanded <- add_partial_residuals(data = simdat, model = fit_gdam, select = NULL)
gdam::smooth_estimates(fit_gdam) %>%
gratia::add_confint() %>%
filter(.smooth == s(x0)) %>%
ggplot() +
geom_point(aes(x = x0, y = `s(x0)`, color = weights), data = simdat_expanded, cex = 1.5) +
geom_line(aes(x = x0, y = .estimate), lwd = 1.2) +
geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci, x = x0), alpha = 0.2) +
scale_color_viridis_c() +
theme_bw()


predict(fit_gdam, se.fit = TRUE)

# Simple check of the in-sample performance.
# Please see the associated manuscript for a more systematic evaluation of
# gdam including against other robust methods for additive models
cbind(simdat$f, fit_mgcv$fitted.values, fit_gdam$fitted_values) %>%
t %>%
dist()




# Checking various selection criteria -- mainly for choosing the gamma tuning parameter
# Personally we prefer using the H-score for this as it is slightly more
# principled in its construction for gamma divergence
fit_gdam00 <- gdam(gamObject = fit_mgcv, gamma_tuning = 1e-6)
fit_gdam001 <- gdam(gamObject = fit_mgcv, gamma_tuning = 0.01)
fit_gdam01 <- gdam(gamObject = fit_mgcv, gamma_tuning = 0.1)
fit_gdam02 <- gdam(gamObject = fit_mgcv, gamma_tuning = 0.2)
fit_gdam03 <- gdam(gamObject = fit_mgcv, gamma_tuning = 0.3)
fit_gdam04 <- gdam(gamObject = fit_mgcv, gamma_tuning = 0.4)
fit_gdam05 <- gdam(gamObject = fit_mgcv, gamma_tuning = 0.5)

fit_gdam00$Hscore
fit_gdam001$Hscore
fit_gdam01$Hscore
fit_gdam02$Hscore
fit_gdam03$Hscore
fit_gdam04$Hscore
fit_gdam05$Hscore

Predictions from a fitted gdam model.

Description

Takes a fitted gam object produced by gdam() and produces predictions given a new set of values for the model covariates or the original values used for the model fit. Predictions can be accompanied by standard errors and uncertainty intervals based on a large sample normality assumption for the distribution of the model coefficients. The routine can optionally return the matrix by which the model coefficients must be premultiplied in order to yield the values of the linear predictor.

Usage

## S3 method for class 'gdam'
predict(
  object,
  newdata,
  type = "response",
  se.fit = FALSE,
  coverage = 0.95,
  na.action = na.pass,
  ...
)

Arguments

Value

A vector of predictions, or a data frame containing the predictions, standard errors, and uncertainty intervals.

Author(s)

Francis K.C. Hui fhui28@gmail.com

Examples

# See the main `gdam` help file for examples.

Print a fitted gdam model.

Description

This function prints out some details of the fitted gdam object in a human-readable format, similar to that of mgcv::print.gam(). This includes the tuning parameter used, the formula, the estimated degrees of freedom, and the gamma divergence. Note the estimated degrees of freedom should be taken with a grain of salt, as they are estimated from a post-hoc GAM fit (we conjecture that they are likely to be marginally underestimated).

Usage

## S3 method for class 'gdam'
print(x, ...)

Arguments

Value

Invisible output object.

Author(s)

Francis K.C. Hui fhui28@gmail.com

Examples

# See the main `gdam` help file for examples.

Residuals from a fitted gdam model.

Description

Computes residuals from a fitted gdam object. The residuals can be of different types: "response" corresponds to the raw residuals, "deviance" corresponds to the deviance residuals (note these are slightly ad-hoc as the concept of deviance itself is not clearly defined in the case of gamma divergence; here the deviance residual is computed assuming normally distributed errors), "pearson" corresponds to the Pearson residuals, "scaled_pearson" corresponds to the Pearson residuals scaled by the square root of the estimated dispersion parameter, and "working" corresponds to the working residuals (i.e., the residuals from a post-hoc GAM fit, though in this case they should correspond to the raw response residuals).

Usage

## S3 method for class 'gdam'
residuals(object, type = "deviance", ...)

Arguments

Value

A vector of residuals.

Author(s)

Francis K.C. Hui fhui28@gmail.com

Examples

# See the main `gdam` help file for examples.

Evaluate smooths from a fitted gdam model at covariate values

Description

This function is basically a wrapper around the gratia::smooth_estimates() function for valuating a smooth at a grid of evenly spaced value over the range of the covariate, but designed to work with gdam objects. Please see the associated help file in the gratia package for more details.

Usage

smooth_estimates(object, ...)

## S3 method for class 'gdam'
smooth_estimates(
  object,
  select = NULL,
  n = 100,
  n_3d = 16,
  n_4d = 4,
  data = NULL,
  dist = NULL,
  partial_match = FALSE,
  ...
)

Arguments

Value

A data frame which is of class smooth_estimates.

Note

Acknowledgments to Gavin Simpson and the gratia package for the original gratia::smooth_estimates() function.

Author(s)

Francis K.C. Hui fhui28@gmail.com

Examples

# See the main `gdam` help file for examples, particularly for constructing partial smooth plots.

Summary for a fitted gdam model.

Description

Takes a fitted gdam object and produces some useful post-hoc summaries of the fitted model. Note all results should be taken with a grain of salt, as they are estimated from a post-hoc GAM fit using mgcv::gam() (we conjecture quantities like degrees of freedom are likely to be marginally underestimated, and hence P-values are likely to be marginally smaller than they should be). Please see mgcv::summary.gam() for details on the output.

Usage

## S3 method for class 'gdam'
summary(
  object,
  digits = max(3, getOption("digits") - 3),
  signif.stars = getOption("show.signif.stars"),
  ...
)

Arguments

Value

A summary of the fitted gdam object, similar to that of mgcv::summary.gam(). Please see the help file for that function for more details.

Note

Warning

The p-values are approximate and neglect smoothing parameter uncertainty. They are likely to be somewhat too low when smoothing parameter estimates are highly uncertain; please see the mgcv::summary.gam() help file for more details.

Author(s)

Francis K.C. Hui fhui28@gmail.com

Examples

# See the main `gdam` help file for examples.