---
title: "Prediction: Foundations and Its Role in Outcome Regression"
author: "Se Yoon Lee"
date: "2026-07-17"
output:
  html_document:
    toc: true
    toc_depth: 3
    toc_float: true
    number_sections: true
    theme: readable
    code_folding: show
---

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

<style type="text/css">
.tutorial-box {
  margin: 1.3em 0 1.5em 0;
  border-radius: 14px;
  overflow: hidden;
  background: #eeeeee;
  box-shadow: 8px 8px 12px rgba(0, 0, 0, 0.18);
}
.tutorial-title {
  background: #cfcfcf;
  padding: 0.55em 0.9em;
  font-weight: 700;
  font-size: 1.05em;
}
.tutorial-title p { margin: 0; }
.tutorial-body { padding: 1.0em 1.25em 1.1em 1.25em; }
.tutorial-body > p:first-child { margin-top: 0; }
.prediction-note {
  background: #DDFBE1;
  border-left: 6px solid #59A14F;
  padding: 0.8em 1.0em;
  margin: 1.1em 0;
}
.causal-note {
  background: #FFF6C7;
  border-left: 6px solid #E3B341;
  padding: 0.8em 1.0em;
  margin: 1.1em 0;
}
.warning-note {
  background: #FFE0E0;
  border-left: 6px solid #E15759;
  padding: 0.8em 1.0em;
  margin: 1.1em 0;
}
</style>

```{r picture-helpers, include=FALSE}
start_picture <- function() {
  par(mar = c(0, 0, 0, 0), xaxs = "i", yaxs = "i")
  plot.new()
  plot.window(xlim = c(0, 1), ylim = c(0, 1))
}

draw_text <- function(x, y, label, cex = 1, font = 1, col = "#222222") {
  lines <- strsplit(label, "\n", fixed = TRUE)[[1]]
  offsets <- seq((length(lines) - 1) / 2, -(length(lines) - 1) / 2) *
    0.058 * cex
  text(rep(x, length(lines)), y + offsets, labels = lines,
       cex = cex, font = font, col = col)
}

draw_box <- function(x, y, w, h, label, fill = "#FFFFFF",
                     border = "#333333", lwd = 2, lty = 1,
                     cex = 1, font = 1) {
  rect(x - w / 2, y - h / 2, x + w / 2, y + h / 2,
       col = fill, border = border, lwd = lwd, lty = lty)
  draw_text(x, y, label, cex = cex, font = font)
}

draw_arrow <- function(x0, y0, x1, y1, col = "#333333",
                       lwd = 2, lty = 1) {
  arrows(x0, y0, x1, y1, length = 0.08, lwd = lwd,
         col = col, lty = lty)
}

bound_probability <- function(p, eps = 1e-8) {
  pmin(pmax(p, eps), 1 - eps)
}

auc_rank <- function(y, p) {
  n1 <- sum(y == 1)
  n0 <- sum(y == 0)
  (sum(rank(p)[y == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}

regression_metrics <- function(y, pred) {
  c(
    MAE = mean(abs(y - pred)),
    RMSE = sqrt(mean((y - pred)^2)),
    R_squared = 1 - sum((y - pred)^2) / sum((y - mean(y))^2)
  )
}
```

# General Prediction

## What Is Prediction?

Prediction uses information that is available now to estimate an outcome that is unknown, unobserved, or not yet observed.

Examples include:

- predicting tomorrow's temperature from recent weather measurements;
- predicting a student's exam score from earlier assignments;
- predicting whether a customer will renew a subscription; and
- predicting which of several image categories an object belongs to.

The general notation is

$$
\widehat{Y} = \hat f(X),
$$

where \(Y\) is the outcome, \(X\) is the available information, and \(\hat f\) is a prediction rule learned from data.

::: {.tutorial-box}
::: {.tutorial-title}
Learning Goals
:::
::: {.tutorial-body}
This note first develops prediction without causal language. We study two frameworks: prediction of a continuous outcome and prediction of a categorical outcome. Only after those foundations are clear do we apply prediction to outcome-regression modeling.
:::
:::

```{r prediction-map-picture, echo=FALSE, fig.width=8.8, fig.height=4.0}
start_picture()

draw_box(0.15, 0.64, 0.23, 0.20,
         "Available information\nX",
         fill = "#F1CE63", cex = 0.72, font = 2)
draw_box(0.50, 0.64, 0.26, 0.20,
         "Prediction rule\nf hat",
         fill = "#A0CBE8", cex = 0.72, font = 2)
draw_box(0.85, 0.64, 0.23, 0.20,
         "Predicted outcome\nY hat",
         fill = "#DDFBE1", cex = 0.68, font = 2)

draw_arrow(0.27, 0.64, 0.36, 0.64)
draw_arrow(0.63, 0.64, 0.73, 0.64)

draw_text(0.50, 0.25,
          "Observed pairs (X, Y) teach the rule how inputs relate to outcomes.",
          cex = 0.75, font = 2)
```

The form of the prediction depends on the type of outcome.

| Outcome type | Example | Natural prediction |
|:--|:--|:--|
| Continuous | Exam score, blood pressure, price | A number on the outcome scale |
| Binary categorical | Pass/fail, disease/no disease | A probability for the event |
| Multicategory | Red/green/blue, diagnosis A/B/C | One probability for each category |

## The Prediction Workflow

Although prediction methods differ, the basic workflow is stable.

1. Define the outcome and the information available at prediction time.
2. Use training data to learn a prediction rule.
3. Apply the rule to new feature values.
4. Compare predictions with observed outcomes.
5. Improve or simplify the rule based on new-data performance.

