Practical Tutorial: Gaussian Linear Regression

Adding a deterministic component to our stochastic model of the mean.

Author

Doug Leasure

Published

11 Jun 2026, 11:35:57

Now that we can design and run a Bayesian model of the mean, let’s add a deterministic linear regression on the mean.

Note

This tutorial assumes that you have already completed the Model of the Mean Practical Tutorial.

Introduction

In this activity, we will run a Bayesian version of a Gaussian linear regression. We will simulate the data, write the model in Stan, and fit the model using MCMC.

Our primary aim is to provide a template for linear regression analysis that we can build on for various types of regression analysis. The power of Bayesian inference lies in its ability to handle bespoke model structures (based on data characteristics and our hypotheses about the system we are modelling) that may include non-Gaussian error structures and non-linear deterministic functions. A Gaussian linear regression is the reference point.

The Model

Our Gaussian linear regression model will take the following form:

\[ \begin{aligned} y_i &\sim Normal(\mu_i, \sigma) \\ \mu_i &= \alpha + \sum_{k=1}^K \beta_k x_{k,i} \end{aligned} \]

with priors, of course:

\[ \begin{aligned} \alpha &\sim Normal(0, 10) \\ \beta_k &\sim Normal(0, 10) \\ \sigma &\sim Uniform(0, 1e4) \end{aligned} \]

Can you draw a DAG for this model? Once you’ve drawn your DAG, you can compare with the DAG prepared below.

Notice that this DAG is very similar to the DAG from the model of the mean. Now we have added deterministic (dashed line) predictors of \(\mu_i\). This allows \(\mu\) to vary among our \(I\) replicate observations based on the effects \(\beta_k\) of some covariates \(x_k\).

Stan Code

Let’s take a look at how to write this model in Stan. This model is setup very similar to our model of the mean from activity 1.

data {
  int<lower=0> n; // sample size
  int<lower=0> K; // number of predictors
  vector[n] y; // response variable
  matrix[n, K] x; // predictor variables
}

parameters {
  real alpha; // intercept
  vector[K] beta; // slope
  real<lower=0> sigma; // residual variation
}

transformed parameters {
  vector[n] mu; // regression expected value
  
  mu = alpha + x * beta; // x * beta is a matrix multiplication
}

model {
  
  // likelihood
  y ~ normal(mu, sigma);
  
  // priors
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma ~ cauchy(0, 10);
}
Important

Save the Stan code into a file named gaussian_linear_regression.stan into your R working directory for this exercise.

Download Stan code

Let’s discuss the differences in each model block compared to our previous work, including introducing an important new parameter block in Stan: transformed parameters {}.

data

The data model block now includes an integer K that defines the number of predictor variables and a matrix \(x_{i,k}\) of values for our predictors. Notice that this matrix has a row for each sample n and a column for each predictor variable \(k\).

parameters

The parameters block now includes our regression coefficients: the intercept alpha (\(\alpha\)) and the effect size beta (\(\beta\)) of each predictor variable.

transformed parameters

This model block was not included in our model of the mean. This is where we can define parameters that are derived deterministically from other parameters and data in the model. In this case, our linear regression determines the mean \(\mu\) of the normal distribution used as our likelihood in the model{} block below.

Note that we are using a matrix multiplication to multiply the matrix \(x_{i,k}\) of predictor variables by our vector of regression coefficients \(\beta_k\). This notation provides flexibility to include more (or fewer) predictor variables in the model by simply modifying the data and without needing to rewrite the model code.

model

This block is very similar to our model of the mean. The difference is only that mu (\(\mu\)) is now derived in the transformed parameters block rather than having a prior in the model block, and we have added priors for alpha (\(\alpha\)) and beta (\(\beta\)).

We have introduced a new distribution as the prior for the standard deviation sigma (\(\sigma\)) in our likelihood. This is the half-cauchy distribution (“half” because sigma is defined to be greater than 0 in the parameters block). The Cauchy distribution is a t-distribution with one degree of freedom. This is similar to a normal distribution with “fat” tails (i.e. a higher probability of outliers). This is commonly used as the prior for standard deviation parameters, although there are other viable alternatives (e.g. uniform, half-normal).

Simulate Data

We provided code below to simulate data for our regression analysis. As usual, let’s start by setting up our R environment.

# cleanup
rm(list = ls())
gc()
cat("\014")
try(dev.off(), silent = T)

# check working directory
getwd()

