---
title: 'Bias and Variance of OLS'
author: "Ramesh Johari"
output:
  pdf_document: default
---

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

### The setup

Consider a population model with outcome variable $Y$ and ten covariates, generated as follows:

* $Y$ is normally distributed with mean zero and variance 1;
* Each of the ten covariates $X_1, \ldots, X_{10}$ are related to $Y$ according to:
\[ X_j = Y + \epsilon_j, \]
where $\epsilon_j$ is normally distributed with mean zero and variance 1; and
* $Y$ and all $\epsilon_1, \ldots, \epsilon_{10}$ are mutually independent.

Suppose we have a dataset that consists of $n$ independent samples from this population.  We assume that the appropriate conditions are satisfied for ordinary least squares regression on this dataset to have a unique solution.

We compare two approaches to model fitting via OLS:

1. (Approach A) OLS using only the first covariate;
2. (Approach B) OLS using all ten covariates.

Our goal is to understand the relative bias and variance of these two OLS approaches.

Regarding bias: Removing covariates that are present in the true population model introduces bias, so **strategy A has nonzero bias relative to B**; that part is unambiguous.

The rest of this notebook shows that **strategy A also has higher variance than B**.

### Technical details

We need the joint distribution of $X_1, \ldots, X_{10}, Y$.  Given the setup, this joint distribution is multivariate normal, with the following properties:

* $Y$ has mean 0 and variance 1.
* Each $X_j$ has mean 0 and variance 2.
* The covariance of $Y$ with each $X_j$ is 1.
* The covariance of $X_j$ with $X_k$ for $j \neq k$ is 1.

We use the preceding observations and standard properties of multivariate normal random vectors to find that the conditional distribution of $Y$ given $X_1, \ldots, X_{10}$ is normal, with mean:
\[ f^*(\vec{X}) = E[Y | X_1, \ldots, X_{10}] = \frac{X_1}{11} + \cdots + \frac{X_{10}}{11}. \]
Further, if we define the error as:
\[ \epsilon = Y - f^*(\vec{X}), \]
then we can show that:
\[ \epsilon = \frac{Y}{11} - \frac{\epsilon_1}{11} - \cdots - \frac{\epsilon_{10}}{11},\]
and from this calculation we can show $\epsilon$ is normally distributed with mean 0 and variance $1/11$, and also that $\epsilon$ has zero covariance with every $X_j$.  Note the variance of $\epsilon$ doesn't depend on $\vec{X}$, i.e., it is homoskedastic.

To summarize, we have shown the conditional distribution of $Y$ given $\vec{X}$ is linear in $\vec{X}$, with error that has zero mean and variance $1/11$, and where the error is uncorrelated with the covariates.  This is an alternative way to specify the population model, and importantly, illustrates that this model satisfies the probabilistic assumptions (A1)-(A3) in the Gauss-Markov Theorem.

### Visualizing bias and variance

We now present some R code to visualize bias and variance.  We draw many "parallel universes" of training data, and make predictions at a fixed point (in this case, the point $X_1 = \cdots = X_{10} = 1$).  We compare our prediction to the conditional expectation of $Y$ given $X_1 = 1, \ldots, X_{10} = 1$, which from the formula above is $f^*(1, \ldots, 1) = 10/11$.  We plot the resulting errors across universes.

```{r view_i}
# Set parameters
num_universes <- 10000
num_samples <- 100

# Initialize vectors to store errors
errors_A <- numeric(num_universes)
errors_B <- numeric(num_universes)

# Simulation loop
for (i in 1:num_universes) {
  # Generate Y
  Y <- rnorm(num_samples, mean = 0, sd = 1)
  
  # Generate Xs
  X <- matrix(nrow = num_samples, ncol = 10)
  for (j in 1:10) {
    epsilon_j <- rnorm(num_samples, mean = 0, sd = 1)
    X[, j] <- Y + epsilon_j
  }
  
  # Linear regression using only X1
  model_A <- lm(Y ~ 0 + X[, 1])
  
  # Linear regression using all Xs
  model_B <- lm(Y ~ 0 + X)
  
  # Prediction at the point X1 = X2 = ... = X10 = 1
  prediction_A <- sum(coef(model_A))
  prediction_B <- sum(coef(model_B))
  
  # Calculate errors
  errors_A[i] <- prediction_A - 10/11 # Correct conditional mean of Y is 10/11 when X is all 1's
  errors_B[i] <- prediction_B - 10/11  
}

# Create data frames for plotting
A_errors_df <- data.frame(Error = errors_A)
B_errors_df <- data.frame(Error = errors_B)

# Add a new column to each data frame to distinguish between them
A_errors_df$type <- 'A: Single Covariate'
B_errors_df$type <- 'B: Multiple Covariates'

# Combine the two data frames
combined_errors_df <- rbind(A_errors_df, B_errors_df)

# Create one density plot with different colors for each type
combined_density_plot <- ggplot(combined_errors_df, aes(x = Error, fill = type)) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(values = c('blue', 'red')) +
  ggtitle('Density of Prediction Errors') +
  xlab('Prediction Error') +
  ylab('Density') +
  scale_x_continuous(limits = c(-0.75,0.25)) +
  theme_minimal()
print(combined_density_plot)
```

From the resulting plot, it is clear that model A has both nonzero bias and higher variance than model B.  (The same would be true even if we made predictions at a point other than $(1,\ldots,1)$.)  This happens because ignoring the additional $X$'s eliminates their contribution to a more precise estimate of $Y$.