Confidence Intervals

If we were interested in the average stopping distance as a function of speed, we can generate a confidence interval using add_ci(). The output is another data frame with three new columns: one for the model predictions, one for the lower confidence bound, and one for the upper confidence bound:

my_data_with_ci <- add_ci(my_data, model, names = c("lcb", "ucb"))
kable(head(my_data_with_ci, n =10), row.names = TRUE)

The data and the model fit can be inspected graphically with ggplot:

my_data_with_ci %>%
    ggplot(aes(x = speed, y = dist)) +
    geom_point(size = 2) +
    geom_line(aes(y = pred), size = 2, color = "maroon") +
    geom_ribbon(aes(ymin = lcb, ymax = ucb), fill =
    "royalblue1", alpha = 0.3) +
    ggtitle("Stopping Distance vs. Car Speed: 95% Confidence Interval") +
    xlab("Car Speed (mph)") +
    ylab("Stopping Distance (ft)")

The red line is the model’s estimated average stopping distance across the different car speeds, and the blue region represents the uncertainty in the average stopping distance. The default confidence level is 95 percent. Users who desire a different confidence level can use the option alpha = * to specify a custom level. For example, if 80 percent intervals are desired, then include alpha = 0.2 in the add_ci() call.

Prediction Intervals

Prediction intervals are similar to confidence intervals, but instead of conveying the uncertainty in the estimated average, prediction intervals convey the uncertainty in a new observation. There is more uncertainty about a single new observation than the average of all new observations, so prediction intervals are wider than confidence intervals. To generate prediction intervals, use add_pi:

my_data_with_pi <- add_pi(my_data, model, names = c("lpb", "upb"))

The data frame is now larger because we have two new columns for lower and upper prediction bounds tacked on to the end:

kable(head(my_data_with_pi, n = 10), row.names = TRUE)

Here is what it looks like when we represent confidence intervals and prediction intervals at the same time:

my_data %>%
  add_ci(model, names = c("lcb", "ucb")) %>%
  add_pi(model, names = c("lpb", "upb")) %>%
  ggplot(aes(x = speed, y = dist)) +
  geom_point(size = 2) +
  geom_line(aes(y = pred), size = 2, color = "maroon") +
  geom_ribbon(aes(ymin = lpb, ymax = upb), fill = "orange2",
              alpha = 0.3) +
  geom_ribbon(aes(ymin = lcb, ymax = ucb), fill =
                                             "royalblue1", alpha = 0.3) +
  ggtitle("Stopping Distance vs. Car Speed: 95% CI and 95% PI") +
  xlab("Car Speed (mph)") +
  ylab("Stopping Distance (ft)")

In the graph above, the blue confidence intervals show the uncertainty in the model fit itself (maroon line), and the orange prediction intervals shows where the model would predict 95% of new responses (Stopping Distances) to fall.

Response Probabilities

Often we want other quantities that depend on the conditional predictive distribution. These include response-level probabilities and response-level quantiles, which are accessed with the functions add_probs and add_quantile respectively.

For example, suppose in the cars data set, I want to know: For each Speed what is the probability that a new Stopping Distance will be less than 70 feet? This may be an important question if, for example, you’re in a car that is hurdling towards a cliff 70 feet away at some speed and you want to know what the probability is that you will be able to stop the car before going over the cliff. Luckily, with add_probs(), you can generate these estimates quickly, allowing you to understand your chance of survival before the problem becomes moot:

my_data %>%
  add_probs(model, q = 70) %>%
  ggplot(aes(x = speed, y = prob_less_than70)) +
  geom_line(aes(y = prob_less_than70), size = 2, color = "maroon") +
  scale_y_continuous(limits = c(0,1)) +
  ggtitle("Probability Stopping Distance is Less Than 70") +
  xlab("Car Speed (mph)") +
  ylab("Pr(Dist < 70)")

The new argument q = * is used to specify the quantile (70 feet in this case) used for computing the probabilities. It’s clear that the probability of surviving declines quickly after 15 mph. (This data set is from the 1920s. It is safer to drive toward cliffs in today’s vehicles).

Another optional argument is comparison, which defaults to "<". In this example, we wanted to know the probability that we’d be able to stop the car before it went over the cliff, which means a stopping distance less than 70 feet. If we were to specify comparison = ">", then add_probs() would return the probability that the stopping distance was greater than 70 feet.

Response Quantiles

On the other hand, suppose my car is hurdling toward a cliff, and I’m comfortable with stopping at whatever distance guarantees about 90% survivability given my speed. These distances are called response quantiles, and what we wish to compute is the 0.9-quantile (or the 90th percentile) of the distibution of Stopping Distances conditional on my car’s speed and the linear model. In ciTools we use the function add_quantile with the argument p = 0.9 to calculate the 90th percentile of the predictive distibution for each row in the data set. These quantiles will be parallel to the bounds of the prediction intervals that we calculated previously.

my_data %>%
  add_pi(model, names = c("lpb", "upb")) %>%
  add_quantile(model, p = 0.9) %>%
  ggplot(aes(x = speed, y = dist)) +
  geom_point(size = 2) +
  geom_line(aes(y = pred), size = 2, color = "maroon") +
  geom_line(aes(y = quantile0.9), size = 2, color = "forestgreen") +
  geom_ribbon(aes(ymin = lpb, ymax = upb), fill = "orange2",
              alpha = 0.3) +
  ggtitle("Stopping Distance vs. Car Speed: 95% PI with 0.9-Quantile") +
  xlab("Car Speed (mph)") +
  ylab("Stopping Distance (ft)")

The 0.9 quantile lies slightly below the upper prediction bound, which in this case is the same as the 0.975-quantile. The red line is the prediction the linear model makes, which is the same as the 0.5-quantile in this case because our model assumes normally distributed errors.

The examples above have taken each piece of analysis one step at a time. However, because ciTools was built to be fully compatible with the tidyverse, these functions can easily be chained or “piped” together in clear, legible code:

data_with_results <- my_data %>%
  add_ci(model) %>%
  add_pi(model) %>%
  add_probs(model, q= 70) %>%
  add_quantile(model, p = 0.9)

kable(head(data_with_results))