Model selection of Bayesian linear mixed models using the marginal likelihood

Software

Multilevel Modeling

Model Specification

Author

Carl F. Falk

Published

June 10, 2026

Overview

Background

Suppose a researcher is interested in estimating a linear mixed effect model (aka, multilevel model or hierarchical linear model) in R. Those who are used to frequentist estimation methods may use lmer() with the lme4 package (Bates et al, 2015). If using a Bayesian approach, brms (Bürkner, 2017) could be used and has similar syntax.

Model selection using Bayesian information criteria

Using Bayesian estimation, there is a small complication if one wants to use something familiar, like information criteria (e.g., AIC or BIC), for model selection. There are a few that can be used in this context:

DIC
WAIC
LOO-CV

However, with multilevel models there is more than one likelihood that can be used to compute these:

Conditional
Marginal¹

What should be used? 🤔

The marginal likelihood, coupled with WAIC or LOO-CV (which are asymptotically equivalent under the marginal likelihood) are sometimes recommended by methodologists (e.g., Merkle et al.,).

The problem

The recommended approach to computing information criteria are rarely, if ever, used by applied researchers. Why?

At the time of this writing, I am unaware of any software package that makes these easily available. 😒

brms produces WAIC and LOO-CV using the conditional likelihood ❌
Blimp appears to be the only software package that allows for computation of WAIC using the marginal likelihood. 👍
blavaan (Merkle & Rosseel) could also be used, but does not accommodate random slopes; given it is also software for structural equation modeling, dealing with many level 1 units and unequal cluster sizes may also be tricky. ❌

So, it may seem a bit strange to hear that the marginal likelihood should be used, if there really aren’t any many software options. 🤷

Purpose of this post

This post is intended to nudge researchers towards being able to use the marginal likelihood to compute WAIC or LOO-CV for linear mixed effects models.

Note that this tutorial might not make it super easy such that anyone could do it, but I am hoping this will help the methodologically-inclined applied researcher. Feel free to get in touch with questions so I can improve this work!

It contains example code related to the following paper (currently a preprint):

Devine, S., Falk, C. F., & Fujimoto, K. A. Comparing three predictive information criteria for linear multilevel Bayesian regression model selection. Preprint: https://psyarxiv.com/p2n8a

In particular, use of a custom rstan file with some helper functions provides lmer-like syntax.

A big thanks to Sean Devine and Ken Fujimoto, without whom this research would have not happened. 🙏

Setup

Prerequisites

I assume:

You are familiar with estimating linear mixed models with brms or lme4
- I will not go into detail regarding the syntax of estimated parameters for these models, nor model equations
You can install stan and related packages, brms and rstan
- Occasionally this is not as straightforward as it is with other R packages
If you want custom modifications to prior distributions…
- I assume you can look at the relevant stan code and change these
- Currently, we wrote it in an attempt to mimic brms default priors
  - LKJ(1) prior on random effects correlation matrix
  - Uniform priors on the fixed effect coefficients
  - half-\(t\) (df = 3; dispersion = 2.5) on variances
You have some familiarity with Bayesian estimation
- If you want a tutorial, this one looks like a good start: Depaoli et al (2017)

Packages and code

fx.R is available in this GitHub repository²
- It also has a number of packages that should be installed: mnormt, MCMCpack, matrixcalc, brms, lme4, dplyr, Matrix, numDeriv, mvtnorm, loo
- This file contains helper functions to set up the model
We will also use this stan file, which is also in the repository
- mlmmarg.stan
- This contains the code to do estimation using the marginal likelihood
If you run this code interactively (i.e., not in an RMarkdown file)
- Probably the rstudioapi package also needs to be installed

# load packages
library(brms)
library(rstan)
library(loo)
library(posterior)

source("fx.R") # custom functions

Model estimation setup. These are used later; if you want to changes the number of chains, processing cores, number of iterations, then modify here — or modify the code directly when used later on.

# This was setup from simulation study...
nchains <- 2 # two independent chains
ncores <- 2 # two cores
iter <- 5000

Load data

The data file example_dataset_77.csv is also in the repository
- It was a simulated dataset from the paper; all predictors are at level 1
- X1 - predictor
- X2 - predictor
- X3 - predictor
- g - cluster index (e.g., level 2 ID)
- y - outcome variable

dat <- read.csv("example_dataset_77.csv")

Setup of models

The paper had the following 5 models. This is analogous to syntax that you would give to lmer(). Recall, this post assumes you are familiar with lmer(), but briefly:

Model A: two predictors, random intercept and one random slope for X1
Model B: two predictors, only a random intercept
Model C: two predictors, random intercept and two random slopes for X1 and X2
Model D: one predictor, random intercept and random slope for X1
Model E: three predictors, random intercept and random slope for X1

Model A was the data generating model (i.e., true model).

modA <- y ~ 1 + X1 + X2 + (X1|g)
modB <- y ~ 1 + X1 + X2 + (1|g)
modC <- y ~ 1 + X1 + X2 + (X1+X2|g)
modD <- y ~ 1 + X1 + (X1|g)
modE <- y ~ 1 + X1 + X2 + X3 + (X1|g)

Custom Marginal Likelihood

`brms`

For sanity, we also estimate all models using brms for comparison.

Estimate Models

Note silent=2 and refresh=0 will suppress output as the model is estimated. You can change these if you like (e.g., remove silent=2 but have refresh=100 if you want updates every 100 iterations).

Note that brms as well as the custom code will need to compile before estimating the models. 🥱

So, you can either twiddle your thumbs and have refresh set to something smallish so that you can stare at it estimating 😐, or you can do something else for a little while (no procrastination or time-wasting).

modA.brms = brm(modA, data=dat, chains=nchains, iter=iter, silent = 2,
                cores=ncores, refresh = 0)
# stancode(modA.brms) # to look at underlying stan code, if you like

modB.brms = brm(modB, data=dat, chains=nchains, iter=iter, silent = 2,
                cores=ncores, refresh = 0)

modC.brms = brm(modC, data=dat, chains=nchains, iter=iter, silent = 2,
                cores=ncores, refresh = 0)

modD.brms = brm(modD, data=dat, chains=nchains, iter=iter, silent = 2,
                cores=ncores, refresh = 0)

modE.brms = brm(modE, data=dat, chains=nchains, iter=iter, silent = 2,
                cores=ncores, refresh = 0)

Custom Marginal Likelihood Code

Here is illustration of estimation of the same 5 models, but using custom stan code.

Estimate Models

mod.mat.mlm is a helper function in fx.R that uses the syntax for the model and the dataset to create model matrices that will be used by the underlying stan code. It might require:
- The data are sorted in order by the cluster index (e.g., g here); you should probably do that before using this code; or we should do something to make the code more foolproof. ⚠️
mlmmarg.stan is then used as the stan code to estimate the model using the stan function

matsA <- mod.mat.mlm(modA, data=dat)
modA.custom <- stan(file='mlmmarg.stan', data = matsA,
                    chains = nchains, seed=5297,
                    iter = iter, refresh = 0, cores = ncores)

matsB <- mod.mat.mlm(modB, data=dat)
modB.custom <- stan(file='mlmmarg.stan', data = matsB,
                    chains = nchains, seed=5297,
                    iter = iter, refresh = 0, cores = ncores)

matsC <- mod.mat.mlm(modC, data=dat)
modC.custom <- stan(file='mlmmarg.stan', data = matsC,
                    chains = nchains, seed=5297,
                    iter = iter, refresh = 0, cores = ncores)

matsD <- mod.mat.mlm(modD, data=dat)
modD.custom <- stan(file='mlmmarg.stan', data = matsD,
                    chains = nchains, seed=5297,
                    iter = iter, refresh = 0, cores = ncores)

matsE <- mod.mat.mlm(modE, data=dat)
modE.custom <- stan(file='mlmmarg.stan', data = matsE,
                    chains = nchains, seed=5297,
                    iter = iter, refresh = 0, cores = ncores)

Look at Diagnostics

You ought to look at some diagnostics before trusting your model. For this, I will defer to tutorials on Bayesian data analysis. But, the following two are often used to check for convergence or that something has not gone completely awry.

Rhats

# error sd, coefficients, random effect covariance matrix
pars <- c("sig", "betas", "Tau")
max(abs(summary(modA.custom, pars=pars)$summary[,"Rhat"]))

[1] 1.000005

