Confidence intervals around predicted values

Usage

get_predicted_ci(x, ...)

# Default S3 method
get_predicted_ci(
  x,
  predictions = NULL,
  data = NULL,
  se = NULL,
  ci = 0.95,
  ci_type = "confidence",
  ci_method = NULL,
  dispersion_method = "sd",
  vcov = NULL,
  vcov_args = NULL,
  verbose = TRUE,
  ...
)

Arguments

x

A statistical model (can also be a data.frame, in which case the second argument has to be a model).

...

Other argument to be passed, for instance to get_predicted_ci().

predictions

A vector of predicted values (as obtained by stats::fitted(), stats::predict() or get_predicted()).

data

An optional data frame in which to look for variables with which to predict. If omitted, the data used to fit the model is used. Visualization matrices can be generated using get_datagrid().

se

Numeric vector of standard error of predicted values. If NULL, standard errors are calculated based on the variance-covariance matrix.

ci

The interval level. Default is NULL, to be fast even for larger models. Set the interval level to an explicit value, e.g. 0.95, for 95% CI).

ci_type

Can be "prediction" or "confidence". Prediction intervals show the range that likely contains the value of a new observation (in what range it would fall), whereas confidence intervals reflect the uncertainty around the estimated parameters (and gives the range of the link; for instance of the regression line in a linear regressions). Prediction intervals account for both the uncertainty in the model's parameters, plus the random variation of the individual values. Thus, prediction intervals are always wider than confidence intervals. Moreover, prediction intervals will not necessarily become narrower as the sample size increases (as they do not reflect only the quality of the fit). This applies mostly for "simple" linear models (like lm), as for other models (e.g., glm), prediction intervals are somewhat useless (for instance, for a binomial model for which the dependent variable is a vector of 1s and 0s, the prediction interval is... [0, 1]).

ci_method

The method for computing p values and confidence intervals. Possible values depend on model type.

NULL uses the default method, which varies based on the model type.
Most frequentist models: "wald" (default), "residual" or "normal".
Bayesian models: "quantile" (default), "hdi", "eti", and "spi".
Mixed effects lme4 models: "wald" (default), "residual", "normal", "satterthwaite", and "kenward-roger".

See get_df() for details.

dispersion_method

Bootstrap dispersion and Bayesian posterior summary: "sd" or "mad".

vcov

Variance-covariance matrix used to compute uncertainty estimates (e.g., for robust standard errors). This argument accepts a covariance matrix, a function which returns a covariance matrix, or a string which identifies the function to be used to compute the covariance matrix.

A covariance matrix
A function which returns a covariance matrix (e.g., stats::vcov())
A string which indicates the kind of uncertainty estimates to return.
- Heteroskedasticity-consistent: "HC", "HC0", "HC1", "HC2", "HC3", "HC4", "HC4m", "HC5". See ?sandwich::vcovHC
- Cluster-robust: "CR", "CR0", "CR1", "CR1p", "CR1S", "CR2", "CR3". See ?clubSandwich::vcovCR
- Bootstrap: "BS", "xy", "residual", "wild", "mammen", "fractional", "jackknife", "norm", "webb". See ?sandwich::vcovBS
- Other sandwich package functions: "HAC", "PC", "CL", "OPG", "PL".
- Kenward-Roger approximation: kenward-roger. See ?pbkrtest::vcovAdj.

One exception are models of class glmgee, which have pre-defined options for the variance-covariance matrix calculation. These are "robust", "df-adjusted", "model", "bias-corrected", and "jackknife". See ?glmtoolbox::vcov.glmgee for details.

vcov_args

List of arguments to be passed to the function identified by the vcov argument. This function is typically supplied by the sandwich or clubSandwich packages. Please refer to their documentation (e.g., ?sandwich::vcovHAC) to see the list of available arguments. If no estimation type (argument type) is given, the default type for "HC" equals the default from the sandwich package; for type "CR", the default is set to "CR3".

verbose

Toggle warnings.

Details

