Un cours en R et Stan avec brms
Ladislas Nalborczyk (LPC, LNC, CNRS, Aix-Marseille Univ)
Cours n°01 : Introduction à l’inférence bayésienne
Cours n°02 : Modèle Beta-Binomial
Cours n°03 : Introduction à brms, modèle de régression linéaire
Cours n°04 : Modèle de régression linéaire (suite)
Cours n°05 : Markov Chain Monte Carlo
Cours n°06 : Modèle linéaire généralisé
Cours n°07 : Comparaison de modèles
Cours n°08 : Modèles multi-niveaux
Cours n°09 : Modèles multi-niveaux généralisés
Cours n°10 : Data Hackathon
\[\newcommand\given[1][]{\:#1\vert\:}\]
Le paquet brms utilise la même syntaxe que les fonctions de base R (comme lm) ou que le paquet lme4.
La partie gauche représente notre variable dépendante (ou “outcome”, i.e., ce qu’on essaye de prédire). Le paquet brms permet également de fitter des modèles multivariés (plusieurs outcomes) en les combinant avec c():
Si l’on veut fitter un modèle sans intercept (why not), il faut le spécifier explicitement comme ci-dessous.
La première partie de la partie droite de la formule représente les effets constants (effets fixes), tandis que la seconde partie (entre parenthèses) représente les effets variables (effets aléatoires).
Lorsqu’on inclut plusieurs effets variables (e.g., un intercept et une pente variables), brms postule qu’on souhaite aussi estimer la corrélation entre ces deux effets. Dans le cas contraire, on peut supprimer cette corrélation (i.e., la fixer à 0) en utilisant ||.
Travailler avec des sujets humains implique un minimum de coopération réciproque. Mais ce n’est pas toujours le cas. Une partie non-négligeable des étudiants qui s’inscrivent pour passer des expériences de psychologie ne se présentent pas le jour prévu… On a voulu estimer la probabilité de présence d’un étudiant inscrit en fonction de l’envoi (ou non) d’un mail de rappel (cet exemple est présenté en détails dans deux blogposts, accessibles ici, et ici).
library(tidyverse)
library(imsb)
# import des données
data <- open_data(absence_multilevel) %>%
mutate(reminder = ifelse(test = reminder == 1, yes = 0.5, no = -0.5) )
# on affiche 12 lignes "au hasard" dans ces données
data %>% sample_frac() %>% head(12) reminder researcher presence absence total
1 -0.5 1 16 86 102
2 0.5 4 92 3 95
3 -0.5 10 23 58 81
4 0.5 6 82 6 88
5 -0.5 3 31 65 96
6 0.5 2 109 3 112
7 -0.5 6 34 54 88
8 0.5 7 64 16 80
9 -0.5 5 34 49 83
10 0.5 8 81 11 92
11 0.5 3 87 9 96
12 -0.5 4 61 34 95\[ \begin{aligned} \color{orangered}{y_{i}} \ &\color{orangered}{\sim \mathrm{Binomial}(n_{i}, p_{i})} \\ \color{black}{\text{logit}(p_{i})} \ &\color{black}{= \alpha_{\text{researcher}_{[i]}} + \beta_{\text{researcher}_{[i]}} \times \text{reminder}_{i}} \\ \color{steelblue}{\begin{bmatrix} \alpha_{\text{researcher}} \\ \beta_{\text{researcher}} \\ \end{bmatrix}} \ & \color{steelblue}{\sim \mathrm{MVNormal}\left(\begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \textbf{S}\right)} \\ \color{steelblue}{\textbf{S}} \ &\color{steelblue}{= \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \\ \end{pmatrix} \textbf{R} \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \\ \end{pmatrix}} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 1)} \\ \color{steelblue}{\beta} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 1)} \\ \color{steelblue}{(\sigma_{\alpha}, \sigma_{\beta})}\ &\color{steelblue}{\sim \mathrm{HalfCauchy}(0, 1)} \\ \color{steelblue}{\textbf{R}} \ &\color{steelblue}{\sim \mathrm{LKJcorr}(2)} \\ \end{aligned} \]
Il s’agit du même modèle de régression logistique vu au Cours n°06, avec une fonction de lien logit, mais cette fois-ci sur plusieurs niveaux.
mod1 <- brm(
formula = presence | trials(total) ~ 1 + reminder + (1 + reminder | researcher),
family = binomial(link = "logit"),
prior = prior1,
data = data,
sample_prior = TRUE,
warmup = 2000, iter = 10000,
chains = 4, cores = parallel::detectCores(),
control = list(adapt_delta = 0.95),
backend = "cmdstanr"
)