Simple Linear Regression (featuring Bayes!) in R, Python, and Julia

$$ \mathbf{y} = \alpha + \mathbf{X \beta} + \varepsilon $$

where:

• \(\mathbf{y}\) is the outcome/dependent variable
• \(\alpha\) is the intercept
• \(\mathbf{X}\) is the feature/covariate matrix
• \(\mathbf{\beta}\) is the coefficient vector
• \(\varepsilon\) is the error term

Linear regression models the mean and variance of a variable. In the case of simple linear regression, the main purpose of this model is to determine how the average value of the continuous outcome \(y\) varies with the value of a single predictor \(x\). The average values of the outcome are assumed to lie on the regression line.

If model assumptions hold, \(\hat{\beta}\) coefficients are unbiased estimates. More specifically, an estimate is considered unbiased if, over many repreated samples drawn from the population, the average value of the estimates based on the different samples would equal the population value of the parameter being estimated.

A major drawback of this approach is that it can be sensitive to outliers.

HERS Dataset

The example used in this post is from the Heart and Estrogen/Progestin Replacement Study (HERS) trial ¹ that I encountered in Regression Methods in Biostatistics by Vittinghoff, Glidden, Shiboski, and McCulloch ². In the figure below, we see the relationship between baseline systolic blood pressure (SBP) by age among HERS trial participants.

hers %>% 
  ggplot(aes(x = age, y = SBP)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "Age (years)", y = "Systolic Blood Pressure") +
  theme_minimal()

The formula for the regression line has interpretable parameters:

$$ \mathbb{E}[y|x] = \alpha + \beta_1 \cdot \text{age} $$

$$ \mathbb{E}[y|x] = \text{avg value of SBP for a given age} $$

\(\mathbb{E}[y|x]\) is the “expectation of \(y\) given \(x\).” In other words, it is the expected or average value of the outcome \(y\) at a given value of the predictor \(x\). In the context of this simple linear regression model, \(\beta_1\) is interpreted as the change in average SBP for a unit change (1-year) increase in age.

Being Bayesian About It

For the Bayesian models we will be using below, we will be using a Gaussian/normal likelihood function. This is represented using the following Bayesian regression model:

In other words, our likelihood function \(P(\mathbf{y} \mid \boldsymbol{\theta})\) is a normal distribution in which \(\mathbf{y}\) depends on the parameters of the model \(\alpha\) and \(\boldsymbol{\beta}\), in addition to having an error \(\sigma\). The remaining 3 outline our “priors” for \(\alpha, \beta, \sigma\).

$$ \text{Prior} \times \text{Likelihood} \propto \text{Posterior} $$

• \(\text{Prior}\): what we already knew
• \(\text{Likelihood}\): what the data tells us
• \(\text{Posterior}\): what we know now (given data)

In this is example, the posterior is \(\Pr(\alpha, \beta, \sigma \ | \ \text{SBP, age})\). We have 3 priors: \(\alpha, \beta, \sigma\) that we define.

R

Linear Regression in R

We can just use the lm model in base-R.

model_lm <- lm(SBP ~ age, data=hers)

summary(model_lm) %>% 
  tidy()

## # A tibble: 2 × 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)  104.       3.60       28.8  6.54e-160
## 2 age            0.472    0.0537      8.79 2.64e- 18

We can interpret the output of this model. The estimate of \(\beta_{\text{age}}\) from the study suggests that among women with heart disease, average SBP increases 0.47 mmHg for each one-year increase in age.

When the two-sided test is statistically significant at \(p < 0.05\), then the corresponding 95% CI will exclude the null parameter value. For example, our \(\beta_1\) for age is statistically significant and (as we would expect) the CI excludes the null value of 0.

confint(model_lm)

##                  2.5 %      97.5 %
## (Intercept) 96.5789056 110.6800709
## age          0.3664636   0.5769925

Centering/Re-scaling for an in interpretable Intercept

An arguably “better” approach to help make the intercept term more interpretable would be to center or re-scale the predictor. We can do so here by subtracting by the sample average age of 67. The centered age variable would have value of 0 for women at age 67, and the new intercept, 135 mmHg, estimates average SBP among women this age. The \(\beta_{\text{age}}\) estimate is unaffected by “centering” the age variable.

