Package 'gamlss2'

Description

Next generation infrastructure for generalized additive models for location, scale, and shape (GAMLSS) and distributional regression more generally. The package provides a fresh reimplementation of the classic 'gamlss' package while being more modular and facilitating the creation of advanced terms and models.

Details

gamlss2 fits an additive predictor for each parameter of a response distribution, using linear terms, smooth terms from mgcv, and user-defined model terms (see specials). Response distributions are specified through gamlss2.family objects, including families with more than four parameters. Fitting algorithms can be selected through gamlss2_control. Fitted models provide methods for summaries, prediction, plotting, and diagnostics. Predictions can be returned as distributions3 objects for evaluating fitted distributions (see predict.gamlss2).

Index of help topics:

ALD                     Asymmetric Laplace Quantile Regression
bamlss2                 Bayesian GAMLSS Wrapper
BS                      Bayesian Sampler for GAMLSS Models
calibration             Calibration Plots for Binary Responses
cf                      Conditional Inference Forest Model Terms
ct                      Conditional Inference Tree Model Terms
cv_gamlss2              Cross Validation for 'gamlss2' Models
deviance.gamlss2        Deviance of GAMLSS Models
discretize_family       Discretize a Distribution Family into Fixed
                        Bins
elm                     Extreme Learning Machine Model Terms
fake_formula            Extended Processing of "Fake" Formulas
family_extensions       Extend Distribution Families
family.distribution     Generate a 'gamlss2' Family from a
                        'distributions3' Object
find_family             Find and Fit GAMLSS Families
fk                      Free-Knot Regression Splines for GAMLSS
gamlss2                 Generalized Additive Models for Location Scale
                        and Shape
GAMLSS2                 Create a GAMLSS2 Distribution
gamlss2_control         Control Parameters
gamlss2_smoothing       Smoothing Parameter Selection
gamlss2_start           Model Initialization and Parameter Constraints
gamlss2-package         GAMLSS Modeling with Advanced Flexible
                        Infrastructures
gamlss2.family          Family Objects in 'gamlss2'
gnet                    Model Terms with 'glmnet'
HarzTraffic             Traffic Counts at Sonnenberg in the Harz Region
im                      Convolutional Image Model Terms for GAMLSS
JR                      Joint REML Estimation
Kumaraswamy             Kumaraswamy Distribution
la                      Lasso Model Terms
lin                     Bundled Linear and Ridge-Penalized Model Terms
lo                      LOESS Smoothers for GAMLSS
logLik.gamlss2          Log-Likelihood of GAMLSS Models
make.link2              Create a Link for Families
marginal_predict        Marginal Prediction Profiles for gamlss2 Models
mcmc                    Run MCMC Sampling for a Fitted 'gamlss2' Model
mixture_family          Construct a Finite Mixture Distribution Family
MN                      Multinomial Logit Family
ms                      Monotonic P-Splines for GAMLSS
n                       Neural Network Model Terms for GAMLSS
new_formula             Extracting a New Formula After Selection
                        Algorithms
OL                      Ordered Logistic Family for Ordinal Regression
pb                      P-Splines for GAMLSS
plot.gamlss2            Plotting GAMLSS
predict.gamlss2         Predictions from gamlss2 Models
prodist.gamlss2         Extracting Fitted or Predicted Probability
                        Distributions from gamlss2 Models
quantile.gamlss2        Quantiles for GAMLSS
re                      Random Effects
residuals.gamlss2       Residuals for GAMLSS Models
response_name           Auxiliary Functions for Formulas and Model
                        Objects
results                 Extract Results from GAMLSS Models
RS                      Rigby and Stasinopoulos (RS) & Cole and Green
                        (CG) Algorithm
Rsq                     GAIC and Generalised (Pseudo) R-squared for
                        GAMLSS Models
select_gamlss2          Smooth Model Term Selection with Additional
                        Penalties
si                      Smooth Image Model Terms for GAMLSS
softplus                Softplus Link Object
special_terms           Special Model Terms for GAMLSS
SpirometryUS            Spirometry Measurements from NHANES 2007-2012
stepwise                Stepwise Model Term Selection Using GAIC
tree                    Regression Tree Model Terms
user                    Reserved User Special Model Term
vcov.gamlss2            Joint Covariance Matrix for gamlss2 Models

Author(s)

Nikolaus Umlauf, Mikis Stasinopoulos, Robert Rigby, Achim Zeileis, Reto Stauffer, Gillian Heller, Fernanda de Bastiani, Andreas Mayr.

Maintainer: Nikolaus Umlauf [email protected].

References

To cite gamlss2 in publications use:

Umlauf N, Stasinopoulos M, De Bastiani F, Heller G, Kneib T, Mayr A, Rigby R, Stauffer R, Zeileis A (2025). “gamlss2: The Next Generation of Distributional Regression”. In O'Sullivan D, Sweeney J, Moghaddam S, Burke K (eds.), Proceedings of the 39th International Workshop on Statistical Modelling, 312–317. Department of Mathematics and Statistics, University of Limerick, Limerick, Ireland. https://iwsm2025.ie/wp-content/uploads/2025/07/IWSM2025_Limerick_Proceedings.pdf

For the underlying methods cite:

Rigby RA, Stasinopoulos DM (2005). “Generalized Additive Models for Location, Scale and Shape (with Discussion)”. Journal of the Royal Statistical Society, Series C (Applied Statistics), 54(3), 507–554. doi:10.1111/j.1467-9876.2005.00510.x

See Also

gamlss2, fake_formula

Description

Asymmetric Laplace family for quantile regression and a wrapper for estimating one or more conditional quantile models.

Usage

ALD(tau = 0.5)

fit_ALD(formula, data, tau, anchor = 0.5, max_step = 0.05,
  restarts = c("always", "on_failure", "never"), start = NULL,
  monotone = FALSE, verbose = interactive(), ...)

Arguments

Details

The asymmetric Laplace density is proportional to $\exp\{-\rho_\tau((y-\mu)/\sigma)\}$, where $\rho_\tau(u) = u\{\tau-I(u<0)\}$. Thus, $\mu$ is the conditional $\tau$-quantile and $\sigma>0$ is a scale parameter.

fit_ALD first estimates the anchor quantile and then proceeds toward the lower and upper requested quantiles. Coefficients and smoothing parameters are transferred between successive fits. Converged candidate fits are preferred and compared using their asymmetric Laplace log-likelihood. If no candidate converges, up to two retries halve and quarter the step length, retaining the smoothing criterion, iteration limits, and convergence tolerances. These retries use cold starts unless restarts = "never". If no candidate converges after these retries, a finite candidate is returned with its convergence flag unchanged; callers should check the individual models or the diagnostics.

Monotone rearrangement affects only the values returned by fitted() and predict(); the individual fitted models are retained unchanged.

The ALD can be used as a working likelihood when the response distribution is not asymmetric Laplace. With constant scale, its data likelihood targets pinball loss for the requested quantile. The location score is nonsmooth, and working-likelihood intervals are not automatically calibrated under misspecification. Ordinary observed-information Gaussian inference for smooth likelihoods should not be interpreted as calibrated quantile confidence intervals here; a suitable calibration or resampling procedure requires separate justification.

Value

ALD returns an object of class "gamlss2.family".

fit_ALD returns an object of class "ALDfit", containing the individual models, fitted values, requested quantiles, and fitting diagnostics.

References

Yu K, Moyeed RA (2001). “Bayesian Quantile Regression”. Statistics & Probability Letters, 54(4), 437–447. doi:10.1016/S0167-7152(01)00124-9

Chernozhukov V, Fernandez-Val I, Galichon A (2010). “Quantile and Probability Curves Without Crossing”. Econometrica, 78(3), 1093–1125. doi:10.3982/ECTA7880

Sriram K (2015). “A Sandwich Likelihood Correction for Bayesian Quantile Regression Based on the Misspecified Asymmetric Laplace Density”. Statistics & Probability Letters, 107, 18–26. doi:10.1016/j.spl.2015.07.035

See Also

gamlss2, predict, mcmc.

Examples

## Not run: ## load data
data("SpirometryUS", package = "gamlss2")

## subset for females
d <- subset(SpirometryUS, gender == "female")

## fit one conditional median model
m50 <- gamlss2(fev1 ~ s(age), data = d, family = ALD(tau = 0.5))

## estimate several conditional quantiles
qu <- c(0.025, 0.5, 0.975)
m <- fit_ALD(fev1 ~ s(age), data = d, tau = qu,
  max_step = 0.1, monotone = TRUE, verbose = FALSE)

## plot fitted quantile curves
fit <- predict(m)
i <- order(d$age)
plot(fev1 ~ age, data = d)
matlines(d$age[i], fit[i, ],
  lty = 1, lwd = 2, col = 4)

## estimate lower limits model using full MCMC
set.seed(1328)
m2.5 <- mcmc(m$models[["tau=0.025"]],
  n.iter = 1200, burnin = 200, thin = 1)

summary(m2.5)

## predict credible intervals
nd <- data.frame(age = seq(6, 80, length = 300))
nd$p <- predict(m2.5, parameter = "mu", newdata = nd,
  FUN = function(x) quantile(x, probs = c(0.025, 0.5, 0.975)))

## plot
plot(fev1 ~ age, data = d, col = adjustcolor(1, 0.3))
matlines(nd$age, nd$p, lty = c(2, 1, 2), lwd = 2, col = 2)