```{r workflow-picture, echo=FALSE, fig.width=8.8, fig.height=3.8}
start_picture()

draw_box(0.09, 0.64, 0.16, 0.18, "Define\nY and X",
         fill = "#E8E8E8", cex = 0.62, font = 2)
draw_box(0.29, 0.64, 0.16, 0.18, "Collect\nexamples",
         fill = "#F1CE63", cex = 0.62, font = 2)
draw_box(0.50, 0.64, 0.16, 0.18, "Learn\nf hat",
         fill = "#A0CBE8", cex = 0.62, font = 2)
draw_box(0.71, 0.64, 0.16, 0.18, "Predict\nnew cases",
         fill = "#DDFBE1", cex = 0.62, font = 2)
draw_box(0.91, 0.64, 0.16, 0.18, "Evaluate\nerrors",
         fill = "#FFE0E0", cex = 0.62, font = 2)

draw_arrow(0.17, 0.64, 0.21, 0.64)
draw_arrow(0.37, 0.64, 0.42, 0.64)
draw_arrow(0.58, 0.64, 0.63, 0.64)
draw_arrow(0.79, 0.64, 0.83, 0.64)

draw_text(0.50, 0.25,
          "The outcome type determines the model output and the error measure.",
          cex = 0.75, font = 2)
```

## Training Data, Test Data, and Generalization

A model should predict observations it has not already seen. We therefore distinguish:

- **training data**, used to estimate the prediction rule; and
- **test data**, held aside and used to evaluate the fitted rule.

Performance measured on the training data is usually optimistic. Test performance better represents **generalization**: how well the learned pattern transfers to new observations from a similar population.

```{r train-test-picture, echo=FALSE, fig.width=8.8, fig.height=4.0}
start_picture()

draw_box(0.14, 0.64, 0.20, 0.19, "Available data",
         fill = "#E8E8E8", cex = 0.72, font = 2)
draw_box(0.45, 0.78, 0.25, 0.17, "Training set\nlearn the rule",
         fill = "#A0CBE8", cex = 0.62, font = 2)
draw_box(0.45, 0.43, 0.25, 0.17, "Test set\nhide outcomes",
         fill = "#F1CE63", cex = 0.62, font = 2)
draw_box(0.81, 0.60, 0.26, 0.21,
         "Test performance\napproximates\nnew-data performance",
         fill = "#DDFBE1", cex = 0.56, font = 2)

draw_arrow(0.24, 0.68, 0.32, 0.75)
draw_arrow(0.24, 0.59, 0.32, 0.47)
draw_arrow(0.58, 0.45, 0.68, 0.55)
```

When the dataset is small, \(K\)-fold cross-validation repeatedly holds out different portions of the data. Each observation is evaluated using a model that was not trained on that observation.

# Framework 1: Continuous-Outcome Prediction

## The Target Is a Number

Suppose the outcome \(Y\) is continuous. Examples include test score, height, blood pressure, cost, or time.

The most common target is the conditional mean:

$$
f^*(x) = E(Y \mid X=x).
$$

For a new case with \(X=x\), the prediction \(\hat f(x)\) is a point on the same scale as \(Y\). If the outcome is an exam score, the prediction is an exam score—not a probability.

## Visual Example: Predicting an Exam Score

The following artificial dataset is used only to explain prediction. It has no causal interpretation. Study hours are the feature \(X\), and exam score is the continuous outcome \(Y\).

```{r simulate-continuous-data}
set.seed(20260717)

n_continuous <- 160
continuous_data <- data.frame(
  study_hours = runif(n_continuous, 0, 10)
)

continuous_data$exam_score <-
  48 + 7 * continuous_data$study_hours -
  0.35 * continuous_data$study_hours^2 +
  rnorm(n_continuous, mean = 0, sd = 5)

training_rows <- sample(
  seq_len(n_continuous),
  size = floor(0.70 * n_continuous),
  replace = FALSE
)

continuous_train <- continuous_data[training_rows, ]
continuous_test <- continuous_data[-training_rows, ]

data.frame(
  Sample = c("Training", "Test"),
  N = c(nrow(continuous_train), nrow(continuous_test)),
  Mean_score = round(c(
    mean(continuous_train$exam_score),
    mean(continuous_test$exam_score)
  ), 1)
)
```

```{r continuous-data-picture, echo=FALSE, fig.width=8.2, fig.height=5.3}
plot(
  continuous_train$study_hours,
  continuous_train$exam_score,
  pch = 19,
  col = adjustcolor("#4E79A7", alpha.f = 0.65),
  xlab = "Study hours",
  ylab = "Exam score",
  main = "Training Data for a Continuous Outcome"
)
```

The points do not lie on a perfect curve. Prediction learns the systematic pattern while recognizing that individual outcomes also vary around that pattern.

## Fit Prediction Rules

A linear regression rule has the form

$$
\hat Y = \hat\beta_0 + \hat\beta_1 X.
$$

A curved pattern can be represented by adding a squared term:

$$
\hat Y = \hat\beta_0 + \hat\beta_1 X + \hat\beta_2 X^2.
$$

```{r fit-continuous-models}
linear_model <- lm(
  exam_score ~ study_hours,
  data = continuous_train
)

curved_model <- lm(
  exam_score ~ study_hours + I(study_hours^2),
  data = continuous_train
)

continuous_test$linear_prediction <- predict(
  linear_model,
  newdata = continuous_test
)

continuous_test$curved_prediction <- predict(
  curved_model,
  newdata = continuous_test
)
```

