{"id":3382,"date":"2024-03-24T23:03:30","date_gmt":"2024-03-24T23:03:30","guid":{"rendered":"http:\/\/optimumsportsperformance.com\/blog\/?p=3382"},"modified":"2024-03-25T12:25:15","modified_gmt":"2024-03-25T12:25:15","slug":"frequentist-bayesian-approaches-to-regression-models-by-hand-compare-and-contrast","status":"publish","type":"post","link":"https:\/\/optimumsportsperformance.com\/blog\/frequentist-bayesian-approaches-to-regression-models-by-hand-compare-and-contrast\/","title":{"rendered":"Frequentist & Bayesian Approaches to Regression Models by Hand – Compare and Contrast"},"content":{"rendered":"

One of the ways I try and learn things is to code them from first principles. It helps me see what is going on under the hood and also allows me wrap my head around how things work. Building regression models in R is incredibly easy using the lm()<\/strong> function and Bayesian regression models can be conveniently built with the same syntax in packages like {rstanarm<\/strong>} and {brms<\/strong>}. However, today I’m going to do everything by hand to get a grasp of how the Bayesian regression model works, at a very basic level. I put together the code using the mathematical presentation of these concepts from William Bolstad’s book, Introduction to Bayesian Statistics<\/a><\/span><\/strong>.<\/p>\n

The entire script is available on my GITHUB page<\/a><\/span><\/strong>.<\/p>\n

Loading Packages & Getting Data<\/strong><\/span><\/p>\n

I’m going to use the data from the {palmerpenguins<\/strong>} package and concentrate on a simple linear regression which will estimate flipper length from bill length.<\/p>\n

\r\n### packages ------------------------------------------------------\r\nlibrary(tidyverse)\r\nlibrary(palmerpenguins)\r\nlibrary(patchwork)\r\n\r\ntheme_set(theme_classic())\r\n\r\n### data ----------------------------------------------------------\r\ndata("penguins")\r\ndat <- penguins %>%\r\n  select(bill_length = bill_length_mm, flipper_length = flipper_length_mm) %>%\r\n  na.omit()\r\n\r\nhead(dat)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
 <\/p>\n
EDA<\/strong><\/span><\/p>\n
These two variables share a relatively large correlation with each other.<\/p>\n
\"\"<\/a> \"\"<\/a><\/p>\n
Ordinary Least Squares Regression<\/strong><\/span><\/p>\n
First, I’ll start by building a regression model under the Frequentist paradigm, specifying no prior values on the parameters (even though in reality the OLS model is using a flat prior, where all values are equally plausible — a discussion for a different time), using the lm()<\/strong> function.<\/p>\n
\r\n### Ordinary Least Squares Regression ------------------------------\r\nfit_ols <- lm(flipper_length ~ I(bill_length - mean(dat$bill_length)), data = dat)\r\nsummary(fit_ols)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Technical Note: <\/strong>Since we don’t observe a bill length of 0 mm in penguins, I chose to transform the bill length data by grand mean centering it — subtracting each bill length from the population mean. This wont change the predictions from the model but does change our interpretation of the coefficients themselves (and the intercept is now directly interpretable, which it wouldn’t be if we left the bill length data in its untrasnformed state).<\/em><\/p>\n
OLS Regression by Hand<\/strong><\/span><\/p>\n
Okay, that was easy. In a single line we constructed a model and we are up and running. But, let’s calculate this by hand so we can see where the coefficients and their corresponding standard errors came from.<\/p>\n
\r\n# Calculate the least squares regression line by hand now\r\ndat_stats <- dat %>%\r\n  summarize(x_bar = mean(bill_length),\r\n            x_sd = sd(bill_length),\r\n            y_bar = mean(flipper_length),\r\n            y_sd = sd(flipper_length),\r\n            r = cor(bill_length, flipper_length),\r\n            x2 = x_bar^2,\r\n            y2 = y_bar^2,\r\n            xy_bar = x_bar * y_bar,\r\n            .groups = "drop")\r\n\r\ndat_stats\r\n<\/pre>\n
In the above code, I’m calculating summary statistics (mean and SD) for both the independent and dependent variables, extracting their correlation coefficient, and getting their squared means and the product of their squared means to be used in downstream calculations of the regression coefficients. Here is what the output looks like:<\/p>\n
\"\"<\/a><\/p>\n
The model intercept is simply the mean of our outcome variable (flipper length).<\/p>\n
\r\nintercept <- dat_stats$y_bar\r\n<\/pre>\n
The coefficient for our independent variable is calculated as the correlation coefficient multiplied by the the SD of the y variable divided by the SD of the x variable.<\/p>\n
\r\nbeta <- with(dat_stats,\r\n             r * (y_sd \/ x_sd))\r\n\r\nbeta\r\n<\/pre>\n
Finally, I’ll store the grand mean of bill length so that we can use it when we need to center bill length and make predictions.<\/p>\n
\r\nx_bar <- dat_stats$x_bar\r\nx_bar\r\n<\/pre>\n
Now that we have the intercept and beta coefficient for bill length calculated by hand we can construct the regression equation:<\/p>\n
\"\"<\/a><\/p>\n
Notice that these are the exact values we obtained from our model using the lm()<\/strong> function.<\/p>\n
We can use the two models to make predictions and show that they are the exact same.<\/p>\n
\r\n## Make predictions with the two models\r\ndat %>%\r\n  mutate(pred_model = predict(fit_ols),\r\n         pred_hand = intercept + beta * (bill_length - x_bar))\r\n<\/pre>\n
\"\"<\/a><\/p>\n
In the model estimated from the lm()<\/strong> function we also get a sigma parameter (residual squared error, RSE), which tells us, on average, how off the model predictions are. We can build this by hand by first calculating the squared error of the observed values and our predictions and then calculating the RSE as the square root of the sum of squared residuals divided by the model degrees of freedom.<\/p>\n
\r\n# Calculate the estimated variance around the line\r\nN_obs <- nrow(dat) dat %>%\r\n  mutate(pred = intercept + beta * (bill_length - x_bar),\r\n         resid = (flipper_length - pred),\r\n         resid2 = resid^2) %>%\r\n  summarize(n_model_params = 2,\r\n            deg_freedom = N_obs - n_model_params,\r\n            model_var = sum(resid2) \/ deg_freedom,\r\n            model_sd = sqrt(model_var))\r\n<\/pre>\n
Again, our by hand calculations equal exactly as what we obtained from the lm() <\/strong>function, 10.6.<\/p>\n
\"\"<\/a> \"\"<\/a><\/p>\n
Bayesian Linear Regression by Hand<\/strong><\/span><\/p>\n
Now let’s turn our attention to building a Bayesian regression model.<\/p>\n
First we need to calculate the sum of squared error for X, which we do with this equation:<\/p>\n
ss_x = N * (mean(mu_x^2) – mu_x^2)<\/strong><\/p><\/blockquote>\n
\r\nN <- nrow(dat)\r\nmean_x2 <- mean(dat$bill_length^2)\r\nmu_x2 <- mean(dat$bill_length)^2\r\n\r\nN\r\nmean_x2\r\nmu_x2\r\n\r\nss_x <- N * (mean_x2 - mu_x2)\r\nss_x\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Next, we need to specify some priors on the model parameters. The priors can come from a number of courses. We could set them subjectively (usually people hate this idea). Or, we could look in the scientific literature and use prior research as a guide for plausible values of the relationship between flipper length and bill length. In this case, I’m just going to specify some priors that are within the range of the data but have enough variance around them to “let the data speak”.<\/em> (NOTE:<\/strong> you could rerun the entire analysis with different priors and smaller or larger variances to see how the model changes. I hope to do a longer blog post about priors in the future).<\/p>\n