## End(Not run)

Description

A wrapper around gamlss2 for fitting Bayesian GAMLSS models by MCMC sampling. Starting values are computed using maxit backfitting iterations before sampling. The controls n.iter, burnin, and thin can be specified directly; weights and offset are passed through to gamlss2.

Usage

bamlss2(formula, n.iter = 1200, burnin = 200, thin = 1, maxit = 2, ...)

Arguments

Value

An object of class "bamlss2" inheriting from "gamlss2" containing the fitted model and posterior samples.

References

Umlauf N, Klein N, Zeileis A (2018). “BAMLSS: Bayesian Additive Models for Location, Scale, and Shape (and Beyond)”. Journal of Computational and Graphical Statistics, 27(3), 612–627. doi:10.1080/10618600.2017.1407325

See Also

gamlss2, mcmc, BS

Examples

## Not run: 
## air quality data
air <- subset(airquality, !is.na(Ozone))

## set the seed
set.seed(123)

## estimate Gamma model with MCMC
m <- bamlss2(
  Ozone ~ s(Temp) + s(Wind) | s(Temp),
  data = air,
  family = GA
)

## posterior summary
summary(m)

## plot estimated effects
plot(m)

## plot samples
plot(m, which = "samples", ask = TRUE)

## predict parameters using samples
pm <- predict(m, type = "parameter", FUN = mean)
print(head(pm))
psd <- predict(m, type = "parameter", FUN = sd)
print(head(psd))

## End(Not run)

Description

The function BS() implements an MCMC sampler for generalized additive models for location, scale, and shape (GAMLSS) as fitted by gamlss2. The sampler uses efficient Metropolis-Hastings updates for regression coefficients of linear and smooth terms.

Usage

BS(x, y, specials, family, offsets,
  weights, start, xterms, sterms, control)

Arguments

Details

Function BS() is typically called via mcmc or via the wrapper function bamlss2. It implements a blocked Metropolis-Hastings sampler using working responses and working weights. Smooth terms are updated using penalized weighted least squares proposals and their smoothing variances are sampled using univariate slice sampling.

The function uses the following control arguments:

• n.iter: integer, the total number of MCMC iterations;

• burnin: integer, the burn-in period;

• thin: integer, the thinning parameter for saved samples;

• trace: logical, should information be printed while the sampler is running;

• flush: logical, whether to use flush.console to display current output in the console;

• nullmodel: logical, should an intercept-only null model be evaluated to compute a null deviance and deviance reduction?;

• binning: logical, should binning be used for smooth terms if available?

Value

A fitted model object containing posterior samples and summary information. Typically it includes:

• samples: posterior draws for coefficients and smoothing variances;

• coefficients: posterior means of coefficients;

• fitted.values: posterior mean fitted values for each distributional parameter;

• edf: estimated degrees of freedom for smooth terms;

• logLik: log-likelihood at the posterior mean predictor.

References

Umlauf N, Klein N, Zeileis A (2018). “BAMLSS: Bayesian Additive Models for Location, Scale, and Shape (and Beyond)”. Journal of Computational and Graphical Statistics, 27(3), 612–627. doi:10.1080/10618600.2017.1407325

See Also

gamlss2, RS, CG, mcmc, bamlss2

Examples

## see ?mcmc

Description

Compute and plot calibration curves for models with 0/1 responses. The function can handle one or several fitted gamlss2 models and compares their calibration in a single plot.

Usage

calibration(..., newdata = NULL,
  y = NULL, parameter = NULL, breaks = seq(0, 1, by = 0.1),
  minn = 20, main = "Calibration plot",
  xlab = "Predicted probability",
  ylab = "Observed proportion", plot = TRUE,
  add_loess = TRUE, smooth_n = 200,
  col = NULL, lty = NULL, legend = TRUE, pos = "topleft",
  xlim = NULL, ylim = NULL)

Arguments

Details

For each fitted model, predicted probabilities are obtained via predict.gamlss2 using type = "parameter" and the selected parameter (all family parameters when NULL). If the family provides a probabilities() function, it converts these predictions to probabilities; otherwise, the predicted parameters are used directly. If the result is a matrix or data frame, only its first column is used, after any conversion. Thus, without probabilities(), multiple predicted parameters are reduced to the first parameter column. The values used for calibration must represent the probability of y = 1.

Predicted probabilities are grouped into bins defined by breaks. Within each bin, the mean predicted probability, the observed proportion of 1s, and the number of observations are computed. Bins with fewer than minn observations are dropped.

Value

Invisibly returns a data frame with the columns interval, probs, y, n, and model. For a single model, the model column is dropped.

References

Van Calster B, McLernon DJ, van Smeden M, Wynants L, Steyerberg EW (2019). “Calibration: The Achilles Heel of Predictive Analytics”. BMC Medicine, 17, 230. doi:10.1186/s12916-019-1466-7

See Also

gamlss2, predict

Examples

## Not run: set.seed(123)
n <- 400
x1 <- runif(n, -3, 3)
x2 <- runif(n, -3, 3)
eta <- sin(x1) + 0.5 * x2
p <- plogis(eta)
y <- rbinom(n, size = 1, prob = p)
d <- data.frame(y = y, x1 = x1, x2 = x2)

m1 <- gamlss2(y ~ x1, family = BI, data = d)
m2 <- gamlss2(y ~ s(x1), family = BI, data = d)
m3 <- gamlss2(y ~ s(x1) + x2, family = BI, data = d)

calibration(m1, m2, m3)
head(calibration(m1, m2, m3, plot = FALSE))

## End(Not run)

Description

Constructs a conditional inference forest model term using partykit.

Usage

cf(formula, ...)

Arguments

Details

The forest is refitted with the current working response and weights during backfitting. The fitted contribution is centered for identifiability.

Value

An object of classes "special" and "cf".

See Also

special_terms

Examples

## Not run: ## air quality data
air <- subset(airquality, !is.na(Ozone))

set.seed(6020)
m <- gamlss2(
  Ozone ~ cf(Temp, minsplit = 2L, minbucket = 2L, alpha = 0.1) | .,
  data = air,
  family = GA
)

plot(m)

## predicted distributions and quantiles
d <- predict(m, newdata = air)
head(d, 3)
pq <- quantile(d, probs = c(0.025, 0.5, 0.975))

plot(Ozone ~ Temp, data = air,
  ylim = range(air$Ozone, pq))
i <- order(air$Temp)
matlines(air$Temp[i], pq[i, ], lty = 1, lwd = 2, col = 4)

## End(Not run)

Description

Constructs a conditional inference tree model term using partykit.

Usage

ct(formula, ...)

Arguments

Details

The tree is refitted with the current working response and weights during backfitting. The fitted contribution is centered for identifiability.

Value

An object of classes "special" and "ct".

See Also

special_terms

Examples

## Not run: 
m <- gamlss2(
  mag ~ ct(depth) | .,
  data = quakes,
  family = truncate_family(NO, lower = 3.99)
)

plot(m)

## predicted distributions and quantiles
d <- predict(m, newdata = quakes)
pq <- quantile(d, probs = c(0.025, 0.5, 0.975))

plot(mag ~ depth, data = quakes,
  ylim = range(quakes$mag, pq),
  col = adjustcolor(1, 0.3))
i <- order(quakes$depth)
matlines(quakes$depth[i], pq[i, ], lty = 1, lwd = 2, col = 4)

## End(Not run)

Description

cv_gamlss2() implements K-fold cross validation for models fitted with gamlss2. Different scoring rules can be supplied via the metric argument. Convenience metric functions (log_pdf_metric(), rqres_metric(), mse_metric()) are provided.

Usage

## k-fold cross-validation
cv_gamlss2(..., data, folds = 5,
  metric = log_pdf_metric, simplify = TRUE, cores = 1L)

## log-pdf for each observation
log_pdf_metric(model, data)

## randomized quantile residuals
rqres_metric(model, data)

## mean squared error
mse_metric(model, data)

Arguments

Details

cv_gamlss2() splits the data into training and test folds. For each fold the model is fitted on the training data, and the chosen metric is evaluated on the held-out test data. By default, the scoring rule is the log predictive density (log_pdf_metric), but other metrics can be used, such as randomized quantile residuals (rqres_metric) or mean squared error of the conditional mean (mse_metric).

The function returns either a list of fold-wise results or, if simplify = TRUE, a named vector or data frame aligned with the original observations.

Set the random seed before calling cv_gamlss2() to reproduce randomly generated folds. Parallel evaluation also uses reproducible random-number streams for stochastic fits or metrics.

Value

If simplify = TRUE and the metric returns scalars, a named numeric vector of fold scores is returned. Otherwise a data frame with fold membership and scores per observation is returned. The convenience metrics return a numeric vector of scores or residuals.

See Also

gamlss2, log_pdf

Examples

## Not run: ## air quality data
air <- subset(airquality, !is.na(Ozone))

## set the seed
set.seed(123)

## cross-validation using the GA distribution
## only model mu parameter with s(Temp)
cv1 <- cv_gamlss2(
  Ozone ~ s(Temp),
  data = air,
  family = GA, maxit = 100
)

## now, add effects
cv2 <- cv_gamlss2(
  Ozone ~ s(Temp) + s(Wind) | s(Temp),
  data = air,
  family = GA,
  folds = cv1$fold, maxit = 100
)

## evaluate log-likelihood
sum(cv1$score)
sum(cv2$score)

