library(tidyverse)        sim <- function(rho){      #Number of samples to draw      N = 50      #Make a covariance matrix      covar = matrix(c(1,rho, rho, 1), byrow = T, nrow = 2)      # Append a column of 1s to N draws from a 2-dimensional        # Gaussian       # With covariance matrix covar      X = cbind(rep(1,N),MASS::mvrnorm(N, mu = c(0,0),                   Sigma = covar))      # True betas for our regression      betas = c(1,2,4)      # Make the outcome      y = X%*%betas + rnorm(N,0,1)      # Fit a linear model      model = lm(y ~ X[,2] + X[,3])      # Return a dataframe of the coefficients      return(tibble(a1 = coef(model)[2], a2 = coef(model)[3]))         }    #Run the function 1000 times and stack the results    zero_covar = rerun(1000, sim(0)) %>%                  bind_rows    #Same as above, but the covariance in covar matrix     #is now non-zero    high_covar = rerun(1000, sim(0.95)) %>% bind_rows    #plot    zero_covar %>%       ggplot(aes(a1,a2)) +      geom_point(data = high_covar, color = 'red') +      geom_point()