```{r continuous-fits-picture, echo=FALSE, fig.width=8.2, fig.height=5.5}
x_grid <- data.frame(study_hours = seq(0, 10, length.out = 300))

plot(
  continuous_train$study_hours,
  continuous_train$exam_score,
  pch = 19,
  col = adjustcolor("#A0CBE8", alpha.f = 0.65),
  xlab = "Study hours",
  ylab = "Exam score",
  main = "Two Prediction Rules Learned from Training Data"
)

lines(
  x_grid$study_hours,
  predict(linear_model, newdata = x_grid),
  col = "#E15759",
  lwd = 3
)

lines(
  x_grid$study_hours,
  predict(curved_model, newdata = x_grid),
  col = "#59A14F",
  lwd = 3
)

legend(
  "bottomright",
  legend = c("Linear rule", "Curved rule"),
  col = c("#E15759", "#59A14F"),
  lwd = 3,
  bty = "n"
)
```

Each fitted line gives one predicted mean score for every value of study hours.

## Prediction Errors for a Continuous Outcome

For observation \(i\), the prediction error or residual is

$$
e_i = Y_i - \hat Y_i.
$$

A positive error means the observed outcome was higher than predicted. A negative error means it was lower than predicted.

```{r continuous-residual-picture, echo=FALSE, fig.width=8.2, fig.height=5.5}
plot(
  continuous_test$study_hours,
  continuous_test$exam_score,
  pch = 19,
  col = "#4E79A7",
  xlab = "Study hours",
  ylab = "Exam score",
  main = "Prediction Errors in the Test Set"
)

segments(
  x0 = continuous_test$study_hours,
  y0 = continuous_test$curved_prediction,
  x1 = continuous_test$study_hours,
  y1 = continuous_test$exam_score,
  col = adjustcolor("#E15759", alpha.f = 0.65),
  lwd = 1.5
)

points(
  continuous_test$study_hours,
  continuous_test$curved_prediction,
  pch = 4,
  cex = 1.1,
  lwd = 2,
  col = "#59A14F"
)

legend(
  "bottomright",
  legend = c("Observed outcome", "Predicted outcome", "Prediction error"),
  pch = c(19, 4, NA),
  lty = c(NA, NA, 1),
  col = c("#4E79A7", "#59A14F", "#E15759"),
  bty = "n"
)
```

Each red vertical segment is one error \(Y_i-\hat Y_i\). Shorter segments mean more accurate predictions.

Common test-set summaries are:

$$
MAE = \frac{1}{n}\sum_{i=1}^n |Y_i-\hat Y_i|,
$$

$$
RMSE =
\sqrt{\frac{1}{n}\sum_{i=1}^n(Y_i-\hat Y_i)^2},
$$

and

$$
R^2 =
1 -
\frac{\sum_i(Y_i-\hat Y_i)^2}
{\sum_i(Y_i-\bar Y)^2}.
$$

MAE and RMSE are on the outcome scale. RMSE penalizes large errors more strongly. \(R^2\) compares the fitted rule with predicting the same test-set mean for everyone. Test-set \(R^2\) can be negative when the fitted rule performs worse than that simple benchmark.

```{r continuous-model-evaluation}
linear_metrics <- regression_metrics(
  continuous_test$exam_score,
  continuous_test$linear_prediction
)

curved_metrics <- regression_metrics(
  continuous_test$exam_score,
  continuous_test$curved_prediction
)

continuous_metrics <- data.frame(
  Model = c("Linear", "Curved"),
  MAE = c(linear_metrics["MAE"], curved_metrics["MAE"]),
  RMSE = c(linear_metrics["RMSE"], curved_metrics["RMSE"]),
  R_squared = c(
    linear_metrics["R_squared"],
    curved_metrics["R_squared"]
  )
)

continuous_metrics[, -1] <- round(continuous_metrics[, -1], 3)
knitr::kable(continuous_metrics)
```

```{r continuous-test-picture, echo=FALSE, fig.width=7.5, fig.height=6.5}
score_plot_limits <- range(
  c(
    continuous_test$exam_score,
    continuous_test$curved_prediction
  )
)
score_plot_padding <- 0.04 * diff(score_plot_limits)
score_plot_limits <- score_plot_limits + c(
  -score_plot_padding,
  score_plot_padding
)
old_par <- par(pty = "s")

plot(
  continuous_test$exam_score,
  continuous_test$curved_prediction,
  pch = 19,
  col = "#4E79A7",
  xlim = score_plot_limits,
  ylim = score_plot_limits,
  xlab = "Observed exam score",
  ylab = "Predicted exam score",
  main = "Predictions in the Held-Out Test Set"
)
abline(0, 1, lty = 2, lwd = 2, col = "#E15759")
par(old_par)
```

Points close to the diagonal have small errors. The remaining scatter reminds us that a predicted conditional mean is not a guarantee for an individual.

## Point Predictions and Prediction Intervals

A point prediction estimates the center of the outcome distribution. A prediction interval describes uncertainty for an individual future outcome.

```{r continuous-prediction-interval}
new_student <- data.frame(study_hours = 6)

predict(
  curved_model,
  newdata = new_student,
  interval = "prediction",
  level = 0.95
)
```

The interval is wider than uncertainty about the mean because individual outcomes vary even when the mean pattern is known.

# Framework 2: Categorical-Outcome Prediction

## The Target Is a Category Probability

For a binary outcome \(Y\in\{0,1\}\), the conditional mean is a probability:

$$
E(Y\mid X=x)=P(Y=1\mid X=x).
$$

The natural prediction is therefore

$$
\hat p(x)=\widehat P(Y=1\mid X=x),
$$

a number between 0 and 1. A probability of 0.70 is not yet a category assignment. It represents estimated uncertainty about the outcome.

## Visual Example: Predicting Pass or Fail

We now use a second artificial dataset. The outcome is pass (\(Y=1\)) or fail (\(Y=0\)), and the feature is a practice score.