# [optional] update/fix working directory
# setwd(file.path('/my/working/directory/model_of_the_mean'))

# directories
wd <- getwd()
dir.create(wd, showWarnings = F, recursive = T)

We’ll load our required Stan libraries.

# load libraries
library(cmdstanr)
library(posterior)
library(bayesplot)

Don’t forget to set a seed for random number generators!

# set seed for random number generators (important for reproducibility)
seed <- round(runif(1, 1, 1e6))
set.seed(seed)

The way the data are simulated will mirror how the model is written in the stan code. The parameter values are randomly drawn from distributions (analogous to defining the priors), the expected value of the response variable is derived deterministically from the linear regression, there is normally distributed random noise added on. We even used the same matrix multiplication trick (x * beta) so that you can convince yourself that this is equivalent to the longer form regression notation (beta1 * x1 + beta2 * x2).

#---- simulate data (Gaussian) ----#

# define the sample size
n <- 1e3

# number of covariates (use K > 1)
K <- 2
if (K < 2) {
  stop("Please define K as 2 or more.")
}

# create covariates
x <- matrix(NA, nrow = n, ncol = K)
for (k in 1:K) {
  x[, k] <- rnorm(n, 0, 1)
}

# define regression parameters (randomly)
alpha <- rnorm(1, 0, 1)
beta <- rnorm(K, 0, 3)
sigma <- runif(1, 0, 2)

# generate response variable
y <- alpha + as.vector(x %*% beta) + rnorm(n, 0, sigma)

Let’s just convince ourselves that the matrix multiplication is doing what is should:

# note: %*% is a matrix multiplication of a matrix of x values times a vector of beta values
# convince yourself that (x %*% beta) is equal to (beta1 * x1 + beta2 * x2)
head(x[, 1:2] %*% beta[1:2])
            [,1]
[1,] -0.07102591
[2,]  0.10026241
[3,]  0.24951907
[4,] -0.47631950
[5,]  0.10065918
[6,]  0.28333667
head(beta[1] * x[, 1] + beta[2] * x[, 2])
[1] -0.07102591  0.10026241  0.24951907 -0.47631950  0.10065918  0.28333667

We can see that we get the exact same results with the matrix multiplication as we do with the element-wise multiplications that are then summed as a second step.

Note

In Stan, the * symbol is a matrix multiplication operator for matrix and vector data types. The .* operator is an elementwise multiplication for two vectors.

Now let’s visualise the relationship between our covariates and our outcome:

# visualise covariate relationships with response variable
plot(y ~ x[, 1])

plot(y ~ x[, 2])

Looks like we’ve got reasonable linear relationships with plenty of random noise to make the analysis interesting.

And review a few summary statistics for our response variable:

# summary statistics of response variable
summary(y)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-3.6170 -0.4169  0.7329  0.6768  1.7115  5.1825 
# density plot of response variable
plot(density(y))
rug(y)

Here we can see that we do indeed have a normally-distributed response variable, and we have confirmed it’s relationships with our simulated covariates.

A few notes

Notice that the regression coefficients (alpha and the betas) are simulated randomly so that they change each time you re-simulate the data. This is useful to see how well the model can recover different combinations of these “true” parameter values in its posterior parameter estimates.

Note that we can easily modify our sample size and the number of predictor variables in our data/model by changing the values in the simulation code. This is useful to explore how the model results might change with small sample sizes and different numbers of covariates/predictors.

Run MCMC

This section is identical to our previous work with the model of the mean. We’re just running MCMC on our model, as usual.

#---- run Bayesian model (CmdStandR) ----#

# format data as a list (required by stan)
md <- list(
  n = n,
  K = K,
  y = as.vector(y),
  x = x,
  alpha_true = alpha,
  beta_true = beta,
  sigma_true = sigma,
  seed = seed # important to set and save seed for reproducibility
)

# save data to disk
saveRDS(
  object = md,
  file = file.path(wd, "md.rds")
)


# define location of the stan model file
model_file <- file.path(wd, "gaussian_linear_regression.stan")

# compile the stan model
mod <- cmdstan_model(model_file)


# function to generate initial values for each parameter in the model
init_generator <- function(md = md, chain_id = 1) {
  result <- list()

  result[["alpha"]] <- runif(1, -10, 10)
  result[["beta"]] <- runif(md$K, -10, 10)
  result[["sigma"]] <- runif(1, 0, 1e2)

  return(result)
}


# MCMC config
chains <- 4
warmup <- 1e3
samples <- 2e3
inits <- lapply(1:chains, function(id) init_generator(md = md, chain_id = id))