Typically, get_predicted() returns confidence intervals based on the standard errors as returned by the predict()-function, assuming normal distribution (+/- 1.96 * SE) resp. a Student's t-distribution (if degrees of freedom are available). If predict() for a certain class does not return standard errors (for example, merMod-objects), these are calculated manually, based on following steps: matrix-multiply X by the parameter vector B to get the predictions, then extract the variance-covariance matrix V of the parameters and compute XVX' to get the variance-covariance matrix of the predictions. The square-root of the diagonal of this matrix represent the standard errors of the predictions, which are then multiplied by the critical test-statistic value (e.g., ~1.96 for normal distribution) for the confidence intervals.

If ci_type = "prediction", prediction intervals are calculated. These are wider than confidence intervals, because they also take into account the uncertainty of the model itself. Before taking the square-root of the diagonal of the variance-covariance matrix, get_predicted_ci() adds the residual variance to these values. For mixed models, get_variance_residual() is used, while get_sigma()^2 is used for non-mixed models.

It is preferred to rely on standard errors returned by get_predicted() (i.e. returned by the predict()-function), because these are more accurate than manually calculated standard errors. Use get_predicted_ci() only if standard errors are not available otherwise. An exception are Bayesian models or bootstrapped predictions, where get_predicted_ci() returns quantiles of the posterior distribution or bootstrapped samples of the predictions. These are actually accurate standard errors resp. confidence (or uncertainty) intervals.

Examples

# Confidence Intervals for Model Predictions
# ------------------------------------------

data(mtcars)

# Linear model
# ------------
x <- lm(mpg ~ cyl + hp, data = mtcars)
predictions <- predict(x)
ci_vals <- get_predicted_ci(x, predictions, ci_type = "prediction")
head(ci_vals)
#>         SE    CI_low  CI_high
#> 1 3.255505 14.558527 27.87504
#> 2 3.255505 14.558527 27.87504
#> 3 3.305931 19.309850 32.83263
#> 4 3.255505 14.558527 27.87504
#> 5 3.303717  8.687625 22.20134
#> 6 3.266957 14.630713 27.99407
ci_vals <- get_predicted_ci(x, predictions, ci_type = "confidence")
head(ci_vals)
#>          SE   CI_low  CI_high
#> 1 0.7281647 19.72752 22.70605
#> 2 0.7281647 19.72752 22.70605
#> 3 0.9279509 24.17337 27.96911
#> 4 0.7281647 19.72752 22.70605
#> 5 0.9200310 13.56281 17.32616
#> 6 0.7777664 19.72168 22.90310
ci_vals <- get_predicted_ci(x, predictions, ci = c(0.8, 0.9, 0.95))
head(ci_vals)
#>                          SE CI_low_0.8 CI_high_0.8 CI_low_0.9 CI_high_0.9
#> Mazda RX4         0.7281647   20.26184    22.17172   19.97954    22.45403
#> Mazda RX4 Wag     0.7281647   20.26184    22.17172   19.97954    22.45403
#> Datsun 710        0.9279509   24.85429    27.28818   24.49453    27.64794
#> Hornet 4 Drive    0.7281647   20.26184    22.17172   19.97954    22.45403
#> Hornet Sportabout 0.9200310   14.23793    16.65104   13.88124    17.00773
#> Valiant           0.7777664   20.29240    22.33238   19.99087    22.63391
#>                   CI_low_0.95 CI_high_0.95
#> Mazda RX4            19.72752     22.70605
#> Mazda RX4 Wag        19.72752     22.70605
#> Datsun 710           24.17337     27.96911
#> Hornet 4 Drive       19.72752     22.70605
#> Hornet Sportabout    13.56281     17.32616
#> Valiant              19.72168     22.90310

