gamlss2

Generalized Additive Models for Location Scale and Shape

Description

Estimation of generalized additive models for location scale and shape (GAMLSS). The model fitting function gamlss2() provides flexible infrastructures to estimate the parameters of a response distribution. The number of distributional parameters is not fixed, see gamlss2.family. Moreover, gamlss2() supports all smooth term constructors from the mgcv package in addition to the classical model terms as provided by gamlss and gamlss.add.

Usage

gamlss2(formula, ...)

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

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

Arguments

formula A GAM-type formula or Formula. All smooth terms of the mgcv package are supported, see also formula.gam. For gamlss.list() formula is a list of formulas.

data A data frame or list or environment containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which gamlss2 is called.

family A gamlss.family or gamlss2.family object used to define distribution and the link functions of the parameters.

subset An optional vector specifying a subset of observations to be used in the fitting process.

na.action NA processing for setting up the model.frame.

weights An optional vector of prior weights to be used in the fitting process. Should be NULL or a numeric vector.

offset This can be used to specify an a priori known components to be included in the linear predictors during fitting. Please note that if only a single numeric vector is provided, the offset will be assigned to the first specified parameter of the distribution. In the case of multiple offsets, a data frame or list must be supplied. Each offset is assigned in the same order as the parameters of the distribution specified in the family object.

start Starting values for estimation algorithms.

control A list of control arguments, see gamlss2_control.

… Arguments passed to gamlss2_control.

Details

The model fitting function gamlss2() provides flexible infrastructures for the estimation of GAMLSS.

Distributional models are specified using family objects, either from the gamlss.dist package or using gamlss2.family objects.
Estimation is carried out through a Newton-Raphson/Fisher scoring algorithm, see function RS. The estimation algorithms can also be exchanged using gamlss2_control. Additionally, if an optimizer is specified by the family object, this optimizer function will be employed for estimation.
The return value is determined by the object returned from the optimizer function, typically an object of class “gamlss2”. Default methods and extractor functions are available for this class. Nevertheless, users have the flexibility to supply their own optimizer function, along with user-specific methods tailored for the returned object.

Value

The return value is determined by the object returned from the optimizer function. By default, the optimization is performed using the RS optimizer function (see gamlss2_control), yielding an object of class “gamlss2”. Default methods and extractor functions are available for this class.

References

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, 507–554. doi:10.1111/j.1467-9876.2005.00510.x

Rigby RA, Stasinopoulos DM, Heller GZ, De Bastiani F (2019). Distributions for Modeling Location, Scale, and Shape: Using GAMLSS in R, Chapman and Hall/CRC. doi:10.1201/9780429298547

Stasinopoulos DM, Rigby RA (2007). “Generalized Additive Models for Location Scale and Shape (GAMLSS) in R.” Journal of Statistical Software, 23(7), 1–46. doi:10.18637/jss.v023.i07

Stasinopoulos DM, Rigby RA, Heller GZ, Voudouris V, De Bastiani F (2017). Flexible Regression and Smoothing: Using GAMLSS in R, Chapman and Hall/CRC. doi:10.1201/b21973

Examples

library("gamlss2")


## load the abdominal circumference data
data("abdom", package = "gamlss.data")

## specify the model Formula
f <- y ~ s(x) | s(x) | s(x) | s(x)

## estimate model
b <- gamlss2(f, data = abdom, family = BCT)

GAMLSS-RS iteration  1: Global Deviance = 4774.5602 eps = 0.534336     
GAMLSS-RS iteration  2: Global Deviance = 4770.2329 eps = 0.000906     
GAMLSS-RS iteration  3: Global Deviance = 4770.1693 eps = 0.000013     
GAMLSS-RS iteration  4: Global Deviance = 4770.1544 eps = 0.000003

## model summary
summary(b)

Call:
gamlss2(formula = f, data = abdom, family = BCT)
---
Family: BCT 
Link function: mu = identity, sigma = log, nu = identity, tau = log
*--------
Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
mu.(Intercept)    226.33309    1.25941  179.71  < 2e-16 ***
sigma.(Intercept)  -2.92261    0.01104 -264.67  < 2e-16 ***
nu.(Intercept)     -0.18026    0.04611   -3.91 0.000103 ***
tau.(Intercept)     2.65014    0.01440  183.98  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
    mu.s(x) sigma.s(x) nu.s(x) tau.s(x)
edf  4.5505     2.5450  1.0015   1.0004
*--------
n = 610 df =  13.1 res.df =  596.9
Deviance = 4770.1544 Null Dev. Red. = 33.39%
AIC = 4796.3494 elapsed =  0.84sec

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

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

## predict parameters
par <- predict(b)

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

## visualize
plot(y ~ x, data = abdom, pch = 19,
  col = rgb(0.1, 0.1, 0.1, alpha = 0.3))
matplot(abdom$x, pq, type = "l", lwd = 2,
  lty = 1, col = 4, add = TRUE)

## use of starting values
m <- gamlss2(f, data = abdom, family = BCT,
  start = c(mu = 200, sigma = 0.1, nu = 0, tau = 10))

GAMLSS-RS iteration  1: Global Deviance = 4816.1239 eps = 0.553570     
GAMLSS-RS iteration  2: Global Deviance = -7383.8761 eps = 2.533157     
GAMLSS-RS iteration  3: Global Deviance = -19583.8761 eps = 1.652248     
GAMLSS-RS iteration  4: Global Deviance = -31783.8761 eps = 0.622961     
GAMLSS-RS iteration  5: Global Deviance = -36347.9794 eps = 0.143598     
GAMLSS-RS iteration  6: Global Deviance = -36347.9794 eps = 0.000000

## fix some parameters
m <- gamlss2(f, data = abdom, family = BCT,
  start = c(mu = 200, sigma = 0.1, nu = 0, tau = 10),
  fixed = c(nu = TRUE, tau = TRUE))

GAMLSS-RS iteration  1: Global Deviance = 5076.3585 eps = 0.529448     
GAMLSS-RS iteration  2: Global Deviance = 5076.3565 eps = 0.000000

## estimated coefficients (intercepts)
coef(m)

   mu.p.(Intercept) sigma.p.(Intercept)    nu.p.(Intercept)   tau.p.(Intercept) 
         226.358484           -2.302585            0.000000            2.302585

## starting values using full predictors
m <- gamlss2(f, data = abdom, family = BCT,
  start = fitted(m))

GAMLSS-RS iteration  1: Global Deviance = 7749.3698 eps = 0.013476     
GAMLSS-RS iteration  2: Global Deviance = 4869.5907 eps = 0.371614     
GAMLSS-RS iteration  3: Global Deviance = 4793.4424 eps = 0.015637     
GAMLSS-RS iteration  4: Global Deviance = 4792.9818 eps = 0.000096     
GAMLSS-RS iteration  5: Global Deviance = 4792.9727 eps = 0.000001

## same with
m <- gamlss2(f, data = abdom, family = BCT,
  start = m)

GAMLSS-RS iteration  1: Global Deviance = 4774.5602 eps = 0.534336     
GAMLSS-RS iteration  2: Global Deviance = 4770.2329 eps = 0.000906     
GAMLSS-RS iteration  3: Global Deviance = 4770.1693 eps = 0.000013     
GAMLSS-RS iteration  4: Global Deviance = 4770.1544 eps = 0.000003