---
title: "Prediction and visualizing uncertainty: O-ring example"
author: "Ramses Llobet"
date: "`r Sys.Date()`"
output:
  pdf_document:
    toc: true
    toc_depth: 4
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

# check if packages are installed; if not, install them
packages <- c("tidyverse", "stargazer", "MASS", "marginaleffects", "cowplot")
not_installed <- setdiff(packages, rownames(installed.packages()))
if (length(not_installed)) install.packages(not_installed)

# load packages
library(MASS)             # always load MASS package before tidyverse due to conflict w/ dplyr
library(stargazer)
library(marginaleffects)
library(cowplot)
library(tidyverse)
```


# Acknowledgments

This illustration was originally created by former TA Brian Leung and was subsequently updated by Ramses Llobet.


# O-ring data


This code imports the O-ring dataset and creates a binary indicator, `damaged_dum`, equal to 1 if at least one O-ring was damaged and 0 otherwise. The temperature variable is also renamed to `temp` for convenience.
 
```{r}
# load data
oring <- read_csv(file="https://www.openintro.org/data/csv/orings.csv")

# some wrangling 
oring <- 
  oring |>
  mutate(damaged_dum = if_else(damaged >= 1, 1, 0)) |>
  rename(temp = temperature)

oring |> print(n=23)
```



# A brief note on logistic regression

We model the probability that an O-ring is damaged as a function of temperature using a logistic regression model:
$$
\text{Pr(Damage|Temp)} = \text{logit}^{-1} (\beta_{0} + \beta_{1} \text{Temp})
$$

In logistic regression, the *link function* that maps the linear predictor $X_i \beta$ to the probability $\pi_i = \Pr(\text{damaged}_i = 1 \mid \text{Temp}_i)$ is the logit function, which is nonlinear. In practice, we typically work with its inverse, which maps real-valued inputs to probabilities in $(0,1)$.

The model can be written equivalently on two scales:

$$
\begin{split}
\text{logit} (\pi_{i}) &= X_{i} \beta \\
\pi_{i} &= \text{logit}^{-1} (X_{i} \beta)
\end{split}
$$ 
The scale on which predictions are made (log-odds versus probabilities) is crucial for interpretation and for constructing predicted values and uncertainty intervals.


\newpage

# Logit model on O-ring data

```{r}
# logit model 
oring_logit <- glm(damaged_dum ~ temp, 
                   data = oring, 
                   family = binomial(link = "logit"))

# summary
summary(oring_logit)
```

```{r, results = "asis"}
# regression table w/ stargazer 
# type = "text" for in-console display
# type = "latex" for knitting PDF 
stargazer(
  oring_logit,
  type = "latex",
  header = FALSE,
  title = "Logistic regression: O-ring damage on temperature",
  dep.var.labels = "Log-odds of O-ring damage",
  covariate.labels = c("Temperature"),
  digits = 2
)
```


## Interpreting log-odds in logistic regression

In a logistic regression, the model is specified as

$$
\operatorname{logit}(\pi_i) = \log\left(\frac{\pi_i}{1 - \pi_i}\right)
= \beta_0 + \beta_1 \text{Temp}_i
$$

where $\pi_i = \Pr(\text{damaged}_i = 1 \mid \text{Temp}_i)$.

The coefficient $\beta_1$ represents the change in the log-odds of O-ring damage associated with a one-unit increase in temperature. A negative value of $\beta_1$ implies that higher temperatures are associated with lower odds of damage.

Because log-odds are not directly interpretable on the probability scale, it is often useful to exponentiate coefficients. The quantity

$$
\exp(\beta_1)
$$

is the odds ratio, which gives the multiplicative change in the odds of damage for a one-unit increase in temperature. Values below 1 indicate a decrease in odds, while values above 1 indicate an increase.

For substantive interpretation and prediction, we typically transform results back to the probability scale using the inverse logit:

$$
\pi_i = \operatorname{logit}^{-1}(\beta_0 + \beta_1 ,\text{Temp}_i).
$$

This step is essential when generating predicted probabilities or visualizing the relationship between temperature and the likelihood of O-ring damage.



\newpage

# Why prediciton and visualization?

Logistic regression, despite its simplicity and widespread use, is often difficult to interpret directly:

-   Logit link function: coefficients are expressed in log-odds.
    -   Exponentiation yields odds ratios, but these remain unintuitive: a one-unit increase in $x_k$ multiplies the odds by $e^{\beta_k}$.
-   Nonlinear link: marginal effects depend on the values of the predictors.
    -   The slope of the logistic curve varies with the baseline probability.
-   Probabilities are more interpretable and substantively meaningful.
-   Prediction requires incorporating uncertainty (e.g., confidence intervals).

These features motivate prediction and visualization on the probability scale rather than interpretation of raw coefficients.

## Prediction with a logit model

```{r}
# create hypothetical values for temperature
temp_hypo <- tibble(temp = 20:90) 

