Confidence Intervals

For any model fit by glm, add_ci() may compute confidence intervals for predictions using either a parametric method or a bootstrap. The parametric method computes confidence intervals on the scale of the linear predictor \(X \beta\) and transforms the intervals to the response level through the inverse link function \(g^{-1}\). Confidence intervals on the linear predictor level are computed using a Normal distribution for Logistic and Poisson regressions or a \(t\) distribution otherwise. (This is consistent with the default behavior for the predict.glm function.) The intervals are given by the following expressions:

\[ \begin{equation} g^{-1}\left(x'\hat{\beta} \pm z_{1 - \alpha/2} \sqrt{\hat{\sigma}^2x'(X'X)^{-1} x}\right) \end{equation} \]

for Binomial and Poisson GLMs or

\[ \begin{equation} \label{eq:glmci} g^{-1}\left(x'\hat{\beta} \pm t_{1 - \alpha/2, n-p-1} \sqrt{\hat{\sigma}^2x'(X'X)^{-1} x}\right) \end{equation} \]

for other generalized linear models. In these expressions, we regard \(X\) as the model matrix from the original fit, and \(x\) as the “new data” matrix. The default method is parametric and is called with add_ci(data, fit, ...). This is the method we generally recommend for constructing confidence intervals for model predictions.

The bootstrap method is called with add_ci(df, fit, type = "boot", ...) and was included originally for making comparisons against the parametric method. There are multiple methods of bootstrap for regression models (resampling cases, resampling residuals, parametric, etc.). The bootstrap method employed by ciTools in add_ci.glm() resamples cases and iteratively refits the model (using the default behavior of boot::boot) to determine confidence intervals. After collecting the bootstrap replicates, a bias-corrected and accelerated (BCa) bootstrap confidence interval is formed for each point in the sample df.

Although there are several methods for computing bootstrap confidence intervals, we don’t provide options to compute all of these types of intervals in ciTools. BCa intervals are slightly larger than parametric intervals, but are less biased than other types of bootstrapped intervals, including percentile based intervals. We may consider adding more types of bootstrap intervals to ciTools in the future.

x <- rnorm(100, mean = 5)
y <- rbinom(n = 100, size = 1, prob = invlogit(-20 + 4*x))
df <- data.frame(x = x, y = y)
fit <- glm(y ~ x, family = binomial)

df1 <- add_ci(df, fit, names = c("lwr", "upr"), alpha = 0.1) %>%
    mutate(type = "parametric")

df2 <- add_ci(df, fit, type = "boot", names = c("lwr", "upr"), alpha = 0.1, nSims = 500) %>%
    mutate(type = "bootstrap")

df <- bind_rows(df1, df2)

ggplot(df, aes(x = x, y = y)) +
    ggtitle("Logistic Regression", subtitle = "Model fit (black line) and confidence intervals (gray)") +
    geom_jitter(height = 0.01) +
    geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.4) +
    geom_line(aes(x =x , y = pred), size = 2) +
    facet_grid(~type)

## New names:
## * x -> x...1
## * y -> y...2
## * pred -> pred...3
## * type -> type...6
## * x -> x...7
## * ...

Prediction Intervals

Generally, parametric prediction intervals for GLMs are not available. The solution ciTools takes is to perform a parametric bootstrap on the model fit, then take quantiles on the bootstrapped data produced for each observation. The procedure is performed via arm::sim. The method of the parametric bootstrap is described by the following algorithm:

1. Fit the GLM, and collect the regression statistics \(\hat{\beta}\) and \(\hat{\mathrm{Cov}}(\hat{\beta})\). Set the number of simulations, \(M\).

2. Simulate \(M\) draws of the regression coefficients, \(\hat{\beta}_{*}\), from \(N(\hat{\beta}, \hat{\mathrm{Cov}}(\hat{\beta}))\), where \(\hat{\mathrm{Cov}}(\hat{\beta}) = \hat{\sigma}^2 (X'X)^{-1}\).

3. Simulate \([y_{*}|x]\) from the response distribution with mean \(g^{-1}(x \hat{\beta}_{*})\) and a variance determined by the response distribution.

4. Determine the \(\alpha/2\) and \(1-\alpha/2\) quantiles of the simulated response \([y_{*}|x]\) for each \(x\).

The parametric bootstrap method propagates the uncertainty in the regression effects \(\hat{\beta}\) into the simulated draws from the predictive distribution.

Generally, there are many different ways to calculate the quantiles of an empirical distribution, but the approach that ciTools takes ensures that estimated quantiles lie in the support set of the response. The choice we make corresponds to setting type = 1 in quantile().

We have seen in our simulations that the parametric bootstrap provides interval estimates with approximately nominal probability coverage. The unfortunate side effect of opting to construct prediction intervals through a parametric bootstrap is that the parameters of the predictive distributions need to hard coded for each model. For this reason, ciTools does not have complete coverage of all models one could fit with glm.

One exception to this scheme are GLMs with Gaussian errors. In this case we may parametrically calculate prediction intervals. Under Gaussian errors, glm() permits the use of the link functions “identity”, “log”, and “inverse”. The corresponding models are given by the following expressions:

\[ \begin{equation} \label{eq:gauss-link} \begin{split} y = X\beta + \epsilon\\ y = \exp(X\beta) + \epsilon \\ y = \frac{1}{X\beta} + \epsilon \end{split} \end{equation} \]

The conditional response distribution in each case may be written \(y \sim N(g^{-1}(X\beta), \sigma^2)\), and prediction intervals may be computed parametrically. Parametric prediction intervals are therefore constructed via

\[ \begin{equation} \label{eq:gauss-pi} g^{-1}(x'\beta) \pm t_{1-\alpha/2, n-p-1} \sqrt{\hat{\sigma}^2 + \hat{\sigma}^2x'(X'X)^{-1}x} \end{equation} \]

for Gaussian GLMs. As in the linear model, \(\hat{\sigma}^2\) estimates the predictive uncertainty and \(\hat{\sigma}^2 x(X'X)^{-1}x\) estimates the inferential uncertainty in the fitted values. Note that \(g^{-1}(X \hat{\beta})\) is a maximum likelihood estimate for the parameter \(g^{-1}(X \beta)\), by the functional invariance of maximum likelihood estimation.

At this point, the only GLMs that we support with add_pi() are Guassian, Poisson, Binomial, Gamma, and Quasipoisson (all of which need to be fit with glm()). The models not supported by add_pi.glm are inverse Gaussian, Quasi, and Quasi-binomial.

Poisson Example

Poisson regression is usually the first line of defense against count data, so we wish to present a complete example of quantifying uncertainty for this type of model with ciTools. For simplicity we fit a model on fake data.

We use rnorm to generate a covariate, but the randomness of \(x\) has no bearing on the model.

x <- rnorm(100, mean = 0)
y <- rpois(n = 100, lambda = exp(1.5 + 0.5*x))
df <- data.frame(x = x, y = y)
fit <- glm(y ~ x , family = poisson(link = "log"))

As seen previously, the commands in ciTools are “pipeable”. Here, we compute confidence and prediction intervals for a model fit at the \(90\%\) level. The warning message only serves to remind the user that precise quantiles cannot be formed for non-continuous distributions.

df_ints <- df %>%
  add_ci(fit, names = c("lcb", "ucb"), alpha = 0.1) %>%
  add_pi(fit, names = c("lpb", "upb"), alpha = 0.1, nSims = 20000) %>%
  head()

## Warning in add_pi.glm(., fit, names = c("lpb", "upb"), alpha = 0.1, nSims =
## 20000): The response is not continuous, so Prediction Intervals are approximate

As with other methods available in ciTools the requested statistics are computed, appended to the data frame, and returned to the user as a data frame.

df_ints %>%
    ggplot(aes(x = x, y = y)) +
    geom_point(size = 2) +
    ggtitle("Poisson Regression", subtitle = "Model fit (black line), with prediction intervals (gray), confidence intervals (dark gray)") +
    geom_line(aes(x = x, y = pred), size = 1.2) +
    geom_ribbon(aes(ymin = lcb, ymax = ucb), alpha = 0.4) +
    geom_ribbon(aes(ymin = lpb, ymax = upb), alpha = 0.2)

Since the response \(y\) is count data, and the method we used to determine the intervals involves simulation, we find that ciTools will produce “jagged” bounds when all the intervals are plotted simultaneously. Increasing the number of simulations using the nSims argument in add_pi can help reduce some of this unsightliness.

We may also wish to compute response-level probabilities and quantiles. ciTools can also handle these estimates with add_probs() and add_quantile() respectively. We use the same parametric bootstrap approach for estimating quantiles and probabilities that we employed for add_pi(). Once again, an error message reminds the user that their support is not continuous.

df %>%
  add_probs(fit, q = 10) %>%
  add_quantile(fit, p = 0.4) %>%
  head()

## Warning in add_probs.glm(., fit, q = 10): The response is not continuous, so
## estimated probabilities are only approximate

## Warning in add_quantile.glm(., fit, p = 0.4): The response is not continuous, so
## estimated quantiles are only approximate

##            x y     pred prob_less_than10 quantile0.4
## 1 -1.2725648 7 2.145159           0.9995           2
## 2 -0.3948852 0 3.492482           0.9965           3
## 3  0.1444020 4 4.711910           0.9755           4
## 4  0.0745726 6 4.532688           0.9780           4
## 5 -1.6400189 3 1.749197           1.0000           1
## 6  0.1952983 4 4.846988           0.9750           4

x <- runif(n = 100, min = 0, max = 2)
mu <- exp(1 + x)
y <- rnegbin(n = 100, mu = mu, theta = mu/(5 - 1))

df <- data.frame(x = x, y = y)
fit <- glm(y ~ x, family = quasipoisson(link = "log"))
summary(fit)$dispersion

## [1] 3.73578

df_ints <- add_ci(df, fit, names = c("lcb", "ucb"), alpha = 0.05) %>%
    add_pi(fit, names = c("lpb", "upb"), alpha = 0.1, nSims = 20000)

## Warning in add_pi.glm(., fit, names = c("lpb", "upb"), alpha = 0.1, nSims =
## 20000): The response is not continuous, so Prediction Intervals are approximate

ggplot(df_ints, aes(x = x, y = y)) +
    ggtitle("Quasipoisson Regression", subtitle = "Model fit (black line), with Prediction intervals (gray), Confidence intervals (dark gray)") +
    geom_point(size = 2) +
    geom_line(aes(x = x, y = pred), size = 1.2) +
    geom_ribbon(aes(ymin = lcb, ymax = ucb), alpha = 0.4) +
    geom_ribbon(aes(ymin = lpb, ymax = upb), alpha = 0.2)

Simulation Study for Prediction Intervals

A simulation study was performed to examine the empirical coverage probabilities of prediction intervals formed using the parametric bootstrap. We focus on these intervals because we could not find results in the literature addressing their performance. Our simulation is not comprehensive, so users of ciTools should exercise care when using these methods. In each simulation, a simple \(y = g^{-1}(mx + b)\) model is fit on a variable number of observations.

New observations were generated from the true model to determine if the empirical coverage was close to the nominal level. The mean interval width was also recorded, as were standard errors of the estimated coverage and mean interval width.

Note that in contrast to a study of confidence intervals, we do not expect interval widths to shrink to \(0\) as sample size tends to infinity. This is due to the predictive error in the conditional response distribution, which is not a factor in the construction of confidence intervals.

We take the same approach in the simulation study of each of the four models described below:

1. Set a sequence of sample sizes e.g. \(n = 20, 50, 100, ...\).
2. For each sample size, set a model matrix
3. Then loop …
   • Generate a response vector from the true model
   • Fit a GLM to the simulated response \(\boldsymbol{y}\) with the fixed model matrix
   • Calculate a prediction interval for \(y_{new}\) given \(x_0\), the midpoint of the range of \(x\) using \(2000\) bootstrap replicates.
   • Store the width of this prediction interval.
   • Generate a response \(y_{new} | x_0\) from the true model and determine if the new response is in the PI.
4. Repeat Step 3 \(M\) times.

pois <- read.csv("pois_pi_results.csv") %>%
  dplyr::select(-total_time) %>%
  rename(nominal_level = level) %>%
  dplyr::select(sample_size, everything())
knitr::kable(pois)

neg_bin <- read.csv("negbin_pi_results.csv") %>%
  dplyr::select(-total_time) %>%
  rename(nominal_level = level) %>%
  dplyr::select(sample_size, everything())
knitr::kable(neg_bin)

gam <- read.csv("gamma_pi_results.csv") %>%
  dplyr::select(-total_time) %>%
  rename(nominal_level = level) %>%
  dplyr::select(sample_size, everything())

knitr::kable(gam)

norm_log <- read.csv("gaussian_pi_loglink_results.csv") %>%
  dplyr::select(-total_time) %>%
  rename(nominal_level = level) %>%
  dplyr::select(sample_size, everything())
knitr::kable(norm_log)

Summary

ciTools is a versatile R package that helps users quantify uncertainty about their generalized linear models. Creating interval estimates that are amenable to plotting is now as simple as fitting a GLM. To date, we provide coverage for many common GLMs used by practitioners. For the models covered by ciTools, our simulations show that our confidence and prediction intervals are trustworthy.

There is still work to be done on this portion of ciTools. We would like to

1. Include more parametric methods for prediction intervals.

2. Add facilities to handle CIs and PIs when offset terms are present in the model fit.

3. Further study interval estimates pertaining to the negative binomial model.

4. Produce a simulation study that compares parametric vs. bootstrap intervals for GLMs.

5. Offer alternative prediction intervals e.g. the shortest intervals that contain \(95\%\) of the simulated data.

6. Include options beyond BCa for creating bootstrap confidence intervals.

sessionInfo()

## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                          
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] arm_1.11-2    lme4_1.1-23   Matrix_1.2-18 MASS_7.3-51.6 knitr_1.29   
## [6] ggplot2_3.3.2 ciTools_0.6.1 dplyr_1.0.1  
## 
## loaded via a namespace (and not attached):
##  [1] statmod_1.4.34      tidyselect_1.1.0    xfun_0.16          
##  [4] purrr_0.3.4         splines_4.0.2       lattice_0.20-41    
##  [7] colorspace_1.4-1    vctrs_0.3.2         generics_0.0.2     
## [10] htmltools_0.5.0     yaml_2.2.1          base64enc_0.1-3    
## [13] utf8_1.1.4          survival_3.1-12     rlang_0.4.7        
## [16] nloptr_1.2.2.2      pillar_1.4.6        withr_2.2.0        
## [19] foreign_0.8-80      glue_1.4.2          RColorBrewer_1.1-2 
## [22] jpeg_0.1-8.1        lifecycle_0.2.0     stringr_1.4.0      
## [25] munsell_0.5.0       gtable_0.3.0        htmlwidgets_1.5.1  
## [28] coda_0.19-3         evaluate_0.14       labeling_0.3       
## [31] latticeExtra_0.6-29 fansi_0.4.1         highr_0.8          
## [34] htmlTable_2.0.1     Rcpp_1.0.5          scales_1.1.1       
## [37] backports_1.1.8     checkmate_2.0.0     Hmisc_4.4-1        
## [40] abind_1.4-5         farver_2.0.3        gridExtra_2.3      
## [43] png_0.1-7           digest_0.6.25       stringi_1.4.6      
## [46] grid_4.0.2          cli_2.0.2           tools_4.0.2        
## [49] magrittr_1.5        tibble_3.0.3        Formula_1.2-3      
## [52] cluster_2.1.0       crayon_1.3.4        pkgconfig_2.0.3    
## [55] ellipsis_0.3.1      data.table_1.13.0   assertthat_0.2.1   
## [58] minqa_1.2.4         rmarkdown_2.3       rstudioapi_0.11    
## [61] R6_2.4.1            boot_1.3-25         rpart_4.1-15       
## [64] nnet_7.3-14         nlme_3.1-148        compiler_4.0.2