High Noise

set.seed(1)
high_noise_tbl <- sim_fcn(noise = 2)

high_noise_tbl %>% plt_fcn()

high_noise_tbl %>% inf_fcn()

## # A tibble: 11 x 6
##    term        estimate std.error statistic p.value `significant?`
##    <chr>          <dbl>     <dbl>     <dbl>   <dbl> <chr>         
##  1 (Intercept)    -0.15      0.1      -1.42    0.16 nope          
##  2 x_1             0.91      0.12      7.42    0    heck_yes      
##  3 x_2             0.05      0.05      0.89    0.37 nope          
##  4 x_3            -0.03      0.1      -0.33    0.74 nope          
##  5 x_4             0.04      0.11      0.34    0.74 nope          
##  6 x_5            -0.03      0.09     -0.36    0.72 nope          
##  7 x_6            -0.01      0.11     -0.1     0.92 nope          
##  8 x_7            -0.01      0.1      -0.09    0.93 nope          
##  9 x_8            -0.03      0.1      -0.34    0.74 nope          
## 10 x_9            -0.13      0.1      -1.3     0.2  nope          
## 11 x_10            0.01      0.1       0.06    0.95 nope

Medium Noise

set.seed(2)
med_noise_tbl <- sim_fcn(noise = 1)

med_noise_tbl %>% plt_fcn()

med_noise_tbl %>% inf_fcn()

## # A tibble: 11 x 6
##    term        estimate std.error statistic p.value `significant?`
##    <chr>          <dbl>     <dbl>     <dbl>   <dbl> <chr>         
##  1 (Intercept)    -0.07      0.11     -0.67    0.51 nope          
##  2 x_1             0.84      0.13      6.55    0    heck_yes      
##  3 x_2             0.12      0.1       1.15    0.25 nope          
##  4 x_3             0.07      0.09      0.7     0.48 nope          
##  5 x_4             0.19      0.11      1.72    0.09 nope          
##  6 x_5             0.06      0.1       0.59    0.56 nope          
##  7 x_6             0.02      0.1       0.21    0.83 nope          
##  8 x_7            -0.05      0.1      -0.49    0.62 nope          
##  9 x_8             0.19      0.1       1.86    0.07 nope          
## 10 x_9             0.27      0.1       2.77    0.01 heck_yes      
## 11 x_10           -0.04      0.1      -0.42    0.68 nope

Low Noise

set.seed(3)
low_noise_tbl <- sim_fcn(noise = .1)

low_noise_tbl %>% plt_fcn()

#low_noise_tbl %>% inf_fcn()

Even in the low noise (large correlation) condition, the model isn’t tempted to see x_2 as important. So is there any cost to pay for leaving all correlated variables in the model? The astute observer might point out that the std. error of the estimate for \(β_1\) seems to trend upwards when the noise is low. Does this mean we are more uncertain about the estimate of the true predictor \(x_1\) when it is correlated with another predictor like \(x_2\)? Or was that just due to chance and the nature of the random draws in the simulation? After all, we only built 1 model for each scenario.

To explore this further we need a bigger experiment - or set of experiments - that looks at a range of correlation strengths between x_1 and x_2 and repeats each setting many times. This will also provide a nice opportunity to practice with nested lists, a workflow that is challenging (for me) but also elegant and powerful. Here’s how we could go about it:

• Set up a preliminary tibble as before that has columns for the betas and the x’s
• Set up a separate tibble with a list of the desired “noise”; this value is effectively the strength of the correlation between x_1 and x_2. When noise is 0, the two variables are perfectly correlated.
• Use crossing() to generate all combinations of noise with the beta’s and x’s

Set up the simulation parameters

prelim_tbl <- betas_tbl %>%
  bind_cols(x_tbl)

error_se_tbl <- tibble(
  error_sd = seq(from = 0.01, to = 3.01, by = .05)
)

setup_tbl <- prelim_tbl %>%
  crossing(error_se_tbl) %>%
  mutate(row_id = row_number()) %>%
  select(row_id, error_sd, everything())

setup_tbl %>%
  gt_preview()

Make Function to Generate Observations