max(abs(summary(modB.custom, pars=pars)$summary[,"Rhat"]))

[1] 1.000446

max(abs(summary(modC.custom, pars=pars)$summary[,"Rhat"]))

[1] 1.000493

max(abs(summary(modD.custom, pars=pars)$summary[,"Rhat"]))

[1] 1.000311

max(abs(summary(modE.custom, pars=pars)$summary[,"Rhat"]))

[1] 0.9999194

Traceplots

rstan::traceplot(modA.custom, pars = pars)

rstan::traceplot(modB.custom, pars = pars)

rstan::traceplot(modC.custom, pars = pars)

rstan::traceplot(modD.custom, pars = pars)

rstan::traceplot(modE.custom, pars = pars)

Compare to `brms`

Parameter Estimates

Some code could automate comparisons between parameters could have been done and probably a helper function could make the output nicer/easier to read, but I did not write any yet. 🫤 Here, I just have output from each model from brms and the custom code printed, one model at a time.

Briefly in custom code,

gssig is a vector of standard deviations for the random effects
- Intercept, slope for X1, etc., in that order
gsmat is a the correlation matrix among the random effects
- Intercept, slope for X1, etc., in that order
betas contains the fixed effects coefficients
- Intercept, slope for X1, X2, etc., in that order
sig is the standard deviation for the error at level 1

As an example for Model A, fixed effects:

Parameter	In `brms`	In custom code
Intercept	Intercept, 1.01	betas[1,1], 1.015
Coefficient for X1	X1, .19	betas[2,1], .187
Coefficient for X2	X2, .40	betas[3,1], .404

As an example for Model A, random effects:

Parameter	In `brms`	In custom code
Random intercept SD	sd(Intercept), .19	gssig[1], .19
Random slope SD	sd(X1), .20	gssig[2], .20
Correlation b/w intercept and slope	cor(Intercept, X1), .25	gsmat[1,2], .248
Level 1 error term	sigma, .85	sig, .85

pars <- c("gssig", "gsmat", "betas","sig")

summary(modA.custom, pars = pars)$summary

                mean      se_mean           sd       2.5%       25%       50%
gssig[1]   0.1918110 3.080043e-04 2.810884e-02  0.1418738 0.1721204 0.1900512
gssig[2]   0.2049995 3.239594e-04 3.030804e-02  0.1513372 0.1834301 0.2033433
gsmat[1,1] 1.0000000          NaN 0.000000e+00  1.0000000 1.0000000 1.0000000
gsmat[1,2] 0.2484226 2.118150e-03 1.850621e-01 -0.1337348 0.1232743 0.2586771
gsmat[2,1] 0.2484226 2.118150e-03 1.850621e-01 -0.1337348 0.1232743 0.2586771
gsmat[2,2] 1.0000000 1.271269e-18 9.195014e-17  1.0000000 1.0000000 1.0000000
betas[1,1] 1.0151835 3.529104e-04 3.154864e-02  0.9532490 0.9946503 1.0144379
betas[2,1] 0.1870510 3.994190e-04 3.396616e-02  0.1195645 0.1645170 0.1869748
betas[3,1] 0.4039262 1.910961e-04 1.706209e-02  0.3702468 0.3922990 0.4040235
sig        0.8520069 1.518920e-04 1.209477e-02  0.8286582 0.8439104 0.8519799
                 75%     97.5%    n_eff      Rhat
gssig[1]   0.2097149 0.2537442 8328.608 1.0000714
gssig[2]   0.2239143 0.2705220 8752.549 0.9996923
gsmat[1,1] 1.0000000 1.0000000      NaN       NaN
gsmat[1,2] 0.3802265 0.5919233 7633.466 0.9998811
gsmat[2,1] 0.3802265 0.5919233 7633.466 0.9998811
gsmat[2,2] 1.0000000 1.0000000 5231.548 0.9995999
betas[1,1] 1.0359007 1.0780544 7991.574 0.9997140
betas[2,1] 0.2102328 0.2535111 7231.616 0.9998718
betas[3,1] 0.4157825 0.4371860 7971.881 0.9998050
sig        0.8599656 0.8762982 6340.527 0.9999245

summary(modA.brms)

 Family: gaussian 
  Links: mu = identity 
Formula: y ~ 1 + X1 + X2 + (X1 | g) 
   Data: dat (Number of observations: 2500) 
  Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 5000

Multilevel Hyperparameters:
~g (Number of levels: 50) 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)         0.19      0.03     0.14     0.25 1.00     2475     2646
sd(X1)                0.20      0.03     0.15     0.27 1.00     2285     2989
cor(Intercept,X1)     0.25      0.18    -0.13     0.59 1.00     1302     2198

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.01      0.03     0.95     1.08 1.00     2466     3561
X1            0.19      0.03     0.12     0.26 1.00     2808     3348
X2            0.40      0.02     0.37     0.44 1.00     9872     3183

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.85      0.01     0.83     0.88 1.00     8218     2922

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).

summary(modB.custom, pars = pars)$summary

                mean      se_mean         sd      2.5%       25%       50%
gssig[1]   0.1923681 0.0003690778 0.02786288 0.1422447 0.1735016 0.1907334
gsmat[1,1] 1.0000000          NaN 0.00000000 1.0000000 1.0000000 1.0000000
betas[1,1] 1.0130818 0.0004613252 0.03263846 0.9496920 0.9914550 1.0130321
betas[2,1] 0.1893812 0.0002292444 0.01786855 0.1541005 0.1772779 0.1890396
betas[3,1] 0.4010111 0.0002196821 0.01749074 0.3664365 0.3892013 0.4013143
sig        0.8736130 0.0001618714 0.01258730 0.8492467 0.8652348 0.8734514
                 75%     97.5%    n_eff      Rhat
gssig[1]   0.2090517 0.2531057 5699.232 1.0002389
gsmat[1,1] 1.0000000 1.0000000      NaN       NaN
betas[1,1] 1.0344937 1.0786197 5005.471 1.0000738
betas[2,1] 0.2014761 0.2239965 6075.491 0.9997895
betas[3,1] 0.4124569 0.4347344 6339.090 0.9998775
sig        0.8818100 0.8991470 6046.789 1.0004460

summary(modB.brms)

 Family: gaussian 
  Links: mu = identity 
Formula: y ~ 1 + X1 + X2 + (1 | g) 
   Data: dat (Number of observations: 2500) 
  Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 5000

Multilevel Hyperparameters:
~g (Number of levels: 50) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.19      0.03     0.14     0.26 1.00     2081     3076

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.01      0.03     0.95     1.08 1.00     3002     3418
X1            0.19      0.02     0.15     0.22 1.00    10349     3228
X2            0.40      0.02     0.37     0.44 1.00     9334     3585

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.87      0.01     0.85     0.90 1.00     9869     3692

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).

summary(modC.custom, pars = pars)$summary

                  mean      se_mean           sd          2.5%          25%
gssig[1]    0.19160290 3.313331e-04 2.897480e-02  0.1401200106  0.171458762
gssig[2]    0.20565102 3.823812e-04 3.000967e-02  0.1515489798  0.184771660
gssig[3]    0.02332682 2.615383e-04 1.765683e-02  0.0009513602  0.009441122
gsmat[1,1]  1.00000000          NaN 0.000000e+00  1.0000000000  1.000000000
gsmat[1,2]  0.23936578 2.267083e-03 1.807565e-01 -0.1299656368  0.115940309
gsmat[1,3] -0.11973518 5.504373e-03 4.735594e-01 -0.9008676107 -0.504083418
gsmat[2,1]  0.23936578 2.267083e-03 1.807565e-01 -0.1299656368  0.115940309
gsmat[2,2]  1.00000000 1.301942e-18 9.161434e-17  1.0000000000  1.000000000
gsmat[2,3]  0.09839165 5.842529e-03 4.712396e-01 -0.8168947751 -0.254230038
gsmat[3,1] -0.11973518 5.504373e-03 4.735594e-01 -0.9008676107 -0.504083418
gsmat[3,2]  0.09839165 5.842529e-03 4.712396e-01 -0.8168947751 -0.254230038
gsmat[3,3]  1.00000000 2.074288e-18 8.544524e-17  1.0000000000  1.000000000
betas[1,1]  1.01449142 4.051494e-04 3.231456e-02  0.9515671688  0.992450492
betas[2,1]  0.18730244 4.294914e-04 3.493284e-02  0.1180532339  0.164203593
betas[3,1]  0.40381828 2.105561e-04 1.738987e-02  0.3692833294  0.392414084
sig         0.85194345 1.459069e-04 1.240023e-02  0.8278588918  0.843579485
                   50%        75%      97.5%    n_eff      Rhat