## evaluate folds with two workers, reusing the same fold assignments
cv_parallel <- cv_gamlss2(
  Ozone ~ s(Temp), data = air, family = GA,
  folds = cv1$fold, maxit = 100, cores = 2L
)
sum(cv_parallel$score)

## End(Not run)

Description

Extracts the deviance from fitted gamlss2 objects. By default, the total deviance is returned. Optionally, observation-wise deviance contributions are returned.

Usage

## S3 method for class 'gamlss2'
deviance(object, ..., newdata = NULL, sum = TRUE)

Arguments

Details

The deviance of a fitted gamlss2 model is defined as minus twice the fitted log-likelihood,

$-2 \log L$

. If sum = TRUE, one deviance value is returned for each supplied model. If more than one model is supplied, the deviances are ordered increasingly.

If sum = FALSE, the function returns the observation-wise deviance contributions,

$-2 \log f(y_i | \hat\theta_i)$

, where $f()$ is the fitted probability density or probability mass function and $\hat\theta_i$ denotes the fitted distribution parameters for observation $i$. These values sum to the total deviance and are useful, for example, in cross-validation.

Value

If sum = TRUE, a numeric vector with the total deviance of the supplied model or models. If sum = FALSE, a numeric vector of observation-wise deviance contributions for one model, or a matrix with one column per model if several models are supplied.

See Also

gamlss2

Examples

## Not run: ## specify a distributional model:
## location and scale vary with x,
## shape parameters remain constant
f <- dist ~ s(speed) | 1

## estimate model
m1 <- gamlss2(f, data = cars, family = GA)

## total deviance
deviance(m1)

## observation-wise deviance contributions
di <- deviance(m1, sum = FALSE)

## the contributions sum to the total deviance
sum(di)
deviance(m1)

## compare two models
m2 <- gamlss2(dist ~ s(speed) | s(speed), data = cars, family = GA)
deviance(m1, m2)

## compare observation-wise deviance contributions
di <- deviance(m1, m2, sum = FALSE)
head(di)

## deviance on new data
nd <- cars[1:5, ]
deviance(m1, newdata = nd)
deviance(m1, newdata = nd, sum = FALSE)

## End(Not run)

Description

Constructs a discrete distribution family by assigning the probability of each fixed interval of a base distribution to a corresponding response value.

Usage

discretize_family(family = NO, breaks,
  values = seq_len(length(breaks) - 1L), right = TRUE)

Arguments

Details

For each interval defined by adjacent elements of breaks, discretize_family() assigns the interval probability under the base distribution to the corresponding element of values. Thus, with right = TRUE, response value $y_j$ has probability

$P(Y \in (b_j, b_{j+1}])$

under the base distribution. With right = FALSE, the intervals are left-closed instead.

The returned family provides probability mass, cumulative distribution, quantile, random-number, mean, and variance functions. Its distributional parameters and link functions are inherited from the base family.

Value

An object of class "gamlss2.family" representing the discretized distribution.

See Also

family_extensions, round_family, gamlss2, gamlss2.family

Examples

## Not run: 
## fixed intervals and their observed values
breaks <- c(-Inf, -1, 0, 1, Inf)
values <- 0:3

## discretize a normal family
fam <- discretize_family(NO, breaks = breaks, values = values)

## simulate data whose location depends on x
set.seed(6020)
n <- 500
x <- sort(runif(n, -3, 3))
par <- data.frame(mu = sin(x), sigma = exp(-0.5 + cos(x)))
y <- fam$random(par, n)
head(y)
df <- data.frame(y = y, x = x)

## plot count data
plot(y ~ x, data = df)

## fitted pdf
fit_family(df$y, family = fam)

## fit a regression model using the discretized family
m <- gamlss2(y ~ s(x) | ., data = df, family = fam)

## inspect the estimated effect and fitted probabilities
plot(m)
d <- predict(m)
head(pdf(d, y))

## predict median
p <- quantile(d, probs = 0.5)
plot(y ~ x, data = df)
lines(p ~ x, lwd = 2, col = 4)

## End(Not run)

Description

Constructor function for Extreme Learning Machine (ELM) model terms for GAMLSS.

Usage

## Model term constructor function.
elm(x, k = 50, a = "tanh", ...)

Arguments

Details

The ELM term constructs a randomized single-hidden-layer representation of the covariate(s) in x. Internally, a design matrix Z is built (including an intercept column), random weights are sampled, and the linear predictors Z %*% W are transformed through the activation function a to obtain the hidden-layer design matrix X. Columns of X are centered before estimation.

For numerical stability and comparability across terms, Z can be scaled internally (scale = TRUE). Intercept columns are left unchanged. If x is a factor, a QR-based group scaling is applied; otherwise, a column-wise center-and-scale normalization is used.

The activation is applied element-wise to Z %*% W. For numerical stability, bounded activations are evaluated on clipped inputs (currently $[-35, 35]$).

By default, hidden-unit weights are generated using a low-discrepancy Halton design in a whitened input space. This tends to produce more stable and less redundant features than purely random draws in many practical applications. With sampler = "pca", directions are additionally aligned with principal component directions of the scaled design. The historical stochastic behavior can be recovered with sampler = "random".

The default sampler = "halton" is deterministic given the input data. Use set.seed only if you explicitly request sampler = "random".

Value

An object representing a specials model term, used internally by gamlss2 during model fitting and prediction.

References

Huang GB, Zhu QY, Siew CK (2006). “Extreme Learning Machine: Theory and Applications”. Neurocomputing, 70(1–3), 489–501. doi:10.1016/j.neucom.2005.12.126

Huang GB, Zhu QY, Siew CK (2004). “Extreme Learning Machine: A New Learning Scheme of Feedforward Neural Networks”. In 2004 IEEE International Joint Conference on Neural Networks, 2, 985–990. IEEE, Budapest, Hungary. doi:10.1109/IJCNN.2004.1380068

See Also

gamlss2, specials.

Examples

## Not run: ## load data
data("SpirometryUS")

## subset for female
d <- subset(SpirometryUS, gender == "female")

## the default sampler is deterministic; set
## a seed only for sampler = "random"
## formula for all 4 parameters
f <- fev1 ~ elm(age,k=100,a="tanh") | . | . | .

## estimate model
m <- gamlss2(f, data = d, family = BCT)

## estimated effects
plot(m)

## predict quantiles
qu <- c(0.025, seq(0.1, 0.9, by = 0.1), 0.975)
fit <- quantile(m, probs = qu)

## plot
plot(fev1 ~ age, data = d)
i <- order(d$age)
matlines(d$age[i], fit[i, ],
  lty = 1, lwd = 2, col = c(2, rep(4, ncol(fit) - 1)))

## main effects and interactions
f <- fev1 ~ s(age) + s(height) + s(weight) +
  elm(~age+height+weight,k=200) | . | . | .

m <- gamlss2(f, data = d, family = BCT)

## summary to inspect effect of interactions
## of the elm()s
summary(m)

## plot main effects
plot(m)

## quantile residuals
plot(m, which = "resid")

## prediction is handled automatically via
## the special term interface
n <- 50
nd <- with(d, expand.grid(
  "age" = seq(min(age), max(age), length = n),
  "weight" = seq(min(weight), max(weight), length = n)
))
nd$height <- mean(d$height)

## compute lower 2.5
nd$fit <- quantile(m, newdata = nd, probs = 0.025)

## visualize
n <- length(unique(nd$age))
age <- sort(unique(nd$age))
weight <- sort(unique(nd$weight))

z <- matrix(nd$fit, n, n)

image(age, weight, z,
      col = hcl.colors(100, "YlOrRd"),
      xlab = "age", ylab = "weight")
contour(age, weight, z, add = TRUE)

## End(Not run)

Description

Create a simplified formula from a formula, a Formula, or a list of formulas. The result retains the variables and transformations needed to build a model.frame while separating out special model terms. This is mainly used internally when gamlss2 prepares formulas for fitting.

Usage

fake_formula(formula, specials = NULL,
  nospecials = FALSE, onlyspecials = FALSE)

Arguments

Details

The function is primarily an internal helper for formula preprocessing in gamlss2, but it can also be useful for inspecting how complex formulas are decomposed before fitting.

Special model terms such as s(), te(), ti(), re(), or package-specific specials are identified and handled separately from ordinary variables. Depending on the control arguments:

• if nospecials = FALSE, variables used inside registered special terms are reintroduced into the returned formula so they remain available in the model frame;

• if nospecials = TRUE, the returned formula excludes the special terms and their variables;

• if onlyspecials = TRUE, the function returns the names of the detected special terms instead of a formula object.

For multi-part or multi-parameter formulas, the output preserves the overall formula structure.

Value

Depending on the input, the function returns a formula or Formula. If onlyspecials = TRUE, a character vector or list of special model term names is returned.

See Also

gamlss2

Examples

## basic formula
f <- y ~ x1 + x2 + log(x3)
ff <- fake_formula(f)
print(ff)

## including special model terms
f <- y ~ x1 + s(x2) + x3 + te(log(x3), x4)
ff <- fake_formula(f)
print(ff)

## multiple parts on the right-hand side
f <- y ~ x1 + s(x2) + x3 + te(log(x3), x4) | x2 + sqrt(x5)
ff <- fake_formula(f)
print(ff)

## collapse all formula parts
print(formula(ff, collapse = TRUE))
print(formula(ff, collapse = TRUE, update = TRUE))