```{r simulate-binary-data}
set.seed(20260718)

n_binary <- 240
binary_data <- data.frame(
  practice_score = runif(n_binary, 0, 10)
)

true_pass_probability <- plogis(
  -4 + 0.75 * binary_data$practice_score
)

binary_data$passed <- rbinom(
  n_binary,
  size = 1,
  prob = true_pass_probability
)

binary_training_rows <- sample(
  seq_len(n_binary),
  size = floor(0.70 * n_binary),
  replace = FALSE
)

binary_train <- binary_data[binary_training_rows, ]
binary_test <- binary_data[-binary_training_rows, ]
```

Logistic regression models the log odds of the event:

$$
\operatorname{logit}\{P(Y=1\mid X)\}
=\beta_0+\beta_1X,
$$

and converts the fitted value back to the probability scale:

$$
\hat p(X)=
\frac{1}{1+\exp\{-(\hat\beta_0+\hat\beta_1X)\}}.
$$

```{r fit-binary-model}
binary_model <- glm(
  passed ~ practice_score,
  data = binary_train,
  family = binomial()
)

binary_test$predicted_probability <- predict(
  binary_model,
  newdata = binary_test,
  type = "response"
)
```

```{r binary-curve-picture, echo=FALSE, fig.width=8.2, fig.height=5.5}
binary_grid <- data.frame(
  practice_score = seq(0, 10, length.out = 300)
)

set.seed(19)
jittered_outcome <- jitter(binary_train$passed, amount = 0.04)

plot(
  binary_train$practice_score,
  jittered_outcome,
  pch = 19,
  col = adjustcolor("#7A7A7A", alpha.f = 0.45),
  xlab = "Practice score",
  ylab = "Pass probability / observed outcome",
  yaxt = "n",
  ylim = c(-0.08, 1.08),
  main = "Binary Outcomes and a Predicted Probability Curve"
)
axis(2, at = c(0, 0.25, 0.50, 0.75, 1))

lines(
  binary_grid$practice_score,
  predict(binary_model, newdata = binary_grid, type = "response"),
  lwd = 3,
  col = "#4E79A7"
)

text(0.5, 0.08, "Fail = 0", pos = 4, font = 2)
text(0.5, 0.92, "Pass = 1", pos = 4, font = 2)
```

The observed outcomes occupy only 0 and 1, but the fitted rule produces a smooth range of probabilities.

## Evaluate Probability Predictions

Probability predictions should be evaluated as probabilities before converting them to categories.

### Brier Score

The Brier score is squared prediction error:

$$
\text{Brier}
=\frac{1}{n}\sum_{i=1}^n(Y_i-\hat p_i)^2.
$$

Smaller values indicate better probability predictions.

### Log Loss

Binary log loss is

$$
\text{Log loss}
=-\frac{1}{n}\sum_{i=1}^n
\left[
Y_i\log(\hat p_i)
+(1-Y_i)\log(1-\hat p_i)
\right].
$$

Smaller is better. Confident but incorrect predictions receive a large penalty.

### Discrimination

Discrimination asks whether cases tend to receive higher predicted probabilities than non-cases. The area under the ROC curve (AUC) is 0.5 for random ranking and 1 for perfect ranking.

### Calibration

Calibration asks whether predicted probabilities agree with observed frequencies. Among observations assigned a probability near 0.70, approximately 70 percent should experience the event in well-calibrated new data.

```{r binary-probability-metrics}
y_binary_test <- binary_test$passed
p_binary_test <- bound_probability(binary_test$predicted_probability)

binary_metrics <- data.frame(
  Metric = c("Brier score", "Log loss", "AUC"),
  Value = c(
    mean((y_binary_test - p_binary_test)^2),
    -mean(
      y_binary_test * log(p_binary_test) +
      (1 - y_binary_test) * log(1 - p_binary_test)
    ),
    auc_rank(y_binary_test, p_binary_test)
  )
)

binary_metrics$Value <- round(binary_metrics$Value, 4)
knitr::kable(binary_metrics)
```

```{r binary-calibration}
calibration_group <- cut(
  rank(p_binary_test, ties.method = "first"),
  breaks = 6,
  labels = FALSE
)

binary_calibration <- aggregate(
  cbind(
    predicted = p_binary_test,
    observed = y_binary_test
  ),
  by = list(Group = calibration_group),
  FUN = mean
)

binary_calibration$N <- as.vector(table(calibration_group))

knitr::kable(
  transform(
    binary_calibration,
    predicted = round(predicted, 3),
    observed = round(observed, 3)
  )
)
```

```{r binary-calibration-picture, echo=FALSE, fig.width=7.5, fig.height=6.5}
old_par <- par(pty = "s")

plot(
  binary_calibration$predicted,
  binary_calibration$observed,
  pch = 19,
  cex = 1.5,
  col = "#4E79A7",
  xlim = c(0, 1),
  ylim = c(0, 1),
  xlab = "Mean predicted probability",
  ylab = "Observed event proportion",
  main = "Test-Set Calibration"
)
abline(0, 1, lty = 2, lwd = 2, col = "#E15759")
text(
  binary_calibration$predicted,
  binary_calibration$observed,
  labels = binary_calibration$Group,
  pos = 3,
  cex = 0.8
)
par(old_par)
```

## From Probability to Category

Classification requires a decision threshold \(c\):

$$
\widehat Y_i =
\begin{cases}
1, & \hat p_i \ge c,\\
0, & \hat p_i < c.
\end{cases}
$$

The threshold is a decision choice, not a fact learned automatically from the outcome. Changing it changes sensitivity, specificity, and the numbers of false positives and false negatives.