gssig[1]    0.18930240 0.21017484 0.25350858 7647.354 0.9998031
gssig[2]    0.20359782 0.22406288 0.27189200 6159.271 0.9999291
gssig[3]    0.01977891 0.03297286 0.06628364 4557.793 1.0000400
gsmat[1,1]  1.00000000 1.00000000 1.00000000      NaN       NaN
gsmat[1,2]  0.24480165 0.36771556 0.57529230 6357.017 0.9999002
gsmat[1,3] -0.15818707 0.23515061 0.82700745 7401.730 1.0001531
gsmat[2,1]  0.24480165 0.36771556 0.57529230 6357.017 0.9999002
gsmat[2,2]  1.00000000 1.00000000 1.00000000 4951.575 0.9995999
gsmat[2,3]  0.12676637 0.47321536 0.89376744 6505.516 1.0001128
gsmat[3,1] -0.15818707 0.23515061 0.82700745 7401.730 1.0001531
gsmat[3,2]  0.12676637 0.47321536 0.89376744 6505.516 1.0001128
gsmat[3,3]  1.00000000 1.00000000 1.00000000 1696.828 0.9995999
betas[1,1]  1.01456233 1.03652346 1.07865477 6361.597 0.9996710
betas[2,1]  0.18756873 0.21026382 0.25797154 6615.440 0.9999144
betas[3,1]  0.40387583 0.41533281 0.43792436 6821.142 1.0000721
sig         0.85170172 0.86023466 0.87635550 7222.842 1.0000742

summary(modC.brms)

 Family: gaussian 
  Links: mu = identity 
Formula: y ~ 1 + X1 + X2 + (X1 + X2 | g) 
   Data: dat (Number of observations: 2500) 
  Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 5000

Multilevel Hyperparameters:
~g (Number of levels: 50) 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)         0.19      0.03     0.14     0.25 1.00     1827     3052
sd(X1)                0.20      0.03     0.15     0.27 1.00     2210     2870
sd(X2)                0.02      0.02     0.00     0.07 1.00     2775     2980
cor(Intercept,X1)     0.25      0.18    -0.12     0.57 1.00     1305     2687
cor(Intercept,X2)    -0.12      0.46    -0.89     0.81 1.00     5760     3624
cor(X1,X2)            0.11      0.47    -0.80     0.88 1.00     6051     3705

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.01      0.03     0.95     1.08 1.00     1596     2842
X1            0.19      0.03     0.12     0.26 1.00     1909     2731
X2            0.40      0.02     0.37     0.44 1.00     7832     3561

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.85      0.01     0.83     0.88 1.00     7283     3611

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).

summary(modD.custom, pars = pars)$summary

                mean      se_mean           sd       2.5%        25%       50%
gssig[1]   0.1935204 4.130046e-04 3.090245e-02  0.1392984 0.17193187 0.1916898
gssig[2]   0.1958051 3.898703e-04 3.075364e-02  0.1418536 0.17419381 0.1937933
gsmat[1,1] 1.0000000          NaN 0.000000e+00  1.0000000 1.00000000 1.0000000
gsmat[1,2] 0.1809829 2.630157e-03 1.985634e-01 -0.2254724 0.04764889 0.1875613
gsmat[2,1] 0.1809829 2.630157e-03 1.985634e-01 -0.2254724 0.04764889 0.1875613
gsmat[2,2] 1.0000000 1.313643e-18 9.123677e-17  1.0000000 1.00000000 1.0000000
betas[1,1] 1.0160345 4.410142e-04 3.303103e-02  0.9523405 0.99330338 1.0160469
betas[2,1] 0.1801277 4.616105e-04 3.381872e-02  0.1137230 0.15815957 0.1800260
sig        0.9424056 1.582922e-04 1.371087e-02  0.9157070 0.93301131 0.9421362
                 75%     97.5%    n_eff      Rhat
gssig[1]   0.2126310 0.2605958 5598.556 0.9998961
gssig[2]   0.2147029 0.2619235 6222.326 1.0002845
gsmat[1,1] 1.0000000 1.0000000      NaN       NaN
gsmat[1,2] 0.3197420 0.5525350 5699.475 0.9997332
gsmat[2,1] 0.3197420 0.5525350 5699.475 0.9997332
gsmat[2,2] 1.0000000 1.0000000 4823.755 0.9995999
betas[1,1] 1.0381456 1.0807423 5609.692 0.9997793
betas[2,1] 0.2023987 0.2469757 5367.389 1.0002954
sig        0.9518284 0.9694847 7502.583 1.0000681

summary(modD.brms)

 Family: gaussian 
  Links: mu = identity 
Formula: y ~ 1 + X1 + (X1 | g) 
   Data: dat (Number of observations: 2500) 
  Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 5000

Multilevel Hyperparameters:
~g (Number of levels: 50) 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)         0.19      0.03     0.14     0.26 1.00     2438     3470
sd(X1)                0.19      0.03     0.14     0.26 1.00     2182     3495
cor(Intercept,X1)     0.18      0.20    -0.23     0.55 1.00     1881     2949

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.02      0.03     0.95     1.08 1.00     3192     3678
X1            0.18      0.03     0.11     0.25 1.00     3150     3804

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.94      0.01     0.92     0.97 1.00    10049     3718

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).

summary(modE.custom, pars = pars)$summary

                   mean      se_mean           sd        2.5%         25%
gssig[1]    0.191584225 3.134064e-04 2.817121e-02  0.14165413  0.17182018
gssig[2]    0.205219533 3.129574e-04 2.940422e-02  0.15388628  0.18452991
gsmat[1,1]  1.000000000          NaN 0.000000e+00  1.00000000  1.00000000
gsmat[1,2]  0.247339082 2.018841e-03 1.812911e-01 -0.12111311  0.12388931
gsmat[2,1]  0.247339082 2.018841e-03 1.812911e-01 -0.12111311  0.12388931
gsmat[2,2]  1.000000000 1.257807e-18 9.165470e-17  1.00000000  1.00000000
betas[1,1]  1.014983726 3.711169e-04 3.258476e-02  0.95299170  0.99325273
betas[2,1]  0.186362834 3.857017e-04 3.374497e-02  0.11819030  0.16419798
betas[3,1]  0.403853970 1.803638e-04 1.727701e-02  0.36920740  0.39259458
betas[4,1] -0.001029122 1.916765e-04 1.729309e-02 -0.03535518 -0.01296458
sig         0.852348714 1.251941e-04 1.239611e-02  0.82807000  0.84424498
                     50%        75%      97.5%    n_eff      Rhat
gssig[1]    0.1898274520 0.20896368 0.25308192 8079.697 0.9997889
gssig[2]    0.2030220708 0.22351691 0.26848477 8827.727 0.9996936
gsmat[1,1]  1.0000000000 1.00000000 1.00000000      NaN       NaN
gsmat[1,2]  0.2555585836 0.37636408 0.58139844 8063.973 0.9997594
gsmat[2,1]  0.2555585836 0.37636408 0.58139844 8063.973 0.9997594
gsmat[2,2]  1.0000000000 1.00000000 1.00000000 5309.840 0.9995999
betas[1,1]  1.0145571484 1.03720499 1.08009665 7709.169 0.9996685
betas[2,1]  0.1866371130 0.20899630 0.25254896 7654.467 0.9999194
betas[3,1]  0.4038066291 0.41541954 0.43780493 9175.678 0.9996951
betas[4,1] -0.0009773671 0.01064513 0.03306578 8139.682 0.9999137
sig         0.8522041555 0.86052695 0.87693732 9803.996 0.9996905

summary(modE.brms)

 Family: gaussian 
  Links: mu = identity 
Formula: y ~ 1 + X1 + X2 + X3 + (X1 | g) 
   Data: dat (Number of observations: 2500) 
  Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 5000