Now we need a function that can look at each row in the setup_tbl and simulate many observations based on those setup parameters. This is basically the same thing that we did before with slice(), only now the noise variable error_sd can vary (which changes the correlation strength between x_1 and x_2.

gen_observations_fcn <- function(sd_row_number = 1, reps = 100) {
  setup_tbl %>%
    mutate(n_row = row_number()) %>%
    filter(n_row == sd_row_number) %>%
    slice(rep(1:n(), each = reps)) %>%
    mutate(n_row = row_number()) %>%
    group_by(n_row) %>%
    mutate(across(x_1:x_10, ~ rnorm(1, mean = 0, sd = 1))) %>%
    mutate(x_2 = x_1 + rnorm(n = 1, mean = 0, sd = error_sd)) %>% # control the strength of the correlation between x_1 and x_2
    mutate(
      bx_1 = b_1 * x_1,
      bx_2 = b_2 * x_2,
      bx_3 = b_3 * x_3,
      bx_4 = b_4 * x_4,
      bx_5 = b_5 * x_5,
      bx_6 = b_6 * x_6,
      bx_7 = b_7 * x_7,
      bx_8 = b_8 * x_8,
      bx_9 = b_9 * x_9,
      bx_10 = b_10 * x_10,
      resid = rnorm(1, mean = 0, sd = 1)
    ) %>%
    ungroup()
}

Use the Function in a Loop to Repeat the Simulation

Sincere apologies to the for-loop haters but I don’t currently know a more elegant way to do this part. First we establish the number of simulations we want - in this case 100 at each noise level will do fine. We then initialize a holder tibble; the data inside this one isn’t important but we want the column names to match what comes out of the simulation loop so that we can bind_rows to build out the superset of simulation data. The final step is to loop the simulation 100 times, generating 100 synthetic observations at a specified correlation each time through. Inside the loop, after the 100 observations are generated at each condition, the model is built and the coefficients, standard error, and p-values for coefficients are extracted and stored. When all is done, for every strength of correlation we should have 100 estimates based on 100 models of 100 data points generated under those specified conditions.

n_sims <- 100

# original_tbl <- tidy(summary(lm(observation ~ ., data = full_wide_tbl))) %>%
#   filter(term == "(Intercept)")

# for (i in 1:n_sims) {
#   unnested_mods_tbl <- setup_tbl %>%
#     rowwise() %>%
#     mutate(observations_tbl = map(row_id, ~ gen_observations_fcn(.x, reps = 100))) %>%
#     ungroup() %>%
#     select(observations_tbl) %>%
#     unnest() %>%
#     mutate(observation = rowSums(select(., matches("bx|resid")))) %>%
#     group_by(error_sd) %>%
#     nest() %>%
#     mutate(model = map(data, function(data) tidy(summary(lm(observation ~ x_1 + x_2 + x_3 + x_4 + x_5 + x_6 + x_7 + x_8 + x_9 + x_10, data = data))))) %>%
#     select(-data) %>%
#     unnest(model) %>%
#     filter(term == "x_1" | term == "x_2")
#
#   original_tbl <- original_tbl %>%
#     bind_rows(unnested_mods_tbl)
#
# }

Here’s what the output looks like:

cln_sim_tbl <- original_tbl %>%
  filter(term != "(Intercept)") %>%
  select(error_sd, everything())

cln_sim_tbl %>%
  gt_preview()

Now we can visualize what happens to the coefficient estimates, their standard errors, and the p-values as the correlation between x_1 and x_2 changes in strength. First, let’s look at the estimates for the coefficients.

cln_sim_tbl %>%
  ggplot(aes(x = estimate, y = error_sd)) +
  geom_jitter(aes(color = term), alpha = .4, height = .05) +
  xlim(c(-4, 4)) +
  geom_vline(xintercept = 1, linetype = 2) +
  geom_vline(xintercept = 0, linetype = 2) +
  coord_flip() +
  theme_bw() +
  scale_color_manual(values = c("purple", "limegreen")) +
  theme(legend.title = element_blank()) +
  labs(
    y = "SD of Random Noise in Relationship Between x_1 and x_2",
    x = "Estimate of Beta Coefficient",
    title = "Variable Correlation Strength vs. Coefficient Estimates",
    subtitle = "True Values are x_1 = 1, x_2 = 0",
    caption = "Lower Noise is Higher Correlation"
  )

The above plot shows that we can indeed trick the model. When the correlation is high enough, the spread of the coefficient estimates get very wide. Even though on average the estimate is still correct, the model is very uncertain about the value for any particular experiment and will occasionally think that x_2 has a large estimated coefficient. And we usually only get to run 1 experiment in real life.

Based on these data, we should expect to see the standard error of the estimate to grow large at high correlation. Let’s see if that’s true:

cln_sim_tbl %>%
  #  mutate(error_sd = as_factor(error_sd)) %>%
  ggplot(aes(x = std.error, y = error_sd)) +
  geom_jitter(aes(color = term), alpha = .4, height = .05) +
  xlim(c(0, 1.2)) +
  coord_flip() +
  theme_bw() +
  scale_color_manual(values = c("purple", "limegreen")) +
  theme(legend.title = element_blank()) +
  labs(
    y = "SD of Random Noise in Relationship Between x_1 and x_2",
    x = "Standard Error of Coefficient Estimate",
    title = "Variable Correlation Strength vs. Standard Error of Coefficient Estimates",
    caption = "Lower Noise is Higher Correlation"
  )

Just as expected - the uncertainty in the estimates spikes as the correlation between x_1 and x_2 approaches 1. But what happens to the p-values of the coefficients for each variable?

cln_sim_tbl %>%
  #  mutate(error_sd = as_factor(error_sd)) %>%
  ggplot(aes(x = p.value, y = error_sd)) +
  geom_jitter(aes(color = term), alpha = .5, height = .05) +
  geom_vline(xintercept = .05, linetype = 2, size = 1.3, alpha = .7) +
  coord_flip() +
  theme_bw() +
  scale_color_manual(values = c("purple", "limegreen")) +
  theme(legend.title = element_blank()) +
  labs(
    y = "SD of Random Noise in Relationship Between x_1 and x_2",
    x = "P-Value for Coefficient Estimate",
    title = "Variable Correlation Strength vs. P-Values for Coefficient Estimates",
    subtitle = "dotted line shows p = .05",
    caption = "Lower Noise is Higher Correlation"
  )

Yep - as noise goes to 0 and the correlation approaches 1, the p-values get very unstable and the true effect of x_1 can no longer be detected consistent relative to the null. But the good news is that as long as the correlation isn’t close to perfect, the model does its job quite well and the p-values for a true effect dive near 0.

Just to drive the point home, we could look at the actual correlation coefficients between x_1 and x_2 at the point that things start to go unstable. From the above plot, it looks like the model starts to get confused at the time the the noise is around 0.6. We can use our helper function from before to quickly make a single dataset at that noise level and then evaluate the correlation of x_1 and x_2 in that set:

unstable_tbl <- sim_fcn(noise = 0.6)

Here’s the correlation plot. We see that in this case, a correlation of around 0.8 is where things start to go unstable with the estimates.

unstable_tbl %>% 
  select(x_1:x_10) %>%
  ggcorr(
  high = "red",
  low = "blue",
  label = TRUE,
  hjust = .75,
  size = 3,
  label_size = 3,
  nbreaks = 4
) +
  labs(
    title = "Correlation Matrix",
    subtitle = "Pearson Method Using Pairwise Obervations"
  )

Here’s the variables plotted against the observation for the above scenario, with x_1 and x_2 correlated with each other at a 0.8 level.

unstable_tbl %>% plt_fcn()

In the end, we found out that we can trick the linear regression model, just not in the way we expected. We were never really able to get the model to think that x_2 was important but we were able to make it uncertain about x_1. The average coefficient estimate remains stable even in the presence of highly correlated predictors, but it’s the uncertainty about the estimates that pays the price for correlated variables. So should we just leave them all in? Like most things in statistics - it depends on the purpose of the model!

Thank you for following along with this post if you’ve made it this far. If you prefer learning concepts from mathematical principles rather than simulation, please refer to this response on Cross Validated for some explanation of the the bahavior we just saw. You could check out variance inflation factor which is the same idea presented a different way.