```{r binary-classification}
threshold <- 0.50
predicted_class <- as.integer(p_binary_test >= threshold)

confusion_matrix <- table(
  Predicted = factor(
    predicted_class,
    levels = 0:1,
    labels = c("Fail", "Pass")
  ),
  Observed = factor(
    y_binary_test,
    levels = 0:1,
    labels = c("Fail", "Pass")
  )
)

confusion_matrix

sensitivity <- sum(predicted_class == 1 & y_binary_test == 1) /
  sum(y_binary_test == 1)
specificity <- sum(predicted_class == 0 & y_binary_test == 0) /
  sum(y_binary_test == 0)
accuracy <- mean(predicted_class == y_binary_test)

data.frame(
  Threshold = threshold,
  Sensitivity = round(sensitivity, 3),
  Specificity = round(specificity, 3),
  Accuracy = round(accuracy, 3)
)
```

Probability prediction and classification are related but distinct:

- probability prediction estimates uncertainty;
- classification turns that uncertainty into an action;
- the best threshold depends on the costs of the two kinds of error.

## More Than Two Categories

If \(Y\) has \(K\) unordered categories, the model predicts a probability vector:

$$
\left\{
\hat p_1(X),\hat p_2(X),\ldots,\hat p_K(X)
\right\},
\qquad
\sum_{k=1}^K \hat p_k(X)=1.
$$

The predicted class is often the category with the largest probability. Multinomial logistic regression, classification trees, and many machine-learning methods can estimate these probabilities.

```{r multicategory-picture, echo=FALSE, fig.width=7.5, fig.height=6.5}
set.seed(20260719)

class_a <- cbind(rnorm(45, -1.5, 0.7), rnorm(45, 0.8, 0.7))
class_b <- cbind(rnorm(45, 1.5, 0.7), rnorm(45, 0.8, 0.7))
class_c <- cbind(rnorm(45, 0.0, 0.7), rnorm(45, -1.3, 0.7))
class_data <- rbind(class_a, class_b, class_c)
class_label <- rep(c("A", "B", "C"), each = 45)
class_color <- c(A = "#4E79A7", B = "#E15759", C = "#59A14F")
class_plot_limits <- range(class_data)
class_plot_padding <- 0.04 * diff(class_plot_limits)
class_plot_limits <- class_plot_limits + c(
  -class_plot_padding,
  class_plot_padding
)
old_par <- par(pty = "s")

plot(
  class_data[, 1],
  class_data[, 2],
  pch = 19,
  col = class_color[class_label],
  xlim = class_plot_limits,
  ylim = class_plot_limits,
  xlab = "Feature 1",
  ylab = "Feature 2",
  main = "A Categorical Outcome with Three Classes"
)

legend(
  "topright",
  legend = names(class_color),
  col = class_color,
  pch = 19,
  title = "Class",
  bty = "n"
)
par(old_par)
```

For ordered categories such as mild, moderate, and severe, models can also use the category ordering.

# Ideas Shared by Both Frameworks

## Continuous and Categorical Prediction Compared

| Feature | Continuous outcome | Categorical outcome |
|:--|:--|:--|
| Example | Exam score | Pass/fail |
| Model output | Predicted number | Predicted category probabilities |
| Common model | Linear regression | Logistic regression |
| Individual error | \(Y-\hat Y\) | Probability loss or classification error |
| Common metrics | MAE, RMSE, \(R^2\) | Brier score, log loss, AUC, calibration |
| Optional decision step | Usually none | Choose a threshold or largest probability |

The common principle is to learn a function from inputs to an outcome-relevant quantity and to judge it on observations not used to learn the function.

## Underfitting and Overfitting

**Underfitting** occurs when a prediction rule is too simple to capture an important pattern. **Overfitting** occurs when a rule follows random details of the training sample that do not repeat in new data.

| Model behavior | Training performance | Test performance |
|:--|:--|:--|
| Underfit | Poor | Poor |
| Useful complexity | Good | Good |
| Overfit | Extremely good | Worse |

```{r complexity-picture, echo=FALSE, fig.width=8.8, fig.height=4.3}
start_picture()

draw_box(0.17, 0.65, 0.25, 0.20,
         "Too simple\nunderfitting",
         fill = "#FFE0E0", cex = 0.68, font = 2)
draw_box(0.50, 0.65, 0.25, 0.20,
         "Useful structure\ngeneralization",
         fill = "#DDFBE1", cex = 0.68, font = 2)
draw_box(0.83, 0.65, 0.25, 0.20,
         "Too flexible\noverfitting",
         fill = "#E4E6FF", cex = 0.68, font = 2)

draw_arrow(0.30, 0.65, 0.37, 0.65)
draw_arrow(0.63, 0.65, 0.70, 0.65)

draw_text(0.50, 0.25,
          "Model complexity should be chosen using new-data performance.",
          cex = 0.75, font = 2)
```

More predictors and more flexible curves can reduce training error without improving test error. Validation, cross-validation, regularization, and subject-matter knowledge help control complexity.

## Prediction Does Not Explain Why

Prediction answers:

> Given the information available, what outcome is likely?

It does not by itself answer:

- Why did the outcome occur?
- What would happen if one input were changed by intervention?
- Is a predictor a cause, consequence, or merely a correlate of the outcome?

This distinction becomes essential when prediction is used inside causal inference.

# Application to Outcome Regression

## Where Prediction Enters Outcome Regression

We now move from general prediction to the causal inference course.

Let:

- \(T\) be a binary treatment or exposure;
- \(X\) be baseline covariates; and
- \(Y\) be the outcome.

An outcome model estimates

$$
m(t,x)=E(Y\mid T=t,X=x).
$$

For a binary outcome, this is a probability:

$$
m(t,x)=P(Y=1\mid T=t,X=x).
$$

Because \(T\) is binary, the prediction function is evaluated at two treatment values:

- **Treatment or exposure set to 1:**

  $$
  m(1,x)
  =
  E(Y\mid T=1,X=x)
  =
  P(Y=1\mid T=1,X=x).
  $$

  For people with baseline covariate profile \(X=x\), \(m(1,x)\) is the conditional mean outcome. When \(Y\) is binary, it is the conditional probability of the event among those with \(T=1\).

- **Treatment or exposure set to 0:**

  $$
  m(0,x)
  =
  E(Y\mid T=0,X=x)
  =
  P(Y=1\mid T=0,X=x).
  $$

  For people with the same baseline covariate profile \(X=x\), \(m(0,x)\) is the conditional mean outcome. For binary \(Y\), it is the conditional probability of the event among those with \(T=0\).

The fitted outcome model estimates these two quantities as \(\hat m(1,x)\) and \(\hat m(0,x)\).

The outcome-regression procedure uses this prediction function in a special way: it predicts each person's outcome under both treatment settings and then averages.

```{r or-roadmap-picture, echo=FALSE, fig.width=8.8, fig.height=4.3}
start_picture()

draw_box(0.11, 0.66, 0.18, 0.19,
         "Observed data\n(T, X, Y)",
         fill = "#E8E8E8", cex = 0.62, font = 2)
draw_box(0.35, 0.66, 0.20, 0.19,
         "Fit outcome rule\nm hat(T, X)",
         fill = "#A0CBE8", cex = 0.60, font = 2)
draw_box(0.62, 0.66, 0.23, 0.19,
         "Predict twice\nm hat(1, X)\nm hat(0, X)",
         fill = "#DDFBE1", cex = 0.55, font = 2)
draw_box(0.89, 0.66, 0.18, 0.19,
         "Average and\ncompare",
         fill = "#FFF6C7", cex = 0.60, font = 2)

draw_arrow(0.20, 0.66, 0.25, 0.66)
draw_arrow(0.45, 0.66, 0.50, 0.66)
draw_arrow(0.74, 0.66, 0.80, 0.66)

draw_text(0.50, 0.24,
  "Prediction supplies m hat. Causal assumptions justify the intervention interpretation.",
  cex = 0.70, font = 2)
```

## Load the Classroom NHANES Data

We use the same complete-case classroom dataset as the other lectures. The example predicts CVD from smoking status and baseline covariates. Survey weights are not used because the purpose is to explain the estimator rather than conduct a full NHANES survey analysis.

```{r load-nhanes-data}
data_candidates <- c(
  file.path("..", "Data", "NHANES_data.csv"),
  file.path("Data", "NHANES_data.csv"),
  "C:/Users/seyoo/OneDrive/GitHub_repository/Causal_Inference_Survival_Analysis/Causal_Inference_Lecture/Data/NHANES_data.csv"
)

data_path <- data_candidates[file.exists(data_candidates)][1]

if (is.na(data_path)) {
  stop("Could not find NHANES_data.csv. Check the data path.")
}

nhanes <- read.csv(data_path, stringsAsFactors = FALSE)

analysis_vars <- c(
  "cvd_indicator",
  "smoker_indicator",
  "age_yr",
  "gender",
  "race",
  "educ_lvl",
  "inc_to_pov_ratio",
  "bmi"
)

or_data <- nhanes[
  complete.cases(nhanes[, analysis_vars]),
  analysis_vars
]

or_data$cvd_indicator <- as.integer(or_data$cvd_indicator)
or_data$smoker_indicator <- as.integer(or_data$smoker_indicator)
or_data$gender <- factor(or_data$gender)
or_data$race <- factor(or_data$race)
or_data$educ_lvl <- factor(or_data$educ_lvl)

data.frame(
  Observations = nrow(or_data),
  CVD_cases = sum(or_data$cvd_indicator),
  CVD_percent = round(100 * mean(or_data$cvd_indicator), 2),
  Smoker_percent = round(100 * mean(or_data$smoker_indicator), 2)
)
```

## Step 1: Fit the Outcome Prediction Model

For the binary CVD outcome, we fit logistic regression:

$$
\operatorname{logit}\{m(T,X)\}
=\beta_0+\beta_TT+\beta_X^\top X.
$$

```{r fit-outcome-model}
outcome_model <- glm(
  cvd_indicator ~ smoker_indicator + age_yr + gender + race +
    educ_lvl + inc_to_pov_ratio + bmi,
  data = or_data,
  family = binomial()
)

model_description <- data.frame(
  Component = c("Outcome", "Exposure", "Baseline covariates", "Model output"),
  Specification = c(
    "CVD indicator (0/1)",
    "Smoking indicator (0/1)",
    "Age, gender, race, education, income-to-poverty ratio, BMI",
    "Predicted conditional CVD probability"
  )
)

knitr::kable(model_description)
```

The model describes the observed conditional CVD probability. Its individual coefficients are not the standardized causal effect; the effect-scale calculation comes after scenario prediction and averaging.

## Step 2: Predict Two Treatment Scenarios

Make two copies of the same analytic sample:

1. set smoking to \(T=1\) for everyone;
2. set smoking to \(T=0\) for everyone;
3. keep age, gender, race, education, income, and BMI unchanged.

```{r predict-two-scenarios}
data_if_smoker <- or_data
data_if_non_smoker <- or_data

data_if_smoker$smoker_indicator <- 1
data_if_non_smoker$smoker_indicator <- 0

m1_hat <- predict(
  outcome_model,
  newdata = data_if_smoker,
  type = "response"
)

m0_hat <- predict(
  outcome_model,
  newdata = data_if_non_smoker,
  type = "response"
)

scenario_examples <- data.frame(
  Age = or_data$age_yr[1:10],
  BMI = round(or_data$bmi[1:10], 1),
  Observed_smoking = or_data$smoker_indicator[1:10],
  Predicted_if_smoker = round(m1_hat[1:10], 3),
  Predicted_if_non_smoker = round(m0_hat[1:10], 3)
)

knitr::kable(scenario_examples)
```