Multilevel Hyperparameters:
~g (Number of levels: 50) 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)         0.19      0.03     0.14     0.25 1.00     2295     3311
sd(X1)                0.20      0.03     0.15     0.27 1.00     2432     3683
cor(Intercept,X1)     0.25      0.19    -0.13     0.60 1.00     1278     2353

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.02      0.03     0.95     1.08 1.00     2536     3397
X1            0.19      0.03     0.12     0.25 1.00     2521     3054
X2            0.40      0.02     0.37     0.44 1.00     8488     3605
X3           -0.00      0.02    -0.03     0.03 1.00    11458     3493

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.85      0.01     0.83     0.88 1.00     9250     3423

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).

Information Criteria

Here is where we have LOO-CV and WAIC for brms and the custom marginal likelihood code. Here, we would not expect these to match since they are using different likelihoods.

Note also that the typical strategy is to pick the model that has the lowest information criteria value (e.g., lowest LOO-CV or lowest WAIC; though you should pick which index ahead of time). Look at looic for LOO-CV or waic for WAIC. The results are also summarized in a Table after this output.

loo(modA.custom)


Computed from 5000 by 50 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3212.0 39.1
p_loo         6.3  0.8
looic      6424.1 78.2
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modA.brms)


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3185.3 36.5
p_loo        73.2  2.8
looic      6370.6 73.0
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.7]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modB.custom)


Computed from 5000 by 50 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3242.3 40.4
p_loo         7.1  0.9
looic      6484.7 80.8
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modB.brms)


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3228.1 37.2
p_loo        39.0  1.2
looic      6456.1 74.3
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 2.6]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modC.custom)


Computed from 5000 by 50 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3212.8 39.0
p_loo         6.8  0.9
looic      6425.5 78.1
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modC.brms)


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3185.8 36.5
p_loo        75.6  2.8
looic      6371.6 72.9
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.2]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modD.custom)


Computed from 5000 by 50 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3454.9 34.9
p_loo         5.4  0.7
looic      6909.8 69.8
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modD.brms)


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3433.6 35.8
p_loo        68.3  2.5
looic      6867.2 71.6
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.5]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modE.custom)


Computed from 5000 by 50 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3213.0 39.1
p_loo         7.2  0.8
looic      6426.0 78.1
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modE.brms)


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3186.0 36.5
p_loo        74.1  2.8
looic      6372.1 72.9
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

waic(extract_log_lik(modA.custom))

Warning: 
1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 50 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3212.0 39.1
p_waic         6.3  0.8
waic        6424.0 78.2

1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(modA.brms)

Warning: 
7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3185.0 36.5
p_waic        72.9  2.8
waic        6370.0 73.0

7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modB.custom))

Warning: 
2 (4.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 50 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3242.3 40.4
p_waic         7.1  0.9
waic        6484.7 80.8

2 (4.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(modB.brms)


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3228.0 37.2
p_waic        38.9  1.2
waic        6456.0 74.3

waic(extract_log_lik(modC.custom))

Warning: 
1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 50 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3212.7 39.0
p_waic         6.8  0.9
waic        6425.5 78.1

1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(modC.brms)

Warning: 
8 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3185.5 36.5
p_waic        75.3  2.8
waic        6371.0 72.9

8 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modD.custom))

Warning: 
2 (4.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 50 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3454.9 34.9
p_waic         5.3  0.7
waic        6909.7 69.8

2 (4.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(modD.brms)

Warning: 
5 (0.2%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3433.4 35.8
p_waic        68.1  2.5
waic        6866.7 71.6

5 (0.2%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modE.custom))

Warning: 
1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 50 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3212.9 39.1
p_waic         7.1  0.8
waic        6425.9 78.1

1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(modE.brms)

Warning: 
7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3185.8 36.5
p_waic        73.8  2.8
waic        6371.6 72.9

7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

Here, it turns out that… Model A would be chosen by both custom code and brms using LOO-CV. Same thing with WAIC, though we get some warnings. It also looks like LOO-CV and WAIC yield almost the same values for the custom code with marginal likelihood (these should be asymptotically equivalent when the marginal likelihood is used).

In the paper, we discuss why picking the lowest value may sometimes lead the researcher astray and some alternatives that can often increase model selection accuracy. Note, for example, how big the standard errors are. Differences in the information criteria among many of these models are negligible given how large these standard errors are.

Model	LOO_custom_marginal	se_custom_marginal	LOO_brms_conditional	se_brms_conditional
A	6424.07	78.19	6370.62	72.98
B	6484.69	80.77	6456.11	74.31
C	6425.54	78.08	6371.62	72.93
D	6909.76	69.84	6867.17	71.57
E	6425.98	78.11	6372.09	72.94

Model	waic_custom_marginal	se_custom_marginal	waic_brms_conditional	se_brms_conditional
A	6424.00	78.19	6370.02	72.96
B	6484.67	80.77	6456.02	74.30
C	6425.48	78.07	6371.04	72.91
D	6909.71	69.84	6866.70	71.56
E	6425.89	78.10	6371.57	72.92

Other Stuff

Caveats

WAIC and LOO-CV are not the only ways to select models or use multiple models
- Bayes factors could be used
- Bayesian model averaging could be used
This is intended as a basic example for others to build upon (respecting the GPL-3 license, of course; we appreciate some credit)
The code is not intended to work with complex set-ups for the random effects, nor parameter constraints; so, the “lmer-like” syntax can’t quite handle all lmer syntax that one could conceivably want to throw at it. Check the model/output before proceeding.

Other Code Surprises

The underlying code also contain things not mentioned in this post:

Utility functions used for simulating some of the data in the manuscript
Code that allows frequentist estimation using the marginal likelihood (as a sanity check)
Estimation of DIC from brms models or the custom stan models (though we don’t necessarily recommend you use this)
Code for conditional estimation of linear mixed effects models using stan
- This is illustrated above as a sanity check
- Under the hood, sometimes the way the random effects are sampled is not transparent with model equations; there are therefore two versions of the stan code for conditional estimation. The explanation can get quite technical, but…
- mlmcond_versionA.stan is how we expect most people would write the model and do sampling
- The other version, mlmcond.stan, is closer to what brms does (sampling each vector of random effects independently from standard normal; transformation to then obtain actual estimates of possibly correlated random effects)

Custom Conditional Likelihood

Custom Conditional Likelihood Code

To avoid clutter with the marginal likelihood part, I moved the code for custom conditional likelihood code here. If you like, compare the results here with what brms obtains for these same models. The information criteria results should line up nicely enough because both brms and the code here are using the conditional likelihood.

matsA <- mod.mat.mlm(modA, data=dat)
modA.custom.cond <- stan(file='mlmcond.stan', data = matsA,
                    chains = nchains, seed=5297,
                    iter = iter, refresh = 0, cores = ncores)

matsB <- mod.mat.mlm(modB, data=dat)
modB.custom.cond <- stan(file='mlmcond.stan', data = matsB,
                    chains = nchains, seed=5297,
                    iter = iter, refresh = 0, cores = ncores)

matsC <- mod.mat.mlm(modC, data=dat)
modC.custom.cond <- stan(file='mlmcond.stan', data = matsC,
                    chains = nchains, seed=5297,
                    iter = iter, refresh = 0, cores = ncores)

matsD <- mod.mat.mlm(modD, data=dat)
modD.custom.cond <- stan(file='mlmcond.stan', data = matsD,
                    chains = nchains, seed=5297,
                    iter = iter, refresh = 0, cores = ncores)

matsE <- mod.mat.mlm(modE, data=dat)
modE.custom.cond <- stan(file='mlmcond.stan', data = matsE,
                    chains = nchains, seed=5297,
                    iter = iter, refresh = 0, cores = ncores)

Look at Diagnostics

Rhats

# error sd, coefficients, random effect covariance matrix
pars <- c("sig", "betas", "Tau")
max(abs(summary(modA.custom.cond, pars=pars)$summary[,"Rhat"]))

[1] 1.002459

max(abs(summary(modB.custom.cond, pars=pars)$summary[,"Rhat"]))

[1] 1.000347

max(abs(summary(modC.custom.cond, pars=pars)$summary[,"Rhat"]))

[1] 1.002137

max(abs(summary(modD.custom.cond, pars=pars)$summary[,"Rhat"]))

[1] 1.00131

max(abs(summary(modE.custom.cond, pars=pars)$summary[,"Rhat"]))

[1] 1.000386

Traceplots

rstan::traceplot(modA.custom.cond, pars = pars)

rstan::traceplot(modB.custom.cond, pars = pars)

rstan::traceplot(modC.custom.cond, pars = pars)

rstan::traceplot(modD.custom.cond, pars = pars)

rstan::traceplot(modE.custom.cond, pars = pars)

Compare to `brms`

Parameter Estimates

Here I compare all 3 approaches (custom marginal likelihood, custom conditional likelihood, brms)

pars <- c("gssig", "gsmat", "betas","sig")

summary(modA.custom, pars = pars)$summary

                mean      se_mean           sd       2.5%       25%       50%
gssig[1]   0.1918110 3.080043e-04 2.810884e-02  0.1418738 0.1721204 0.1900512
gssig[2]   0.2049995 3.239594e-04 3.030804e-02  0.1513372 0.1834301 0.2033433
gsmat[1,1] 1.0000000          NaN 0.000000e+00  1.0000000 1.0000000 1.0000000
gsmat[1,2] 0.2484226 2.118150e-03 1.850621e-01 -0.1337348 0.1232743 0.2586771
gsmat[2,1] 0.2484226 2.118150e-03 1.850621e-01 -0.1337348 0.1232743 0.2586771
gsmat[2,2] 1.0000000 1.271269e-18 9.195014e-17  1.0000000 1.0000000 1.0000000
betas[1,1] 1.0151835 3.529104e-04 3.154864e-02  0.9532490 0.9946503 1.0144379
betas[2,1] 0.1870510 3.994190e-04 3.396616e-02  0.1195645 0.1645170 0.1869748
betas[3,1] 0.4039262 1.910961e-04 1.706209e-02  0.3702468 0.3922990 0.4040235
sig        0.8520069 1.518920e-04 1.209477e-02  0.8286582 0.8439104 0.8519799
                 75%     97.5%    n_eff      Rhat
gssig[1]   0.2097149 0.2537442 8328.608 1.0000714
gssig[2]   0.2239143 0.2705220 8752.549 0.9996923
gsmat[1,1] 1.0000000 1.0000000      NaN       NaN
gsmat[1,2] 0.3802265 0.5919233 7633.466 0.9998811
gsmat[2,1] 0.3802265 0.5919233 7633.466 0.9998811
gsmat[2,2] 1.0000000 1.0000000 5231.548 0.9995999
betas[1,1] 1.0359007 1.0780544 7991.574 0.9997140
betas[2,1] 0.2102328 0.2535111 7231.616 0.9998718
betas[3,1] 0.4157825 0.4371860 7971.881 0.9998050
sig        0.8599656 0.8762982 6340.527 0.9999245

summary(modA.custom.cond, pars = pars)$summary

                mean      se_mean           sd       2.5%       25%       50%
gssig[1]   0.1917614 5.708049e-04 2.783647e-02  0.1422070 0.1718994 0.1904965
gssig[2]   0.2050996 6.471360e-04 2.977677e-02  0.1523955 0.1843445 0.2032047
gsmat[1,1] 1.0000000          NaN 0.000000e+00  1.0000000 1.0000000 1.0000000
gsmat[1,2] 0.2615447 5.032553e-03 1.826167e-01 -0.1154052 0.1387897 0.2692926
gsmat[2,1] 0.2615447 5.032553e-03 1.826167e-01 -0.1154052 0.1387897 0.2692926
gsmat[2,2] 1.0000000 1.351782e-18 9.104739e-17  1.0000000 1.0000000 1.0000000
betas[1,1] 1.0159249 7.210667e-04 3.230296e-02  0.9541794 0.9943023 1.0152046
betas[2,1] 0.1873588 7.601275e-04 3.409776e-02  0.1210169 0.1643743 0.1869247
betas[3,1] 0.4036757 1.787176e-04 1.733937e-02  0.3679687 0.3922464 0.4036457
sig        0.8522053 1.397209e-04 1.240695e-02  0.8283015 0.8439205 0.8519354
                 75%     97.5%    n_eff      Rhat
gssig[1]   0.2092188 0.2523644 2378.225 0.9997012
gssig[2]   0.2241427 0.2688319 2117.211 1.0005606
gsmat[1,1] 1.0000000 1.0000000      NaN       NaN
gsmat[1,2] 0.3933291 0.5922740 1316.753 0.9998251
gsmat[2,1] 0.3933291 0.5922740 1316.753 0.9998251
gsmat[2,2] 1.0000000 1.0000000 4536.512 0.9995999
betas[1,1] 1.0377332 1.0800986 2006.937 1.0024589
betas[2,1] 0.2097332 0.2571362 2012.235 0.9996596
betas[3,1] 0.4153987 0.4374996 9413.083 0.9997361
sig        0.8603520 0.8765842 7885.099 0.9996684

summary(modA.brms)

 Family: gaussian 
  Links: mu = identity 
Formula: y ~ 1 + X1 + X2 + (X1 | g) 
   Data: dat (Number of observations: 2500) 
  Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 5000

Multilevel Hyperparameters:
~g (Number of levels: 50) 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)         0.19      0.03     0.14     0.25 1.00     2475     2646
sd(X1)                0.20      0.03     0.15     0.27 1.00     2285     2989
cor(Intercept,X1)     0.25      0.18    -0.13     0.59 1.00     1302     2198

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.01      0.03     0.95     1.08 1.00     2466     3561
X1            0.19      0.03     0.12     0.26 1.00     2808     3348
X2            0.40      0.02     0.37     0.44 1.00     9872     3183

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.85      0.01     0.83     0.88 1.00     8218     2922

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).