## list of formulas
f <- list(
  y ~ x1 + s(x2) + x3 + te(log(x3), x4),
    ~ x2 + sqrt(x5),
    ~ z2 + x1 + exp(x3)
)
ff <- fake_formula(f)
print(ff)

## extract separate parts on the right-hand side
formula(ff, rhs = 1)
formula(ff, rhs = 2)
formula(ff, rhs = 3)

## formula with multiple responses and multiple parts
f <- y1 | y2 | y3 ~ x1 + s(x2) + x3 + te(log(x3), x4) | x2 + ti(x5)
ff <- fake_formula(f)
print(ff)

## list of formulas with multiple responses
f <- list(
  y1 ~ x1 + s(x2) + x3 + te(log(x3), x4),
  y2 ~ x2 + sqrt(x5),
  y3 ~ z2 + x1 + exp(x3) + s(x10)
)
ff <- fake_formula(f)

## extract only without special terms
ff <- fake_formula(f, nospecials = TRUE)
print(ff)

## extract only special terms
ff <- fake_formula(f, onlyspecials = TRUE)
print(ff)

Description

Extends an existing distribution family by transforming, conditioning, coarsening, or adding point masses. The resulting "gamlss2.family" objects can be used directly in gamlss2 models, and family extensions can be nested to compose several operations.

Usage

transform_family(family = NO, transform, inverse, log_jacobian,
  increasing = TRUE, support = c(-Inf, Inf), valid.response = NULL,
  name = NULL)

logit_family(family = NO)

truncate_family(family = NO, lower = -Inf, upper = Inf,
  include = c(FALSE, TRUE))

censor_family(family = NO, lower = -Inf, upper = Inf)

inflate_family(family = NO, at = 0)

hurdle_family(family = PO, at = 0)

round_family(family = NO, digits = 0L)

Arguments

Details

For a fixed monotone transformation $Y = g(X)$, transform_family uses

$\log f_Y(y) = \log f_X(g^{-1}(y)) + \log |d g^{-1}(y) / dy|.$

Increasing and decreasing transformations are supported. The convenience function logit_family constructs the logistic push-forward, which has strict support $0 < Y < 1$. Exact zeros and ones are not valid observations for this continuous transformation.

The function truncate_family conditions the base distribution on the interval selected by lower, upper, and include. Its density or probability mass is divided by the probability of the retained interval. With at least one finite boundary, censor_family observes

$Y = \min(\max(X, lower), upper)$

and consequently creates point masses at finite censoring boundaries. With both boundaries left infinite, it instead accepts a right-censored survival::Surv(time, status) response. An observed event contributes the base density or mass at time; a censored observation contributes the base survival probability beyond time.

For inflation points $c_j$, inflate_family constructs

$F(y) = \pi_0 F_0(y) + \sum_j \pi_j I(c_j \le y).$

Natural base-distribution mass at an inflation point is retained. The appended parameters pi1, pi2, and so on have log links and represent positive odds $w_j$ relative to the base component, not probabilities. Thus,

$\pi_0 = \frac{1}{1 + \sum_j w_j}, \qquad \pi_j = \frac{w_j}{1 + \sum_\ell w_\ell}.$

This is the baseline-odds parameterization also used by mixture_family. Quantiles account for every probability jump.

The function hurdle_family instead constructs

$F(y) = \pi_1 I(c \le y) + \pi_0 F_0(y \mid X \ne c).$

Here pi1 is also the odds $w_1$, so the hurdle probability is $w_1/(1 + w_1)$. This is the complete mass at $c$, including when the base distribution already assigns positive mass to that point.

The function round_family uses intervals centered on a decimal grid and is restricted to continuous base distributions. Its quantile and random-number functions follow R's round() convention. For fixed, user-defined intervals, see discretize_family.

Inflated, hurdle, and fixed-bin families provide exact means and variances when the required base-family moments are available. Censored and rounded moments requiring unavailable partial moments report an error.

Continuous distributions with added point masses provide randomized quantile residuals. Bounds, atom locations, breaks, and transformations are fixed when the family is constructed. Other extensions requiring numerical normalization or derivatives during fitting are not provided.

Value

An object of class "gamlss2.family". Existing base-family parameters retain their names, order, and links. Inflation and hurdle families append positive-odds parameters named pi1, pi2, and so on.

See Also

gamlss2, gamlss2.family, discretize_family, mixture_family, distributions3_family

Examples

## Not run: ## simulate data
set.seed(6020)
n <- 300
x <- sort(runif(n, -3, 3))

## logistic transformation of a normal distribution
fam <- logit_family(NO)
y <- fam$random(list(mu = sin(x), sigma = exp(-1 + cos(x))), n)
hist(y)
m <- gamlss2(y ~ s(x) | s(x), family = fam)
plot(m)

## zero-truncated Poisson distribution
fam <- truncate_family(PO, lower = 0)
y <- fam$random(list(mu = exp(1/2 + sin(x))), n)
m <- gamlss2(y ~ s(x), family = fam)
plot(x, y)

## small helper for plotting fitted quantiles
fq <- function(model) {
  matlines(x,
    quantile(model, probs = c(0.05, 0.5, 0.95)),
      lty = 1, lwd = 2, col = 4
  )
}

fq(m)

## normal distribution censored at both boundaries
fam <- censor_family(NO, lower = -1, upper = 1)
y <- fam$random(list(mu = sin(x), sigma = exp(-1 + cos(x))), n)
hist(y)
m <- gamlss2(y ~ s(x) | s(x), family = fam)
plot(m)
plot(x, y)
fq(m)

## zero-inflated Poisson distribution
fam <- inflate_family(PO, at = 0)
## pi1 is the added-mass odds; its mixture probability is plogis(cos(x))
y <- fam$random(list(mu = exp(1/2 + sin(x)), pi1 = exp(-1 + cos(x))), n)
barplot(table(y))
m <- gamlss2(y ~ s(x) | s(x), family = fam)
plot(m)
plot(x, y)
fq(m)

## zero-hurdle Poisson distribution
fam <- hurdle_family(PO, at = 0)
## pi1 is the hurdle odds; the probability of zero is plogis(cos(x))
y <- fam$random(list(mu = exp(1/2 + sin(x)), pi1 = exp(cos(x))), n)
barplot(table(y))
m <- gamlss2(y ~ s(x) | s(x), family = fam)
plot(m)
plot(x, y)
fq(m)

## decimal rounding
fam <- round_family(NO, digits = 1)
y <- fam$random(list(mu = sin(x), sigma = exp(-1 + cos(x))), n)
m <- gamlss2(y ~ s(x) | s(x), family = fam)
plot(m)
plot(m, which = "resid")

## End(Not run)

Description

Set up a gamlss2.family object based on an object inheriting from "distribution" (as provided by the distributions3 package).

Usage

## S3 method for class 'distribution'
family(
  object,
  links,
  score = TRUE,
  hessian = FALSE,
  update = FALSE,
  expected = FALSE,
  ...
)

## S3 method for class 'GAMLSS'
family(object, ...)

## S3 method for class 'gamlss.family'
family(object, score = TRUE, hessian = FALSE, update = FALSE, ...)

distributions3_family(
  distribution,
  links,
  score = TRUE,
  hessian = FALSE,
  update = FALSE,
  expected = FALSE,
  create_distribution = "structure",
  type = NULL,
  valid.response = NULL,
  initialize = NULL,
  ...
)

Arguments

Details

distributions3_family maps fitted parameter lists to distributions3 objects and uses their methods for fitting and assessing the model. The family methods call this function and infer some family settings automatically.

Value

A list inheriting from gamlss2.family.

Examples

library("distributions3")
f <- family(Poisson(), links = c(lambda = "log"))
str(f)

Description

These functions provide useful infrastructures for finding suitable GAMLSS families for a response variable.

Usage

## list of available families from gamlss.dist package
available_families(type = c("continuous", "discrete"), families = NULL, ...)

## find suitable response distribution
find_family(y, families = NULL, k = 2, verbose = TRUE, ...)

## fit distribution parameters
fit_family(y, family = NO, plot = TRUE, ...)

## find gamlss2 model
find_gamlss2(formula, families = NULL, k = 2,
  select = FALSE, verbose = TRUE, ...)

Arguments

Details

The function find_family() employs gamlss2 to estimate intercept-only models for each specified family object in the families argument. Note that model estimation occurs within a try block with warnings suppressed. Additionally, the function calculates the GAIC for each family whenever feasible and returns the values sorted in decreasing order of the selected fit criterion, so smaller values indicate a better fit.

Function fit_family() fits a single intercept-only model using the specified family and creates a plot of the fitted density.

Value

Function find_family() returns a vector of GAIC values for the different fitted families. The best model is the model with the lowest GAIC.

Function fit_family() returns the fitted intercept-only model.

Function find_gamlss2() returns the final selected model.

See Also

gamlss2.

Examples

## Not run: ## find a suitable response to the response
ic <- find_family(airquality$Ozone)
print(ic)

## fit parameters using the WEI family
fit_family(airquality$Ozone, family = WEI)

## count data,
## search best count model
ic <- find_family(InsectSprays$count, k = 0,
  families = available_families(type = "discrete"))
print(ic)

## fit parameters using the ZANBI family
fit_family(InsectSprays$count, family = ZANBI)

## search complete model
## m <- find_gamlss2(Ozone ~ s(Temp) + s(Wind) | . | . | .,
##   data = airquality, select = TRUE, mu.link = "log")

