gamlss2

Extracting Fitted or Predicted Parameters or Terms from gamlss2 Models

Description

Methods for gamlss2 model objects that extract fitted values or predictions for distribution parameters, additive predictors, individual terms, or the mean of the response distribution.

Usage

## S3 method for class 'gamlss2'
predict(object, model = NULL, newdata = NULL,
  type = c("parameter", "link", "response", "terms"), terms = NULL,
  se.fit = FALSE, drop = TRUE, ...)

Arguments

object A model object of class gamlss2.

model Character. Which distribution parameter(s) should be predicted? This can be one or more of “mu”, “sigma”, etc. By default, predictions are returned for all model parts.

newdata Data frame. Optionally, a new data frame in which to look for variables with which to predict. If omitted, the original observations are used.

type Character. Which type of prediction should be computed? Use “link” for additive predictors on the link scale, “parameter” for distribution parameters on their natural scale, “terms” for individual additive terms, and “response” for the mean of the fitted response distribution.

terms Character. Which of the terms in the additive predictor(s) should be included? By default all terms are included.

se.fit Logical. Should standard errors for the predictions be included? Standard errors are computed by simulating from the approximate multivariate normal distribution of the maximum likelihood estimates. The number of simulations is controlled by the argument R, which defaults to R = 200, and can be passed via ….

drop Logical. Should the predictions be simplified to a vector if possible (TRUE) or always returned as a data frame (FALSE)?

… Further control arguments. Supported aliases are what and parameter for model. Additional arguments are used for simulation-based summaries, including FUN, R, seed, burnin, and nogrep.

Details

Prediction for a gamlss2 object proceeds in stages. First, the required model terms are evaluated for the original data or for newdata. These terms are then added to form the full additive predictor for each distribution parameter. Applying the inverse link function yields the predicted distribution parameters. If requested, the mean of the corresponding response distribution is computed in a final step.

If se.fit = TRUE, or if a summary function is supplied via FUN, predictions are based on simulation draws from the approximate multivariate normal distribution of the fitted coefficients, unless the object already contains stored draws in samples. For objects of class “bamlss2”, stored posterior draws are used directly.

For type = “terms”, the terms argument is matched against linear and special-term labels. By default, partial matching is used; set nogrep = TRUE in … to require exact matching.

See also prodist.gamlss2 for setting up a full distributions3 object from which moments, probabilities, quantiles, or random numbers can be obtained.

Value

If drop = FALSE, a data frame is returned. If drop = TRUE (the default), the result may be simplified to a numeric vector when possible. For type = “terms”, a data frame is returned with one column per selected term. For se.fit = TRUE, the returned object contains fitted values and corresponding simulation-based standard errors.

Examples

library("gamlss2")

## fit heteroscedastic normal GAMLSS model
## stopping distance (ft) explained by speed (mph)
data("cars", package = "datasets")
m <- gamlss2(dist ~ s(speed) | s(speed), data = cars, family = NO)

GAMLSS-RS iteration  1: Global Deviance = 405.737 eps = 0.128969     
GAMLSS-RS iteration  2: Global Deviance = 405.7368 eps = 0.000000

## new data for predictions
nd <- data.frame(speed = c(10, 20, 30))

## default: model parameter(s) for all model parts
predict(m, newdata = nd)

        mu    sigma
1 23.05774 10.06075
2 58.98303 18.52920
3 96.14450 34.11606

## additive predictors on the link scale
predict(m, newdata = nd, type = "link")

        mu    sigma
1 23.05774 2.308642
2 58.98303 2.919348
3 96.14450 3.529768

## mean of the response distribution
predict(m, newdata = nd, type = "response")

       1        2        3 
23.05774 58.98303 96.14450

## model parameter(s)
predict(m, newdata = nd, model = "sigma")

       1        2        3 
10.06075 18.52920 34.11606

predict(m, newdata = nd, model = "sigma", drop = FALSE)

## individual terms in additive predictor(s)
predict(m, newdata = nd, type = "terms", model = "sigma")

  (Intercept)   s(speed)
1     2.63843 -0.3297884
2     2.63843  0.2809179
3     2.63843  0.8913382

predict(m, newdata = nd, type = "terms", model = "sigma", terms = "s(speed)")

    s(speed)
1 -0.3297884
2  0.2809179
3  0.8913382

## predict quantiles
quantile(m, newdata = nd, probs = c(0.1, 0.5, 0.9))

       10%      50%       90%
1 10.16437 23.05774  35.95111
2 35.23691 58.98303  82.72916
3 52.42302 96.14450 139.86599

## standard errors
predict(m, newdata = nd, se.fit = TRUE, R = 200)

    mu.fit    mu.se sigma.fit  sigma.se
1 23.05774 1.908156  10.18556  1.619731
2 58.98303 3.459757  18.77165  3.062022
3 96.14450 8.981025  36.24582 12.975020