Special Model Terms

Special model terms extend a gamlss2 formula with a model component that is fitted inside the backfitting algorithm. Ordinary linear terms and mgcv smooths cover many models; specials are useful when the fitting method needs its own estimator, prediction method, state, or tuning parameters.

The standard fitting workflow is described in First Steps. This vignette focuses on extending the predictor; Family Objects instead extends the response distribution.

1 Built-in special terms

The package already provides several specials. For example, lo() fits a weighted local-polynomial smoother, tree(), ct(), and cf() fit tree-based terms, n() fits a neural-network term, and lin(), re(), random(), and gnet() provide structured linear or grouped effects. The mgcv terms s(), te(), ti(), and t2() are also supported.

air <- subset(airquality, !is.na(Ozone))

m <- gamlss2(
  Ozone ~ lo(~ Temp, span = 0.7) + s(Wind) | s(Temp),
  data = air,
  family = GA,
  trace = FALSE
)

summary(m)

Call:
gamlss2(formula = Ozone ~ lo(~Temp, span = 0.7) + s(Wind) | s(Temp), 
    data = air, family = GA, trace = FALSE)
---
Family: GA(log(mu), log(sigma)) 
*--------
Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
mu.(Intercept)     3.54163    0.02714  130.48   <2e-16 ***
sigma.(Intercept) -0.90192    0.06343  -14.22   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
                         edf
mu.s(Wind)            3.3464
mu.lo(~Temp,span=0.7) 5.7708
sigma.s(Temp)         4.3414
*--------
n = 116 df =  15.46 res.df =  100.54
Deviance = 921.4427 Null Dev. Red. = 14.92%
AIC = 952.3597 elapsed =  0.74sec

Specials can occur in any distributional parameter formula. They are refitted from the current working response and weights at each backfitting iteration, then centered before being added to the additive predictor. Centering keeps the intercept and the special term identifiable.

2 Defining a custom special

A custom special consists of three pieces:

1. a constructor used in the formula;
2. a special_fit() method that fits the term to the current working response and weights; and
3. a special_predict() method for fitted values on new data.

The constructor should return a list of class c("special", "user") containing the formula, variables, data, controls, and anything else needed by fitting and prediction. The reserved user name is recognized by fake_formula(), so it can be used without changing package internals. For a package-level special, register its name in fake_formula() and use a dedicated class instead.

A fitting method receives at least x, z, and w. The framework can also supply y, eta, j, family, and control; accepting ... keeps a method compatible with these additional arguments. The returned list must contain numeric fitted.values. It should usually also contain:

• edf, the effective degrees of freedom;
• model, the fitted underlying model;
• shift, the centering constant; and
• optionally transfer, for state passed to the next backfitting iteration.

A prediction method receives the fitted object and data, and returns a numeric vector. If it supports se.fit = TRUE, it should return a data frame with a fit column and, optionally, interval columns.

3 A custom weighted local smoother

The following example implements a small user() term around stats::loess(). It is deliberately compact: a production term should validate its controls and handle missing values and extrapolation explicitly.

user <- function(formula, span = 0.75, degree = 2) {
  if(!inherits(formula, "formula"))
    formula <- as.formula(paste("~", deparse1(substitute(formula))))
  st <- list(
    formula = formula,
    data = model.frame(formula),
    control = list(span = span, degree = degree),
    term = all.vars(formula),
    label = paste0("user(", paste(all.vars(formula), collapse = "+"), ")")
  )
  class(st) <- c("special", "user")
  st
}

special_fit.user <- function(x, z, w, control, ...) {
  fit_formula <- update(x$formula, response_z ~ .)
  dat <- x$data
  dat$response_z <- z
  dat$weights_w <- w
  fit <- loess(
    fit_formula,
    data = dat,
    weights = weights_w,
    span = x$control$span,
    degree = x$control$degree
  )
  fitted <- as.numeric(predict(fit, newdata = dat))
  shift <- mean(fitted)
  out <- list(
    model = fit,
    fitted.values = fitted - shift,
    shift = shift,
    edf = fit$trace.hat
  )
  class(out) <- "user.fitted"
  out
}

special_predict.user.fitted <- function(x, data, se.fit = FALSE, ...) {
  fit <- as.numeric(predict(x$model, newdata = data)) - x$shift
  if(se.fit) return(data.frame(fit = fit))
  fit
}

The user constructor is available while the formula is parsed because it is one of the reserved special names. The corresponding methods are found by S3 dispatch during fitting and prediction.

m <- gamlss2(
  Ozone ~ user(~ Temp, span = 0.7) | s(Temp),
  data = air,
  family = GA,
  trace = FALSE
)

summary(m)

Call:
gamlss2(formula = Ozone ~ user(~Temp, span = 0.7) | s(Temp), 
    data = air, family = GA, trace = FALSE)
---
Family: GA(log(mu), log(sigma)) 
*--------
Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
mu.(Intercept)     3.55870    0.02891  123.08   <2e-16 ***
sigma.(Intercept) -0.82894    0.06315  -13.13   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
                           edf
mu.user(~Temp,span=0.7) 5.8084
sigma.s(Temp)           4.3942
*--------
n = 116 df =  12.2 res.df =  103.8
Deviance = 939.5092 Null Dev. Red. = 13.26%
AIC = 963.9144 elapsed =  0.04sec

## predicted distributions
predict(m, newdata = air[1:5, , drop = FALSE])

                                       1 
"GAMLSS2 GA(mu = 18.91, sigma = 0.4630)" 
                                       2 
"GAMLSS2 GA(mu = 19.02, sigma = 0.4831)" 
                                       3 
"GAMLSS2 GA(mu = 19.89, sigma = 0.5001)" 
                                       4 
"GAMLSS2 GA(mu = 16.22, sigma = 0.5930)" 
                                       6 
"GAMLSS2 GA(mu = 18.83, sigma = 0.4703)" 

4 Practical guidance

Use an existing special when it already provides the required estimator and prediction behavior. Write a custom special when the term has a distinct fitting algorithm, needs custom state, or wraps an external modelling method. Keep the constructor lightweight, store all information needed for prediction, return centered fitted contributions, and test both fitting and newdata predictions. Prediction support alone does not imply support for joint coefficient uncertainty: the Gaussian intervals in Prediction and Uncertainty require a coefficient-linear design and a verifiable quadratic penalty. A custom fitting method may need its own uncertainty procedure. The help topic ?special_terms documents the method arguments and return-value contract; ?fake_formula documents special-name registration.