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()
## # A tibble: 11 x 6
##    term        estimate std.error statistic p.value `significant?`
##    <chr>          <dbl>     <dbl>     <dbl>   <dbl> <chr>         
##  1 (Intercept)    -0.27      0.09     -2.88    0    heck_yes      
##  2 x_1             1.14      0.79      1.43    0.16 nope          
##  3 x_2            -0.02      0.77     -0.02    0.98 nope          
##  4 x_3             0.08      0.09      0.8     0.42 nope          
##  5 x_4            -0.11      0.09     -1.28    0.2  nope          
##  6 x_5             0.06      0.08      0.75    0.46 nope          
##  7 x_6             0.02      0.11      0.21    0.83 nope          
##  8 x_7            -0.12      0.1      -1.16    0.25 nope          
##  9 x_8            -0.13      0.09     -1.41    0.16 nope          
## 10 x_9             0.1       0.1       1.04    0.3  nope          
## 11 x_10            0.05      0.09      0.62    0.54 nope

Even in the high noise 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 change 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()) %>%
  dplyr::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 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() %>%
#     dplyr::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 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 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 = "Standard Error of 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.

In the end, we found out that we can trick the linear regression model, just not in the way we expected. The average coefficient estimate remains stable even in the presence of highly correlated predictors, but it’s the uncertainty about this estimate that pays the price.

Thank you for following along with this post if you’ve made it this far. If you are interested in the math behind what we showed in the above simulations, please refer to this response on Cross Validated: https://stats.stackexchange.com/a/366579/31189