Distribution of the sample slope

1 Loading packages and creating in the data set

We’ll start by loading the tidyverse:

Since the variance of the slope, \(b_1\), depends on \(S_{XX}\):

\[S_{XX} = \sum_{i=1}^n(X_i-\bar{X})^2\]

We need to keep the \(X_i\) the same throughout the simulation for each iteration. We’ll keep it simple and let \(x_i = 2i\) with \(i = 1, ..., 10\)

X <- 2*(1:10)

X

##  [1]  2  4  6  8 10 12 14 16 18 20

The \(S_{XX}\) is:

# S_XX = sum (X_i - X_bar)^2
S_XX <- sum((X - mean(X))^2)
S_XX

## [1] 330

When performing a simulation, we need to create a model to generate our data. The population model for our simulation is:

\[Y_i = 20 + 40X_i + \varepsilon_i \\ \varepsilon_i \sim N(0, 660)\]

For each simulation we run, we’ll generate \(Y_i\) using the model above

2 Theoretical distribution for \(b_1\)

According to the slides we saw in class, as long as

\[Y_i \sim N(\beta_0 + \beta_1 X_i, \sigma^2)\]

are are independent, then the distribution of \(b_1\) is:

\[b_1 \sim N\left(\beta_1, \frac{\sigma^2}{S_{XX}}\right)\]

Since we control the population, we know the distribution \(b_1\) for our example should follow:

\[b_1 \sim N\left(40, \frac{660}{330}\right)\]

3 Simulating 10,000 samples of n = 10

We’ll generate 10,000 samples of 10 \(Y_i\) values each based on the model stated above, which is theoretically akin to taking a random sample from a population.

From each sample, we’ll save the following info:

1. The sample
2. \(b_0\)
3. \(b_1\)
4. \(SSE\)

# Simulation and sample sizes
m <- 1e5
sample_size <- 10

# Setting the parameters of the model
beta0 <- 20; beta1 <- 40; sigma2_e <- 660

# Data frame to save the result:
sim_df <- 
  data.frame(
    sample = 1:m,
    b0 = -100,
    b1 = -100,
    SSE = -1
  )

Before we conduct the full simulation, let’s generate one sample and calculate the stats needed as an example:

### Let's generate a random sample for sake of showing how:
set.seed(4210)
# The response from the model:
Y_loop <- beta0 + beta1 * X + 
          # Random error term
          rnorm(sample_size, mean = 0, sd = sqrt(sigma2_e))

# Find b1:
b1_loop <- 
  # S_XY
  sum((X - mean(X))*(Y_loop - mean(Y_loop))) / S_XX

# Finding b0:
b0_loop <- mean(Y_loop) - b1_loop * mean(X)

# Calculating y_hat
Y_hat_loop <- b0_loop + b1_loop * X

# Residuals:
res_loop <- Y_loop - Y_hat_loop

# SSE:
SSE_loop <- sum(res_loop^2)

# Now we'll save the results in the first row
sim_df[1, 2:ncol(sim_df)] <- c(b0_loop, b1_loop, SSE_loop)


tibble(sim_df)

## # A tibble: 100,000 × 4
##    sample     b0     b1   SSE
##     <int>  <dbl>  <dbl> <dbl>
##  1      1   23.5   39.8 1657.
##  2      2 -100   -100     -1 
##  3      3 -100   -100     -1 
##  4      4 -100   -100     -1 
##  5      5 -100   -100     -1 
##  6      6 -100   -100     -1 
##  7      7 -100   -100     -1 
##  8      8 -100   -100     -1 
##  9      9 -100   -100     -1 
## 10     10 -100   -100     -1 
## # ℹ 99,990 more rows

Now we’ll add the remaining sample stats to the data frame:

# Performing the simulation:
for (i in 2:m){
  # The response from the model:
  Y_loop <- beta0 + beta1 * X + 
            # Random error term
            rnorm(sample_size, mean = 0, sd = sqrt(sigma2_e))
  
  # Find b1:
  b1_loop <- sum((X - mean(X))*(Y_loop - mean(Y_loop))) / S_XX
  
  # Finding b0:
  b0_loop <- mean(Y_loop) - b1_loop * mean(X)
  
  # Calculating y_hat
  Y_hat_loop <- b0_loop + b1_loop * X
  
  # Residuals:
  res_loop <- Y_loop - Y_hat_loop
  
  # SSE:
  SSE_loop <- sum(res_loop^2)
  
  # Now we'll save the results in the first row
  sim_df[i, 2:ncol(sim_df)] <- c(b0_loop, b1_loop, SSE_loop)
}