## End(Not run)

Description

Constructs a univariate regression spline whose interior knot locations, and optionally the number of interior knots, are selected from the current weighted working response.

Usage

fk(x, k = 10L, nk = NULL, degree = 1L,
  criterion = "bic", knots = NULL, K = NULL,
  reselect = FALSE, optim.control = list(), ...)

Arguments

Details

For each candidate number of interior knots, fk() uses a B-spline basis with fixed, repeated boundary knots. Interior knots are represented by positive ordered gaps. This removes permutation ambiguity and prevents knot crossings during optimization.

Conditional on the knot locations, spline coefficients are profiled out by weighted least squares. Knot optimization therefore only operates on the interior locations. The fitted model term is centered.

When nk = NULL, models with one through k interior knots are compared using criterion. Effective degrees of freedom include both the centered spline coefficients and the estimated interior knot locations. Fixed locations supplied through knots do not add knot-location degrees of freedom.

Value

The constructor returns an object of classes "special" and "fk".

See Also

gamlss2, special_terms

Examples

## Not run: ## load data
data("SpirometryUS")

## subset for female
d <- subset(SpirometryUS, gender == "male")

## Automatically select one free interior knot
m <- gamlss2(fev1 ~ fk(age, nk = 1) | ., data = d, family = NO)

## estimated effects
plot(m)

## predict quantiles
qu <- c(0.025, seq(0.1, 0.9, by = 0.1), 0.975)
fit <- quantile(m, probs = qu)

## plot
plot(fev1 ~ age, data = d)
i <- order(d$age)
matlines(d$age[i], fit[i, ],
  lty = 1, lwd = 2, col = c(2, rep(4, ncol(fit) - 1)))

## End(Not run)

Description

The function gamlss2() fits generalized additive models for location, scale, and shape (GAMLSS). It estimates one additive predictor for each parameter of the chosen response distribution, such as mu, sigma, nu, or tau. The number and meaning of these parameters are determined by the selected family object, see gamlss2.family.

In addition to ordinary linear terms, gamlss2() supports smooth terms from mgcv as well special user-defined model-terms, see specials.

Usage

gamlss2(formula, ...)

## S3 method for class 'formula'
gamlss2(formula, data, family = NO,
  subset, na.action, weights, offset, start = NULL,
  knots = NULL, control = gamlss2_control(...), ...)

## S3 method for class 'list'
gamlss2(formula, ...)

Arguments

Details

The function gamlss2() is the main entry point for fitting distributional regression models in gamlss2.

• The response distribution is specified through a family object, either from gamlss.dist or created as a gamlss2.family object.

• By default, estimation is carried out with the RS algorithm. The optimizer can be changed through gamlss2_control. A family object may also provide its own optimizer.

• The returned object depends on the chosen optimizer. In the standard case, this is an object of class "gamlss2", for which summary, plotting, prediction, and extractor methods are available.

Value

The object returned by the selected optimizer. By default, this is an object of class "gamlss2" produced by RS. Standard methods and extractor functions are available for this class.

See Also

RS, gamlss2_control, gamlss2.family, gamlss2_start

Examples

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

## positive skewed response
hist(air$Ozone)

## specify a distributional model
f <- Ozone ~ s(Temp) + s(Wind) | s(Temp)

## estimate model
m1 <- gamlss2(f, data = air, family = GA)

## model summary
summary(m1)

## plot parameter-specific effects
plot(m1, which = "effects")

## plot diagnostics
plot(m1, which = "resid")

## predict distribution parameters
par <- predict(m1, parameter = c("mu", "sigma"))
print(head(par))

## create new data:
## vary temperature while holding wind constant
nd <- with(air, data.frame(
  Temp = seq(min(Temp), max(Temp), length.out = 100),
  Wind = median(Wind)
))

## fitted conditional quantiles
pq <- quantile(
  m1,
  newdata = nd,
  probs = c(0.05, 0.5, 0.95)
)

## visualize
plot(Ozone ~ Temp, data = air, pch = 19,
  col = rgb(0.1, 0.1, 0.1, alpha = 0.3),
  xlab = "Temperature (degrees F)",
  ylab = "Ozone (ppb)")
matlines(nd$Temp, pq,
  lwd = 2,
  lty = c(2, 1, 2),
  col = 4)

## exceedance probabilities
nd5 <- with(air, data.frame(
  Temp = seq(min(Temp), max(Temp), length.out = 5),
  Wind = min(Wind)
))

## use distributions3 methods
d <- predict(m1, newdata = nd5)
nd5$Probs <- 1 - cdf(d, 90)

barplot(
  nd5$Probs,
  names.arg = round(nd5$Temp),
  xlab = "Temperature (degrees F)",
  ylab = "Pr(Ozone > 90 ppb)",
  ylim = c(0, 1)
)

## count response example using HarzTraffic data
data("HarzTraffic", package = "gamlss2")

## seasonal count model for motorcycles
m2 <- gamlss2(
  bikes ~ s(yday, bs = "cc") | s(yday, bs = "cc"),
  data = HarzTraffic,
  family = NBI
)

## summary
summary(m2)

## effect plots
plot(m2)

## residual diagnostics
plot(m2, which = "resid")

## quantiles over the day of the year
nd <- data.frame(yday = 1:365)
d <- predict(m2, newdata = nd)
pq <- quantile(d, probs = c(0.05, 0.5, 0.95))

## visualize fitted quantiles
plot(bikes ~ yday, data = HarzTraffic, pch = 19,
  col = gray(0.1, alpha = 0.3))
matlines(nd$yday, pq, lty = c(2, 1, 2),
  lwd = 2, col = 4)

Description

The GAMLSS2 family is a unified class with corresponding methods that represent all distributional families supported by the gamlss2 package. It enables seamless integration with the distributions3 workflow and provides a consistent interface for model fitting and distributional computations.

Usage

GAMLSS2(family, parameters)

Arguments

Details

The S3 class GAMLSS2 is a slightly more general implementation of the S3 class GAMLSS tailored for gamlss2. For details please see the documentation of GAMLSS

Value

A GAMLSS2 object, inheriting from distribution.

References

Zeileis A, Lang MN, Hayes A (2022). “distributions3: From Basic Probability to Probabilistic Regression”. Presented at useR! 2022 - The R User Conference. Slides, video, vignette, code at https://www.zeileis.org/news/user2022/.

See Also

gamlss2.family

Examples

## package and random seed
library("distributions3")
set.seed(6020)

## one normal distribution
X <- GAMLSS2("NO", c(mu = 1, sigma = 2))
X

## two normal distributions
X <- GAMLSS2("NO", cbind(c(1, 1.5), c(0.6, 1.2)))
X

## three Weibull distributions
X <- GAMLSS2("WEI", list(mu = c(1, 1, 2), sigma = c(1, 2, 2)))
X

## see ?gamlss.dist::GAMLSS for the remainder of this example

Description

Control parameters for fitting GAMLSS models with gamlss2.

Usage

gamlss2_control(optimizer = RS, trace = TRUE,
  flush = TRUE, light = FALSE, expand = TRUE,
  model = TRUE, x = TRUE, y = TRUE,
  fixed = FALSE, ...)

Arguments

Details

The set of control parameters can be extended. For example, if a different optimizer is used, newly specified control parameters are passed on to that optimizer automatically.

Optimizer-specific controls and algorithm details are documented in RS and CG, and in JR for joint REML estimation. Smoothing parameter selection is described in gamlss2_smoothing.

Value

A list with the arguments specified.

See Also

RS, CG, JR, gamlss2, gamlss2_start, gamlss2_smoothing

Examples

## Not run: ## specify the model formula
f <- dist ~ s(speed) | s(speed)

## specify fitting controls
ctrl <- gamlss2_control(trace = FALSE, maxit = c(50, 100))
m <- gamlss2(f, data = cars, family = GA, control = ctrl)

## End(Not run)

Description

Smoothing parameter selection for smooth terms in models fitted with gamlss2.

Details

For mgcv-based smooth terms, smoothing parameters are selected by local REML updates by default (criterion = "reml"). These updates use the current working response and weights. With one penalty matrix this is the usual Schall update; with multiple penalty matrices all smoothing parameters are updated together using the extended Fellner–Schall trace equations. This is a local working-model REML update. Explicit choices such as "gcv", "aic", "aicc", "bic", and "ncv" retain their respective selection methods.

These local updates are used by the RS and CG optimizers. With optimizer = JR, coefficients and smoothing parameters are estimated jointly by Laplace approximate REML; see JR.

Neighbourhood cross-validation can be selected with criterion = "ncv". If ncv is omitted, singleton neighbourhoods are used, corresponding to leave-one-observation-out cross-validation. For ordered data, a numeric ncv = k uses the observation together with up to k neighbours on either side. Larger neighbourhoods generally provide more smoothing and can account for short-range dependence, at increased computational cost. For ordered data, small values such as ncv = 1 or ncv = 2 are often preferable to the singleton default. Custom neighbourhoods can also be supplied as a list through a smooth term's xt argument, for example s(x, xt = list(ncv = neighbourhoods)), where neighbourhoods has one integer index vector per observation and each vector contains its target observation.

All methods return coefficients and covariance at the reported smoothing parameters.

