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

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 = 909.151 eps = 0.173776     
GAMLSS-RS iteration  2: Global Deviance = 901.5501 eps = 0.008360     
GAMLSS-RS iteration  3: Global Deviance = 901.1446 eps = 0.000449     
GAMLSS-RS iteration  4: Global Deviance = 901.1033 eps = 0.000045     
GAMLSS-RS iteration  5: Global Deviance = 901.0673 eps = 0.000039     
GAMLSS-RS iteration  6: Global Deviance = 901.036 eps = 0.000034     
GAMLSS-RS iteration  7: Global Deviance = 901.0085 eps = 0.000030     
GAMLSS-RS iteration  8: Global Deviance = 900.9843 eps = 0.000026     
GAMLSS-RS iteration  9: Global Deviance = 900.963 eps = 0.000023     
GAMLSS-RS iteration 10: Global Deviance = 900.9442 eps = 0.000020     
GAMLSS-RS iteration 11: Global Deviance = 900.9275 eps = 0.000018     
GAMLSS-RS iteration 12: Global Deviance = 900.9126 eps = 0.000016     
GAMLSS-RS iteration 13: Global Deviance = 900.8994 eps = 0.000014     
GAMLSS-RS iteration 14: Global Deviance = 900.8875 eps = 0.000013     
GAMLSS-RS iteration 15: Global Deviance = 900.8768 eps = 0.000011     
GAMLSS-RS iteration 16: Global Deviance = 900.8672 eps = 0.000010     
GAMLSS-RS iteration 17: Global Deviance = 900.8585 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.52962    0.02099  168.19   <2e-16 ***
sigma.(Intercept) -0.98746    0.06363  -15.52   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---
Smooth terms:
    mu.s(Temp) mu.s(Wind) sigma.s(Temp)
edf     7.9120     3.1873        7.5602
*--------
n = 116 df =  20.66 res.df =  95.34
Deviance = 900.8585 Null Dev. Red. = 16.82%
AIC = 942.1774 elapsed =  0.23sec

## 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.60531 0.3717324
2 18.98143 0.4718580
3 13.92848 0.3693678
4 14.93042 0.5548871
6 19.23528 0.3498466
7 23.06716 0.3551156

## 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:
    mu.s(yday) sigma.s(yday)
edf     6.3823        6.2095
*--------
n = 1057 df =  14.59 res.df =  1042.41
Deviance = 10150.8818 Null Dev. Red. = 12.62%
AIC = 10180.0654 elapsed =  0.08sec

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