These are two model predictions for each person. They are not two observed outcomes: only one smoking status was actually observed.

```{r two-scenario-picture, echo=FALSE, fig.width=7.5, fig.height=6.5}
scenario_plot_limit <- min(1, 1.05 * max(c(m0_hat, m1_hat)))
example_person <- which.min(abs(m0_hat - 0.30))
example_m0 <- m0_hat[example_person]
example_m1 <- m1_hat[example_person]
example_risk_contrast <- example_m1 - example_m0
old_par <- par(pty = "s")

plot(
  m0_hat,
  m1_hat,
  pch = 19,
  col = adjustcolor("#4E79A7", alpha.f = 0.35),
  xlim = c(0, scenario_plot_limit),
  ylim = c(0, scenario_plot_limit),
  xlab = expression(
    paste("Predicted CVD risk if non-smoker,  ", hat(m)(0, X[i]))
  ),
  ylab = expression(
    paste("Predicted CVD risk if smoker,  ", hat(m)(1, X[i]))
  ),
  main = "Two Predictions for Every Person"
)
abline(0, 1, lty = 2, lwd = 2, col = "#E15759")

segments(
  x0 = example_m0,
  y0 = 0,
  x1 = example_m0,
  y1 = example_m1,
  lty = 3,
  lwd = 2,
  col = "#F28E2B"
)

segments(
  x0 = 0,
  y0 = example_m1,
  x1 = example_m0,
  y1 = example_m1,
  lty = 3,
  lwd = 2,
  col = "#F28E2B"
)

points(
  example_m0,
  example_m1,
  pch = 21,
  cex = 1.8,
  lwd = 2,
  col = "#9C2F21",
  bg = "#F28E2B"
)

text(
  example_m0 - 0.01,
  example_m1 + 0.035,
  labels = sprintf(
    "Example person j: (%.2f, %.2f)",
    example_m0,
    example_m1
  ),
  pos = 2,
  offset = 0.5,
  cex = 0.82,
  font = 2,
  col = "#9C2F21"
)

par(old_par)
```

<div class="prediction-note">
**Interpreting one ordered pair.** The highlighted point is an example person \(j\). Reading the horizontal coordinate first and the vertical coordinate second, this person's ordered pair is

$$
\left(\hat m(0,X_j),\hat m(1,X_j)\right)
=
\left(
`r sprintf("%.3f", example_m0)`,
`r sprintf("%.3f", example_m1)`
\right).
$$

- The first coordinate, \(\hat m(0,X_j)=\) `r sprintf("%.3f", example_m0)`, is the predicted CVD risk when this person's smoking status is set to \(T=0\). In percentage terms, that is about `r sprintf("%.1f", 100 * example_m0)` percent.
- The second coordinate, \(\hat m(1,X_j)=\) `r sprintf("%.3f", example_m1)`, is the predicted CVD risk for the same person, with the same baseline covariates \(X_j\), when smoking is set to \(T=1\). That is about `r sprintf("%.1f", 100 * example_m1)` percent.
- The model-based scenario contrast is \(\hat m(1,X_j)-\hat m(0,X_j)=\) `r sprintf("%.3f", example_risk_contrast)`, or approximately `r sprintf("%.1f", 100 * example_risk_contrast)` percentage points higher in the smoking scenario.

The orange point is above the dashed equality line because its second coordinate is larger than its first. This ordered pair does **not** represent two outcomes observed for person \(j\), and its difference is not a directly observed individual causal effect. It contains two conditional mean predictions from the fitted model; only one smoking state and one outcome were actually observed for this person.

**Causal interpretation under the identification assumptions.** If consistency, conditional exchangeability given \(X\), and positivity hold—and the outcome model is adequately specified—the two fitted predictions can identify conditional potential-outcome risks:

$$
\hat m(0,X_j)
\approx
P\{Y(0)=1\mid X=X_j\}
=
`r sprintf("%.3f", example_m0)`,
$$

and

$$
\hat m(1,X_j)
\approx
P\{Y(1)=1\mid X=X_j\}
=
`r sprintf("%.3f", example_m1)`.
$$

Therefore,

$$
\hat m(1,X_j)-\hat m(0,X_j)
\approx
E\{Y(1)-Y(0)\mid X=X_j\}
=
`r sprintf("%.3f", example_risk_contrast)`.
$$

For people with this baseline profile, the model estimates that setting smoking to \(1\), rather than \(0\), raises CVD risk by about `r sprintf("%.1f", 100 * example_risk_contrast)` percentage points. This is a **conditional average causal risk difference** for the covariate profile \(X_j\), not the unobservable individual causal effect \(Y_j(1)-Y_j(0)\). The point lies above the equality line because this estimated conditional causal contrast is positive.
</div>

## Step 3: Standardize the Predictions

Average each set of predicted probabilities over the same covariate distribution:

$$
\hat\psi_1^{OR}
=\frac{1}{n}\sum_{i=1}^n\hat m(1,X_i),
$$

$$
\hat\psi_0^{OR}
=\frac{1}{n}\sum_{i=1}^n\hat m(0,X_i).
$$

Then compare the standardized risks:

$$
\widehat{RD}^{OR}
=\hat\psi_1^{OR}-\hat\psi_0^{OR},
$$

$$
\widehat{RR}^{OR}
=\frac{\hat\psi_1^{OR}}{\hat\psi_0^{OR}}.
$$