$$ \mu_i = \alpha + \beta (\text{age}_i - \overline{\text{age}}_i) $$

where \(\overline{\text{age}}_i\) is the mean value of \(\text{age}_i\). Now \(\alpha\) can be interpreted as the value of \(\mu\) when \(\text{age}_i - \overline{\text{age}}_i = 0\), i.e. \(\alpha\) estimates the average SBP when age is \(\overline{\text{age}}_i\), 67 in this case. In other words, \(\alpha\) is the expected value of SBP when an individual is at the average age.

hers <- hers %>% 
  mutate(centered_age = age - mean(age))

Now let’s run our simple regression model again, but this time where \(x_i = \text{age}_i - \overline{\text{age}}\), which is captured by the centered_age variable.

model_lm_centered <- lm(SBP ~ centered_age, data = hers) 

model_lm_centered %>% 
  tidy()

## # A tibble: 2 × 5
##   term         estimate std.error statistic  p.value
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   135.       0.357     378.   0       
## 2 centered_age    0.472    0.0537      8.79 2.64e-18

Comparing the original model against the model using the “centered” values, we can confirm that \(\beta_{\text{age}}\) remains the same, but only the Intercept value changes.

huxtable::huxreg(
  "Original model" = model_lm, 
  "Centered model" = model_lm_centered
  )

Simple regression slope = Pearson correlation coefficient

For simple linear regression models, the slope coefficient ($\beta_1$) is systematically related to the Pearson correlation coefficient ($r$ or \(\rho\)).

Also, \(R^2 = \rho^2 = \frac{MSS}{TSS}\)

$$ \rho = \beta_1 \sigma_x / \sigma_y $$

where \(\sigma_x\) and \(\sigma_y\) are the standard deviations of the predictor and outcome variables, respectively.

We can verify this here:

broom::glance(model_lm)

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic  p.value    df  logLik    AIC    BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>   <dbl>  <dbl>  <dbl>
## 1    0.0272        0.0269  18.8      77.2 2.64e-18     1 -12021. 24049. 24067.
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

cor(hers$SBP, hers$age)^2

## [1] 0.02720518

Both \(R^2\) and \(\rho^2\) are equal to \(0.027\).

Bayesian Regression in R

brms

library(brms)

get_prior(SBP ~ 1 + age,
          data = hers,
          family = gaussian()
         )

##                    prior     class coef group resp dpar nlpar bound
##                   (flat)         b                                 
##                   (flat)         b  age                            
##  student_t(3, 134, 19.3) Intercept                                 
##    student_t(3, 0, 19.3)     sigma                                 
##        source
##       default
##  (vectorized)
##       default
##       default

We can use the default priors listed above from get_prior like so:

model_brms_defaultpriors <- brm(
  SBP ~ 1 + age, 
  data = hers, 
  family = gaussian()
  )

summary(model_brms_defaultpriors)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: SBP ~ 1 + age 
##    Data: hers (Number of observations: 2763) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   103.69      3.60    96.72   110.87 1.00     3701     2630
## age           0.47      0.05     0.36     0.57 1.00     3622     2590
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    18.78      0.24    18.31    19.26 1.00     4208     2661
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

We can specify our desired priors laid out below in a list to the prior argument within brm():

model_brms <- 
  brm(data = hers, family = gaussian,
      SBP ~ 1 + age,
      prior = c(prior(normal(120, 30), class = Intercept),
                prior(normal(0, 1), class = b),
                prior(cauchy(0, 1), class = sigma)),
      sample_prior = TRUE,
      iter = 2000,
      warmup = 1000,
      chains = 4,
      cores = 2,
      seed = 4
      )

summary(model_brms)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: SBP ~ 1 + age 
##    Data: hers (Number of observations: 2763) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   103.81      3.61    96.75   111.07 1.00     4578     2991
## age           0.47      0.05     0.36     0.57 1.00     4568     3086
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    18.77      0.26    18.28    19.28 1.00     3998     2648
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Rhat - \(\hat{R}\) - measures the ratio of the average variance of samples within each chain to the variance of the pooled samples across chains. \(\hat{R}\) values provide information about the convergence of the algorithm. \(\hat{R}\) values close to 1 suggest that the model has converged. This makes sense as \(\hat{R} \approx 1\) implies that all chains are at equilibrium.

