Predictions from gamlss2 Models

Description

Predict response distributions, parameters, additive predictors, terms, or response means from a gamlss2 model.

Usage

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

Arguments

object	A gamlss2 model.
parameter	Distribution parameter name(s) or index(es), such as “mu” or “sigma”; defaults to all parameters.
newdata	Data frame for predictions; defaults to the fitted data.
type	“distribution” (default) returns distributions3 objects; “parameter” returns parameters on their natural scale; “link” returns additive predictors; “response” returns response means; “terms” returns individual additive terms. Specifying parameter without type selects “parameter”.
terms	Additive terms to include; defaults to all. Names are partially matched unless nogrep = TRUE is passed in ….
se.fit	For ordinary predictions, TRUE computes standard errors from coefficient draws; R in … sets the number of draws (default 200).
drop	Simplify the result when possible.
level	Numeric confidence level(s) between zero and one; supplying levels selects pointwise Wald intervals.
interval	“none” (default) returns ordinary predictions; “wald” returns analytic pointwise confidence intervals, at the 95% confidence level by default.
unconditional	If TRUE, include approximate smoothing-parameter uncertainty in Wald intervals and coefficient draws.
…	Further controls, including FUN, R, seed, and burnin for simulation.

Details

Coefficient draws stored in the fit are used when available; otherwise, draws use the full joint coefficient covariance. They approximate the coefficient sampling distribution and may give unreliable tails after nonlinear transformations. Simulation stops if a marginal 95% Gaussian interval for a log smoothing parameter spans a factor of 100 on the smoothing-parameter scale.

Wald intervals are approximate pointwise confidence intervals for fitted values, not prediction intervals for future observations. They are computed without simulation and may be asymmetric on transformed scales. They condition on fitted smoothing parameters by default and exclude model-selection uncertainty; see vcov.gamlss2.

Wald intervals support maximum-likelihood fits with linear terms and coefficient-linear mgcv smooths, but not nonlinear custom special terms. They cannot be combined with FUN, R, seed, or burnin. An indefinite joint information matrix or non-estimable coefficient causes an error rather than a working-information approximation.

Value

Ordinary predictions return a distribution object or values in the selected form. Term predictions are data frames or, for multiple parameters, lists of matrices. With se.fit = TRUE, fitted values and standard errors are returned. With FUN = identity, simulated distribution predictions form a list of distribution objects.

Wald predictions return fit, se.fit, lower, and upper. For multiple levels, lower and upper are lists named by confidence percentage (for example, “80%” and “95%”). fit and each bound retain the structure selected by parameter, terms, and drop; se.fit is shared across levels.

See Also

predict, marginal_predict, prodist.gamlss2, quantile.gamlss2

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.0352 eps = 0.130476     
GAMLSS-RS iteration  2: Global Deviance = 405.6152 eps = 0.001432     
GAMLSS-RS iteration  3: Global Deviance = 405.657 eps = 0.000102     
GAMLSS-RS iteration  4: Global Deviance = 405.6596 eps = 0.000006     

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

## default: a "distribution" object
d <- predict(m, newdata = nd)
print(d)

                                      1                                       2 
"GAMLSS2 NO(mu = 23.06, sigma = 10.03)" "GAMLSS2 NO(mu = 58.89, sigma = 18.54)" 
                                      3 
"GAMLSS2 NO(mu = 96.20, sigma = 34.27)" 

## can be used with distributions3 methods
pdf(d, 20)

           1            2            3 
0.0379452246 0.0023848152 0.0009821864 

cdf(d, 20)

        1         2         3 
0.3800749 0.0179784 0.0130812 

quantile(d, 0.8)

        1         2         3 
 31.50914  74.50001 125.04054 

mean(d)

       1        2        3 
23.06351 58.89341 96.20125 

variance(d)

        1         2         3 
 100.7002  343.8611 1174.1825 

random(d, 3)

       r_1      r_2       r_3
1 34.95525 12.84773  18.33270
2 67.67621 70.83299  63.16037
3 57.31322 25.59377 169.09065

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

        mu    sigma
1 23.06351 2.306074
2 58.89341 2.920119
3 96.20125 3.534164

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

       1        2        3 
23.06351 58.89341 96.20125 

## same as
mean(d)

       1        2        3 
23.06351 58.89341 96.20125 

## all model parameters
predict(m, newdata = nd, type = "parameter")

        mu    sigma
1 23.06351 10.03495
2 58.89341 18.54349
3 96.20125 34.26635

## parameter selects parameter predictions if type is omitted
predict(m, newdata = nd, parameter = "sigma")

       1        2        3 
10.03495 18.54349 34.26635 

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

     sigma
1 10.03495
2 18.54349
3 34.26635

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

  (Intercept)   s(speed)
1    2.637658 -0.3315843
2    2.637658  0.2824607
3    2.637658  0.8965055

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

    s(speed)
1 -0.3315843
2  0.2824607
3  0.8965055

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

       10%      50%       90%
1 10.20321 23.06351  35.92381
2 35.12897 58.89341  82.65786
3 52.28716 96.20125 140.11535

## same as
quantile(d, probs = c(0.1, 0.5, 0.9), elementwise = FALSE)

     q_0.1    q_0.5     q_0.9
1 10.20321 23.06351  35.92381
2 35.12897 58.89341  82.65786
3 52.28716 96.20125 140.11535

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

    mu.fit    mu.se sigma.fit  sigma.se
1 23.06351 1.917939  10.15932  1.614735
2 58.89341 3.470156  18.78608  3.064083
3 96.20125 9.341785  36.40538 13.041566

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

        mu    sigma
1 20.47978  8.16792
2 54.40545 15.36917
3 84.96474 21.78597

bands$upper[["95%"]]

         mu    sigma
1  27.01499 13.74822
2  65.75715 24.71144
3 113.38602 68.49845

## prediction uncertainty using multiple draws
d <- predict(m, newdata = cars, R = 1000, FUN = identity)

## predict 1000 quantiles
q <- sapply(d, quantile, probs = 0.95)

i95 <- function(x) {
  setNames(quantile(x, c(0.025, 0.5, 0.975), names = FALSE),
    c("lower", "fit", "upper"))
}

## quantiles of quantile samples
q <- t(apply(q, 1, i95))

## plot quantile with uncertainty bands
plot(dist ~ speed, data = cars,
  ylim = range(cars$dist, q))
matlines(cars$speed, q,
  lty = c(2, 1, 2), col = 4, lwd = 2)