```{r standardize-scenario-predictions}
psi1_hat <- mean(m1_hat)
psi0_hat <- mean(m0_hat)
rd_hat <- psi1_hat - psi0_hat
rr_hat <- psi1_hat / psi0_hat

or_results <- data.frame(
  Quantity = c(
    "Standardized risk if everyone smoked",
    "Standardized risk if no one smoked",
    "Standardized risk difference",
    "Standardized risk ratio"
  ),
  Value = c(psi1_hat, psi0_hat, rd_hat, rr_hat)
)

or_results$Value <- round(or_results$Value, 4)
knitr::kable(or_results)
```

```{r standardized-risks-picture, echo=FALSE, fig.width=7.4, fig.height=5.1}
standardized_risks <- 100 * c(
  "If everyone were a non-smoker" = psi0_hat,
  "If everyone were a smoker" = psi1_hat
)

barplot(
  standardized_risks,
  col = c("#A0CBE8", "#F28E2B"),
  ylim = c(0, max(standardized_risks) * 1.25),
  ylab = "Standardized CVD risk (%)",
  main = "Outcome-Regression Standardized Risks"
)

text(
  x = c(0.7, 1.9),
  y = standardized_risks,
  labels = paste0(round(standardized_risks, 1), "%"),
  pos = 3,
  font = 2
)
```

The key operation is **average after predicting**. We do not classify people as CVD or no CVD. Outcome regression retains the predicted probabilities and averages them.

## Ordinary Prediction Versus Outcome Regression

| Ordinary outcome prediction | Outcome-regression standardization |
|:--|:--|
| Insert each person's observed \(T\) and \(X\) | Set \(T=1\) and \(T=0\) for every person |
| Produce one fitted risk per person | Produce two scenario-specific risks per person |
| Evaluate prediction of observed \(Y\) | Average each scenario over the same \(X\) distribution |
| Target observed conditional risk | Target marginal potential-outcome risk under assumptions |

The fitted function is a prediction model in both columns. The data operation and interpretation differ.

## Why Good Prediction Is Not Enough for Causality

Prediction quality cannot establish that a contrast is causal. A model may predict observed CVD well and still produce a biased intervention comparison.

For a causal interpretation, outcome regression requires:

1. **Consistency:** the observed outcome under the observed treatment equals the corresponding potential outcome.
2. **Conditional exchangeability:** after conditioning on \(X\), treatment is independent of the potential outcomes.
3. **Positivity:** each relevant covariate pattern has a positive probability of each treatment.
4. **Adequate outcome modeling:** \(\hat m(t,X)\) represents the relevant conditional outcome risks sufficiently well.

<div class="warning-note">
Variables that improve prediction are not automatically appropriate causal adjustment variables. A variable caused by treatment may predict the outcome well but create post-treatment adjustment bias. A collider may also improve prediction while distorting a causal comparison.
</div>

## Continuous Outcomes in Outcome Regression

The same standardization logic works when \(Y\) is continuous. For example, if \(Y\) were blood pressure, a regression model would estimate

$$
m(t,x)=E(Y\mid T=t,X=x).
$$

We would predict a blood-pressure value under \(T=1\) and \(T=0\) for every person and average:

$$
\hat\psi_t^{OR}
=\frac{1}{n}\sum_{i=1}^n\hat m(t,X_i).
$$

The only change is the scale of the prediction:

- binary outcome: predicted probability;
- continuous outcome: predicted conditional mean.

# Summary

## Prediction Checklist

Before interpreting a prediction analysis, ask:

1. What is the outcome?
2. Is it continuous, binary, or multicategory?
3. What information is available at prediction time?
4. What quantity does the model output?
5. Was performance evaluated outside the training data?
6. Which metric matches the outcome and intended use?
7. Is the model calibrated as well as discriminative?
8. If categories are assigned, how was the decision threshold chosen?

For causal outcome regression, also ask:

9. Are the adjustment variables appropriate baseline confounders?
10. Are consistency, exchangeability, and positivity plausible?
11. Does the model predict reasonably under both treatment settings?
12. Does the uncertainty calculation repeat the full fitting and standardization procedure?

## Key Takeaways

1. Prediction learns a rule that maps available information to an unknown outcome.
2. Continuous-outcome prediction returns a number on the outcome scale and is commonly evaluated with MAE, RMSE, and \(R^2\).
3. Categorical-outcome prediction returns category probabilities and is evaluated with probability, calibration, discrimination, and classification measures.
4. Probability prediction and category assignment are separate steps.
5. Test data or cross-validation are needed to evaluate generalization.
6. Outcome regression uses prediction as a building block: predict under both treatment settings, average, and compare.
7. Predictive performance alone does not give a causal interpretation; causal assumptions and appropriate variable selection are still required.

## Appendix: Connection to the G-Formula

Under consistency, conditional exchangeability, and positivity,

$$
E\{Y(t)\}
=E_X[E(Y\mid T=t,X)].
$$

The inner expectation is the prediction function

$$
m(t,X)=E(Y\mid T=t,X),
$$

and the outer expectation standardizes it over the covariate distribution. Its sample analogue is

$$
\hat\psi_t^{OR}
=\frac{1}{n}\sum_{i=1}^n\hat m(t,X_i).
$$

Outcome-regression standardization is therefore also called parametric g-computation. Prediction provides \(\hat m\); identification assumptions connect the observed conditional mean to the potential-outcome mean.

## References

Harrell, F. E. (2015). *Regression Modeling Strategies* (2nd ed.). Springer.

Hernan, M. A., & Robins, J. M. (2020). *Causal Inference: What If*. Chapman & Hall/CRC.

James, G., Witten, D., Hastie, T., Tibshirani, R., & Taylor, J. (2023). *An Introduction to Statistical Learning* (2nd ed.). Springer.

Steyerberg, E. W. (2019). *Clinical Prediction Models* (2nd ed.). Springer.