We can also visually inspect the convergence using a plot:

plot(model_brms)

Posterior Predictive Check

Run a “posterior predictive check” using pp_check (docs):

# Overall
pp_check(model_brms, ndraws=100)

# How well did we capture the mean? (can check out other stats similarly)
pp_check(model_brms, type = "stat", stat = "mean")

# Scatter plot
pp_check(model_brms, type = "scatter_avg")

This plot shows the predictions from our prior distributions in blue and the actual data in black. This serves as a nice “check” to see if our specified prior distributions are appropriate.

Prior-Posterior Update Plots

Another way to ensure that our model represents a good fit to the data is to plot prior-posterior update plots. These plots show how our model updates from our priors after seeing the data.

To plot these data, we’ll use as_draws_df() to sample from the prior and posterior distributions for the relevant parameters:

#Overview of model parameters:
variables(model_brms)

## [1] "b_Intercept"     "b_age"           "sigma"           "prior_Intercept"
## [5] "prior_b"         "prior_sigma"     "lp__"

#Sample the parameters of interest:
posterior_m <- as_draws_df(model_brms)

#Plot the prior-posterior update plot for the *Intercept*:
ggplot(posterior_m) +
  geom_density(aes(prior_Intercept), fill="steelblue", color="black",alpha=0.6) +
  geom_density(aes(b_Intercept), fill="#FC4E07", color="black",alpha=0.6) + 
  theme_classic()

#Plot the prior-posterior update plot for *age*:
ggplot(posterior_m) +
  geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
  geom_density(aes(b_age), fill="#FC4E07", color="black",alpha=0.6) + 
  theme_classic()

The posterior is shown in blue and the prior shown in orange.

rstanarm

Like brms, we can use the default priors in rstanarm:

library(rstanarm)

# `rstanarm` can be expressed like our `lm` model
model_rstanarm_defaultpriors <- stan_glm(SBP ~ age, data = hers)

summary(model_rstanarm_defaultpriors)

## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      SBP ~ age
##  algorithm:    sampling
##  sample:       4000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 2763
##  predictors:   2
## 
## Estimates:
##               mean   sd    10%   50%   90%
## (Intercept) 103.7    3.7  98.8 103.8 108.4
## age           0.5    0.1   0.4   0.5   0.5
## sigma        18.8    0.3  18.5  18.8  19.1
## 
## Fit Diagnostics:
##            mean   sd    10%   50%   90%
## mean_PPD 135.1    0.5 134.4 135.1 135.7
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.1  1.0  3852 
## age           0.0  1.0  3885 
## sigma         0.0  1.0  3986 
## mean_PPD      0.0  1.0  3759 
## log-posterior 0.0  1.0  1883 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Use the prior_summary function to see the default priors

# You can see what priors are used using `prior_summary`
prior_summary(model_rstanarm_defaultpriors)

## Priors for model 'model_rstanarm_defaultpriors' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 135, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 135, scale = 48)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = 0, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0, scale = 7.1)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.053)
## ------
## See help('prior_summary.stanreg') for more details

We can also explicitly define our priors using rstanarm as follows:

model_rstanarm <- stan_glm(SBP ~ age,
                            family = gaussian,
                            prior_intercept = normal(120, 30),
                            prior = normal(0, 1),
                            prior_aux = cauchy(0, 1),
                            data = hers
                           )

summary(model_rstanarm)

## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      SBP ~ age
##  algorithm:    sampling
##  sample:       4000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 2763
##  predictors:   2
## 
## Estimates:
##               mean   sd    10%   50%   90%
## (Intercept) 103.8    3.5  99.3 103.7 108.2
## age           0.5    0.1   0.4   0.5   0.5
## sigma        18.8    0.3  18.5  18.8  19.1
## 
## Fit Diagnostics:
##            mean   sd    10%   50%   90%
## mean_PPD 135.1    0.5 134.4 135.1 135.7
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.1  1.0  4264 
## age           0.0  1.0  4237 
## sigma         0.0  1.0  3839 
## mean_PPD      0.0  1.0  3942 
## log-posterior 0.0  1.0  1875 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