summary(modB.custom, pars = pars)$summary

                mean      se_mean         sd      2.5%       25%       50%
gssig[1]   0.1923681 0.0003690778 0.02786288 0.1422447 0.1735016 0.1907334
gsmat[1,1] 1.0000000          NaN 0.00000000 1.0000000 1.0000000 1.0000000
betas[1,1] 1.0130818 0.0004613252 0.03263846 0.9496920 0.9914550 1.0130321
betas[2,1] 0.1893812 0.0002292444 0.01786855 0.1541005 0.1772779 0.1890396
betas[3,1] 0.4010111 0.0002196821 0.01749074 0.3664365 0.3892013 0.4013143
sig        0.8736130 0.0001618714 0.01258730 0.8492467 0.8652348 0.8734514
                 75%     97.5%    n_eff      Rhat
gssig[1]   0.2090517 0.2531057 5699.232 1.0002389
gsmat[1,1] 1.0000000 1.0000000      NaN       NaN
betas[1,1] 1.0344937 1.0786197 5005.471 1.0000738
betas[2,1] 0.2014761 0.2239965 6075.491 0.9997895
betas[3,1] 0.4124569 0.4347344 6339.090 0.9998775
sig        0.8818100 0.8991470 6046.789 1.0004460

summary(modB.custom.cond, pars = pars)$summary

                mean      se_mean         sd      2.5%       25%       50%
gssig[1]   0.1922894 0.0006734092 0.02791011 0.1425325 0.1729303 0.1905776
gsmat[1,1] 1.0000000          NaN 0.00000000 1.0000000 1.0000000 1.0000000
betas[1,1] 1.0121623 0.0006811161 0.03305573 0.9467564 0.9902695 1.0121083
betas[2,1] 0.1894188 0.0002134748 0.01801588 0.1540873 0.1769402 0.1896585
betas[3,1] 0.4008553 0.0002147105 0.01782444 0.3659216 0.3888636 0.4007410
sig        0.8734694 0.0001510031 0.01244701 0.8500514 0.8649370 0.8732035
                 75%     97.5%    n_eff      Rhat
gssig[1]   0.2101924 0.2517036 1717.770 1.0000110
gsmat[1,1] 1.0000000 1.0000000      NaN       NaN
betas[1,1] 1.0341945 1.0764556 2355.327 0.9996959
betas[2,1] 0.2019651 0.2240854 7122.264 1.0003473
betas[3,1] 0.4126570 0.4363902 6891.678 0.9999336
sig        0.8817166 0.8992709 6794.511 1.0001442

summary(modB.brms)

 Family: gaussian 
  Links: mu = identity 
Formula: y ~ 1 + X1 + X2 + (1 | g) 
   Data: dat (Number of observations: 2500) 
  Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 5000

Multilevel Hyperparameters:
~g (Number of levels: 50) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.19      0.03     0.14     0.26 1.00     2081     3076

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.01      0.03     0.95     1.08 1.00     3002     3418
X1            0.19      0.02     0.15     0.22 1.00    10349     3228
X2            0.40      0.02     0.37     0.44 1.00     9334     3585

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.87      0.01     0.85     0.90 1.00     9869     3692

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).