For the RS and CG optimizers, argument demmler.reinsch controls the use of a weighted Demmler–Reinsch reparameterization when selecting the smoothing parameter of an eligible single-penalty smooth, including multivariate smooths. The default, "auto", uses the native reparameterization for local REML updates, and for quadratic criteria after the initial outer iteration. Set demmler.reinsch = TRUE to force it for eligible updates, or FALSE to disable it. For quadratic criteria, the initial outer iteration continues to use the direct solver; with numerical gradients or native.wfit = FALSE, automatic reparameterization is delayed until the current search reaches the edge of its one-decade search window.

Within an RS or CG fit, the reparameterization is cached while the design matrix, weights, penalty, and binning are unchanged. A new working response only requires updating its projection into the cached basis. The declared penalty rank keeps the null space unpenalized even at large smoothing parameters.

Analytic gradients are used by default for quadratic smoothing criteria (GCV, AIC, GAIC, AICc, and BIC), including smooths with multiple penalties. Set analytic.gradient = FALSE to use numerical gradients. Criteria using logLik = TRUE continue to use numerical gradients, with native direct solves and family likelihood evaluation in R. Nearly exact fits use observation-space residuals to avoid cancellation in the quadratic RSS.

Smooths with fixed smoothing parameters and multiple-penalty smooths use the direct solver. For local REML, the common range of the penalty matrices is cached and the coupled update is evaluated in native code. Single-penalty updates use the constrained penalty null space and the Demmler–Reinsch basis when available.

References

Demmler A, Reinsch C (1975). “Oscillation Matrices with Spline Smoothing”. Numerische Mathematik, 24(5), 375–382. doi:10.1007/BF01437406

Wood SN (in press). “On Neighbourhood Cross Validation”. Statistical Science. Preprint: doi:10.48550/arXiv.2404.16490

See Also

gamlss2, gamlss2_control, RS, CG, JR, pb

Examples

## Not run: ## specify the model formula
f <- dist ~ s(speed) | s(speed)

## estimate models with different smoothing criteria
m0 <- gamlss2(f, data = cars, family = GA, criterion = "reml")
m1 <- gamlss2(f, data = cars, family = GA, criterion = "aic")
m2 <- gamlss2(f, data = cars, family = GA, criterion = "bic")
m3 <- gamlss2(f, data = cars, family = GA, criterion = "gaic")
m4 <- gamlss2(f, data = cars, family = GA, criterion = "gaic", K = 4)
m5 <- gamlss2(f, data = cars, family = GA, criterion = "ncv")
m6 <- gamlss2(f, data = cars, family = GA, criterion = "ncv", ncv = 2)

plot(c(m0, m1, m2, m3, m4, m5, m6))

## End(Not run)

Description

Guidance on supplying starting values through argument start in gamlss2 and on fixing distribution parameters via fixed in gamlss2_control.

Details

The start argument can be used to initialize estimation in several ways. Supported inputs include:

• a named numeric vector or list giving one starting value per distribution parameter, such as mu, sigma, nu, or tau;

• a data frame or matrix with one column per distribution parameter and one row per observation, providing starting values for the full additive predictors on the link scale;

• a fitted "gamlss2" object, in which case fitted() values are used as starting values for all distribution parameters;

• a coef.gamlss2 object, allowing coefficient-wise initialization of linear and smooth terms.

Parameters can be fixed during estimation using the fixed control argument. This is useful, for example, when shape parameters should remain constant while the remaining parameters are estimated.

Starting values from a previous fit can reduce optimizer iterations. For coefficient-wise initialization, coef(object, full = TRUE) includes smooth coefficients; lambdas = TRUE also includes smoothing parameters. Compare iteration counts using the same model, observations, and controls.

See Also

gamlss2, gamlss2_control, RS

Examples

## Not run: air <- subset(airquality, !is.na(Ozone))

## hold the smoothing parameter fixed for comparable refits
f <- Ozone ~ s(Wind, sp = 1) | 1
ctrl <- gamlss2_control(maxit = 100, eps = 1e-7)

## scalar starting values for the distribution parameters
m1 <- gamlss2(f, data = air, family = GA,
  start = c(mu = 200, sigma = 0.1), control = ctrl)

## reuse fitted additive predictors on the link scale
m2 <- gamlss2(f, data = air, family = GA,
  start = fitted(m1, type = "link"), control = ctrl)

## coef() without full = TRUE initializes only the linear coefficients
m3 <- gamlss2(f, data = air, family = GA,
  start = coef(m1), control = ctrl)

## reuse all linear and smooth coefficients from the previous fit
m4 <- gamlss2(f, data = air, family = GA,
  start = coef(m1, full = TRUE, lambdas = TRUE, dropall = FALSE),
  control = ctrl)

## full coefficient starts reduce iterations in this example
## check convergence and compare the final log-likelihoods as well
fits <- list(scalar = m1, predictors = m2, linear = m3, full = m4)
data.frame(
  iterations = sapply(fits, function(m) m$iterations),
  converged = sapply(fits, function(m) m$converged),
  logLik = sapply(fits, function(m) as.numeric(logLik(m)))
)

## fix BCT shape parameters at their supplied starting values
m5 <- gamlss2(Ozone ~ s(Wind, sp = 1) | 1 | 1 | 1,
  data = air, family = BCT,
  start = c(mu = 200, sigma = 0.1, nu = 1, tau = 10),
  fixed = c(nu = TRUE, tau = TRUE), maxit = 100)

## fixed flags must also be supplied when reusing coefficients
m6 <- gamlss2(Ozone ~ s(Wind, sp = 1) | 1 | 1 | 1,
  data = air, family = BCT,
  start = coef(m5, full = TRUE, dropall = FALSE),
  fixed = c(nu = TRUE, tau = TRUE), maxit = 100)
unique(fitted(m6, parameter = c("nu", "tau")))

## End(Not run)

Description

A "gamlss2.family" object describes a response distribution for use with gamlss2. It tells the fitting algorithm which distributional parameters exist, which link functions they use, and which distributional functions are available, such as densities, quantiles, random-number generators, or score functions.

Usage

## S3 method for class 'gamlss2'
family(object, ...)
## S3 method for class 'gamlss2.family'
print(x, full = TRUE, ...)

Arguments

Details

A gamlss2 family object is a list of class "gamlss2.family". It defines the response distribution and the functions needed to fit and use that distribution in a model.

The minimum required elements are:

• family: the name of the distribution (character string).

• names: a character vector of parameter names (e.g., c("mu", "sigma")).

• links: a named character vector specifying the link function for each parameter (e.g., c(mu = "identity", sigma = "log")), or a list of link functions, see softplus.

• pdf(par, y, log = FALSE, ...): a function to evaluate the (log-)density.

The pdf() function must accept the response y, a named list par of evaluated parameter values (e.g., par$mu, par$sigma), a logical log, and optional additional arguments.

Further elements can be added depending on how fully the family should support fitting, prediction, diagnostics, and simulation. Common optional elements are:

• score: a named list of functions (one per parameter), each computing the first derivative of the log-likelihood with respect to the linear predictor: score[[param]](par, y, ...).

• hessian: a named list of functions computing second derivatives (the negative Hessian). For parameters mu and sigma, this includes: hessian[["mu"]](), hessian[["sigma"]](), and optionally cross derivatives like hessian[["mu:sigma"]]().

• log_likelihood(par, y, ...): a function computing the total log-likelihood.

• cdf(par, y, ...): cumulative distribution function.

• quantile(par, p, ...): quantile function.

• support(par, ...): lower and upper response support. It must return a numeric matrix with columns min and max, with one row per parameter combination. A fixed numeric vector c(min, max) can be supplied instead. If support is omitted but a quantile function is available, it is inferred from the quantiles at probabilities zero and one. Without either, a discrete family is assumed to be a count distribution on the nonnegative integers. For a continuous family with a density, a parameter-dependent support is inferred heuristically by probing the log-density. Explicit support is recommended for shifted or truncated counts and for unusual, disconnected, or numerically extreme continuous supports.

• random(par, n, ...): random number generator.

• mean(par, ...): mean function.

• variance(par, ...): variance function.

• skewness(par, ...), kurtosis(par, ...): optional higher-order moment functions.

• initialize: a named list of initialization functions, one for each parameter (e.g., initialize$mu(y, ...)), used to generate starting values.

• valid.response(x): a function that checks whether the response is valid (e.g., numeric, non-factor).

• optimizer(): an optional function defining a custom optimization method for use with gamlss2, see also RS.

If analytical score or Hessian functions are not provided, they are approximated numerically. If cdf() is omitted for a family with an available density and support, it is constructed by numerical integration for a continuous family or by summing the probability mass over the integer support for a discrete family. Analytical CDFs are retained unchanged. Numerical CDFs accept lower.tail, log.p, rel.tol, and abs.tol; continuous CDFs additionally accept subdivisions, while discrete CDFs accept max.terms.

If quantile() is omitted but a CDF and support are available, it is constructed by numerical inversion. Continuous CDFs are inverted by bracketing and bisection. For discrete families, the generated function returns the smallest integer whose CDF is at least the requested probability. Numerical quantiles accept lower.tail, log.p, tol, and maxiter; quantiles generated directly from a count PMF additionally accept max.terms. Analytical quantile functions are retained unchanged.

If mean() or variance() is omitted for a family with a density and support, the missing functions are constructed numerically. Continuous moments use adaptive integration; the variance is integrated as a centered second moment. Discrete moments are obtained by batched summation over the integer support. Numerical moments accept rel.tol and abs.tol; continuous moments additionally accept subdivisions, while discrete moments accept max.terms. A failure to establish a finite moment is reported as an error. Analytical moment functions are retained unchanged.