We can use the bayesplot package to visually inspect the model.

See tables here for more details on explicitly defining priors using rstanarm.

Python

import pandas as pd

# Load data
hers = pd.read_csv('../data/Chapter3/hersdata.csv')

X = hers.age.values
y = hers.SBP.values

Linear Regression in Python

scikit-learn

# Fit linear regression using scikit-learn
from sklearn.linear_model import LinearRegression

linreg = LinearRegression()
linreg.fit(X.reshape(-1, 1), y)

print(linreg.intercept_)
print(linreg.coef_)

103.62948825976596
[0.47172807]

statsmodels

# Fit linear regression using statsmodels
import statsmodels.api as sm

# Add constant to X
x_w_int = sm.add_constant(X) # adds a column of ones to the left of X

# Fit regression
linreg_sm = sm.OLS(y, x_w_int).fit()

print(linreg_sm.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.027
Model:                            OLS   Adj. R-squared:                  0.027
Method:                 Least Squares   F-statistic:                     77.21
Date:                Wed, 12 Jan 2022   Prob (F-statistic):           2.64e-18
Time:                        09:48:06   Log-Likelihood:                -12021.
No. Observations:                2763   AIC:                         2.405e+04
Df Residuals:                    2761   BIC:                         2.406e+04
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        103.6295      3.596     28.820      0.000      96.579     110.680
x1             0.4717      0.054      8.787      0.000       0.366       0.577
==============================================================================
Omnibus:                       75.882   Durbin-Watson:                   1.788
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               81.523
Skew:                           0.414   Prob(JB):                     1.98e-18
Kurtosis:                       3.150   Cond. No.                         675.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Bayesian Regression using PyMC3

# PyMC3 model
import pymc3 as pm

with pm.Model() as pymc3_model:
    # Define priors
    alpha = pm.Normal('alpha', mu=120, sd=30)
    beta = pm.Normal('beta', mu=0, sd=1)
    sigma = pm.HalfCauchy('sigma', 1)
    mu = alpha + beta * X

    # Define likelihood
    y_pred = pm.Normal('y_pred', mu=mu, sd=sigma, observed=y)
    
    # Sample from posterior
    posterior = pm.sample(2000, tune=1000)

# Print summary
pm.summary(posterior).round(3)

/var/folders/5x/pv9gwtw1041bc41pbp6j4j5m0000gn/T/ipykernel_19067/1390012775.py:15: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.
  posterior = pm.sample(2000, tune=1000)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, beta, alpha]

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 43 seconds.
The acceptance probability does not match the target. It is 0.909142303380742, but should be close to 0.8. Try to increase the number of tuning steps.
/Users/sami/miniconda3/envs/myenv/lib/python3.9/site-packages/arviz/data/io_pymc3.py:96: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
  warnings.warn(

Plot the posterior distributions

# Plot posterior distributions
pm.traceplot(posterior)

Julia

using CSV
using DataFrames
using Random:seed!
seed!(123)

# Read CSV into DataFrame object
hers = CSV.read("../data/Chapter3/hersdata.csv", DataFrame)

Linear Regression in Julia

using GLM

# Run a simple linear regression model
model_lm = lm(@formula(SBP ~ age), hers)

StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}}}}, Matrix{Float64}}

SBP ~ 1 + age

Coefficients:
───────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      t  Pr(>|t|)  Lower 95%   Upper 95%
───────────────────────────────────────────────────────────────────────────
(Intercept)  103.629      3.59573    28.82    <1e-99  96.5789    110.68
age            0.471728   0.0536838   8.79    <1e-17   0.366464    0.576993
───────────────────────────────────────────────────────────────────────────

Bayesian Regression using Turing.jl

using Turing

# Run a simple Bayesian linear regression model
@model function model_turing(X, y)
    #Priors
    # Intercept prior
    α ~ Normal(120, 30)
    # β coefficient prior
    β ~ Normal(0, 1)
    # Variance prior
    σ ~ Cauchy(0, 1)

    μ = α .+ X * β
    y ~ MvNormal(μ, σ)