# Bootstrapped
# ------------
predictions <- get_predicted(x, iterations = 500)
get_predicted_ci(x, predictions)
#>           SE   CI_low  CI_high
#> 1  0.6778386 20.06367 22.73426
#> 2  0.6778386 20.06367 22.73426
#> 3  1.1213280 24.07299 28.19569
#> 4  0.6778386 20.06367 22.73426
#> 5  0.7137549 13.96361 16.74841
#> 6  0.7196399 20.11708 22.92205
#> 7  0.9616939 11.54478 15.32004
#> 8  1.1099003 24.70988 29.07236
#> 9  1.1281547 24.00743 28.15526
#> 10 0.5968210 19.90476 22.18073
#> 11 0.5968210 19.90476 22.18073
#> 12 0.6867313 13.94831 16.57341
#> 13 0.6867313 13.94831 16.57341
#> 14 0.6867313 13.94831 16.57341
#> 15 0.6619519 13.42214 15.93349
#> 16 0.7053541 13.08691 15.77129
#> 17 0.8151358 12.40637 15.55107
#> 18 1.1012410 24.63705 28.94998
#> 19 1.1440413 24.87906 29.34752
#> 20 1.1031298 24.65526 28.99487
#> 21 1.1356722 23.94017 28.10764
#> 22 0.9243855 14.17452 17.75612
#> 23 0.9243855 14.17452 17.75612
#> 24 0.9616939 11.54478 15.32004
#> 25 0.7137549 13.96361 16.74841
#> 26 1.1012410 24.63705 28.94998
#> 27 1.1152049 24.14757 28.24896
#> 28 1.2189054 23.41192 27.86701
#> 29 1.1785211 10.35989 15.04183
#> 30 0.8033095 18.06171 21.24781
#> 31 2.1096285  6.00274 14.50050
#> 32 1.1944721 23.58991 27.93128

ci_vals <- get_predicted_ci(x, predictions, ci = c(0.80, 0.95))
head(ci_vals)
#>          SE CI_low_0.8 CI_high_0.8 CI_low_0.95 CI_high_0.95
#> 1 0.6778386   20.45157    22.16835    20.06367     22.73426
#> 2 0.6778386   20.45157    22.16835    20.06367     22.73426
#> 3 1.1213280   24.67816    27.54677    24.07299     28.19569
#> 4 0.6778386   20.45157    22.16835    20.06367     22.73426
#> 5 0.7137549   14.54509    16.32089    13.96361     16.74841
#> 6 0.7196399   20.53000    22.33227    20.11708     22.92205
datawizard::reshape_ci(ci_vals)
#>           SE   CI    CI_low  CI_high
#> 1  0.6778386 0.80 20.451569 22.16835
#> 2  0.6778386 0.95 20.063672 22.73426
#> 3  0.6778386 0.80 20.451569 22.16835
#> 4  0.6778386 0.95 20.063672 22.73426
#> 5  1.1213280 0.80 24.678158 27.54677
#> 6  1.1213280 0.95 24.072992 28.19569
#> 7  0.6778386 0.80 20.451569 22.16835
#> 8  0.6778386 0.95 20.063672 22.73426
#> 9  0.7137549 0.80 14.545092 16.32089
#> 10 0.7137549 0.95 13.963615 16.74841
#> 11 0.7196399 0.80 20.530002 22.33227
#> 12 0.7196399 0.95 20.117084 22.92205
#> 13 0.9616939 0.80 12.655150 14.95435
#> 14 0.9616939 0.95 11.544778 15.32004
#> 15 1.1099003 0.80 25.382182 28.18554
#> 16 1.1099003 0.95 24.709884 29.07236
#> 17 1.1281547 0.80 24.601358 27.52583
#> 18 1.1281547 0.95 24.007430 28.15526
#> 19 0.5968210 0.80 20.254035 21.81961
#> 20 0.5968210 0.95 19.904758 22.18073
#> 21 0.5968210 0.80 20.254035 21.81961
#> 22 0.5968210 0.95 19.904758 22.18073
#> 23 0.6867313 0.80 14.457044 16.21674
#> 24 0.6867313 0.95 13.948306 16.57341
#> 25 0.6867313 0.80 14.457044 16.21674
#> 26 0.6867313 0.95 13.948306 16.57341
#> 27 0.6867313 0.80 14.457044 16.21674
#> 28 0.6867313 0.95 13.948306 16.57341
#> 29 0.6619519 0.80 13.896709 15.57043
#> 30 0.6619519 0.95 13.422141 15.93349
#> 31 0.7053541 0.80 13.603118 15.38838
#> 32 0.7053541 0.95 13.086914 15.77129
#> 33 0.8151358 0.80 13.135102 15.14540
#> 34 0.8151358 0.95 12.406367 15.55107
#> 35 1.1012410 0.80 25.282352 28.10660
#> 36 1.1012410 0.95 24.637053 28.94998
#> 37 1.1440413 0.80 25.552262 28.39162
#> 38 1.1440413 0.95 24.879062 29.34752
#> 39 1.1031298 0.80 25.307084 28.11418
#> 40 1.1031298 0.95 24.655261 28.99487
#> 41 1.1356722 0.80 24.536732 27.48013
#> 42 1.1356722 0.95 23.940169 28.10764
#> 43 0.9243855 0.80 14.860147 17.07831
#> 44 0.9243855 0.95 14.174521 17.75612
#> 45 0.9243855 0.80 14.860147 17.07831
#> 46 0.9243855 0.95 14.174521 17.75612
#> 47 0.9616939 0.80 12.655150 14.95435
#> 48 0.9616939 0.95 11.544778 15.32004
#> 49 0.7137549 0.80 14.545092 16.32089
#> 50 0.7137549 0.95 13.963615 16.74841
#> 51 1.1012410 0.80 25.282352 28.10660
#> 52 1.1012410 0.95 24.637053 28.94998
#> 53 1.1152049 0.80 24.740830 27.58820
#> 54 1.1152049 0.95 24.147568 28.24896
#> 55 1.2189054 0.80 24.021410 27.27151
#> 56 1.2189054 0.95 23.411922 27.86701
#> 57 1.1785211 0.80 12.017648 14.69466
#> 58 1.1785211 0.95 10.359894 15.04183
#> 59 0.8033095 0.80 18.845838 20.82505
#> 60 0.8033095 0.95 18.061714 21.24781
#> 61 2.1096285 0.80  9.139553 13.95768
#> 62 2.1096285 0.95  6.002740 14.50050
#> 63 1.1944721 0.80 24.160657 27.32296
#> 64 1.1944721 0.95 23.589908 27.93128