Let’s just take a quick look at our initials to make sure we understand the values being generated by our init_generator() function. Remember, this is provided a starting value for the first MCMC iteration for each parameter, and we need set of indpendent initials for each of our four chains.

inits
[[1]]
[[1]]$alpha
[1] -9.130684

[[1]]$beta
[1]  4.343447 -4.150267

[[1]]$sigma
[1] 87.4989


[[2]]
[[2]]$alpha
[1] -3.098959

[[2]]$beta
[1] -4.8577970 -0.7693757

[[2]]$sigma
[1] 59.7614


[[3]]
[[3]]$alpha
[1] 8.963895

[[3]]$beta
[1] -5.592822  7.971530

[[3]]$sigma
[1] 71.36906


[[4]]
[[4]]$alpha
[1] 0.1216526

[[4]]$beta
[1] -5.563133  5.253487

[[4]]$sigma
[1] 2.697108
# run MCMC
fit <- mod$sample(
  data = md,
  parallel_chains = chains,
  init = inits,
  iter_sampling = samples,
  iter_warmup = warmup,
  save_warmup = TRUE,
  seed = md$seed
)

# save model to disk
fit$save_object(file = file.path(wd, "fit.rds"))

Did Stan return any warnings at the end of the model run? If so, follow the links provided by Stan to learn more about the issues.

Model Diagnostics

Let’s start by looking at our traceplots.

# traceplots to visually inspect MCMC chains
bayesplot::mcmc_trace(
  fit$draws(),
  pars = c(
    "alpha",
    paste0("beta[", 1:md$K, "]"),
    "sigma"
  )
)

What do you think? Do we have convergence?

If the traceplots visually appear converged and there are no obvious issues, let’s formally confirm convergence for each parameter with the r-hat statistics (i.e. we’re looking for r-hat < 1.01 for publication quality convergence or r-hat < 1.1 for a reliable exploratory model).

# Check r-hat for convergence
posterior::rhat(fit$draws("alpha"))
[1] 1.000466
posterior::rhat(fit$draws("beta[1]"))
[1] 1.000008
posterior::rhat(fit$draws("beta[2]"))
[1] 0.9999318
posterior::rhat(fit$draws("sigma"))
[1] 1.000144

Are all of our parameters converged?

Compare Results to Simulated “Truth”

Let’s remind ourselves of the true population-level parameter values that we defined in our simulation.

# parameters from the true unobserved population (i.e. from simulation)
alpha
[1] 0.6818281
beta
[1] -0.2068608 -0.3324488
sigma
[1] 1.486319

Now, let’s compare them to the summary statistics from the posterior distribution for these parameters based on our model and simulated sample data.

# summary of model results
pars <- c("alpha", "beta", "sigma")
fit$summary(variables = pars)
# A tibble: 4 × 10
  variable   mean median     sd    mad     q5     q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl> <dbl>    <dbl>    <dbl>
1 alpha     0.694  0.694 0.0468 0.0455  0.617  0.771  1.000    9789.    6144.
2 beta[1]  -0.159 -0.159 0.0430 0.0431 -0.230 -0.0883 1.00     8478.    6431.
3 beta[2]  -0.378 -0.378 0.0471 0.0478 -0.456 -0.301  1.00     9025.    5743.
4 sigma     1.47   1.47  0.0329 0.0326  1.42   1.52   1.00     9195.    6020.

Do our credible intervals include the true parameter values? Should the true parameter values always be inside the credible intervals?

Are the credible intervals reflected in columns q5 and q95 giving us 95% CI or 90% CI?

Now let’s just get a nice visual of the marginal posteriors for our parameters using density plots.

# visualise marginal posterior distributions for parameters as density plots
pars <- c("alpha", paste0("beta[", 1:md$K, "]"), "sigma")
bayesplot::mcmc_areas(fit$draws(), pars = pars, prob = 0.95)

# ... as histograms
bayesplot::mcmc_hist(fit$draws(), pars = pars)
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Full Analysis Scripts

Stan Model

gaussian_linear_regression.stan

data {
  int<lower=0> n; // sample size
  int<lower=0> K; // number of predictors
  vector[n] y; // response variable
  matrix[n, K] x; // predictor variables
}

parameters {
  real alpha; // intercept
  vector[K] beta; // slope
  real<lower=0> sigma; // residual variation
}