end

model_turing (generic function with 1 method)

# Specify the model
model = model_turing(hers.age, hers.SBP)

# Run the model
chain = sample(model, NUTS(), 2_000)

┌ Info: Found initial step size
│   ϵ = 0.000390625
└ @ Turing.Inference /Users/sami/.julia/packages/Turing/uAz5c/src/inference/hmc.jl:195
[32mSampling: 100%|█████████████████████████████████████████| Time: 0:00:16[39m





Chains MCMC chain (2000×15×1 Array{Float64, 3}):

Iterations        = 1:1:2000
Number of chains  = 1
Samples per chain = 2000
parameters        = α, β, σ
internals         = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth

Summary Statistics
 [1m parameters [0m [1m     mean [0m [1m     std [0m [1m naive_se [0m [1m    mcse [0m [1m       ess [0m [1m    rhat [0m
 [90m     Symbol [0m [90m  Float64 [0m [90m Float64 [0m [90m  Float64 [0m [90m Float64 [0m [90m   Float64 [0m [90m Float64 [0m

           α   103.6575    3.5464     0.0793    0.1244    673.2840    0.9995
           β     0.4714    0.0529     0.0012    0.0018    664.6015    0.9995
           σ    18.7731    0.2625     0.0059    0.0049   1053.6198    1.0001

Quantiles
 [1m parameters [0m [1m    2.5% [0m [1m    25.0% [0m [1m    50.0% [0m [1m    75.0% [0m [1m    97.5% [0m
 [90m     Symbol [0m [90m Float64 [0m [90m  Float64 [0m [90m  Float64 [0m [90m  Float64 [0m [90m  Float64 [0m

           α   96.3406   101.2145   103.7226   106.0859   110.3310
           β    0.3743     0.4348     0.4703     0.5076     0.5810
           σ   18.2740    18.5973    18.7583    18.9489    19.3000




┌ Warning: timestamp of type Missing unknown
└ @ MCMCChains /Users/sami/.julia/packages/MCMCChains/IKF6o/src/chains.jl:364
┌ Warning: timestamp of type Missing unknown
└ @ MCMCChains /Users/sami/.julia/packages/MCMCChains/IKF6o/src/chains.jl:364

# Print summary statistics
summarystats(chain)

Summary Statistics
 [1m parameters [0m [1m     mean [0m [1m     std [0m [1m naive_se [0m [1m    mcse [0m [1m       ess [0m [1m    rhat [0m
 [90m     Symbol [0m [90m  Float64 [0m [90m Float64 [0m [90m  Float64 [0m [90m Float64 [0m [90m   Float64 [0m [90m Float64 [0m

           α   103.6575    3.5464     0.0793    0.1244    673.2840    0.9995
           β     0.4714    0.0529     0.0012    0.0018    664.6015    0.9995
           σ    18.7731    0.2625     0.0059    0.0049   1053.6198    1.0001

# Visualize the posterior
using MCMCChains, Plots, StatsPlots

plot(chain)

References

1. Vittinghoff E, ed. Regression Methods in Biostatistics: Linear, Logistic, Survival, and Repeated Measures Models. 2nd ed. Springer; 2012.
2. Hulley S, Grady D, Bush T, et al. Randomized trial of estrogen plus progestin for secondary prevention of coronary heart disease in postmenopausal women. JAMA. 1998;280(7):605-613. doi:10.1001/jama.280.7.605
3. McElreath R. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. 2nd ed. Taylor and Francis, CRC Press; 2020.
4. Kurz AS. Statistical Rethinking with Brms, Ggplot2, and the Tidyverse: Second Edition. https://bookdown.org/content/4857/
5. Clark M. Easy Bayes with Rstanarm and Brms. Accessed January 12, 2022. https://m-clark.github.io/easy-bayes/
6. Fusaroli R & Cox C. Workshop on Bayesian Inference: Priors and Workflow. https://4ccoxau.github.io/PriorsWorkshop/
7. Storopoli J. Bayesian Statistics with Julia and Turing. https://storopoli.io/Bayesian-Julia.