 B-Coefficients Simulation

# How does Collinearity Influence Linear Regressions?

This is a short simulation study trying to figure out the impact of collinearity on linear regressions.

Overview:

## Packages

# install pacman once if not avaible on your machine
# install.packages("pacman")

pacman::p_load(arm, purrr, MASS, broom, ggthemes, tidyverse, ecodist, viridis, gridExtra, grid, lm.beta, tidyr, ggrepel)


## Simulation Function

First, I write a little function to simulate collinearity.

generate_multi <- function(n, cor_seq){
set.seed(2017)
x <- runif(n, 1, 10)
models <- list()
std.models <- list()
for (jj in seq_along(cor_seq)) {
## generate correlated variable x2
dat <- data.frame(corgen(x = x, r = cor_seq[jj],  epsilon = 0))
colnames(dat) <- c("x1", "x2")
## generate y variable
dat$y <- 0.5 * dat$x1 + 0.5 * dat$x2 + rnorm(n, sd = 10) ## modelling and tidy dataframe models[[jj]] <- tidy(lm(y ~ x1 + x2, data = dat)) ## get standardized betas std.models[[jj]] <- data.frame(lm.beta::coef.lm.beta((lm.beta::lm.beta(lm(y ~ x1 + x2, data = dat))))) colnames(std.models[[jj]]) <- c("std.estimate") ## bind it together models[[jj]] <- std.models[[jj]] %>% bind_cols(models[[jj]]) models[[jj]]$cors <- cor_seq[jj]
}
sim_dat <- bind_rows(models)
sim_dat$col <- n return(sim_dat) } draw.data <- function(cor_seq = NULL, step_seq = NULL){ sim.list <- list() for(jj in seq_along(step_seq)) { sim.list[[jj]] <- generate_multi(n = step_seq[jj], cor_seq) sim.list[[jj]]$n <- step_seq[jj]

cat(paste0("Batch: ", jj, "\t"))

}
sim_data <- bind_rows(sim.list)
return(sim_data)
}


## Simulate and Save Data

Draw data from function and save it.

sim_data <- draw.data(cor_seq = seq(0,.99,0.01), step_seq = seq(50, 10000, by = 50))

if(!dir.exists("data")) dir.create("data")

save(sim_data, file = "data/sim_data.Rdata")


## Visualizing the Influence of Collinearity

load("data/sim_data.Rdata")


Now, consider the following linear regression: y ~ x1 + x2.

\begin{equation} Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \epsilon \end{equation}

For this simulation, I consistently increase the correlation between x1 and x2 (from 0 to .99) and the sample size (from 50 to 10.000) and estimate separate models for each of these combinations,

### Standard Errors

First, I take a look at the impact that sample size and collinearity have on the standard error of x1. On the x-axis, you can see the increasing collinearity between x1 and x2 and the individual lines are colored in by sample size. There are a few things we can observe here: the standard error diminishes with greater sample size, correlations around .75 increase the standard error quite drastically and this effect is weaker the higher your sample size is.

get_smooths <- function(smooth_dat, n_val, y) {
smooth_dat <- filter(smooth_dat, n == n_val & term == "x1")
fm <- paste0(y," ~ cors")
smooth_vals <- predict(loess(fm, smooth_dat), smooth_dat\$cors)

smooth_dat %>%
mutate(smooth = smooth_vals) %>%
group_by(n) %>%
summarise(max_smooth = max(smooth),
min_smooth = min(smooth)
)
}

smooth_dat <- c(50, 100, 150, 200, 10000) %>%
map_df(~get_smooths(sim_data, n_val = .x, y = "std.error")) %>%
mutate(n_lab = ifelse(n == 50, "Sample Size: 50", n))

sim_data %>%
filter(term == "x1") %>%
ggplot(aes(cors, std.error, colour = n, group = n)) +
geom_smooth(method = "loess", se = F, size = 1, alpha = 0.5) +
xlab(expression("Pearson's"~r~correlation~between~x~and~x)) +
ylab(expression(x~Standard~Error)) +
theme_hc() +
scale_color_viridis("Sample Size", direction = -1,
breaks = seq(1000, 10000, 3000),
labels = seq(1000, 10000, 3000)) +
ggtitle("Sample Size and Collinearity Influence on Standard Error") +
geom_point(data = smooth_dat, aes(x = .99, y = max_smooth)) +
geom_text_repel(data = smooth_dat, aes(x = .99, y = max_smooth, label = n_lab),
nudge_y = 0.07, nudge_x = 0.03) +
guides(colour = guide_colourbar(barwidth = 20, label.position = "bottom")) ggsave(filename = "images/std_static.png", width = 10, height = 7)