tibble(sim_df)

## # A tibble: 100,000 × 4
##    sample    b0    b1   SSE
##     <int> <dbl> <dbl> <dbl>
##  1      1 23.5   39.8 1657.
##  2      2 16.1   41.0 3318.
##  3      3 28.9   38.7 5290.
##  4      4 57.0   37.9 2793.
##  5      5 24.6   40.2 2391.
##  6      6 19.8   39.6 4262.
##  7      7 13.7   40.4 2323.
##  8      8 40.0   38.4 3262.
##  9      9 10.3   40.1 3662.
## 10     10 -6.51  41.2 4166.
## # ℹ 99,990 more rows

4 Distribution of \(b_1\)

Let’s create a histogram of our simulated slopes and superimpose a Normal curve on it.

Reminder: If the theory is correct, then:

\[b_1 \sim N\left(40, \frac{660}{330}\right)\]

ggplot(
  data = sim_df,
  mapping = aes(
    x = b1
  )
) + 
  geom_histogram(
    # Need to swap from count to density on the y-axis to place a normal curve
    mapping = aes(y = after_stat(density)),
    bins = 40,
    fill = 'steelblue',
    color = 'white'
  ) + 
  # Adding the Normal curve:
  stat_function(
    fun = dnorm,
    args = list(mean = beta1, sd = sqrt(sigma2_e/S_XX)),
    color = 'red',
    linewidth = 1,
    alpha = 0.5
  ) + 
  theme_bw()  +
  labs(
    title = 'Histogram of simulated slopes'
  ) + 
  scale_y_continuous(
    expand = c(0, 0, 0.05, 0)
  )

Histogram looks pretty accurate!

What about the test statistics and p-values? Do they ‘line up’ with what the theory says?

Result: If \(H_0\) is true, then the p-value is a random variable that follows a uniform distribution between 0 and 1:

\[\text{p-val} | H_0 \sim \text{Uniform}(0, 1)\]

4.1 Type I error rate: \(\sigma^2\) known

So if we calculate the test statistic and p-value for all of our simulated samples, then 5% of our samples should give us a p-value below 0.05 (or whatever significance level you choose)

\[P(\text{p-val} < \alpha) = \alpha\]

# Finding the p-values:
sim_df <- 
  sim_df |>
  mutate(
    # Test statistic
    test_stat = (b1 - beta1) / sqrt(sigma2_e/S_XX),
    # Now p-value for Beta1 =\= 40
    p_val_sigma = 2*pnorm(abs(test_stat), lower.tail = F)
  )

sim_df |>
  dplyr::select(sample, b1, test_stat, p_val_sigma) |>
  tibble()

## # A tibble: 100,000 × 4
##    sample    b1 test_stat p_val_sigma
##     <int> <dbl>     <dbl>       <dbl>
##  1      1  39.8   -0.138        0.890
##  2      2  41.0    0.693        0.488
##  3      3  38.7   -0.937        0.349
##  4      4  37.9   -1.50         0.134
##  5      5  40.2    0.164        0.870
##  6      6  39.6   -0.259        0.795
##  7      7  40.4    0.255        0.799
##  8      8  38.4   -1.17         0.244
##  9      9  40.1    0.0506       0.960
## 10     10  41.2    0.849        0.396
## # ℹ 99,990 more rows

Let’s plot them and calculate the probability of rejecting a true \(H_0\):

# Calculating the probability of a type 1 error
prob_type1_alpha_5 <- mean(sim_df$p_val_sigma < 0.05)

# How wide to make each bin
binwidth = 0.01

ggplot(
  data = sim_df,
  mapping = aes(
    x = p_val_sigma,
    fill = if_else(p_val_sigma < 0.05, 'tomato', 'steelblue')
  )
) + 
  geom_histogram(
    # Need to swap from count to density on the y-axis to place a normal curve
    mapping = aes(y = after_stat(count) / (m * binwidth)),
    breaks = seq(0, 1, binwidth),
    #fill = 'steelblue',
    color = 'white'
  ) + 
  # Adding the Uniform 'curve':
  stat_function(
    fun = dunif,
    args = list(min = 0, max = 1),
    color = 'red',
    linewidth = 1,
    alpha = 0.5
  ) + 
  # Adding the probability of Type 1
  labs(
    title = 'Distribution of p-values when Ho is true and sigma is known',
    subtitle = paste('When alpha = 5% the chance of Type I error is', 
                     scales::percent(prob_type1_alpha_5)),
    caption = paste('Results from a simulation of', scales::comma(m), 
                    'simulations each with a sample size of', sample_size),
    x = 'P-value',
    y = 'Density'
  ) +
  theme_bw()  +
  scale_fill_identity() + 
  scale_y_continuous(
    expand = c(0, 0, 0.05, 0)
  )