summary(modC.custom, pars = pars)$summary

                  mean      se_mean           sd          2.5%          25%
gssig[1]    0.19160290 3.313331e-04 2.897480e-02  0.1401200106  0.171458762
gssig[2]    0.20565102 3.823812e-04 3.000967e-02  0.1515489798  0.184771660
gssig[3]    0.02332682 2.615383e-04 1.765683e-02  0.0009513602  0.009441122
gsmat[1,1]  1.00000000          NaN 0.000000e+00  1.0000000000  1.000000000
gsmat[1,2]  0.23936578 2.267083e-03 1.807565e-01 -0.1299656368  0.115940309
gsmat[1,3] -0.11973518 5.504373e-03 4.735594e-01 -0.9008676107 -0.504083418
gsmat[2,1]  0.23936578 2.267083e-03 1.807565e-01 -0.1299656368  0.115940309
gsmat[2,2]  1.00000000 1.301942e-18 9.161434e-17  1.0000000000  1.000000000
gsmat[2,3]  0.09839165 5.842529e-03 4.712396e-01 -0.8168947751 -0.254230038
gsmat[3,1] -0.11973518 5.504373e-03 4.735594e-01 -0.9008676107 -0.504083418
gsmat[3,2]  0.09839165 5.842529e-03 4.712396e-01 -0.8168947751 -0.254230038
gsmat[3,3]  1.00000000 2.074288e-18 8.544524e-17  1.0000000000  1.000000000
betas[1,1]  1.01449142 4.051494e-04 3.231456e-02  0.9515671688  0.992450492
betas[2,1]  0.18730244 4.294914e-04 3.493284e-02  0.1180532339  0.164203593
betas[3,1]  0.40381828 2.105561e-04 1.738987e-02  0.3692833294  0.392414084
sig         0.85194345 1.459069e-04 1.240023e-02  0.8278588918  0.843579485
                   50%        75%      97.5%    n_eff      Rhat
gssig[1]    0.18930240 0.21017484 0.25350858 7647.354 0.9998031
gssig[2]    0.20359782 0.22406288 0.27189200 6159.271 0.9999291
gssig[3]    0.01977891 0.03297286 0.06628364 4557.793 1.0000400
gsmat[1,1]  1.00000000 1.00000000 1.00000000      NaN       NaN
gsmat[1,2]  0.24480165 0.36771556 0.57529230 6357.017 0.9999002
gsmat[1,3] -0.15818707 0.23515061 0.82700745 7401.730 1.0001531
gsmat[2,1]  0.24480165 0.36771556 0.57529230 6357.017 0.9999002
gsmat[2,2]  1.00000000 1.00000000 1.00000000 4951.575 0.9995999
gsmat[2,3]  0.12676637 0.47321536 0.89376744 6505.516 1.0001128
gsmat[3,1] -0.15818707 0.23515061 0.82700745 7401.730 1.0001531
gsmat[3,2]  0.12676637 0.47321536 0.89376744 6505.516 1.0001128
gsmat[3,3]  1.00000000 1.00000000 1.00000000 1696.828 0.9995999
betas[1,1]  1.01456233 1.03652346 1.07865477 6361.597 0.9996710
betas[2,1]  0.18756873 0.21026382 0.25797154 6615.440 0.9999144
betas[3,1]  0.40387583 0.41533281 0.43792436 6821.142 1.0000721
sig         0.85170172 0.86023466 0.87635550 7222.842 1.0000742

summary(modC.custom.cond, pars = pars)$summary

                  mean      se_mean           sd         2.5%         25%
gssig[1]    0.19214914 6.730741e-04 2.907744e-02  0.142298262  0.17167740
gssig[2]    0.20445448 6.803969e-04 2.950764e-02  0.152877679  0.18347172
gssig[3]    0.02367172 3.489747e-04 1.785130e-02  0.000837383  0.00943531
gsmat[1,1]  1.00000000          NaN 0.000000e+00  1.000000000  1.00000000
gsmat[1,2]  0.24195410 5.044851e-03 1.807986e-01 -0.131923742  0.11723144
gsmat[1,3] -0.12034850 6.507067e-03 4.597356e-01 -0.878255710 -0.47990076
gsmat[2,1]  0.24195410 5.044851e-03 1.807986e-01 -0.131923742  0.11723144
gsmat[2,2]  1.00000000 1.269140e-18 9.108801e-17  1.000000000  1.00000000
gsmat[2,3]  0.10299380 6.994985e-03 4.729334e-01 -0.828721219 -0.24591166
gsmat[3,1] -0.12034850 6.507067e-03 4.597356e-01 -0.878255710 -0.47990076
gsmat[3,2]  0.10299380 6.994985e-03 4.729334e-01 -0.828721219 -0.24591166
gsmat[3,3]  1.00000000 8.972476e-19 6.065311e-17  1.000000000  1.00000000
betas[1,1]  1.01433997 8.295802e-04 3.254999e-02  0.951528582  0.99280894
betas[2,1]  0.18657261 8.449635e-04 3.392476e-02  0.119475336  0.16359012
betas[3,1]  0.40392262 2.178400e-04 1.792986e-02  0.368567611  0.39140375
sig         0.85210301 1.606109e-04 1.242186e-02  0.828115581  0.84370768
                   50%        75%      97.5%    n_eff      Rhat
gssig[1]    0.18950242 0.21051196 0.25497265 1866.322 1.0009214
gssig[2]    0.20192920 0.22321936 0.26696559 1880.807 1.0008805
gssig[3]    0.02027037 0.03382843 0.06503983 2616.687 1.0007433
gsmat[1,1]  1.00000000 1.00000000 1.00000000      NaN       NaN
gsmat[1,2]  0.25111754 0.37284176 0.56568546 1284.380 0.9999281
gsmat[1,3] -0.16190907 0.20802457 0.79833055 4991.669 1.0008697
gsmat[2,1]  0.25111754 0.37284176 0.56568546 1284.380 0.9999281
gsmat[2,2]  1.00000000 1.00000000 1.00000000 5151.138 0.9995999
gsmat[2,3]  0.13935178 0.47998879 0.88808065 4571.159 0.9997500
gsmat[3,1] -0.16190907 0.20802457 0.79833055 4991.669 1.0008697
gsmat[3,2]  0.13935178 0.47998879 0.88808065 4571.159 0.9997500
gsmat[3,3]  1.00000000 1.00000000 1.00000000 4569.635 0.9995999
betas[1,1]  1.01406604 1.03565631 1.07829590 1539.519 1.0006541
betas[2,1]  0.18631058 0.20921032 0.25298182 1611.972 1.0021369
betas[3,1]  0.40400965 0.41583230 0.43813461 6774.525 0.9997167
sig         0.85187168 0.86006417 0.87721235 5981.682 0.9999569

summary(modC.brms)

 Family: gaussian 
  Links: mu = identity 
Formula: y ~ 1 + X1 + X2 + (X1 + X2 | g) 
   Data: dat (Number of observations: 2500) 
  Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 5000

Multilevel Hyperparameters:
~g (Number of levels: 50) 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)         0.19      0.03     0.14     0.25 1.00     1827     3052
sd(X1)                0.20      0.03     0.15     0.27 1.00     2210     2870
sd(X2)                0.02      0.02     0.00     0.07 1.00     2775     2980
cor(Intercept,X1)     0.25      0.18    -0.12     0.57 1.00     1305     2687
cor(Intercept,X2)    -0.12      0.46    -0.89     0.81 1.00     5760     3624
cor(X1,X2)            0.11      0.47    -0.80     0.88 1.00     6051     3705

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.01      0.03     0.95     1.08 1.00     1596     2842
X1            0.19      0.03     0.12     0.26 1.00     1909     2731
X2            0.40      0.02     0.37     0.44 1.00     7832     3561

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.85      0.01     0.83     0.88 1.00     7283     3611

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).

summary(modD.custom, pars = pars)$summary

                mean      se_mean           sd       2.5%        25%       50%