ci_vals <- get_predicted_ci(x,
  predictions,
  dispersion_method = "MAD",
  ci_method = "HDI"
)
head(ci_vals)
#>          SE   CI_low  CI_high
#> 1 0.6511853 20.08775 22.75861
#> 2 0.6511853 20.08775 22.75861
#> 3 1.1886593 24.12011 28.22016
#> 4 0.6511853 20.08775 22.75861
#> 5 0.6922916 13.92620 16.69426
#> 6 0.7042060 20.11663 22.92984


# Logistic model
# --------------
x <- glm(vs ~ wt, data = mtcars, family = "binomial")
predictions <- predict(x, type = "link")
ci_vals <- get_predicted_ci(x, predictions, ci_type = "prediction")
head(ci_vals)
#>                   CI_low CI_high
#> Mazda RX4           -Inf     Inf
#> Mazda RX4 Wag       -Inf     Inf
#> Datsun 710          -Inf     Inf
#> Hornet 4 Drive      -Inf     Inf
#> Hornet Sportabout   -Inf     Inf
#> Valiant             -Inf     Inf
ci_vals <- get_predicted_ci(x, predictions, ci_type = "confidence")
head(ci_vals)
#>          SE     CI_low   CI_high
#> 1 0.5623444 -0.3931282 1.8112213
#> 2 0.4690190 -0.6974034 1.1411172
#> 3 0.7195076 -0.1279982 2.6924199
#> 4 0.4459072 -1.3016913 0.4462326
#> 5 0.5021936 -1.8418839 0.1266787
#> 6 0.5094490 -1.8943152 0.1026881