# predicted probabilities (response scale)
damaged_prob <- predict(oring_logit, 
                        newdata = temp_hypo, 
                        type = "response")

# linear predictor (log-odds scale)
damaged_link <- predict(oring_logit, 
                        newdata = temp_hypo, 
                        type = "link")

# inverse-logit relationship
inv.logit <- plogis

all.equal(damaged_prob, inv.logit(damaged_link))


# combine predictions with hypothetical values
damaged_pred <- bind_cols(temp_hypo,
                          damaged_prob = damaged_prob)

# visualize predicted probabilities
ggplot(damaged_pred, aes(x = temp, y = damaged_prob)) +
  geom_line()
```

What is missing from the graph?

\newpage

### Confidence intervals

#### Linear regression

```{r}
# use linear regression instead
oring_lm <- lm(damaged_dum ~ temp, data = oring)



damaged_prob_lm <- predict(oring_lm, newdata = temp_hypo, 
                           interval = "confidence", 
                           level = 0.95)

damaged_pred_lm <- bind_cols(temp_hypo, damaged_prob_lm)

# visualize 
ggplot(damaged_pred_lm, aes(x = temp, y = fit, ymin = lwr, ymax = upr)) +
  geom_line() +
  geom_ribbon(alpha = 0.2)
```

\newpage

#### logit regression (or other GLMs)

```{r}
# predict doesn't work with glm objects in terms of calculating CIs
class(oring_logit)

predict(oring_logit, 
        newdata = temp_hypo, 
        interval = "confidence", 
        level = 0.95)
```

\newpage

#### Computing CIs for logit: inverse link function

```{r}
# prediction on logit scale w/ standard errors
link_pred <- predict(oring_logit, 
                     newdata = temp_hypo,
                     type = "link", 
                     se = TRUE)

# some wrangling
link_pred <- tibble(fit = as.numeric(link_pred$fit),
                    sefit = as.numeric(link_pred$se.fit))

# critical values for 95% and 67% CIs (normal approximation)
# ignore problem of small-n for simplicity
z_95 <- qnorm(p = (1 - 0.95)/2, lower.tail = FALSE) # ~ 1.96
z_67 <- qnorm(p = (1 - 0.67)/2, lower.tail = FALSE) # ~0.97


# manually compute CIs: transform linear predictor back to probability via inverse logit 
link_pred_CIs <- 
  link_pred |>
  mutate(
    pred_prob = plogis(fit),
    lwr_95 = plogis(fit - z_95 * sefit),
    upr_95 = plogis(fit + z_95 * sefit),
    lwr_67 = plogis(fit - z_67 * sefit),
    upr_67 = plogis(fit + z_67 * sefit)
    ) |>
  bind_cols(temp_hypo)