### T-Statistic and P-Values

Next, let’s take a look at the t-statistic and p-values. In line with what you would expect with greater standard errors, statistical significance also suffers from collinearity. Again, greater sample sizes seem to partly remedy this but the shift is clearly observable.

sim_data  %>%
filter(term == "x1") %>%
ggplot(aes(cors, statistic, colour = n, group = n)) +
geom_smooth(method = "loess", se = F, size = 1, alpha = 0.5) +
xlab(expression("Pearson's"~r~correlation~between~x~and~x)) +
ylab(expression(x~t-Statistic)) +
theme_hc() +
scale_color_viridis("Sample Size", direction = -1,
breaks = seq(1000, 10000, 3000),
labels = seq(1000, 10000, 3000)) +
geom_hline(yintercept = 1.96, linetype = "dashed", alpha = 0.9) +
annotate(geom = "text", x = 0, y = 2.3, label = "t = 1.96") +
ggtitle("Sample Size and Collinearity Influence on t-Statistic") +
guides(colour = guide_colourbar(barwidth = 20, label.position = "bottom")) ggsave(filename = "images/t_static.png", width = 10, height = 7)

sim_data  %>%
filter(term=="x1") %>%
ggplot(aes(cors, p.value, colour = n, group = n)) +
geom_smooth(method = "loess", se = F, size = 1, alpha = 0.5) +
xlab(expression("Pearson's"~r~correlation~between~x~and~x)) +
ylab(expression(x~p-value)) +
theme_hc() +
scale_color_viridis("Sample Size", direction = -1,
breaks = seq(1000, 10000, 3000),
labels = seq(1000, 10000, 3000)) +
geom_hline(yintercept = 0.05, linetype = "dashed", alpha = 0.9) +
annotate(geom = "text", x = 0, y = 0.08, label = "p = 0.05") +
ggtitle("Sample Size and Collinearity Influence on p-values") +
guides(colour = guide_colourbar(barwidth = 20, label.position = "bottom")) ggsave(filename = "images/p_static.png", width = 10, height = 7)


### B-Coefficients

This one I find most interesting. In addition to statistical signif. being affected, you can also see how the coefficients change with increasing collinearity. The most drastic impact is for small sample sizes (~1000) where a coefficient of .5 can even become negative.

sim_data  %>%
filter(term=="x1") %>%
filter(n>200) %>%
ggplot(aes(cors, estimate, colour = n, group = n)) +
geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.9) +
geom_line(alpha = 0.5) +
xlab(expression("Pearson's"~r~correlation~between~x~and~x)) +
ylab(expression(x~b-coefficient)) +
theme_hc() +
scale_color_viridis("Sample Size", direction = -1,
breaks = seq(1000, 10000, 3000),
labels = seq(1000, 10000, 3000)) +
annotate(geom = "text", x = 0, y = 0.08, label = "b = 0.5") +
ggtitle("Sample Size and Collinearity Influence on b-coefficients") +
guides(colour = guide_colourbar(barwidth = 20, label.position = "bottom")) ##### Animation
library(gganimate)

sim_data_sub <- sim_data  %>%
filter(term == "x1") %>%
filter(n > 200) %>%
filter(n %in% c(300, 1000, 10000)) %>%
mutate(estimate_lab = round(estimate, 2) %>% as.character) %>%
mutate(n = as.factor(n))

