I was teaching varying-intercept, varying-slope multilevel models, and . . . I can get them to fit in Stan, but the code is kinda ugly, so I was struggling to clean it up, with no success. This will be a real research project, to add appropriate functions and possibly expand the Stan language so that these models can be written at a higher, more intuitive level.
Varying-intercept models aren’t so bad. In lme4 or blme or rstanarm or brms, you write something like:
y ~ 1 | group + x + z + x:z
and that transfers pretty directly into Stan. Just create the X matrix and go from there. Indeed, you can add as many batches of varying coefficients and it’s no biggie to code it up.
But once you get to varying intercepts and slopes, it all changes. In lme4 or blme or rstanarm or brms, you can just write things like:
y ~ (1 + z | group) + x + z + x:z
But if you want to program this directly in Stan, once you have varying intercepts and slopes, you have to deal with covariance-matrix decompositions and arrays of coefficient vectors, and it’s all a hairy mess.
What to do? For this semester’s class, Imma just gonna go with lme4/blme/rstanarm when fitting varying-intercept, varying-slope models. All this Stan coding is a rabbit hole that’s getting us away from the goal, which is to be able to fit, use, and evaluate statistical models for measurement and variation.
I would like to be able to more easily fit these in Stan, though. Why, you might ask? If we can fit them in lme4, or blme for more stability, or rstanarm for including more uncertainty in the inferences, then why bother coding directly in Stan?
The answer for why we want to code directly in Stan is that we’re often wanting to expand our models, for example adding mixture components, measurement error terms, time series or spatial dependence, etc.
For that reason, you will want to be able to code varying-intercept, varying-slope models in Stan—even if I won’t be teaching that in class this semester.
The good news is that I did some googling and found this tutorial by Will Hipson on programming hierarchical regressions in Stan. It’s from 2020 and I have not looked at every line of code there, but it all looks reasonable and there’s lots of explanation of the workflow. So maybe this is the best place to start, if you want to go in this direction, as you should!
I keep trying to learn Stan but always find I can do what I need with brms.
Agreed. I often think something along the lines of ‘if I was doing analysis X I would need to learn how to program that in Stan’, but in practise it doesn’t happen.
Also, Aki pointed out these recent features that allow you to build extensions of a model within brms:
Define Custom Response Distributions with brms
Estimating Non-Linear Models with brms
Hi Andrew, I hope you don’t mind that I point people to some relevant work we are doing. Hierarchical models are bread and butter stuff for psycholinguists, so we are trying hard to make Stan/brms mainstream through various means. Teaching this stuff feels like the most important work I am doing right now, more important even than the scientific side of things.
We have chapters on hierarchical modeling in our book (to be published soon with CRC Press), we use both brms and Stan:
https://vasishth.github.io/bayescogsci/book/
The online version will remain available for free. Comments/corrections are welcome; one can open issues: https://github.com/vasishth/bayescogsci/issues
This summer, I am teaching an intro to Bayes using brms/Stan, with a focus on hierarchical modeling, especially directed at researchers in linguistics who do experimental work. This is in Gent, Belgium:
https://www.mils.ugent.be/courses/module-9-bayesian-data-analysis/
Plus, at Potsdam, for the last seven years I have been running an annual summer school on stats for linguistics and psych, where our focus is on hierarchical modeling using Stan/brms:
https://vasishth.github.io/smlp2024/
Here, we teach both frequentist and Bayesian approaches to hierarchical modeling. This year, Doug Bates will be teaching his Julia version of lme4 with Phillip Alday and Reinhold Kliegl.
Does anyone know of a *comprehensive* overview of Rs (fitting) expression language? I’ve see many, but they all seem
to end before such expressions as (1+ z|group) are explained, and operator precedence is never mentioned.
An R -> stan translation would be great!
I don’t know of a comprehensive guide for R, but here’s one for Julia’s Formula.jl:
https://juliastats.org/StatsModels.jl/stable/formula/
And the extensions for MixedModels.jl: https://juliastats.org/MixedModels.jl/stable/constructors/#Model-constructors
In addition to using the @formula macro, you can construct the terms programmatically as they’re represented by data types like ConstantTerm, Term, and FunctionTerm etc.
Julia as a programming language is far far far superior to R, though it builds on computer science / programming language historical ideas which people might not know if they only learned R or Python, so there’s some stuff to be learned to take advantage of it all.
The vast majority of what you’d do in R is doable already in convenient form in Julia with DataFrames.jl, DataFramesMeta.jl, CSV.jl, GLM.jl, MixedModels.jl and you have access to both Turing.jl and Stan.jl for Bayesian models. There are database packages for MySQL, Postgres, SQLite, DuckDB, and others. Julia has the truly most powerful differential equations solvers of any language out there today.
I’ve heard good things about Arviz.jl but haven’t had a chance to use it yet.
I gave up on R in about 2019 and never looked back. I do all my analyses in Julia and if I had to call some special purpose pre-made R package to do special analysis I’d just use RCall.jl
One of the major advantages is that Julia doesn’t have Fexpr / “nonstandard evaluation”. It only has explicit macros. The advantage is that it makes Julia way way easier to reason about (in general you can’t rely on R to mean anything, nonstandard evaluation could always apply without you knowing about it). The disadvantage for people who are pure users of libraries and don’t code much themselves is that if you rely on some big convenience packages from Hadley you might need to learn some more explicit programming constructs.
I also personally like the community at discourse.julialang.org a lot.