Prediction intervals for GLMs

Overview

Our implementation of prediction intervals (when simulate_pi = FALSE) follows that described by Gavin Simpson in two posts on his blog. Whilst what follows is a brief overview, more detail, including discussion on whether or not it makes sense to calculate these intervals, can be found in the following post.

Confidence interval

To calculate prediction intervals we first calculate the confidence interval on the scale of the linear predictor. The upper and lower bounds of this interval, are then fed in to the inverse link function which in turn gives us a confidence interval on the expected response.

Prediction interval

Once we have calculated the confidence interval on the response we feed the upper and lower bounds, in to the quantile function associated with the relevant distribution. The maximum and minimum values of the output are then used as the upper and lower bounds of our prediction interval.

Comparison to a bootstrap approach

Below we compare the prediction intervals from trending with those generated by the ciTools package. ciTools uses a parametric bootstrap approach so the expectation is that trending will produce a more conservative (wider) interval when we allow for uncertainty around the estimate, and a less conservative (narrower) interval when uncertainty is ignored.

The following examples build on those discussed in the ciTools glm vignette:

Example 1 - Poisson

library(ciTools)
#> ciTools version 0.6.1 (C) Institute for Defense Analyses
library(trending)
library(ggplot2)
library(patchwork)
library(MASS)
#> 
#> Attaching package: 'MASS'
#> The following object is masked from 'package:patchwork':
#> 
#>     area
#> The following object is masked from 'package:dplyr':
#> 
#>     select


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

# use ciTools to add prediction interval
dat1 <- add_pi(dat, fit, names = c("lpb", "upb"), alpha = 0.1, nsims = 20000)
#> Warning in add_pi.glm(dat, fit, names = c("lpb", "upb"), alpha = 0.1, nsims =
#> 20000): The response is not continuous, so Prediction Intervals are approximate
head(dat1)
#>             x y     pred lpb upb
#> 1  0.56810730 5 5.756049   2  10
#> 2  0.80218142 6 6.449702   3  11
#> 3 -0.01838925 4 4.328226   1   8
#> 4 -0.89913899 3 2.820823   0   6
#> 5  0.74591569 6 6.275689   2  11
#> 6 -1.57154246 1 2.034349   0   5

# add intervals with trending (no uncertainty in parameters)
poisson_model <- glm_model(y ~ x, family = "poisson")
fitted_model <- fit(poisson_model, dat)
dat2 <- predict(fitted_model, simulate_pi = FALSE, uncertain = FALSE, alpha = 0.1)
head(dat2)
#>   y           x estimate lower_ci upper_ci lower_pi upper_pi
#> 1 5  0.56810730 5.756049 5.347189 6.196173        2       10
#> 2 6  0.80218142 6.449702 5.972608 6.964906        3       11
#> 3 4 -0.01838925 4.328226 3.978685 4.708474        1        8
#> 4 3 -0.89913899 2.820823 2.473551 3.216850        0        6
#> 5 6  0.74591569 6.275689 5.818104 6.769264        2       11
#> 6 1 -1.57154246 2.034349 1.704848 2.427533        0        5

# add intervals with trending (uncertainty in parameters)
dat3 <- predict(fitted_model, simulate_pi = FALSE, alpha = 0.1)
head(dat3)
#>   y           x estimate lower_ci upper_ci lower_pi upper_pi
#> 1 5  0.56810730 5.756049 5.347189 6.196173        2       11
#> 2 6  0.80218142 6.449702 5.972608 6.964906        2       12
#> 3 4 -0.01838925 4.328226 3.978685 4.708474        1        9
#> 4 3 -0.89913899 2.820823 2.473551 3.216850        0        6
#> 5 6  0.74591569 6.275689 5.818104 6.769264        2       11
#> 6 1 -1.57154246 2.034349 1.704848 2.427533        0        5

# plots
p1 <- ggplot(dat1, aes(x, y)) +
  geom_point(size = 1) +
  geom_line(aes(y = pred), size = 1.2) +
  geom_ribbon(aes(ymin = lpb, ymax = upb), alpha = 0.2) +
  geom_ribbon(aes(ymin = `lower_pi`, ymax = `upper_pi`), data = dat2, alpha = 0.4) +
  ggtitle("Poisson regression with prediction intervals and no uncertainty in parameters", 
          subtitle = "Model fit (black line), with bootstrap intervals (gray), parametric intervals (dark gray)") +
  coord_cartesian(ylim=c(0, 30))