transformed parameters {
  vector[n] mu; // regression expected value
  
  mu = alpha + x * beta; // x * beta is a matrix multiplication
}

model {
  
  // likelihood
  y ~ normal(mu, sigma);
  
  // priors
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma ~ cauchy(0, 10);
}

R Analysis

gaussian_linear_regression.R

# cleanup
rm(list = ls())
gc()
cat("\014")
try(dev.off(), silent = T)

# check working directory
getwd()

# [optional] update/fix working directory
# setwd(file.path('/my/working/directory/model_of_the_mean'))

# directories
wd <- getwd()
dir.create(wd, showWarnings = F, recursive = T)

# load libraries
library(cmdstanr)
library(posterior)
library(bayesplot)

# set seed for random number generators (important for reproducibility)
seed <- round(runif(1, 1, 1e6))
set.seed(seed)

#---- simulate data (Gaussian) ----#

# define the sample size
n <- 1e3

# number of covariates (use K > 1)
K <- 2
if (K < 2) {
  stop("Please define K as 2 or more.")
}

# create covariates
x <- matrix(NA, nrow = n, ncol = K)
for (k in 1:K) {
  x[, k] <- rnorm(n, 0, 1)
}

# define regression parameters (randomly)
alpha <- rnorm(1, 0, 1)
beta <- rnorm(K, 0, 3)
sigma <- runif(1, 0, 2)

# generate response variable
y <- alpha + as.vector(x %*% beta) + rnorm(n, 0, sigma)

# note: %*% is a matrix multiplication of a matrix of x values times a vector of beta values
# convince yourself that (x %*% beta) is equal to (beta1 * x1 + beta2 * x2)
head(x[, 1:2] %*% beta[1:2])
head(beta[1] * x[, 1] + beta[2] * x[, 2])

# visualise covariate relationships with response variable
plot(y ~ x[, 1])
plot(y ~ x[, 2])

# summary statistics of response variable
summary(y)

# density plot of response variable
plot(density(y))
rug(y)

#---- run Bayesian model (CmdStandR) ----#

# format data as a list (required by stan)
md <- list(
  n = n,
  K = K,
  y = as.vector(y),
  x = x,
  alpha_true = alpha,
  beta_true = beta,
  sigma_true = sigma,
  seed = seed # important to set and save seed for reproducibility
)

# save data to disk
saveRDS(
  object = md,
  file = file.path(wd, "md.rds")
)


# define location of the stan model file
model_file <- file.path(wd, "gaussian_linear_regression.stan")

# compile the stan model
mod <- cmdstan_model(model_file)


# function to generate initial values for each parameter in the model
init_generator <- function(md = md, chain_id = 1) {
  result <- list()

  result[["alpha"]] <- runif(1, -10, 10)
  result[["beta"]] <- runif(md$K, -10, 10)
  result[["sigma"]] <- runif(1, 0, 1e2)

  return(result)
}


# MCMC config
chains <- 4
warmup <- 1e3
samples <- 2e3
inits <- lapply(1:chains, function(id) init_generator(md = md, chain_id = id))

# run MCMC
fit <- mod$sample(
  data = md,
  parallel_chains = chains,
  init = inits,
  iter_sampling = samples,
  iter_warmup = warmup,
  save_warmup = TRUE,
  seed = md$seed
)

# save model to disk
fit$save_object(file = file.path(wd, "fit.rds"))

# traceplots to visually inspect MCMC chains
bayesplot::mcmc_trace(
  fit$draws(),
  pars = c(
    "alpha",
    paste0("beta[", 1:md$K, "]"),
    "sigma"
  )
)

# Check r-hat for convergence
posterior::rhat(fit$draws("alpha"))
posterior::rhat(fit$draws("beta[1]"))
posterior::rhat(fit$draws("beta[2]"))
posterior::rhat(fit$draws("sigma"))

# parameters from the true unobserved population (i.e. from simulation)
alpha
beta
sigma

# summary of model results
pars <- c("alpha", "beta", "sigma")
fit$summary(variables = pars)

# visualise marginal posterior distributions for parameters as density plots
pars <- c("alpha", paste0("beta[", 1:md$K, "]"), "sigma")
bayesplot::mcmc_areas(fit$draws(), pars = pars, prob = 0.95)

# ... as histograms
bayesplot::mcmc_hist(fit$draws(), pars = pars)
© 2026 Real Good Research.
All rights reserved.
Source code for this website
is available on GitHub.
Open Source Content: CC-BY 4.0
Open Source Code: MIT