# Multilevel Models (a.k.a. Hierarchical Models)

Bayesian hierarchical models (also called multilevel models) are a statistical model written at multiple levels (hierarchical form) that estimates the parameters of the posterior distribution using the Bayesian approach. The sub-models combine to form the hierarchical model, and Bayes' theorem is used to integrate them with the observed data and to account for all the uncertainty that is present. The result of this integration is the posterior distribution, also known as an updated probability estimate, as additional evidence of the likelihood function is integrated together with the prior distribution of the parameters.

Hierarchical modeling is used when information is available at several different levels of observation units. The hierarchical form of analysis and organization helps to understand multiparameter problems and also plays an important role in the development of computational strategies.

Hierarchical models are mathematical statements that involve several parameters, so that the estimates of some parameters depend significantly on the values of other parameters. The figure below shows a hierarchical model in which there is a $$\phi$$ hyperparameter that parameterizes the parameters $$\theta_1, \theta_2, \dots, \theta_N$$ that are finally used to infer the posterior density of some variable of interest $$\mathbf{y} = y_1, y_2, \dots, y_N$$. Hierarchical Model

## When to use Multilevel Models?

Multilevel models are particularly suitable for research projects where participant data is organized at more than one level, i.e. nested data. Units of analysis are usually individuals (at a lower level) that are nested in contextual/aggregate units (at a higher level). An example is when we are measuring the performance of individuals and we have additional information about belonging to different groups such as sex, age group, hierarchical level, educational level or housing status.

There is a main assumption that cannot be violated in multilevel models which is exchangeability (de Finetti, 1974; Nau, 2001). Yes, this is the same assumption that we discussed in 2. What is Bayesian Statistics?. This assumption assumes that groups are exchangeable. The figure below shows a graphical representation of the exchangeability. The groups shown as "cups" that contain observations shown as "balls". If in the model's inferences, this assumption is violated, then multilevel models are not appropriate. This means that, since there is no theoretical justification to support exchangeability, the inferences of the multilevel model are not robust and the model can suffer from several pathologies and should not be used for any scientific or applied analysis.  Exchangeability – Images from Michael Betancourt

## Hyperpriors