If random() is omitted but a quantile function is available, random values are generated by inverse transform sampling. Its ... arguments are forwarded to the quantile function. Existing random-number generators are retained unchanged.

Value

For family.gamlss2(), the family object stored in a fitted "gamlss2" model is returned. A "gamlss2.family" object itself is a list describing the response distribution, its parameters, their link functions, and the distribution-specific methods used during fitting, prediction, diagnostics, and simulation.

See Also

gamlss2

Examples

## Not run: ## new family object for the normal distribution
Normal <- function(...) {
  fam <- list(
    "family" = "Normal",
    "names" = c("mu", "sigma"),
    "links" = c("mu" = "identity", "sigma" = "log"),
    "score" = list(
      "mu" = function(par, y, ...) {
        (y - par$mu) / (par$sigma^2)
      },
      "sigma" = function(par, y, ...) {
        -1 + (y - par$mu)^2 / (par$sigma^2)
      }
    ),
    "hessian" = list(
      "mu" = function(par, y, ...) {
        1 / (par$sigma^2)
      },
      "sigma" = function(par, y, ...) {
        rep(2, length(y))
      },
      "mu:sigma" = function(par, y, ...) {
        rep(0, length(y))
      }
    ),
    "log_likelihood" = function(par, y, ...) {
      sum(dnorm(y, par$mu, par$sigma, log = TRUE))
    },
    "pdf" = function(par, y, log = FALSE) {
      dnorm(y, mean = par$mu, sd = par$sigma, log = log)
    },
    "cdf" = function(par, y, ...) {
      pnorm(y, mean = par$mu, sd = par$sigma, ...)
    },
    "random" = function(par, n) {
      rnorm(n, mean = par$mu, sd = par$sigma)
    },
    "quantile" = function(par, p) {
      qnorm(p, mean = par$mu, sd = par$sigma)
    },
    "support" = c(-Inf, Inf),
    "initialize" = list(
      "mu"    = function(y, ...) { (y + mean(y)) / 2 },
      "sigma" = function(y, ...) { rep(sd(y), length(y)) }
    ),
    "mean"      = function(par) par$mu,
    "variance"  = function(par) par$sigma^2,
    "valid.response" = function(x) {
      if(is.factor(x) | is.character(x))
        stop("the response should be numeric!")
      return(TRUE)
    }
  )

  class(fam) <- "gamlss2.family"

  return(fam)
}

## specify the model formula
f <- Ozone ~ s(Wind) | s(Wind)

## estimate model
m1 <- gamlss2(f, data = airquality, family = Normal)

## plot estimated effects
plot(m1, which = "effects")

## plot diagnostics
plot(m1, which = "resid")

## predict distributions
nd <- model.frame(m1)
nd <- nd[order(nd$Wind), ]
d <- predict(m1, newdata = nd)

## predict quantiles
pq <- quantile(d, probs = c(0.05, 0.5, 0.95), elementwise = FALSE)

## visualize
plot(Ozone ~ Wind, data = nd,
  pch = 19, col = rgb(0.1, 0.1, 0.1, alpha = 0.3),
  ylim = range(nd$Ozone, pq))
matlines(nd$Wind, pq, lwd = 2,
  lty = 1, col = 4)

## another example that defines only the density
## function; in this case all derivatives are
## approximated numerically. For residual diagnostics,
## the $cdf() and $quantile() functions are needed as well
Gamma <- function(...) {
  fam <- list(
    "family" = "Gamma",
    "names" = c("mu", "sigma"),
    "links" = c("mu" = "log", "sigma" = "log"),
    "pdf" = function(par, y, log = FALSE, ...) {
      shape <- par$sigma
      scale <- par$mu/par$sigma
      dgamma(y, shape = shape, scale = scale, log = log)
    },
    "cdf" = function(par, y, lower.tail = TRUE, log.p = FALSE) {
      shape <- par$sigma
      scale <- par$mu/par$sigma
      pgamma(y, shape = shape, scale = scale,
        lower.tail = lower.tail, log.p = log.p)
    },
    "quantile" = function(par, p, lower.tail = TRUE, log.p = FALSE) {
      shape <- par$sigma
      scale <- par$mu/par$sigma
       qgamma(p, shape = shape, scale = scale,
         lower.tail = lower.tail, log.p = log.p)
    }
  )

  class(fam) <- "gamlss2.family"

  return(fam)
}

## model formula
f <- Ozone ~ s(Temp) + s(Wind) | .

## estimate model
m2 <- gamlss2(f, data = airquality, family = Gamma)

## visualize estimated effects
plot(m2, which = "effects")

## diagnostics, needs the $cdf() and $quantile() function!
plot(m2, which = "resid")

## End(Not run)

Description

Constructor for estimating penalized model terms using the glmnet package.

Usage

gnet(formula, ...)

Arguments

Details

gnet() constructs a model.matrix from the formula, removes its intercept column, and fits the term using glmnet within backfitting. The penalty lambda minimizes the selected criterion. For criterion = "gaic", the multiplier K is set in gamlss2_control (default: 2).

Value

A special model-term object used internally by gamlss2.

See Also

la, gamlss2, specials.

Examples

## Not run: ## transform to categorical
air <- subset(airquality, !is.na(Ozone))

air$Tc <- cut(air$Temp, breaks = seq(50, 100, by = 10),
  include.lowest = TRUE)
air$Wc <- cut(air$Wind, breaks = seq(2, 22, by = 4),
  include.lowest = TRUE)

## model formula, for each parameter we use
## the same gnet() model term
f <- Ozone ~ gnet(~Tc + Wc) | .

## estimate model
m <- gamlss2(f, data = air, family = NO)

## summary with edf
summary(m)

## plot coefficient paths
sp <- specials(m, parameter = "mu", terms = "gnet(~Tc+Wc)", elements = "model")
print(class(sp))
plot(sp)

## End(Not run)

Description

This dataset contains daily traffic counts near Sonnenberg in the Harz region of Germany. It covers nearly three years, from 2021-01-01 to 2023-11-30, and combines traffic counts with weather covariates. The data are useful for illustrating seasonal count regression models, for example with cyclic smoothers over the day of the year.

Usage

data("HarzTraffic", package = "gamlss2")

Format

A data frame containing 1057 observations on 16 variables.

date

Date, the date of the record.

yday

Integer, the day of the year extracted from date. This is convenient for modeling within-year seasonality.

bikes

Integer, the number of motorcycles on that day.

cars

Integer, the number of cars on that day.

trucks

Integer, the number of trucks on that day.

others

Integer, the number of other vehicles on that day.

tempmin

Numeric, minimum temperature in $^{\circ}C$.

tempmax

Numeric, maximum temperature in $^{\circ}C$.

temp

Numeric, mean temperature in $^{\circ}C$.

humidity

Numeric, mean relative humidity in percent.

tempdew

Numeric, average dewpoint temperature in $^{\circ}C$.

cloudiness

Numeric, average cloud cover in percent.

rain

Numeric, amount of precipitation in mm (snow and rain).

sunshine

Numeric, sunshine duration in minutes.

wind

Numeric, mean wind speed in m/s.

windmax

Numeric, maximum wind speed in m/s.

Details

The traffic variables record daily counts by vehicle type. The weather variables summarize daily conditions measured at a nearby station. In package examples, bikes is often modeled as a seasonal response because motorcycle traffic shows a pronounced annual cycle and substantial variation across weather conditions.

Source

Weather Data:

Data Source:

Deutscher Wetterdienst (DWD), Climate Data Center (CDC).

Licence:

CC BY 4.0

URL:

https://opendata.dwd.de/climate_environment/CDC/

Station:

Wernigerode (5490; Sachsen-Anhalt)

Position:

10.7686/51.8454/233 (lon, lat, alt, EPSG 4326)

Traffic Data:

Data Source:

Bundesanstalt für Strassenwesen (BASt)

Licence:

CC BY 4.0

URL:

https://www.bast.de

Examples

## seasonal variation of motorcycle counts at Sonnenberg/Harz
data("HarzTraffic", package = "gamlss2")
plot(bikes ~ yday, data = HarzTraffic)

## count distribution
barplot(table(HarzTraffic$bikes))

## negative binomial seasonal model using cyclic splines
m <- gamlss2(bikes ~ s(yday, bs = "cc") | s(yday, bs = "cc"),
  data = HarzTraffic, family = NBI)

## visualize effects
plot(m)

## residual diagnostics
plot(m, which = "resid")

## fitted parameters for each day of the year
nd <- data.frame(yday = 1:365)

## corresponding quantiles
p <- quantile(m, newdata = nd, probs = c(0.05, 0.5, 0.95))

## visualization
plot(bikes ~ yday, data = HarzTraffic, pch = 19, col = gray(0.1, alpha = 0.3))
matlines(nd$yday, p, lty = c(2, 1, 2), lwd = 2, col = 4)

Description

Constructs a convolutional neural-network image model term for use with gamlss2. The term learns convolutional filters, applies a ReLU activation and spatial-pyramid pooling, and fits a penalized linear output layer to the current working response.

Usage

im(x, dim = NULL, ...)

Arguments

Details

The image term is intended for images stored as one observation per row. For a flattened matrix, use dim to specify the image dimensions; for an array or image list, the dimensions are checked and retained. The same image representation and dimensions must be supplied when predicting on new data.