gssig[1]   0.1935204 4.130046e-04 3.090245e-02  0.1392984 0.17193187 0.1916898
gssig[2]   0.1958051 3.898703e-04 3.075364e-02  0.1418536 0.17419381 0.1937933
gsmat[1,1] 1.0000000          NaN 0.000000e+00  1.0000000 1.00000000 1.0000000
gsmat[1,2] 0.1809829 2.630157e-03 1.985634e-01 -0.2254724 0.04764889 0.1875613
gsmat[2,1] 0.1809829 2.630157e-03 1.985634e-01 -0.2254724 0.04764889 0.1875613
gsmat[2,2] 1.0000000 1.313643e-18 9.123677e-17  1.0000000 1.00000000 1.0000000
betas[1,1] 1.0160345 4.410142e-04 3.303103e-02  0.9523405 0.99330338 1.0160469
betas[2,1] 0.1801277 4.616105e-04 3.381872e-02  0.1137230 0.15815957 0.1800260
sig        0.9424056 1.582922e-04 1.371087e-02  0.9157070 0.93301131 0.9421362
                 75%     97.5%    n_eff      Rhat
gssig[1]   0.2126310 0.2605958 5598.556 0.9998961
gssig[2]   0.2147029 0.2619235 6222.326 1.0002845
gsmat[1,1] 1.0000000 1.0000000      NaN       NaN
gsmat[1,2] 0.3197420 0.5525350 5699.475 0.9997332
gsmat[2,1] 0.3197420 0.5525350 5699.475 0.9997332
gsmat[2,2] 1.0000000 1.0000000 4823.755 0.9995999
betas[1,1] 1.0381456 1.0807423 5609.692 0.9997793
betas[2,1] 0.2023987 0.2469757 5367.389 1.0002954
sig        0.9518284 0.9694847 7502.583 1.0000681

summary(modD.custom.cond, pars = pars)$summary

                mean      se_mean           sd       2.5%        25%       50%
gssig[1]   0.1927912 6.183227e-04 3.006365e-02  0.1385658 0.17205747 0.1910379
gssig[2]   0.1959378 6.386475e-04 3.067985e-02  0.1403043 0.17471799 0.1944403
gsmat[1,1] 1.0000000          NaN 0.000000e+00  1.0000000 1.00000000 1.0000000
gsmat[1,2] 0.1929143 4.892967e-03 1.986138e-01 -0.2147432 0.05866578 0.1997407
gsmat[2,1] 0.1929143 4.892967e-03 1.986138e-01 -0.2147432 0.05866578 0.1997407
gsmat[2,2] 1.0000000 1.283064e-18 9.073542e-17  1.0000000 1.00000000 1.0000000
betas[1,1] 1.0155747 6.859032e-04 3.340000e-02  0.9483559 0.99357811 1.0151354
betas[2,1] 0.1809350 5.776368e-04 3.400896e-02  0.1147238 0.15816750 0.1811421
sig        0.9421652 1.459965e-04 1.393512e-02  0.9157080 0.93225269 0.9422138
                 75%     97.5%    n_eff      Rhat
gssig[1]   0.2122468 0.2557692 2364.030 1.0003723
gssig[2]   0.2155238 0.2606005 2307.725 1.0014114
gsmat[1,1] 1.0000000 1.0000000      NaN       NaN
gsmat[1,2] 0.3358380 0.5584623 1647.685 1.0001474
gsmat[2,1] 0.3358380 0.5584623 1647.685 1.0001474
gsmat[2,2] 1.0000000 1.0000000 5001.006 0.9995999
betas[1,1] 1.0379409 1.0813420 2371.196 1.0011947
betas[2,1] 0.2030142 0.2480804 3466.387 1.0004660
sig        0.9515949 0.9691680 9110.387 0.9997178

summary(modD.brms)

 Family: gaussian 
  Links: mu = identity 
Formula: y ~ 1 + X1 + (X1 | g) 
   Data: dat (Number of observations: 2500) 
  Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 5000

Multilevel Hyperparameters:
~g (Number of levels: 50) 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)         0.19      0.03     0.14     0.26 1.00     2438     3470
sd(X1)                0.19      0.03     0.14     0.26 1.00     2182     3495
cor(Intercept,X1)     0.18      0.20    -0.23     0.55 1.00     1881     2949

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.02      0.03     0.95     1.08 1.00     3192     3678
X1            0.18      0.03     0.11     0.25 1.00     3150     3804

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.94      0.01     0.92     0.97 1.00    10049     3718

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).

summary(modE.custom, pars = pars)$summary

                   mean      se_mean           sd        2.5%         25%
gssig[1]    0.191584225 3.134064e-04 2.817121e-02  0.14165413  0.17182018
gssig[2]    0.205219533 3.129574e-04 2.940422e-02  0.15388628  0.18452991
gsmat[1,1]  1.000000000          NaN 0.000000e+00  1.00000000  1.00000000
gsmat[1,2]  0.247339082 2.018841e-03 1.812911e-01 -0.12111311  0.12388931
gsmat[2,1]  0.247339082 2.018841e-03 1.812911e-01 -0.12111311  0.12388931
gsmat[2,2]  1.000000000 1.257807e-18 9.165470e-17  1.00000000  1.00000000
betas[1,1]  1.014983726 3.711169e-04 3.258476e-02  0.95299170  0.99325273
betas[2,1]  0.186362834 3.857017e-04 3.374497e-02  0.11819030  0.16419798
betas[3,1]  0.403853970 1.803638e-04 1.727701e-02  0.36920740  0.39259458
betas[4,1] -0.001029122 1.916765e-04 1.729309e-02 -0.03535518 -0.01296458
sig         0.852348714 1.251941e-04 1.239611e-02  0.82807000  0.84424498
                     50%        75%      97.5%    n_eff      Rhat
gssig[1]    0.1898274520 0.20896368 0.25308192 8079.697 0.9997889
gssig[2]    0.2030220708 0.22351691 0.26848477 8827.727 0.9996936
gsmat[1,1]  1.0000000000 1.00000000 1.00000000      NaN       NaN
gsmat[1,2]  0.2555585836 0.37636408 0.58139844 8063.973 0.9997594
gsmat[2,1]  0.2555585836 0.37636408 0.58139844 8063.973 0.9997594
gsmat[2,2]  1.0000000000 1.00000000 1.00000000 5309.840 0.9995999
betas[1,1]  1.0145571484 1.03720499 1.08009665 7709.169 0.9996685
betas[2,1]  0.1866371130 0.20899630 0.25254896 7654.467 0.9999194
betas[3,1]  0.4038066291 0.41541954 0.43780493 9175.678 0.9996951
betas[4,1] -0.0009773671 0.01064513 0.03306578 8139.682 0.9999137
sig         0.8522041555 0.86052695 0.87693732 9803.996 0.9996905

summary(modE.custom.cond, pars = pars)$summary

                    mean      se_mean           sd        2.5%         25%
gssig[1]    0.1922497909 5.878509e-04 2.853238e-02  0.14087655  0.17260126
gssig[2]    0.2054511746 6.728583e-04 3.006053e-02  0.15307039  0.18444785
gsmat[1,1]  1.0000000000          NaN 0.000000e+00  1.00000000  1.00000000
gsmat[1,2]  0.2520563312 4.591635e-03 1.860165e-01 -0.13136699  0.12775814
gsmat[2,1]  0.2520563312 4.591635e-03 1.860165e-01 -0.13136699  0.12775814
gsmat[2,2]  1.0000000000 1.381783e-18 9.104739e-17  1.00000000  1.00000000
betas[1,1]  1.0155099452 6.790466e-04 3.306841e-02  0.95057689  0.99320109
betas[2,1]  0.1882400793 6.519954e-04 3.399363e-02  0.12263950  0.16623551
betas[3,1]  0.4034706732 1.692261e-04 1.740901e-02  0.36905682  0.39228388
betas[4,1] -0.0007688253 1.689383e-04 1.691220e-02 -0.03453865 -0.01223937
sig         0.8522343439 1.209930e-04 1.221243e-02  0.82933382  0.84381928
                     50%        75%      97.5%     n_eff      Rhat