As the priors of the parameters are sampled from another prior of the hyperparameter (upper-level's parameter), which are called hyperpriors. This makes one group's estimates help the model to better estimate the other groups by providing more robust and stable estimates.

We call the global parameters as population effects (or population-level effects, also sometimes called fixed effects) and the parameters of each group as group effects (or group-level effects, also sometimes called random effects). That is why multilevel models are also known as mixed models in which we have both fixed effects and random effects.

## Three Approaches to Multilevel Models

Multilevel models generally fall into three approaches:

1. Random-intercept model: each group receives a different intercept in addition to the global intercept.

2. Random-slope model: each group receives different coefficients for each (or a subset of) independent variable(s) in addition to a global intercept.

3. Random-intercept-slope model: each group receives both a different intercept and different coefficients for each independent variable in addition to a global intercept.

### Random-Intercept Model

The first approach is the random-intercept model in which we specify a different intercept for each group, in addition to the global intercept. These group-level intercepts are sampled from a hyperprior.

To illustrate a multilevel model, I will use the linear regression example with a Gaussian/normal likelihood function. Mathematically a Bayesian multilevel random-slope linear regression model is:

\begin{aligned} \mathbf{y} &\sim \text{Normal}\left( \alpha + \alpha_j + \mathbf{X} \cdot \boldsymbol{\beta}, \sigma \right) \\ \alpha &\sim \text{Normal}(\mu_\alpha, \sigma_\alpha) \\ \alpha_j &\sim \text{Normal}(0, \tau) \\ \boldsymbol{\beta} &\sim \text{Normal}(\mu_{\boldsymbol{\beta}}, \sigma_{\boldsymbol{\beta}}) \\ \tau &\sim \text{Cauchy}^+(0, \psi_{\alpha})\\ \sigma &\sim \text{Exponential}(\lambda_\sigma) \end{aligned}

The priors on the global intercept $$\alpha$$, global coefficients $$\boldsymbol{\beta}$$ and error $$\sigma$$, along with the Gaussian/normal likelihood on $$\mathbf{y}$$ are the same as in the linear regression model. But now we have new parameters. The first are the group intercepts prior $$\alpha_j$$ that denotes that every group $$1, 2, \dots, J$$ has its own intercept sampled from a normal distribution centered on 0 with a standard deviation $$\psi_\alpha$$. This group intercept is added to the linear predictor inside the Gaussian/normal likelihood function. The group intercepts' standard deviation $$\tau$$ have a hyperprior (being a prior of a prior) which is sampled from a positive-constrained Cauchy distribution (a special case of the Student-$$t$$ distribution with degrees of freedom $$\nu = 1$$) with mean 0 and standard deviation $$\sigma_\alpha$$. This makes the group-level intercept's dispersions being sampled from the same parameter $$\tau$$ which allows the model to use information from one group intercept to infer robust information regarding another group's intercept dispersion and so on.

This is easily accomplished with Turing:

using Turing
using Statistics: mean, std
using Random:seed!
seed!(123)

@model varying_intercept(X, idx, y; n_gr=length(unique(idx)), predictors=size(X, 2)) = begin
#priors
α ~ Normal(mean(y), 2.5 * std(y))       # population-level intercept
β ~ filldist(Normal(0, 2), predictors)  # population-level coefficients
σ ~ Exponential(1 / std(y))             # residual SD
#prior for variance of random intercepts
#usually requires thoughtful specification
τ ~ truncated(Cauchy(0, 2), 0, Inf)     # group-level SDs intercepts
αⱼ ~ filldist(Normal(0, τ), n_gr)       # group-level intercepts

#likelihood
ŷ = α .+ X * β .+ αⱼ[idx]
y ~ MvNormal(ŷ, σ)
end;

### Random-Slope Model

The second approach is the random-slope model in which we specify a different slope for each group, in addition to the global intercept. These group-level slopes are sampled from a hyperprior.

To illustrate a multilevel model, I will use the linear regression example with a Gaussian/normal likelihood function. Mathematically a Bayesian multilevel random-slope linear regression model is:

\begin{aligned} \mathbf{y} &\sim \text{Normal}\left( \alpha + \mathbf{X} \cdot \boldsymbol{\beta}_j \cdot \boldsymbol{\tau}, \sigma \right) \\ \alpha &\sim \text{Normal}(\mu_\alpha, \sigma_\alpha) \\ \boldsymbol{\beta}_j &\sim \text{Normal}(0, 1) \\ \boldsymbol{\tau} &\sim \text{Cauchy}^+(0, \psi_{\boldsymbol{\beta}})\\ \sigma &\sim \text{Exponential}(\lambda_\sigma) \end{aligned}

Here we have a similar situation from before with the same hyperprior, but now it is a hyperprior for the the group coefficients' standard deviation prior: $$\boldsymbol{\beta}_j$$. This makes the group-level coefficients's dispersions being sampled from the same parameter $$\tau$$ which allows the model to use information from one group coefficients to infer robust information regarding another group's coefficients dispersion and so on.

In Turing we can accomplish this as:

@model varying_slope(X, idx, y; n_gr=length(unique(idx)), predictors=size(X, 2)) = begin
#priors
α ~ Normal(mean(y), 2.5 * std(y))                   # population-level intercept
σ ~ Exponential(1 / std(y))                         # residual SD
#prior for variance of random slopes
#usually requires thoughtful specification
τ ~ filldist(truncated(Cauchy(0, 2), 0, Inf), n_gr) # group-level slopes SDs
βⱼ ~ filldist(Normal(0, 1), predictors, n_gr)       # group-level standard normal slopes

#likelihood
ŷ = α .+ X * βⱼ * τ
y ~ MvNormal(ŷ, σ)
end;

### Random-Intercept-Slope Model

The third approach is the random-intercept-slope model in which we specify a different intercept and slope for each group, in addition to the global intercept. These group-level intercepts and slopes are sampled from hyperpriors.

To illustrate a multilevel model, I will use the linear regression example with a Gaussian/normal likelihood function. Mathematically a Bayesian multilevel random-intercept-slope linear regression model is:

\begin{aligned} \mathbf{y} &\sim \text{Normal}\left( \alpha + \alpha_j + \mathbf{X} \cdot \boldsymbol{\beta}_j \cdot \boldsymbol{\tau}_{\boldsymbol{\beta}}, \sigma \right) \\ \alpha &\sim \text{Normal}(\mu_\alpha, \sigma_\alpha) \\ \alpha_j &\sim \text{Normal}(0, \tau_{\alpha}) \\ \boldsymbol{\beta}_j &\sim \text{Normal}(0, 1) \\ \tau_{\alpha} &\sim \text{Cauchy}^+(0, \psi_{\alpha})\\ \boldsymbol{\tau}_{\boldsymbol{\beta}} &\sim \text{Cauchy}^+(0, \psi_{\boldsymbol{\beta}})\\ \sigma &\sim \text{Exponential}(\lambda_\sigma) \end{aligned}

Here we have a similar situation from before with the same hyperpriors, but now we fused both random-intercept and random-slope together.

In Turing we can accomplish this as:

@model varying_intercept_slope(X, idx, y; n_gr=length(unique(idx)), predictors=size(X, 2)) = begin
#priors
α ~ Normal(mean(y), 2.5 * std(y))                    # population-level intercept
σ ~ Exponential(1 / std(y))                          # residual SD
#prior for variance of random intercepts and slopes
#usually requires thoughtful specification
τₐ ~ truncated(Cauchy(0, 2), 0, Inf)                 # group-level SDs intercepts
τᵦ ~ filldist(truncated(Cauchy(0, 2), 0, Inf), n_gr) # group-level slopes SDs
αⱼ ~ filldist(Normal(0, τₐ), n_gr)                   # group-level intercepts
βⱼ ~ filldist(Normal(0, 1), predictors, n_gr)        # group-level standard normal slopes

#likelihood
ŷ = α .+ αⱼ[idx] .+ X * βⱼ * τᵦ
y ~ MvNormal(ŷ, σ)
end;

## Example - Cheese Ratings

For our example, I will use a famous dataset called cheese (Boatwright, McCulloch & Rossi, 1999), which is data from cheese ratings. A group of 10 rural and 10 urban raters rated 4 types of different cheeses (A, B, C and D) in two samples. So we have $$4 \cdot 20 \cdot2 = 160$$ observations and 4 variables:

• cheese: type of cheese from A to D

• rater: id of the rater from 1 to 10

• background: type of rater, either rural or urban

• y: rating of the cheese

Ok let's read our data with CSV.jl and output into a DataFrame from DataFrames.jl:

using DataFrames, CSV, HTTP

url = "https://raw.githubusercontent.com/storopoli/Bayesian-Julia/master/datasets/cheese.csv"
describe(cheese)
4×7 DataFrame
Row │ variable    mean     min    median  max    nmissing  eltype
│ Symbol      Union…   Any    Union…  Any    Int64     DataType
─────┼───────────────────────────────────────────────────────────────
1 │ cheese               A              D             0  String1
2 │ rater       5.5      1      5.5     10            0  Int64
3 │ background           rural          urban         0  String7
4 │ y           70.8438  33     71.5    91            0  Int64

As you can see from the describe() output, the mean cheese ratings is around 70 ranging from 33 to 91.

In order to prepare the data for Turing, I will convert the Strings in variables cheese and background to Ints. Regarding cheese, I will create 4 dummy variables one for each cheese type; and background will be converted to integer data taking two values: one for each background type. My intent is to model background as a group both for intercept and coefficients. Take a look at how the data will look like for the first 5 observations:

for c in unique(cheese[:, :cheese])
cheese[:, "cheese_\$c"] = ifelse.(cheese[:, :cheese] .== c, 1, 0)
end

cheese[:, :background_int] = map(cheese[:, :background]) do b
b == "rural" ? 1 :
b == "urban" ? 2 : missing
end

first(cheese, 5)
5×9 DataFrame
Row │ cheese    rater  background  y      cheese_A  cheese_B  cheese_C  cheese_D  background_int
│ String1…  Int64  String7…    Int64  Int64     Int64     Int64     Int64     Int64
─────┼────────────────────────────────────────────────────────────────────────────────────────────
1 │ A             1  rural          67         1         0         0         0               1
2 │ A             1  rural          66         1         0         0         0               1
3 │ B             1  rural          51         0         1         0         0               1
4 │ B             1  rural          53         0         1         0         0               1
5 │ C             1  rural          75         0         0         1         0               1

Now let's us instantiate our model with the data. Here, I will specify a vector of Ints named idx to represent the different observations' group memberships. This will be used by Turing when we index a parameter with the idx, e.g. αⱼ[idx].

X = Matrix(select(cheese, Between(:cheese_A, :cheese_D)));
y = cheese[:, :y];
idx = cheese[:, :background_int];

The first model is the varying_intercept:

model_intercept = varying_intercept(X, idx, y)
chain_intercept = sample(model_intercept, NUTS(), MCMCThreads(), 2_000, 4)
summarystats(chain_intercept)  |> DataFrame  |> println
9×8 DataFrame
Row │ parameters  mean       std       naive_se    mcse        ess      rhat     ess_per_sec
│ Symbol      Float64    Float64   Float64     Float64     Float64  Float64  Float64
─────┼────────────────────────────────────────────────────────────────────────────────────────
1 │ α            70.9202   5.33286   0.0596232   0.127038    1557.59  0.99995      33.9841
2 │ β          3.23581  1.25415   0.0140218   0.023304    3224.24  1.00135      70.3476
3 │ β        -11.6113   1.24938   0.0139685   0.0228865   2929.34  1.00231      63.9134
4 │ β          7.20121  1.26355   0.0141269   0.0228724   3419.05  1.00111      74.598
5 │ β          1.23636  1.25291   0.0140079   0.0202693   3546.94  1.00106      77.3884
6 │ σ             5.99912  0.269684  0.00301516  0.00323454  5381.43  0.99992     117.414
7 │ τ             6.60777  6.68737   0.0747671   0.148843    1579.61  1.00167      34.4645
8 │ αⱼ        -3.67621  5.22431   0.0584096   0.125214    1520.63  1.00025      33.1775
9 │ αⱼ         3.49673  5.21951   0.0583559   0.125978    1529.61  1.00016      33.3735


Here we can see that the model has a population-level intercept α along with population-level coefficients βs for each cheese dummy variable. But notice that we have also group-level intercepts for each of the groups αⱼs. Specifically, αⱼ are the rural raters and αⱼ are the urban raters.

Now let's go to the second model, varying_slope:

model_slope = varying_slope(X, idx, y)
chain_slope = sample(model_slope, NUTS(), MCMCThreads(), 2_000, 4)
summarystats(chain_slope)  |> DataFrame  |> println
12×8 DataFrame
Row │ parameters  mean        std       naive_se    mcse        ess      rhat     ess_per_sec
│ Symbol      Float64     Float64   Float64     Float64     Float64  Float64  Float64
─────┼─────────────────────────────────────────────────────────────────────────────────────────
1 │ α           70.8114     4.79924   0.0536571   0.101646    2101.85  1.00087     17.8627
2 │ σ            6.54192    0.285938  0.00319688  0.00416645  5383.29  1.00037     45.7502
3 │ τ         6.05703    4.91366   0.0549364   0.137627    1281.8   1.00115     10.8934
4 │ τ         6.05027    5.00599   0.0559686   0.147636    1153.73  1.00123      9.80502
5 │ βⱼ[1,1]      0.243822   0.811422  0.00907197  0.0156275   3187.9   1.0015      27.0926
6 │ βⱼ[2,1]     -0.927608   1.03722   0.0115965   0.0264846   1644.87  1.0001      13.979
7 │ βⱼ[3,1]      0.562864   0.877836  0.0098145   0.0199541   2372.51  1.0006      20.1629
8 │ βⱼ[4,1]      0.0969922  0.799195  0.00893528  0.0136793   3056.19  1.00085     25.9732
9 │ βⱼ[1,2]      0.258359   0.807021  0.00902277  0.0156428   2795.54  1.00193     23.7581
10 │ βⱼ[2,2]     -0.917479   1.05359   0.0117794   0.024764    1797.64  1.00027     15.2774
11 │ βⱼ[3,2]      0.567299   0.902078  0.0100855   0.0175936   2502.98  1.0006      21.2717
12 │ βⱼ[4,2]      0.107948   0.797702  0.00891858  0.0132952   3685.84  1.00135     31.3244


Here we can see that the model has still a population-level intercept α. But now our population-level coefficients βs are replaced by group-level coefficients βⱼs along with their standard deviation τᵦs. Specifically βⱼ's first index denotes the 4 dummy cheese variables' and the second index are the group membership. So, for example βⱼ[1,1] is the coefficient for cheese_A and rural raters (group 1).

Now let's go to the third model, varying_intercept_slope:

model_intercept_slope = varying_intercept_slope(X, idx, y)
chain_intercept_slope = sample(model_intercept_slope, NUTS(), MCMCThreads(), 2_000, 4)
summarystats(chain_intercept_slope)  |> DataFrame  |> println
15×8 DataFrame
Row │ parameters  mean        std       naive_se    mcse        ess       rhat      ess_per_sec
│ Symbol      Float64     Float64   Float64     Float64     Float64   Float64   Float64
─────┼───────────────────────────────────────────────────────────────────────────────────────────
1 │ α           71.0626     7.83209   0.0875655   0.240281    1103.11   1.00238       4.80392
2 │ σ            5.87618    0.265073  0.00296361  0.00345513  7032.87   1.00004      30.6274
3 │ τₐ           6.70437    7.51895   0.0840645   0.224019    1372.05   1.0018        5.97511
4 │ τᵦ        6.36806    5.68283   0.063536    0.135849    1508.82   1.00161       6.57073
5 │ τᵦ        6.37695    5.73307   0.0640977   0.176432     963.377  1.00279       4.1954
6 │ αⱼ       -3.89769    5.94114   0.066424    0.219604     799.028  1.0036        3.47968
7 │ αⱼ        3.28275    5.92363   0.0662282   0.216349     811.2    1.00343       3.53268
8 │ βⱼ[1,1]      0.246034   0.831311  0.00929434  0.0176387   2208.44   1.0013        9.61752
9 │ βⱼ[2,1]     -0.888226   1.05528   0.0117984   0.0231345   2157.27   1.00124       9.39467
10 │ βⱼ[3,1]      0.556325   0.921074  0.0102979   0.0189355   2003.7    1.00009       8.72587
11 │ βⱼ[4,1]      0.0953832  0.779033  0.00870985  0.013333    3880.57   1.00038      16.8994
12 │ βⱼ[1,2]      0.264749   0.806497  0.00901691  0.0122747   3961.52   1.00036      17.252
13 │ βⱼ[2,2]     -0.889554   1.0405    0.0116331   0.0203335   2363.47   1.00235      10.2926
14 │ βⱼ[3,2]      0.550742   0.897587  0.0100353   0.0174808   2898.1    1.0006       12.6209
15 │ βⱼ[4,2]      0.103328   0.791032  0.00884401  0.0115532   4253.39   0.999846     18.523


Now we have fused the previous model in one. We still have a population-level intercept α. But now we have in the same model group-level intercepts for each of the groups αⱼs and group-level along with their standard deviation τₐ. We also have the coefficients βⱼs with their standard deviation τᵦs. The parameters are interpreted exactly as the previous cases.

Boatwright, P., McCulloch, R., & Rossi, P. (1999). Account-level modeling for trade promotion: An application of a constrained parameter hierarchical model. Journal of the American Statistical Association, 94(448), 1063–1073.

de Finetti, B. (1974). Theory of Probability (Volume 1). New York: John Wiley & Sons.

Nau, R. F. (2001). De Finetti was Right: Probability Does Not Exist. Theory and Decision, 51(2), 89–124. https://doi.org/10.1023/A:1015525808214