The model applies a same-padded convolution, a rectified-linear activation, and spatial-pyramid pooling. The pooled features are mapped to the current working response by a penalized linear output layer. Convolutional parameters are initialized on the first update and warm-started on later backfitting updates. The fitted contribution is centered before it is added to the distributional predictor.

This term is a CNN-based alternative to si, which represents images through a penalized tensor-product coefficient surface. The fitted model can be inspected with plot.gamlss2 and used for new images with predict.gamlss2.

Value

An object of classes "special" and "im". During fitting, the corresponding special_fit.im() method returns the fitted CNN, centered fitted values, transfer information for warm starts, and an effective degrees-of-freedom estimate.

See Also

si, special_terms, gamlss2

Examples

## Not run: 
set.seed(1328)
n <- 120
nr <- 16
nc <- 16

## images with a three-by-three pattern at random locations
images <- array(rnorm(n * nr * nc, sd = 0.25), c(n, nr, nc))
locations <- data.frame(
  row = sample(2:(nr - 2), n, replace = TRUE),
  col = sample(2:(nc - 2), n, replace = TRUE)
)
amplitude <- runif(n, 0.5, 1.5)
for(i in seq_len(n)) {
  images[i, locations$row[i]:(locations$row[i] + 2),
    locations$col[i]:(locations$col[i] + 2)] <-
    images[i, locations$row[i]:(locations$row[i] + 2),
      locations$col[i]:(locations$col[i] + 2)] + amplitude[i]
}

## known spatial response surface
truth <- outer(seq_len(nr), seq_len(nc), function(row, col)
  exp(-((row - 5) / 3)^2 - ((col - 11) / 4)^2))
truth <- truth / max(truth)
d <- data.frame(y = vapply(seq_len(n), function(i) {
  amplitude[i] * mean(truth[locations$row[i]:(locations$row[i] + 2),
    locations$col[i]:(locations$col[i] + 2)])
}, 0) + rnorm(n, sd = 0.2))
image(images[1, , ], col = hcl.colors(20, "Inferno"),
  main = "Example image")
d$image <- I(matrix(images, nrow = n))

m <- gamlss2(
  y ~ im(image, dim = c(nr, nc), filters = 3, epochs = 15,
    warm_epochs = 2, seed = 1904) | 1,
  data = d,
  family = NO
)

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

## predicted response map for a moving image pattern
map.grid <- expand.grid(row = 2:(nr - 2), col = 2:(nc - 2))
map.images <- array(0, c(nrow(map.grid), nr, nc))
for(i in seq_len(nrow(map.grid))) {
  map.images[i, map.grid$row[i]:(map.grid$row[i] + 2),
    map.grid$col[i]:(map.grid$col[i] + 2)] <- 1
}
map.data <- data.frame(y = rep(0, nrow(map.grid)))
map.data$image <- I(matrix(map.images, nrow = nrow(map.grid)))
map.pred <- predict(m, parameter = "mu", newdata = map.data)
true.map <- vapply(seq_len(nrow(map.grid)), function(i) {
  sum(truth[map.grid$row[i]:(map.grid$row[i] + 2),
    map.grid$col[i]:(map.grid$col[i] + 2)]) / 9
}, numeric(1))
true.map <- matrix(true.map,
  nrow = length(unique(map.grid$row)))
pred.map <- matrix(map.pred,
  nrow = length(unique(map.grid$row)))

op <- par(mfrow = c(1, 2))
image(unique(map.grid$row), unique(map.grid$col), true.map,
  col = hcl.colors(20, "Inferno"),
  xlab = "image row", ylab = "image column",
  main = "True response map")
image(unique(map.grid$row), unique(map.grid$col), pred.map,
  col = hcl.colors(20, "Inferno"),
  zlim = range(true.map, pred.map),
  xlab = "image row", ylab = "image column",
  main = "Predicted response map")
par(op)

## End(Not run)

Description

Fit a GAMLSS model by joint Laplace approximate REML.

Usage

## optimizer function
JR(x, y, specials, family, offsets,
  weights, start, xterms, sterms, control)

## refit function
jr(object, trace = TRUE, ...)

Arguments

Details

Use optimizer = JR to fit a model or jr(fit) to continue from an existing fit. The direct solver jointly estimates coefficients and smoothing parameters from cached designs and quadratic penalties. It requires a family density, a predictor-to-parameter map, and an identifiable finite joint mode. Selected penalties and fixed distributional parameters are unsupported. Duplicate unpenalized directions in overlapping smooths are removed from the joint solve. The stored coefficients and covariance use the minimum-norm representation in the original coefficient space, and effective degrees of freedom exclude these redundant directions. If joint REML fails, the input fit is retained and reml.failure records the reason. Check isTRUE(fit$jr$converged) before using joint uncertainty estimates. Power-exponential families may have a non-smooth REML profile when a fitted location coincides with an observation. Such fits retain RS/CG and report that smoothing-parameter curvature could not be verified. Fully penalized smooths that start nearly suppressed are also checked at less penalized starting values, to avoid stopping in a flat profile tail. If the outer optimizer reports false convergence, a stricter independent profile Newton check may verify the fit. In that case fit$jr$outer.verified is true; the original optimizer status is retained in fit$jr$outer. Iteration and evaluation limits still retain the initial fit.

The following controls can be supplied through gamlss2_control or as named arguments to jr():

• jr.inner: Function, RS or CG, used to obtain the initial fit. Default is RS.

• jr.outer.maxit: Integer, the maximum number of outer smoothing parameter iterations. Default is 60.

• jr.outer.maxeval: Integer, the maximum number of outer objective evaluations. Default is 150.

• jr.eps: Positive number, the relative coefficient step tolerance in the joint mode solve. Default is 1e-8. The solver may refine further to stabilize power-exponential curvature.

• jr.maxit: Positive integer, the maximum number of coefficient iterations per objective evaluation. Default is 100.

The initial RS/CG fit uses eps and maxit, as documented in RS. Fully and partially fixed smoothing parameters and unpenalized fx = TRUE smooths are supported. At extreme smoothing parameters, penalty roots and an equilibrated solve preserve accuracy in the unpenalized coefficient directions.

For local REML updates with RS/CG and other smoothing parameter selection methods, see gamlss2_smoothing.

Value

A fitted gamlss2 model with joint REML results in fit$jr.

See Also

RS, gamlss2_control, gamlss2_smoothing, vcov.gamlss2

Examples

## Not run: 
m <- gamlss2(dist ~ s(speed), data = cars, family = NO)
m <- jr(m)
vcov(m, unconditional = TRUE)

## End(Not run)

Description

The two-parameter Kumaraswamy family for continuous responses strictly between 0 and 1, together with a shifted logarithmic link specification.

Usage

## the Kumaraswamy family
Kumaraswamy(a.link = shiftlog, b.link = shiftlog, ...)
KS(a.link = shiftlog, b.link = shiftlog, ...)

## the log(mu - shift) link specification
shiftlog(shift = 1)

Arguments

Details

The Kumaraswamy distribution is continuous on $(0, 1)$ with strictly positive shape parameters $a$ and $b$. Its probability density function is

$f(y; a, b) = a b y^{a - 1} (1 - y^a)^{b - 1},$

and its cumulative distribution function is

$F(y; a, b) = 1 - (1 - y^a)^b.$

Consequently, the quantile function is

$F^{-1}(p; a, b) = \left(1 - (1 - p)^{1/b}\right)^{1/a}.$

For a parameter $\mu$, shiftlog(shift = s) defines the link and inverse link

$g(\mu) = \log(\mu - s), \qquad g^{-1}(\eta) = \exp(\eta) + s,$

with domain $\mu > s$. The default shift = 1 therefore restricts both shape parameters to be greater than 1. In that case the distribution has the unique interior mode

$\left(\frac{a - 1}{a b - 1}\right)^{1/a}.$

When alternative links allow $a \leq 1$ or $b \leq 1$, the family returns NA for the interior mode.

Value

Kumaraswamy() and KS() return the same "gamlss2.family" object. It provides probability density, cumulative distribution, quantile, and random-generation functions, as well as the mean, variance, and interior mode. Its random generator follows the family convention random(par, n).

shiftlog() returns a link specification of class "link-glm".

References

Kumaraswamy P (1980). “A Generalized Probability Density Function for Double-Bounded Random Processes”. Journal of Hydrology, 46(1–2), 79–88. doi:10.1016/0022-1694(80)90036-0

See Also

gamlss2.family, gamlss2

Examples

## create family object with different link specifications
fam <- Kumaraswamy(a.link = shiftlog, b.link = "log")

## simulate data
set.seed(1328)
n <- 1000
d <- data.frame("x" = sort(runif(n, -pi, pi)))

## true parameters
par <- data.frame(
  "a" = exp(1.2 + sin(d$x)) + 1,
  "b" = 1
)

## sample one response for every parameter row
d$y <- fam$random(par, 1)

## estimate model using the Kumaraswamy family
m <- gamlss2(y ~ s(x), data = d, family = fam)

## plot estimated effect
plot(m)

## plot residual diagnostics
plot(m, which = "resid")

## predict quantiles
p <- quantile(m, probs = seq(0.05, 0.95, by = 0.05))

## plot
plot(d, pch = 19, col = adjustcolor(1, 0.3))
matlines(d$x, p, lty = 1, col = 4, lwd = 2)