There are a number of assumptions that underpin linear regression models. Simulation can be a useful way of exploring these assumptions and understanding how violating these assumptions can lead to bias, large variance in the regression coefficients, and\/or poor predictions.<\/p>\n
Some typical assumptions include:<\/p>\n
- \n
- Homoskedasticity<\/li>\n
- Multicollinearity of independent variables<\/li>\n
- Measurement Error<\/li>\n
- Serial correlation<\/li>\n<\/ol>\n
Today, we will explore the assumption of homoskedasticity.<\/p>\n
As always, all code is freely available in the Github repository<\/a><\/strong><\/span>.<\/p>\n
Creating the baseline simulation<\/strong><\/span><\/p>\n
Before exploring how violations of the homoskedasticity assumption influence a regression model, we need a baseline model to compare it against. So, we will begin by simulating a simple linear regression with 1 predictor. Our model will look like this:<\/p>\n
y = 2 + 5*x + e<\/em><\/p>\n
Where e<\/strong><\/em> will be random error from a normal distribution with a mean of 0 and standard deviation of 1.<\/p>\n
The code below should look familiar as we’ve been building up simulations like this in the previous 4 tutorials. We specify the intercept to be 2 and the slope to be 5. The independent variable, x<\/strong><\/em>, <\/em>is drawn from a uniform distribution between -1 and 1. With each of the 500 iterations of the for()<\/strong> loop we store the simulated intercept, slope, and their corresponding standard errors, which we calculate using the variance-covariance matrix (which we discussed in the previous tutorial). Finally, we also store the residual standard error (RSE) of each of the simulated models.<\/p>\n
\r\nlibrary(tidymodels)\r\nlibrary(patchwork)\r\n\r\n## set seed for reproducibility\r\nset.seed(58968)\r\n\r\n## create a data frame to store intercept values, slope values, their standard errors, and the model residual standard error, for each simulation\r\nsim_params <- data.frame(intercept = NA,\r\n slope = NA,\r\n intercept_se = NA,\r\n slope_se = NA,\r\n model_rse = NA)\r\n\r\n# true intercept value\r\nintercept <- 2\r\n\r\n# true slope value\r\nslope <- 5\r\n\r\n## Number of simulations to run\r\nn <- 500\r\n\r\n# random draw from a uniform distribution to simulate the independent variable\r\nX <- runif(n = n, min = -1, max = 1)\r\n\r\n## loop for regression model\r\nfor(i in 1:n){\r\n \r\n # create dependent variable, Y\r\n Y <- intercept + slope*X + rnorm(n = n, mean = 0, sd = 1)\r\n \r\n # build model\r\n model <- lm(Y ~ X)\r\n \r\n ## store predictions\r\n fitted_vals <- model$fitted.values\r\n\r\n ## store residuals\r\n # output_df[i, 2] <- model$residuals\r\n \r\n # variance-covariance matrix for the model\r\n vcv <- vcov(model)\r\n \r\n # estimates for the intercept\r\n sim_params[i, 1] <- model$coef[1]\r\n \r\n # estimates for the slope\r\n sim_params[i, 2] <- model$coef[2]\r\n \r\n # SE for the intercept\r\n sim_params[i, 3] <- sqrt(diag(vcv)[1])\r\n \r\n # SE for the slope\r\n sim_params[i, 4] <- sqrt(diag(vcv)[2])\r\n \r\n # model RSE\r\n sim_params[i, 5] <- summary(model)$sigma\r\n \r\n}\r\n\r\nhead(sim_params)\r\n<\/pre>\n![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.30.25-PM.png\")
<\/a><\/p>\nNow we summarize the data to see if we have values close to the specified model parameters.<\/p>\n
\r\nsim_params %>%\r\n summarize(across(.cols = everything(),\r\n ~mean(.x)))\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.36.22-PM.png\")
<\/a><\/p>\nThe average intercept and slope of the 500 simulated models are pretty much identical to the specified intercept and slope of 2 and 5, respectively.<\/p>\n
The final model of the 500 iterations is also stored from our for loop and we can look directly at it and create plots of the model fit.<\/p>\n
\r\n# model summary\r\nsummary(model)\r\n\r\n# model fit plots\r\npar(mfrow = c(2,2))\r\nplot(model)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.38.36-PM.png\")
<\/a> ![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.38.48-PM-1024x862.png\")
<\/a><\/p>\nWe can also create a function that lets us evaluate how often the 95% confidence interval of our simulated beta coefficients cover the true beta coefficients that we specified for the simulation. From there, we can get a coverage probability and a 95% probability coverage interval.<\/p>\n
<\/p>\n
\r\n### Create a coverage probability function\r\ncoverage_interval95 <- function(beta_coef, se_beta_coef, true_beta_val, model_df){\r\n \r\n level95 <- 1 - (1 - 0.95) \/ 2\r\n \r\n # lower 95\r\n lower95 <- beta_coef - qt(level95, df = model_df)*se_beta_coef\r\n \r\n # upper 95\r\n upper95 <- beta_coef + qt(level95, df = model_df)*se_beta_coef\r\n \r\n # what rate did we cover the true value (hits and misses)\r\n hits <- ifelse(true_beta_val >= lower95 & true_beta_val <= upper95, 1, 0)\r\n prob_cover <- mean(hits)\r\n \r\n # create the probability coverage intervals\r\n low_coverage_interval <- prob_cover - 1.96 * sqrt((prob_cover * (1 - prob_cover)) \/ length(beta_coef))\r\n \r\n upper_coverage_interval <- prob_cover + 1.96 * sqrt((prob_cover * (1 - prob_cover)) \/ length(beta_coef))\r\n \r\n # results in a list\r\n return(list('Probability of Covering the True Value' = prob_cover,\r\n '95% Prob ability Coverage Intervals' = c(low_coverage_interval, upper_coverage_interval)))\r\n \r\n}\r\n<\/pre>\nLet’s apply the function to the intercept.<\/p>\n
\r\ncoverage_interval95(beta_coef = sim_params$intercept,\r\n se_beta_coef = sim_params$intercept_se,\r\n true_beta = intercept,\r\n model_df = model$df.residual)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.42.36-PM.png\")
<\/a><\/p>\nNow apply the function to the slope.<\/p>\n
\r\ncoverage_interval95(beta_coef = sim_params$slope,\r\n se_beta_coef = sim_params$slope_se,\r\n true_beta = slope,\r\n model_df = model$df.residual)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.43.20-PM.png\")
<\/a><\/p>\nIn both cases we are covering the true betas around 95% of the time, with relatively small intervals.<\/p>\n
Homoskedasticity<\/strong><\/span><\/p>\n
Linear models make an assumption that the variance of the residuals remain constant across the predicted values (homoskedastic). We can see what this looks like by plotting the fitted values relative to the residuals, which was the first plot in the model check plots we created for the 500th simulation above. We can see that the residuals exhibit relatively the same amount of variance across the fitted values.<\/p>\n
\r\nplot(model, which = 1)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.47.13-PM-1024x870.png\")
<\/a><\/p>\nLet’s simulate a model with heteroskedastic residuals and see what it looks like. We will keep the same intercept and slope parameters as above. The only thing will we do is add an exponential parameter to the error term\u00a0 of the model to create a heteroskedastic outcome in the residuals.<\/p>\n
\r\n## parameter for heteroskedasticity \r\nheteroskedasticity_param <- 2\r\n\r\n## set seed for reproducibility\r\nset.seed(22)\r\n\r\n## data frame for results\r\nheteroskedastic_sim_params <- data.frame(intercept = NA,\r\n slope = NA,\r\n intercept_se = NA,\r\n slope_se = NA,\r\n model_rse = NA)\r\n\r\n## for loop\r\nfor(i in 1:n ){\r\n \r\n # the error variance of Y is a function of X plus some random noise\r\n Y <- intercept + slope*X + rnorm(n = n, mean = 0, sd = exp(X*heteroskedasticity_param))\r\n \r\n # model\r\n heteroskedastic_model <- lm(Y ~ X)\r\n \r\n \r\n # variance-covariance matrix\r\n vcv <- vcov(heteroskedastic_model)\r\n \r\n # estimates for the intercept\r\n heteroskedastic_sim_params[i, 1] <- heteroskedastic_model$coef[1]\r\n \r\n # estimates for the slope\r\n heteroskedastic_sim_params[i, 2] <- heteroskedastic_model$coef[2]\r\n \r\n # SE for the intercept\r\n heteroskedastic_sim_params[i, 3] <- sqrt(diag(vcv)[1])\r\n \r\n # SE for the slope\r\n heteroskedastic_sim_params[i, 4] <- sqrt(diag(vcv)[2])\r\n \r\n # model RSE\r\n heteroskedastic_sim_params[i, 5] <- summary(heteroskedastic_model)$sigma\r\n \r\n}\r\n\r\nhead(heteroskedastic_sim_params)\r\n\r\n\r\nplot(X, Y, pch = 19)\r\n<\/pre>\n![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.49.54-PM-1024x825.png\")
<\/a><\/p>\nThe relationship between X and Y certainly looks weird given how it starts very tightly on the left side and then fans out on the right side.<\/p>\n
Let’s take the average across all 500 simulations for each coefficient and their corresponding standard errors.<\/p>\n
\r\nheteroskedastic_sim_params %>%\r\n summarize(across(.cols = everything(),\r\n ~mean(.x)))\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.50.53-PM.png\")
<\/a><\/p>\nThe coefficients of 2.0 for the intercept and 5 for the slope are exactly what we set them as for the simulation. However, notice how much larger the standard errors are for the intercept and slope compared to the original model above. Additionally, notice that the model residual standard error has increased substantially compared to the previous model.<\/p>\n
Let’s get the 500th model again and check out the fitted vs residual plot.<\/p>\n
\r\n# fitted vs residuals\r\nplot(heteroskedastic_model, which = 1)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.52.16-PM-1024x868.png\")
<\/a><\/p>\nThat looks like a large amount of heteroskedasticity as the residual variance is no longer homogenous across the range of fitted values. Notice the large fanning out towards the right side of the plot. As the predictions get larger so two does the variability in residuals, which we noticed when we plotted Y and X above.<\/p>\n
What we’ve learned is that the estimate of intercept and slope is unbiased for both the heteroskedastic and homoskedastic models, as they both are centered on the parameters that we specified for the simulation (intercept = 2, slope = 5). However, the heteroskedastic model creates greater variance in our coefficients. We can visualize how much uncertainty there is under the heteroskedastic model relative to the homoskedastic model by visualizing the density of the coefficient estimates from our two model simulations.<\/p>\n
\r\nplt_intercept <- sim_params %>%\r\n mutate(model = 'homoskedastic model') %>%\r\n bind_rows(\r\n heteroskedastic_sim_params %>%\r\n mutate(model = 'heteroskedastic model')\r\n ) %>%\r\n ggplot(aes(x = intercept, fill = model)) +\r\n geom_density(alpha = 0.6) +\r\n theme_classic() +\r\n theme(legend.position = "top")\r\n\r\nplt_slope <- sim_params %>%\r\n mutate(model = 'homoskedastic model') %>%\r\n bind_rows(\r\n heteroskedastic_sim_params %>%\r\n mutate(model = 'heteroskedastic model')\r\n ) %>%\r\n ggplot(aes(x = slope, fill = model)) +\r\n geom_density(alpha = 0.6) +\r\n theme_classic() +\r\n theme(legend.position = "none")\r\n\r\nplt_intercept | plt_slope\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.56.26-PM-1024x742.png\")
<\/a><\/p>\nFinally, let’s see how often the 95% coverage interval is covering the true intercept and slope in the heteroskedastic model.<\/p>\n
\r\ncoverage_interval95(beta_coef = heteroskedastic_sim_params$intercept,\r\n se_beta_coef = heteroskedastic_sim_params$intercept_se,\r\n true_beta = intercept,\r\n model_df = model$df.residual)\r\n\r\n\r\ncoverage_interval95(beta_coef = heteroskedastic_sim_params$slope,\r\n se_beta_coef = heteroskedastic_sim_params$slope_se,\r\n true_beta = slope,\r\n model_df = model$df.residual)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/08\/Screenshot-2023-08-11-at-7.57.17-PM-1024x656.png\")
<\/a><\/p>\nNotice that we are no longer covering the true model values at the 95% level.<\/p>\n
Wrapping Up<\/strong><\/span><\/p>\n
Simulations can be useful for evaluating model assumptions and seeing how a model may be have under different circumstances. In this tutorial we saw that a model suffering from severe heteroskedasticity is still able to return the true values for the intercept and slope. However, the variance around those values is large, meaning that we have a lot of uncertainty about what those true parameters actually are. In upcoming tutorials we will explore other linear regression assumptions. If this is a topic you enjoy, a book that I recommend (which heavily influenced today’s tutorial) is Monte Carlo Simulation and Resampling Methods for Social Sciences by Carsey & Harden<\/a><\/span><\/strong>.<\/p>\n