Sums of Squares and Hypothesis Testing

Problem 1

Given the following data:

\[\mathbf{y} = \begin{bmatrix}10\\13\\14\\16\\18\end{bmatrix}\]

\[\mathbf{X} = \begin{bmatrix}1&1\\1&2\\1&3\\1&4\\1&5\end{bmatrix}\]

a) Calculate TSS by hand using both the definition and the matrix form.

b) Fit the linear model and calculate the SSE and \(\textrm{SS}_\textrm{Reg}\). Verify that TSS = SSE + SS_Reg.

c) Calculate R² and interpret its meaning in context.

Problem 2

Prove algebraically that the F-statistic can be written in terms of R² as:

\(F = \frac{R^2/(p-1)}{(1-R^2)/(n-p)}\)

Show each algebraic step clearly.

Problem 3

In this problem, you’ll verify through simulation that \(\frac{\hat{\beta}_j - \beta_j}{\text{se}(\hat{\beta}_j)} \sim t_{n-p}\).

Introduction to purrr

Before we start, we’ll use the purrr package, which provides tools for functional programming in R. The key function we’ll use is map_dbl(), which applies a function to each element of a list or vector and outputs a numeric vector. Think of map_dbl() as a more elegant alternative to writing loops when you expect a numeric output.

For example: - map_dbl(1:5, sqrt) applies the square root function to numbers 1 through 5 - map_dbl(1:1000, my_simulation_function) runs your simulation function 1000 times

This approach is cleaner than loops and encourages you to think about your simulation as a function that gets repeated many times.

Step-by-step simulation guide:

Step 1: Set up simulation parameters and load libraries

library(purrr)
library(ggplot2)
set.seed(1)

# Simulation parameters
n <- 50          # sample size
p <- 3           # number of parameters (intercept + 2 predictors)
n_sims <- 10000  # number of simulations
beta_true <- c(2, 1.5, -0.8)  # true coefficients [intercept, x1, x2]
sigma_true <- 2  # true error standard deviation

# Create design matrix (stays fixed across simulations)
X <- cbind(1, rnorm(n, 0, 1), rnorm(n, 2, 1.5))

Step 2: Create a simulation function

Fill in the function below. This function should: 1. Generate random errors 2. Create the response variable using the true model 3. Fit the linear model 4. Extract the coefficient and standard error for x1 (the second coefficient) 5. Calculate and return the t-statistic

simulate_t_stat <- function(sim_number) {
  # Generate random errors with a mean = 0 and sd = 1 
  epsilon <- rnorm(n, 0, 1)
  
  # Create response variable: y = X * beta_true + epsilon
  y <- # YOUR CODE HERE
  
  # Calculate coefficient and standard error for x1
  beta_hat <- # YOUR CODE HERE (coefficient for x1)
  se_beta <- # YOUR CODE HERE (standard error for x1)
  
  # Calculate t-statistic: (beta_hat - true_beta) / se_beta
  t_stat <- # YOUR CODE HERE
  
  return(t_stat)
}

Step 3: Run simulation using purrr::map_dbl()

# Use map to run simulation n_sims times
t_stats <- map_dbl(1:n_sims, simulate_t_stat)

Step 4: Analyze results

# Calculate degrees of freedom
dft <- n - p

# Create data frame for ggplot
sim_data <- data.frame(t_statistic = t_stats)

# Create comparison plot
ggplot(sim_data, aes(x = t_statistic)) +
  geom_histogram(aes(y = after_stat(density)), bins = 50, fill = "lightblue") +
  stat_function(fun = dt, args = list(df = dft), color = "red", linewidth = 1.2) +
  labs(title = 
         paste("Simulated t-statistics vs Theoretical t-distribution (df =", dft, ")"),
       subtitle = paste("Based on", n_sims, "simulations"),
       x = "t-statistic", 
       y = "Density") +
  theme_minimal()

Questions to answer:

• What are the degrees of freedom for this t-distribution?
  
• Do the simulated t-statistics follow the expected distribution?

Problem 4

Now verify that the overall F-statistic follows an F-distribution under the null hypothesis.

Simulation setup:

Step 1: Set up parameters for testing under the null

set.seed(1)

# Simulation parameters
n <- 30
p <- 4  # intercept + 3 predictors
n_sims <- 5000
sigma_true <- 1.5

# Under null hypothesis: all slope coefficients are 0
beta_null <- c(5, 0, 0, 0)  # only intercept is non-zero

# Create design matrix (fixed across simulations)
X <- cbind(1,
           rnorm(n, 0, 1),
           rnorm(n, 1, 2),
           rnorm(n, -1, 1.2))

Step 2: Create your F-statistic simulation function

Fill in the function below. This function should: 1. Generate random errors 2. Create response variable under the null hypothesis 3. Fit the linear model 4. Extract the overall F-statistic 5. Return the F-statistic

simulate_f_stat <- function(sim_number) {
  # Generate random errors
  epsilon <- # YOUR CODE HERE
  
  # Create response under null: y = X * beta_null + epsilon
  y <- # YOUR CODE HERE
  
  # Fit the model (be careful about how you specify the predictors)
  fit <- # YOUR CODE HERE
  
  # Extract overall F-statistic from model summary
  # Hint: summary(fit)$fstatistic[1] gives the F-statistic
  f_stat <- # YOUR CODE HERE
  
  return(f_stat)
}

Step 3: Run simulation using purrr

# Use map_dbl to get a numeric vector of F-statistics
f_stats <- # YOUR CODE HERE (use map_dbl and your function)

Step 4: Compare to theoretical F-distribution with improved visualizations

# Calculate degrees of freedom
df1 <- p - 1  # numerator degrees of freedom (number of predictors excluding intercept)
df2 <- n - p  # denominator degrees of freedom

# Create data frame for ggplot
sim_data_f <- data.frame(f_statistic = f_stats)

# Create comparison plot
ggplot(sim_data_f, aes(x = f_statistic)) +
  geom_histogram(aes(y = after_stat(density)), bins = 50, fill = "lightgreen") +
  stat_function(fun = df, args = list(df1 = df1, df2 = df2), color = "red", linewidth = 1.2) +
  xlim(0, quantile(f_stats, 0.99)) +  # Trim extreme values 
  labs(title = paste("Simulated F-statistics vs Theoretical F-distribution"),
       subtitle = paste("df1 =", df1, ", df2 =", df2, "| Based on", n_sims, "simulations"),
       x = "F-statistic", 
       y = "Density") +
  theme_minimal()

Questions to answer: - What are the numerator and denominator degrees of freedom?
- Do the simulated F-statistics match the theoretical distribution?
- What would happen if the null hypothesis were false (i.e., if some slope coefficients were non-zero)?

Problem 5

Consider the model: \(y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + \varepsilon\)

a) Set up the contrast matrix \(\mathbf{C}\) and vector \(\mathbf{d}\) to test the hypothesis: \(H_0: \beta_1 + \beta_2 = 0 \text{ and } \beta_3 = 2\beta_4\)

b) How many degrees of freedom would the F-test have? Explain your reasoning.

Problem 6

You are given the following information from a regression analysis:

• n = 25 (sample size)
  
• p = 5 (number of parameters including intercept)
  
• TSS = 500
  
• R² = 0.72

Complete the ANOVA table:

Source	df	Sum of Squares	Mean Square	F
Regression
Error
Total