The distribution above assumes we know what \(\sigma^2\) is, which we usually don’t. What happens if we swap \(\sigma^2\) with the MSE?

4.2 Type I error rate: \(\sigma^2\) unknown

Instead of using the standard deviation of \(b_1\), we’ll have to use the standard error:

\[SE(b_1) = \sqrt{\frac{MSE}{S_{XX}}}\]

We can calculate the \(MSE\) for each sample by:

\[MSE = \frac{SSE}{n-2}\]

So let’s find the Standard error, test stat, and p-value using the \(MSE\):

\[z = \frac{b_1 - \beta_{1,0}}{\sqrt{\frac{MSE}{S_{XX}}}}\]

# Finding the p-values:
sim_df <- 
  sim_df |>
  mutate(
    MSE = SSE / (sample_size - 2),
    # Test statistic using mse
    test_stat_mse = (b1 - beta1) / sqrt(MSE/S_XX),
    # Now p-value for Beta1 =\= 40
    p_val_mse = 2*pnorm(abs(test_stat_mse), lower.tail = F)
  )

sim_df |>
  dplyr::select(sample, b1, MSE, test_stat_mse, p_val_mse) |>
  tibble()

## # A tibble: 100,000 × 5
##    sample    b1   MSE test_stat_mse p_val_mse
##     <int> <dbl> <dbl>         <dbl>     <dbl>
##  1      1  39.8  207.       -0.246     0.806 
##  2      2  41.0  415.        0.874     0.382 
##  3      3  38.7  661.       -0.936     0.349 
##  4      4  37.9  349.       -2.06      0.0396
##  5      5  40.2  299.        0.244     0.807 
##  6      6  39.6  533.       -0.289     0.773 
##  7      7  40.4  290.        0.384     0.701 
##  8      8  38.4  408.       -1.48      0.138 
##  9      9  40.1  458.        0.0607    0.952 
## 10     10  41.2  521.        0.956     0.339 
## # ℹ 99,990 more rows

4.2.1 Plot of the test statistic using MSE instead of \(\sigma^2\)

The histogram of the test statistic is shown below. Does it look like the Standard Normal distribution fit the histogram?

gg_test_stat_mse_hist <- 
  ggplot(
    # There are a few extreme values, so we'll only look at -5 < z < 5
    data = sim_df |>
      filter(between(test_stat_mse, -5, 5)),
    mapping = aes(
      x = test_stat_mse
    )
  ) + 
  geom_histogram(
    # Need to swap from count to density on the y-axis to place a normal curve
    mapping = aes(y = after_stat(density)),
    bins = 40,
    fill = 'steelblue',
    color = 'white'
  ) + 
  theme_bw()  +
  labs(
    title = paste('Histogram of', scales::comma(m), 'test statistics for the slope'),
    subtitle = 'Calculated using the MSE',
    x = 'Test statistics'
  ) + 
  scale_y_continuous(
    expand = c(0, 0, 0.05, 0)
  ) + 
  scale_x_continuous(
    breaks = -5:5
  )


gg_test_stat_mse_hist + 
  # Adding the Normal curve:
  stat_function(
    fun = dnorm,
    args = list(mean = 0, sd = 1),
    color = 'red',
    linewidth = 1,
    alpha = 0.5
  ) + 
  labs(
    caption = 'Standard Normal Distribution Overlayed on the Histogram'
  )

The Normal curve looks close, but the middle bars are consistently below the line while the tails are above the line. That typically indicates that a wider distribution is needed instead of a Normal distribution.

Let’s check the p-values as well!

4.2.2 Plot of the p-values using MSE and Normal distribution

Again, plotting them and calculate the probability of rejecting a true \(H_0\) assuming the test statistic follows a Standard Normal:

# Calculating the probability of a type 1 error
prob_type1_alpha_5_mse <- mean(sim_df$p_val_mse < 0.05)

