Predict Method for a VGLM fit

Predicted values based on a vector generalized linear model (VGLM) object.

Usage

predictvglm(object, newdata = NULL,
            type = c("link", "response", "terms"),
            se.fit = FALSE, deriv = 0, dispersion = NULL,
            untransform = FALSE,
            type.fitted = NULL, percentiles = NULL, ...)

Arguments

object

Object of class inheriting from "vlm", e.g., vglm.

newdata

An optional data frame in which to look for variables with which to predict. If omitted, the fitted linear predictors are used.

type

The value of this argument can be abbreviated. The type of prediction required. The default is the first one, meaning on the scale of the linear predictors. This should be a \(n \times M\) matrix.

The alternative "response" is on the scale of the response variable, and depending on the family function, this may or may not be the mean. Often this is the fitted value, e.g., fitted(vglmObject) (see fittedvlm). Note that the response is output from the @linkinv slot, where the eta argument is the \(n \times M\) matrix of linear predictors.

The "terms" option returns a matrix giving the fitted values of each term in the model formula on the linear predictor scale. The terms have been centered.

se.fit

logical: return standard errors?

deriv

Non-negative integer. Currently this must be zero. Later, this may be implemented for general values.

dispersion

Dispersion parameter. This may be inputted at this stage, but the default is to use the dispersion parameter of the fitted model.

type.fitted

Some VGAM family functions have an argument by the same name. If so, then one can obtain fitted values by setting type = "response" and choosing a value of type.fitted from what's available. If type.fitted = "quantiles" is available then the percentiles argument can be used to specify what quantile values are requested.

percentiles

Used only if type.fitted = "quantiles" is available and is selected.

untransform

Logical. Reverses any parameter link function. This argument only works if type = "link", se.fit = FALSE, deriv = 0. Setting untransform = TRUE does not work for all VGAM family functions; only ones where there is a one-to-one correspondence between a simple link function and a simple parameter might work.

...

Arguments passed into predictvlm.

Details

Obtains predictions and optionally estimates standard errors of those predictions from a fitted vglm object. By default, each row of the matrix returned can be written as \(\eta_i^T\), comprising of \(M\) components or linear predictors. If there are any offsets, these are included.

This code implements smart prediction (see smartpred).

Value

If se.fit = FALSE, a vector or matrix of predictions. If se.fit = TRUE, a list with components

fitted.values: Predictions
se.fit: Estimated standard errors
df: Degrees of freedom
sigma: The square root of the dispersion parameter (but these are being phased out in the package)

References

Yee, T. W. (2015). Vector Generalized Linear and Additive Models: With an Implementation in R. New York, USA: Springer.

Yee, T. W. and Hastie, T. J. (2003). Reduced-rank vector generalized linear models. Statistical Modelling, 3, 15–41.

Author

Thomas W. Yee

Note

Setting se.fit = TRUE and type = "response" will generate an error.

The arguments type.fitted and percentiles are provided in this function to give more convenience than modifying the extra slot directly.

Warning

This function may change in the future.

Examples

# Illustrates smart prediction
pneumo <- transform(pneumo, let = log(exposure.time))
fit <- vglm(cbind(normal, mild, severe) ~ poly(c(scale(let)), 2),
            propodds, pneumo, trace = TRUE, x.arg = FALSE)
#> Iteration 1: deviance = 4.162527
#> Iteration 2: deviance = 3.951629
#> Iteration 3: deviance = 3.946558
#> Iteration 4: deviance = 3.94655
#> Iteration 5: deviance = 3.94655
class(fit)
#> [1] "vglm"
#> attr(,"package")
#> [1] "VGAM"

(q0 <- head(predict(fit)))
#>   logitlink(P[Y>=2]) logitlink(P[Y>=3])
#> 1         -6.6420717          -7.540193
#> 2         -2.7470610          -3.645182
#> 3         -1.6175447          -2.515666
#> 4         -0.9547998          -1.852921
#> 5         -0.4876278          -1.385749
#> 6         -0.1414232          -1.039544
(q1 <- predict(fit, newdata = head(pneumo)))
#>   logitlink(P[Y>=2]) logitlink(P[Y>=3])
#> 1       -5.773141107         -6.6712621
#> 2       -2.059406372         -2.9575274
#> 3       -1.017291252         -1.9154122
#> 4       -0.420223629         -1.3183446
#> 5       -0.009188218         -0.9073092
#> 6        0.287785773         -0.6103352
(q2 <- predict(fit, newdata = head(pneumo)))
#>   logitlink(P[Y>=2]) logitlink(P[Y>=3])
#> 1       -5.773141107         -6.6712621
#> 2       -2.059406372         -2.9575274
#> 3       -1.017291252         -1.9154122
#> 4       -0.420223629         -1.3183446
#> 5       -0.009188218         -0.9073092
#> 6        0.287785773         -0.6103352
all.equal(q0, q1)  # Should be TRUE
#> [1] "Mean relative difference: 0.2354654"
all.equal(q1, q2)  # Should be TRUE
#> [1] TRUE

head(predict(fit))
#>   logitlink(P[Y>=2]) logitlink(P[Y>=3])
#> 1         -6.6420717          -7.540193
#> 2         -2.7470610          -3.645182
#> 3         -1.6175447          -2.515666
#> 4         -0.9547998          -1.852921
#> 5         -0.4876278          -1.385749
#> 6         -0.1414232          -1.039544
head(predict(fit, untransform = TRUE))
#>       P[Y>=2]      P[Y>=3]
#> 1 0.001302623 0.0005310131
#> 2 0.060252851 0.0254519377
#> 3 0.165543769 0.0747672299
#> 4 0.277920572 0.1355303291
#> 5 0.380452564 0.2000873101
#> 6 0.464703016 0.2612379561

p0 <- head(predict(fit, type = "response"))
p1 <- head(predict(fit, type = "response", newdata = pneumo))
p2 <- head(predict(fit, type = "response", newdata = pneumo))
p3 <- head(fitted(fit))
all.equal(p0, p1)  # Should be TRUE
#> [1] TRUE
all.equal(p1, p2)  # Should be TRUE
#> [1] TRUE
all.equal(p2, p3)  # Should be TRUE
#> [1] TRUE

predict(fit, type = "terms", se = TRUE)
#> $fitted.values
#>   poly(c(scale(let)), 2):1 poly(c(scale(let)), 2):2
#> 1               -5.1276963               -5.1276963
#> 2               -1.2326855               -1.2326855
#> 3               -0.1031692               -0.1031692
#> 4                0.5595757                0.5595757
#> 5                1.0267477                1.0267477
#> 6                1.3729523                1.3729523
#> 7                1.6576137                1.6576137
#> 8                1.8466615                1.8466615
#> attr(,"vterm.assign")
#> attr(,"vterm.assign")$`poly(c(scale(let)), 2)`
#> [1] 1 2
#> 
#> 
#> $se.fit
#>   poly(c(scale(let)), 2):1 poly(c(scale(let)), 2):2
#> 1                1.7242629                1.7242629
#> 2                0.2327739                0.2327739
#> 3                0.3311413                0.3311413
#> 4                0.3768132                0.3768132
#> 5                0.3653147                0.3653147
#> 6                0.3317412                0.3317412
#> 7                0.3044709                0.3044709
#> 8                0.3113166                0.3113166
#> attr(,"vterm.assign")
#> attr(,"vterm.assign")$`poly(c(scale(let)), 2)`
#> [1] 1 2
#> 
#> 
#> $df
#> [1] 12
#> 
#> $sigma
#> [1] 1
#>