p2 <- ggplot(dat1, aes(x, y)) +
  geom_point(size = 1) +
  geom_line(aes(y = pred), size = 1.2) +
  geom_ribbon(aes(ymin = lpb, ymax = upb), alpha = 0.4) +
  geom_ribbon(aes(ymin = `lower_pi`, ymax = `upper_pi`), data = dat3, alpha = 0.2) +
  ggtitle("Poisson regression with prediction intervals and uncertainty in parameters", 
          subtitle = "Model fit (black line), with parametric intervals (gray), bootstrap intervals (dark gray)") +
  coord_cartesian(ylim=c(0, 30))

p1 / p2

Example 2 - Quassipoisson

# generate data
x <- runif(n = 100, min = 0, max = 2)
mu <- exp(1 + x)
y <- rnegbin(n = 100, mu = mu, theta = mu/(5 - 1))
dat <- data.frame(x = x, y = y)
fit <- glm(y ~ x, family = quasipoisson(link = "log"))

# use ciTools to add prediction interval
dat1 <- add_pi(dat, fit, names = c("lpb", "upb"), alpha = 0.1, nsims = 20000)
#> Warning in add_pi.glm(dat, fit, names = c("lpb", "upb"), alpha = 0.1, nsims =
#> 20000): The response is not continuous, so Prediction Intervals are approximate
head(dat1)
#>           x  y      pred lpb upb
#> 1 0.4486188  5  3.763902   0  13
#> 2 1.7909775 11 17.557779   5  36
#> 3 1.6927022  8 15.685697   4  32
#> 4 1.2761543 13  9.726543   1  23
#> 5 1.9332004 18 20.669608   7  39
#> 6 1.7737271  8 17.213712   5  35

# add intervals with trending (no uncertainty in parameters)
quasipoisson_model <- glm_model(y ~ x, family = quasipoisson(link = "log"))
fitted_model <- fit(quasipoisson_model, dat)
dat2 <- predict(fitted_model, simulate_pi = FALSE,  uncertain = FALSE, alpha = 0.1)
head(dat2)
#>    y         x  estimate  lower_ci  upper_ci lower_pi upper_pi
#> 1  5 0.4486188  3.763902  2.883352  4.913365        0       12
#> 2 11 1.7909775 17.557779 15.181288 20.306287        5       35
#> 3  8 1.6927022 15.685697 13.732269 17.917002        4       32
#> 4 13 1.2761543  9.726543  8.588198 11.015771        1       23
#> 5 18 1.9332004 20.669608 17.476957 24.445484        7       39
#> 6  8 1.7737271 17.213712 14.919380 19.860872        5       34

# add intervals with trending (uncertainty in parameters)
dat3 <- predict(fitted_model, simulate_pi = FALSE, alpha = 0.1)
head(dat3)
#>    y         x  estimate  lower_ci  upper_ci lower_pi upper_pi
#> 1  5 0.4486188  3.763902  2.883352  4.913365        0       15
#> 2 11 1.7909775 17.557779 15.181288 20.306287        4       39
#> 3  8 1.6927022 15.685697 13.732269 17.917002        3       35
#> 4 13 1.2761543  9.726543  8.588198 11.015771        1       25
#> 5 18 1.9332004 20.669608 17.476957 24.445484        5       45
#> 6  8 1.7737271 17.213712 14.919380 19.860872        4       38

# plots
p3 <- ggplot(dat1, aes(x, y)) +
  geom_point(size = 1) +
  geom_line(aes(y = pred), size = 1.2) +
  geom_ribbon(aes(ymin = lpb, ymax = upb), alpha = 0.2) +
  geom_ribbon(aes(ymin = `lower_pi`, ymax = `upper_pi`), data = dat2, alpha = 0.4) +
  ggtitle("Quasipoisson regression with prediction intervals and no uncertainty in parameters", 
          subtitle = "Model fit (black line), with bootstrap intervals (gray), parametric intervals (dark gray)") +
  coord_cartesian(ylim=c(0, 30))

p4 <- ggplot(dat1, aes(x, y)) +
  geom_point(size = 1) +
  geom_line(aes(y = pred), size = 1.2) +
  geom_ribbon(aes(ymin = lpb, ymax = upb), alpha = 0.4) +
  geom_ribbon(aes(ymin = `lower_pi`, ymax = `upper_pi`), data = dat3, alpha = 0.2) +
  ggtitle("Quasipoisson regression with prediction intervals and uncertainty in parameters", 
          subtitle = "Model fit (black line), with parametric intervals (gray), bootstrap intervals (dark gray)") +
  coord_cartesian(ylim=c(0, 30))

p3 / p4