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 gamlss.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.

See Also

RS, gamlss2_control, gamlss2.family, gamlss2_start

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 = 906.1982 eps = 0.176460     
GAMLSS-RS iteration  2: Global Deviance = 899.9294 eps = 0.006917     
GAMLSS-RS iteration  3: Global Deviance = 899.8272 eps = 0.000113     
GAMLSS-RS iteration  4: Global Deviance = 899.7914 eps = 0.000039     
GAMLSS-RS iteration  5: Global Deviance = 899.7641 eps = 0.000030     
GAMLSS-RS iteration  6: Global Deviance = 899.7434 eps = 0.000022     
GAMLSS-RS iteration  7: Global Deviance = 899.7278 eps = 0.000017     
GAMLSS-RS iteration  8: Global Deviance = 899.7161 eps = 0.000013     
GAMLSS-RS iteration  9: Global Deviance = 899.7073 eps = 0.000009     

## model summary
summary(m1)

Call:
gamlss2(formula = f, data = air, family = GA)
---
Family: GA 
Link functions: log, log
*--------
Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
mu.(Intercept)     3.53026    0.02105   167.7   <2e-16 ***
sigma.(Intercept) -0.99253    0.06362   -15.6   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
                 edf
mu.s(Temp)    7.8814
mu.s(Wind)    3.2562
sigma.s(Temp) 7.4967
*--------
n = 116 df =  20.63 res.df =  95.37
Deviance = 899.7073 Null Dev. Red. = 16.93%
AIC = 940.9758 elapsed =  0.27sec

## 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 28.35701 0.3717591
2 18.74042 0.4561988
3 13.83960 0.3632633
4 14.84861 0.5606800
6 19.37133 0.3533476
7 23.09536 0.3606976

## 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)
))
par <- predict(m1, newdata = nd5)
nd5$Probs <- 1 - family(m1)$cdf(par, 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 = 10151.1432 eps = 0.149402     
GAMLSS-RS iteration  2: Global Deviance = 10150.8948 eps = 0.000024     
GAMLSS-RS iteration  3: Global Deviance = 10150.8818 eps = 0.000001     

## summary
summary(m2)

Call:
gamlss2(formula = bikes ~ s(yday, bs = "cc") | s(yday, bs = "cc"), 
    data = HarzTraffic, family = NBI)
---
Family: NBI 
Link functions: log, log
*--------
Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
mu.(Intercept)     3.99347    0.03062  130.43   <2e-16 ***
sigma.(Intercept)  0.46777    0.04416   10.59   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
                 edf
mu.s(yday)    6.3823
sigma.s(yday) 6.2095
*--------
n = 1057 df =  14.59 res.df =  1042.41
Deviance = 10150.8818 Null Dev. Red. = 12.62%
AIC = 10180.0654 elapsed =  0.12sec

## effect plots
plot(m2)

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

## quantiles over the day of the year
nd <- data.frame(yday = 1:365)
pq <- quantile(m2, newdata = nd, 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)