# visualize w/ ggplot2 
link_pred_vis <- 
 ggplot(link_pred_CIs, aes(x = temp, y = pred_prob, ymin = lwr_67, ymax = upr_67)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_line(aes(y = upr_95), linetype = 2) +
  geom_line(aes(y = lwr_95), linetype = 2)

link_pred_vis
```

\newpage

#### Computing CIs for logit: simulation method via MASS::mvrnorm()

Consult Chris's lecture on [Maximum Likelihood](https://faculty.washington.edu/cadolph/mle/topic3.p.pdf) and [King et al (2000)](https://gking.harvard.edu/files/making.pdf) for reference.

Let's take a step back and think about `predict()` function: how does it calculate the predicted probability?

You can do it by hand. Quick example:

Let's say we want to know the probability of damage given that temperature is 50 degree. We know that the intercept coefficient is 15.0429 and the temperature coefficient is -0.2322.

$$
\begin{split}
\pi_{t=50} &= \text{logit}^{-1} (15.0429 \times 1 + -0.2322 \times 50) \\ 
&\approx 0.9687
\end{split}
$$

We can check the result with `predict()`

```{r}
predict(oring_logit, newdata = data.frame(temp = 50), type = "response")
```

But the problem is that we treat the estimated coefficients as certain and fail to *propagate uncertainty* from our estimation

How can we propagate uncertainty to our prediction? Counterfactual simulation!

Basic logic of counterfactual simulation:

1.  Choose a set of counterfactual value for $x_{c}$
2.  Estimate the model and obtain the parameter vector, $\hat{\beta}$, and its variance covariance matrix, $\hat{V}(\hat{\beta})$
3.  Draw $\tilde{\beta}$ from the multivariate normal $f_{MVN}(\hat{\beta}, \hat{V}(\hat{\beta}))$
4.  Calculate $\tilde{\pi_{c}} = \text{logit}^{-1} (x_{c} \tilde{\beta})$
5.  Repeat the procedure many times, summarizing this vector to get expected values and confidence intervals

```{r}
# point estimate of the parameters
pe <- coef(oring_logit)

# variance covariance of the parameters
vc <- vcov(oring_logit)
 
# set N of simulations
sims <- 100000

# simulate many betas 
sim_beta <- MASS::mvrnorm(n=sims, mu=pe, Sigma=vc)

dim(sim_beta)
```

Each row represents one trial in the simulation; there are 1,000 simulations, hence 1,000 rows.

Each column represents one simulated $\tilde{\beta}$; there are two parameters, hence 2 columns.

They encapsulate the uncertainties in our estimation.

Now we can calculate $\tilde{\pi_{c}}$ with matrix multiplication. To see this:

$$
\underbrace{
  \begin{bmatrix}
  \tilde{\pi}_{t=20, n=1} & \tilde{\pi}_{t=20, n=2} &\dots & \tilde{\pi}_{t=20, n=1,000} \\
  \tilde{\pi}_{t=21, n=1} & \tilde{\pi}_{t=21, n=2} &\dots & \tilde{\pi}_{t=21, n=1,000} \\
  \vdots & \vdots & \ddots & \vdots \\
  \tilde{\pi}_{t=90, n=1} & \tilde{\pi}_{t=90, n=2} &\dots & \tilde{\pi}_{t=90, n=1,000} \\
  \end{bmatrix} 
}_\text{71 counterfactuals; 1,000 simulations}
=
\underbrace{
  \begin{bmatrix}
  1 & x_{t=20} \\
  1 & x_{t=21} \\
  \vdots & \vdots \\ 
  1 & x_{t=90} \\
  \end{bmatrix}
}_\text{71 counterfactural; 2 predictors}
\times
\underbrace{
  \begin{bmatrix}
  \tilde{\beta}_{0, n=1} & \tilde{\beta}_{0, n=2} & \dots &\tilde{\beta}_{0, n=1,000}  \\
  \tilde{\beta}_{1, n=1} & \tilde{\beta}_{1, n=2} & \dots &\tilde{\beta}_{1, n=1,000}  \\
  \end{bmatrix}
}_\text{2 parameters; 1,000 simulations}
$$

Intuitively, let's imagine there are $n = 1,000$ parallel universes, each of which is one simulation where $\tilde{\beta}$ exhibits some particular value (from a random MVN draw).

In each parallel universe (simulation), you calculate the particular $\tilde{\pi}$ for each and every counterfactual temperature ${= \{20, 21, 22, \dots, 90\}}$. Essentially, you're repeating the manual calculation we've done above for $k = 71$ times.

Then, you repeat the procedure for each and every parallel universe (simulation).

You should get $n \times k = 1,000 \times 71$ different $\tilde{\pi}$

```{r}
# create hypothetical values for temp; plus constant
hypo_temp <- cbind(1, 20:90)

# check dimensions: 71 rows, 2 columns
dim(hypo_temp)

# check dimensions for simulated betas
dim(sim_beta)
sim_beta <- t(sim_beta)

# matrix multiplication 
sim_prob <- hypo_temp %*% sim_beta

# check dimensions for simulated probabilities
dim(sim_prob)
```

```{r}
# calculate expected values via mean()
expected_values <- apply(sim_prob, 1, mean)

# calculate confidence intervals via quantile()
CIs_95 <- apply(sim_prob, 1, quantile, prob = c(0.025, 0.975))
CIs_95 <- t(CIs_95)

CIs_67 <- apply(sim_prob, 1, quantile, prob = c(0.165, 0.835))
CIs_67 <- t(CIs_67)

# put everything together 
mvrnorm_sim_CIs <- 
  bind_cols(expected_values = expected_values, 
            CIs_95, 
            CIs_67) |>
  mutate_all(inv.logit) |>
  mutate(hypo_temp = 20:90)

# visualize w/ ggplot2
mvrnorm_sim_vis <- 
  ggplot(mvrnorm_sim_CIs, aes(x = hypo_temp, y = expected_values, ymin = `16.5%`, ymax = `83.5%`)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_line(aes(y = `97.5%`), linetype = 2) +
  geom_line(aes(y = `2.5%`), linetype = 2)

print(mvrnorm_sim_vis)
```

\newpage

#### Computing CIs for logit: marginaleffects package

```{r}
# use predictions() function from marginaleffects package
margin_pred_95 <- predictions(oring_logit, newdata = datagrid(temp = 20:90), conf_level = 0.95)
margin_pred_67 <- predictions(oring_logit, newdata = datagrid(temp = 20:90), conf_level = 0.67)

# some wrangling
margin_pred_95 <- 
  margin_pred_95 |>
  as_tibble() |>
  select(temp, estimate, conf.low, conf.high) |>
  rename(conf_low_95 = conf.low,
         conf_high_95 = conf.high)

margin_pred_67 <- 
  margin_pred_67 |>
  as_tibble() |>
  select(conf.low, conf.high) |>
  rename(conf_low_67 = conf.low,
         conf_high_67 = conf.high)

# put everything together 
margin_pred_CIs <- bind_cols(margin_pred_95, margin_pred_67)

# visualize w/ ggplot2
margin_pred_vis <- 
  ggplot(margin_pred_CIs, aes(x = temp, y = estimate, ymin = conf_low_67, ymax = conf_high_67)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_line(aes(y = conf_low_95), linetype = 2) +
  geom_line(aes(y = conf_high_95), linetype = 2)

print(margin_pred_vis)
```

\newpage

#### Computing CIs for logit: compare all three methods

```{r}
# use plot_grid() function from cowplot
plot_grid(link_pred_vis, mvrnorm_sim_vis, margin_pred_vis, 
          labels = c("Inverse link function", "Mvrnorm simulation", "marginaleffects package"))
```

\newpage

# Final remarks

-   It's reassuring that all three methods produce equivalent results

-   Hazards of over-relying on off-the-shelf functions or packages: opaque computation can produce unintended, or often wrong, results

    -   Especially when your models become more complex and with more variables

-   Manual simulations can be flexible: e.g. computing first difference and its uncertainties

    -   Given a 10 degree increase in temperature, what is the change in probabilities in damage (and its uncertainties)
    -   Also, a great conceptual check on your fundamental understanding of regression

-   We didn't talk about how to improve the graphs visually

    -   Ugly defaults; no annotation
    -   Also, there are more to the inner working of `ggplot2`
    -   After the lectures have covered more on scientific principles on visual displays, we'll return to this example

# Knitting PDF

You have to install `tinytex` before you can knit a PDF file. Run the following code. We'll talk about LaTeX next week.

```{r}
# install.packages("tinytex")
# tinytex::install_tinytex()
```
