gamlss2

Generalized Additive Models for Location Scale and Shape

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

formula A GAM-type formula or Formula describing the additive predictors for the distribution parameters. All smooth terms from mgcv are supported, see also formula.gam. For gamlss2.list(), formula is a list of formulas.

data A data frame, 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 defining the response distribution and the link functions of its 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 Optional offset terms to be included in the linear predictors during fitting. If a single numeric vector is supplied, it is assigned to the first parameter of the distribution. If offsets are needed for several parameters, provide a data frame or list in the same order as the parameter names in the family object.

start Optional starting values for the estimation algorithm, see gamlss2_start.

knots An optional list of user-specified knots, see smoothCon.

control A list of control arguments, see gamlss2_control.

… Arguments passed to gamlss2_control.

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.

Examples

library("gamlss2")

## 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)

GAMLSS-RS iteration  1: Global Deviance = 934.9808 eps = 0.150303     
GAMLSS-RS iteration  2: Global Deviance = 927.3127 eps = 0.008201     
GAMLSS-RS iteration  3: Global Deviance = 924.8517 eps = 0.002653     
GAMLSS-RS iteration  4: Global Deviance = 925.5146 eps = 0.000716     
GAMLSS-RS iteration  5: Global Deviance = 926.9917 eps = 0.001595     
GAMLSS-RS iteration  6: Global Deviance = 928.0621 eps = 0.001154     
GAMLSS-RS iteration  7: Global Deviance = 928.6793 eps = 0.000665     
GAMLSS-RS iteration  8: Global Deviance = 928.9916 eps = 0.000336     
GAMLSS-RS iteration  9: Global Deviance = 929.1382 eps = 0.000157     
GAMLSS-RS iteration 10: Global Deviance = 929.2117 eps = 0.000079     
GAMLSS-RS iteration 11: Global Deviance = 929.2443 eps = 0.000035     
GAMLSS-RS iteration 12: Global Deviance = 929.2588 eps = 0.000015     
GAMLSS-RS iteration 13: Global Deviance = 929.2643 eps = 0.000005

## model summary
summary(m1)

Call:
gamlss2(formula = f, data = air, family = GA)
---
Family: GA(log(mu), log(sigma)) 
*--------
Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
mu.(Intercept)     3.55382    0.02682  132.49   <2e-16 ***
sigma.(Intercept) -0.87025    0.06331  -13.75   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
                 edf
mu.s(Temp)    4.4762
mu.s(Wind)    1.0718
sigma.s(Temp) 3.9513
*--------
n = 116 df =  11.5 res.df =  104.5
Deviance = 929.2643 Null Dev. Red. = 14.2%
AIC = 952.263 elapsed =  0.15sec

## 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))

        mu     sigma
1 22.86957 0.4849417
2 23.62326 0.5090063
3 20.71022 0.5219236
4 16.65926 0.6021705
6 16.22035 0.4917132
7 20.78846 0.5063877

## 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
)

GAMLSS-RS iteration  1: Global Deviance = 10163.0634 eps = 0.148403     
GAMLSS-RS iteration  2: Global Deviance = 10151.0415 eps = 0.001182     
GAMLSS-RS iteration  3: Global Deviance = 10150.8407 eps = 0.000019     
GAMLSS-RS iteration  4: Global Deviance = 10150.7454 eps = 0.000009

## summary
summary(m2)

Call:
gamlss2(formula = bikes ~ s(yday, bs = "cc") | s(yday, bs = "cc"), 
    data = HarzTraffic, family = NBI)
---
Family: NBI(log(mu), log(sigma)) 
*--------
Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
mu.(Intercept)     3.99888    0.03065 130.452   <2e-16 ***
sigma.(Intercept)  0.46596    0.04913   9.485   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
                 edf
mu.s(yday)    7.1494
sigma.s(yday) 5.9745
*--------
n = 1057 df =  15.12 res.df =  1041.88
Deviance = 10150.7454 Null Dev. Red. = 12.62%
AIC = 10180.9933 elapsed =  0.11sec

## 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)