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, ...,
  level = NULL, interval.cache = NULL)

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

level NULL, or a numeric vector of confidence levels for Wald bands, with every value strictly between zero and one.

interval.cache Optional information factor obtained from the “interval.cache” attribute of a previous Wald prediction for the same, unchanged model. Reusing it avoids rebuilding the joint information matrix for subsequent predictions.

… 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

Predictions can be computed for the original observations or for newdata. Depending on type, results are returned on the link scale, as distribution parameters, as individual model terms, or as the fitted response mean.

If se.fit = TRUE, standard errors are obtained from coefficient draws. Existing draws stored in the fitted object are used when available; otherwise, draws are generated from the approximate distribution of the fitted coefficients.

If level is specified, approximate pointwise confidence intervals and standard errors are computed analytically, without simulation. Intervals on transformed parameter or response scales are transformed as needed and may be asymmetric. These are confidence intervals for the fitted model, not prediction intervals for future observations, and they do not include smoothing-parameter or model-selection uncertainty.

Wald intervals currently support maximum-likelihood fits with linear terms and standard mgcv smooths. They are not available for Bayesian fits, custom special terms, or automatic term-selection penalties. Wald intervals cannot be combined with simulation controls such as FUN, R, seed, or burnin. The optional interval.cache can be reused for repeated predictions from the same unchanged model.

For type = “terms”, term names are partially matched by default. Use nogrep = TRUE in … for exact matching.

See also prodist.gamlss2 for creating a full distributions3 object for moments, probabilities, quantiles, and random numbers.

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. Multi-parameter models return a list of term matrices. For se.fit = TRUE, the returned object contains fitted values and corresponding simulation-based standard errors.

When level is specified, a list with elements fit, se.fit, lower, and upper is returned. For one confidence level, lower and upper follow the structure of the ordinary point prediction, including model, terms, and drop selection. For multiple confidence levels, each is a list named by confidence percentage (for example, “80%” and “95%”), whose elements have that same prediction structure. The fit element is the ordinary, non-simulated prediction and se.fit is shared by all intervals. Attributes record the confidence levels, method, and reusable interval.cache. The cache retains a reference to the fitted model; reuse it only with that unchanged model.

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.4242 eps = 0.129640     
GAMLSS-RS iteration  2: Global Deviance = 405.2294 eps = 0.000480     
GAMLSS-RS iteration  3: Global Deviance = 405.2069 eps = 0.000055     
GAMLSS-RS iteration  4: Global Deviance = 405.1868 eps = 0.000049     
GAMLSS-RS iteration  5: Global Deviance = 405.17 eps = 0.000041     
GAMLSS-RS iteration  6: Global Deviance = 405.1562 eps = 0.000033     
GAMLSS-RS iteration  7: Global Deviance = 405.1452 eps = 0.000027     
GAMLSS-RS iteration  8: Global Deviance = 405.1365 eps = 0.000021     
GAMLSS-RS iteration  9: Global Deviance = 405.1297 eps = 0.000016     
GAMLSS-RS iteration 10: Global Deviance = 405.1244 eps = 0.000012     
GAMLSS-RS iteration 11: Global Deviance = 405.1206 eps = 0.000009

## 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.02272  9.956954
2 59.18651 18.423672
3 97.22282 31.930725

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

        mu    sigma
1 23.02272 2.298271
2 59.18651 2.913636
3 97.22282 3.463569

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

       1        2        3 
23.02272 59.18651 97.22282

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

        1         2         3 
 9.956954 18.423672 31.930725

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

      sigma
1  9.956954
2 18.423672
3 31.930725

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

  (Intercept)   s(speed)
1    2.632278 -0.3340064
2    2.632278  0.2813587
3    2.632278  0.8312911

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

    s(speed)
1 -0.3340064
2  0.2813587
3  0.8312911

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

       10%      50%       90%
1 10.26237 23.02272  35.78307
2 35.57562 59.18651  82.79740
3 56.30195 97.22282 138.14369

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

    mu.fit    mu.se sigma.fit  sigma.se
1 23.02272 1.907417  10.07766  1.583013
2 59.18651 3.447267  18.67225  3.095756
3 97.22282 9.561225  34.39733 13.684037

## fast pointwise confidence bands, without sampling
bands <- predict(m, newdata = nd, level = c(0.8, 0.95))
bands$lower[["80%"]]

        mu     sigma
1 20.45409  8.094526
2 54.72243 15.273264
3 85.75127 19.109395

bands$upper[["95%"]]

         mu    sigma
1  26.95111 13.66691
2  66.01373 24.54339
3 114.76705 70.01636