# Creating the graph
ggplot(
  data = sim_df,
  mapping = aes(
    x = p_val_mse,
    fill = if_else(p_val_mse < 0.05, 'tomato', 'steelblue'),
   # group = 1
  )
) + 
  geom_histogram(
    # Need to swap from count to density on the y-axis to place a normal curve
    mapping = aes(y = after_stat(count) / (m * binwidth)),
    breaks = seq(0, 1, binwidth),
    color = 'white'
  ) + 
  # Adding the Uniform 'curve':
  stat_function(
    fun = dunif,
    args = list(min = 0, max = 1),
    color = 'red',
    linewidth = 1,
    alpha = 0.5
  ) + 
  # Adding the probability of Type 1
  labs(
    title = 'Distribution of p-values when H0 is true, using a Normal with sigma unknown',
    subtitle = paste('alpha = 5% and Type I probability is', 
                     scales::percent(prob_type1_alpha_5_mse)),
    caption = paste('Results from a simulation of', scales::comma(m), 
                    'simulations each with a sample size of', sample_size),
    x = 'P-value',
    y = 'Density'
  ) +
  theme_bw()  +
  scale_fill_identity() + 
  scale_y_continuous(
    expand = c(0, 0, 0.05, 0)
  )

Our chance of a Type I error is almost twice the significance level! Why?

4.3 P-value using a t-distribution

Like what you did waaaay back in your introductory stats class, when you replace \(\sigma^2\) with its sample statistic equivalent, the test statistic is no longer Normally distributed.

Instead, the test stat will follow a **t-distribution* with degrees of freedom equal to the denominator of the formula to calculate \(\hat{\sigma}^2\):

\[\hat{\sigma}^2 = \sqrt{\frac{SSE}{n - 2}}\]

So if we use \(MSE\) in place of \(\sigma^2\), then the test statistic will follow a t-distribution with \(n - 2\) degrees of freedom!

4.3.1 Histogram of the test statistic with a t-distribution

gg_test_stat_mse_hist + 
  # Adding the Normal curve:
  stat_function(
    fun = dt,
    args = list(df = sample_size - 2),
    color = 'red',
    linewidth = 1,
    alpha = 0.5
  ) + 
  labs(
    caption = paste('t-distribution with', sample_size - 2, 
                    'df overlayed on the histogram')
  )

The t-distribution fits great!

4.3.2 Histogram of p-values using MSE and t-distribution

The code chunk below will calculate the p-value using the t-distribution with df = 8:

# Finding the p-values:
sim_df <- 
  sim_df |>
  mutate(
    # Now p-value for Beta1 =\= 40 using a t-distribution
    t_p_val = 2*pt(abs(test_stat_mse), df = sample_size - 2, lower.tail = F)
  )

sim_df |>
  dplyr::select(sample, b1, MSE, test_stat_mse, t_p_val) |>
  tibble()

## # A tibble: 100,000 × 5
##    sample    b1   MSE test_stat_mse t_p_val
##     <int> <dbl> <dbl>         <dbl>   <dbl>
##  1      1  39.8  207.       -0.246   0.812 
##  2      2  41.0  415.        0.874   0.407 
##  3      3  38.7  661.       -0.936   0.377 
##  4      4  37.9  349.       -2.06    0.0735
##  5      5  40.2  299.        0.244   0.813 
##  6      6  39.6  533.       -0.289   0.780 
##  7      7  40.4  290.        0.384   0.711 
##  8      8  38.4  408.       -1.48    0.176 
##  9      9  40.1  458.        0.0607  0.953 
## 10     10  41.2  521.        0.956   0.367 
## # ℹ 99,990 more rows

Again, plotting the p-values calculated above and calculate the probability of rejecting a true \(H_0\) assuming the test statistic follows a t-distribution:

# Calculating the probability of a type 1 error
prob_type1_alpha_5_t <- mean(sim_df$t_p_val < 0.05)

# Graph
ggplot(
  data = sim_df,
  mapping = aes(
    x = t_p_val,
    fill = if_else(t_p_val < 0.05, 'tomato', 'steelblue'),
   # group = 1
  )
) + 
  geom_histogram(
    # Need to swap from count to density on the y-axis to place a normal curve
    mapping = aes(y = after_stat(count) / (m * binwidth)),
    breaks = seq(0, 1, binwidth),
    color = 'white'
  ) + 
  # Adding the Uniform 'curve':
  stat_function(
    fun = dunif,
    args = list(min = 0, max = 1),
    color = 'red',
    linewidth = 1,
    alpha = 0.5
  ) + 
  # Adding the probability of Type 1
  labs(
    title = 'Distribution of p-values when Ho is true, using MSE, and t-distribution',
    subtitle = paste('alpha = 5% and Type I probability is', 
                     scales::percent(prob_type1_alpha_5_t)),
    caption = paste('Results from a simulation of', scales::comma(m), 
                    'simulations each with a sample size of', sample_size),
    x = 'P-value',
    y = 'Density'
  ) +
  theme_bw()  +
  scale_fill_identity() + 
  scale_y_continuous(
    expand = c(0, 0, 0.05, 0)
  )

Now it matches!