gssig[1]    0.1901421430 0.20988566 0.25326256  2355.817 1.0001158
gssig[2]    0.2031737811 0.22409647 0.26995233  1995.933 1.0004933
gsmat[1,1]  1.0000000000 1.00000000 1.00000000       NaN       NaN
gsmat[1,2]  0.2617509046 0.38430188 0.58476762  1641.226 0.9997707
gsmat[2,1]  0.2617509046 0.38430188 0.58476762  1641.226 0.9997707
gsmat[2,2]  1.0000000000 1.00000000 1.00000000  4341.656 0.9995999
betas[1,1]  1.0153087408 1.03722314 1.07973751  2371.523 1.0002524
betas[2,1]  0.1874117749 0.21065831 0.25522049  2718.354 1.0001173
betas[3,1]  0.4031768067 0.41533354 0.43701037 10583.120 0.9996807
betas[4,1] -0.0007379099 0.01080002 0.03116623 10021.763 0.9996639
sig         0.8520701542 0.86061600 0.87659987 10187.880 0.9996518

summary(modE.brms)

 Family: gaussian 
  Links: mu = identity 
Formula: y ~ 1 + X1 + X2 + X3 + (X1 | g) 
   Data: dat (Number of observations: 2500) 
  Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 5000

Multilevel Hyperparameters:
~g (Number of levels: 50) 
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)         0.19      0.03     0.14     0.25 1.00     2295     3311
sd(X1)                0.20      0.03     0.15     0.27 1.00     2432     3683
cor(Intercept,X1)     0.25      0.19    -0.13     0.60 1.00     1278     2353

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.02      0.03     0.95     1.08 1.00     2536     3397
X1            0.19      0.03     0.12     0.25 1.00     2521     3054
X2            0.40      0.02     0.37     0.44 1.00     8488     3605
X3           -0.00      0.02    -0.03     0.03 1.00    11458     3493

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.85      0.01     0.83     0.88 1.00     9250     3423

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).

Information Criteria

Here, the custom conditional likelihood code should yield indices close to brms. But, recall, these are not expected to be close to the custom marginal likelihood code.

loo(modA.custom) # custom code, marginal likelihood


Computed from 5000 by 50 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3212.0 39.1
p_loo         6.3  0.8
looic      6424.1 78.2
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modA.custom.cond) # custom code, conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3185.0 36.5
p_loo        73.2  2.8
looic      6370.0 73.0
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.5]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modA.brms) # conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3185.3 36.5
p_loo        73.2  2.8
looic      6370.6 73.0
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.7]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modB.custom) # custom code, marginal likelihood


Computed from 5000 by 50 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3242.3 40.4
p_loo         7.1  0.9
looic      6484.7 80.8
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modB.custom.cond) # custom code, conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3228.6 37.2
p_loo        39.2  1.2
looic      6457.2 74.3
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 2.0]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modB.brms) # conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3228.1 37.2
p_loo        39.0  1.2
looic      6456.1 74.3
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 2.6]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modC.custom) # custom code, marginal likelihood


Computed from 5000 by 50 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3212.8 39.0
p_loo         6.8  0.9
looic      6425.5 78.1
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modC.custom.cond) # custom code, conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3186.0 36.5
p_loo        75.6  2.9
looic      6372.0 73.0
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modC.brms) # conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3185.8 36.5
p_loo        75.6  2.8
looic      6371.6 72.9
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.2]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modD.custom) # custom code, marginal likelihood


Computed from 5000 by 50 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3454.9 34.9
p_loo         5.4  0.7
looic      6909.8 69.8
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modD.custom.cond) # custom code, conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3433.1 35.8
p_loo        68.0  2.5
looic      6866.2 71.6
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.6]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modD.brms) # conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3433.6 35.8
p_loo        68.3  2.5
looic      6867.2 71.6
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.5]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modE.custom) # custom code, marginal likelihood


Computed from 5000 by 50 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3213.0 39.1
p_loo         7.2  0.8
looic      6426.0 78.1
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modE.custom.cond) # custom code, conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3185.2 36.5
p_loo        73.4  2.8
looic      6370.3 72.9
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.4]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo(modE.brms) # conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

         Estimate   SE
elpd_loo  -3186.0 36.5
p_loo        74.1  2.8
looic      6372.1 72.9
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

waic(extract_log_lik(modA.custom)) # custom code, marginal likelihood

Warning: 
1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 50 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3212.0 39.1
p_waic         6.3  0.8
waic        6424.0 78.2

1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modA.custom.cond)) # custom code, conditional likelihood

Warning: 
7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3184.7 36.5
p_waic        72.9  2.8
waic        6369.5 72.9

7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(modA.brms) # conditional likelihood

Warning: 
7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3185.0 36.5
p_waic        72.9  2.8
waic        6370.0 73.0

7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modB.custom)) # custom code, marginal likelihood

Warning: 
2 (4.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 50 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3242.3 40.4
p_waic         7.1  0.9
waic        6484.7 80.8

2 (4.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modB.custom.cond)) # custom code, conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3228.5 37.2
p_waic        39.2  1.2
waic        6457.1 74.3

waic(modB.brms) # conditional likelihood


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3228.0 37.2
p_waic        38.9  1.2
waic        6456.0 74.3

waic(extract_log_lik(modC.custom)) # custom code, marginal likelihood

Warning: 
1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 50 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3212.7 39.0
p_waic         6.8  0.9
waic        6425.5 78.1

1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modC.custom.cond)) # custom code, conditional likelihood

Warning: 
7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3185.7 36.5
p_waic        75.4  2.8
waic        6371.4 72.9

7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(modC.brms) # conditional likelihood

Warning: 
8 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3185.5 36.5
p_waic        75.3  2.8
waic        6371.0 72.9

8 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modD.custom)) # custom code, marginal likelihood

Warning: 
2 (4.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 50 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3454.9 34.9
p_waic         5.3  0.7
waic        6909.7 69.8

2 (4.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modD.custom.cond)) # custom code, conditional likelihood

Warning: 
5 (0.2%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3432.9 35.8
p_waic        67.7  2.5
waic        6865.8 71.5

5 (0.2%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(modD.brms) # conditional likelihood

Warning: 
5 (0.2%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3433.4 35.8
p_waic        68.1  2.5
waic        6866.7 71.6

5 (0.2%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modE.custom)) # custom code, marginal likelihood

Warning: 
1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 50 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3212.9 39.1
p_waic         7.1  0.8
waic        6425.9 78.1

1 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(extract_log_lik(modE.custom.cond)) # custom code, conditional likelihood

Warning: 
7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3184.9 36.5
p_waic        73.1  2.7
waic        6369.8 72.9

7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic(modE.brms) # conditional likelihood

Warning: 
7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.


Computed from 5000 by 2500 log-likelihood matrix.

          Estimate   SE
elpd_waic  -3185.8 36.5
p_waic        73.8  2.8
waic        6371.6 72.9

7 (0.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

Model	LOO_custom_marginal	se_custom_marginal	LOO_brms_conditional	se_brms_conditional	LOO_custom_conditional	se_custom_conditional
A	6424.07	78.19	6370.62	72.98	6369.98	72.96
B	6484.69	80.77	6456.11	74.31	6457.16	74.30
C	6425.54	78.08	6371.62	72.93	6371.97	72.96
D	6909.76	69.84	6867.17	71.57	6866.22	71.55
E	6425.98	78.11	6372.09	72.94	6370.30	72.94

Model	waic_custom_marginal	se_custom_marginal	waic_brms_conditional	se_brms_conditional	WAIC_custom_conditional	se_custom_conditional
A	6424.00	78.19	6370.02	72.96	6369.47	72.95
B	6484.67	80.77	6456.02	74.30	6457.06	74.30
C	6425.48	78.07	6371.04	72.91	6371.38	72.93
D	6909.71	69.84	6866.70	71.56	6865.80	71.54
E	6425.89	78.10	6371.57	72.92	6369.76	72.92

Footnotes

Don’t confuse these with the marginal likelihood that is computed for Bayes factors; not the same thing.↩︎
You should be able to click on the file and then click a download button; those of you who know git can go ahead and clone the repo.↩︎

Overview

Background

Model selection using Bayesian information criteria

What should be used? 🤔

The problem

Purpose of this post

Setup

Prerequisites

Packages and code

Load data

Setup of models

Custom Marginal Likelihood

brms

Estimate Models

Custom Marginal Likelihood Code

Estimate Models

Look at Diagnostics

Rhats

Traceplots

Compare to brms

Parameter Estimates

Information Criteria

Other Stuff

Caveats

Other Code Surprises

Custom Conditional Likelihood

Custom Conditional Likelihood Code

Look at Diagnostics

Rhats

Traceplots

Compare to brms

Parameter Estimates

Information Criteria

Footnotes

`brms`

Compare to `brms`

Compare to `brms`