anim1 <-  sim_data_sub  %>%
# filter(term == "x1") %>%
# filter(n > 200) %>%
ggplot(aes(cors, estimate, colour = n, group = n)) +
geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.9) +
geom_line() +
geom_segment(aes(xend = 1, yend = estimate), linetype = 2, colour = 'grey') +
geom_point(size = 2) +
geom_text(aes(x = 1, label = n), hjust = 0, size = 4, fontface = "bold") +
geom_text(aes(x = 0.15, y = 1.8, label = paste0("Correlation: ", cors)),
hjust = 1, size = 5, color = "black") +
geom_text(aes(label = estimate_lab), hjust = 0, size = 3, fontface = "bold", nudge_y = 0.1) +
xlab("Pearson's r correlation between x1 and x2") +
ylab("x1 b-coefficient") +
coord_cartesian(clip = 'off') +
theme_hc() +
scale_color_viridis("Sample Size",
direction = -1,
discrete = T,
begin = 0.3,
# limits = seq(1000, 10000, 3000),
breaks = seq(1000, 10000, 3000),
labels = seq(1000, 10000, 3000)) +
guides(colour = F) +
theme(title = element_text(size = 15, face = "bold"),
axis.text.x = element_text(size = 14, face = "bold"),
axis.text.y = element_text(size = 10, face = "italic")) +
ggtitle("Sample Size and Collinearity Influence on b-coefficients (Sample Sizes: 300, 1000 and 10.000)") +
# guides(colour = guide_colourbar(barwidth = 20, label.position = "bottom")) +
# Here comes the gganimate code
transition_reveal(n, cors)

anim1 %>% animate(
nframes = 500, fps = 15, width = 1000, height = 600, detail = 3
)

anim_save("images/b_anim.gif") ggsave(filename = "images/b_static.png", width = 10, height = 7)

##### Standardized

Same spiel, just this time with standardized coefficients.

sim_data  %>%
filter(term=="x1") %>%
filter(n>200) %>%
ggplot(aes(cors, std.estimate, colour = n, group = n)) +
geom_line(alpha = 0.5) +
xlab(expression("Pearson's"~r~correlation~between~x~and~x)) +
ylab(expression(x~b-coefficient)) +
theme_hc() +
scale_color_viridis("Sample Size", direction = -1,
breaks = seq(1000, 10000, 3000),
labels = seq(1000, 10000, 3000)) +
ggtitle("Sample Size and Collinearity Influence on standardized b-coefficients") +
guides(colour = guide_colourbar(barwidth = 20, label.position = "bottom")) ggsave(filename = "images/b_standardized_static.png", width = 10, height = 7)

##### Session Info
sessionInfo()

## R version 3.5.0 (2018-04-23)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17134)
##
## Matrix products: default
##
## locale:
##  LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252
##  LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
##  LC_TIME=German_Germany.1252
##
## attached base packages:
##  grid      stats     graphics  grDevices utils     datasets  methods
##  base
##
## other attached packages:
##   bindrcpp_0.2.2     ggrepel_0.8.0      lm.beta_1.5-1
##   gridExtra_2.3      viridis_0.5.1      viridisLite_0.3.0
##   ecodist_2.0.1      forcats_0.3.0      stringr_1.3.0
##  tibble_1.4.2       ggplot2_3.0.0.9000 tidyverse_1.2.1
##  ggthemes_4.0.0     broom_0.4.4        purrr_0.2.4
##  arm_1.10-1         lme4_1.1-17        Matrix_1.2-14
##  MASS_7.3-49
##
## loaded via a namespace (and not attached):
##   Rcpp_0.12.18     lubridate_1.7.4  lattice_0.20-35  assertthat_0.2.0
##   rprojroot_1.3-2  digest_0.6.15    psych_1.8.3.3    R6_2.2.2
##   cellranger_1.1.0 plyr_1.8.4       backports_1.1.2  evaluate_0.10.1
##  coda_0.19-1      httr_1.3.1       pillar_1.2.1     rlang_0.2.1
##  lazyeval_0.2.1   readxl_1.1.0     minqa_1.2.4      rstudioapi_0.7
##  nloptr_1.0.4     rmarkdown_1.9    labeling_0.3     splines_3.5.0
##  foreign_0.8-70   munsell_0.4.3    compiler_3.5.0   modelr_0.1.1
##  pkgconfig_2.0.1  mnormt_1.5-5     htmltools_0.3.6  tidyselect_0.2.4
##  crayon_1.3.4     withr_2.1.2      nlme_3.1-137     jsonlite_1.5
##  gtable_0.2.0     pacman_0.4.6     magrittr_1.5     scales_0.5.0
##  cli_1.0.0        stringi_1.1.7    reshape2_1.4.3   xml2_1.2.0
##  tools_3.5.0      glue_1.3.0       hms_0.4.2        abind_1.4-5
##  parallel_3.5.0   yaml_2.1.19      colorspace_1.4-0 rvest_0.3.2
##  knitr_1.20       bindr_0.1.1      haven_1.1.2