Introduction à la modélisation statistique bayésienne

Un cours en R et Stan avec brms

Ladislas Nalborczyk (LPC, LNC, CNRS, Aix-Marseille Univ)

Planning

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\:}\]

Rappels de syntaxe

Le paquet brms utilise la même syntaxe que les fonctions de base R (comme lm) ou que le paquet lme4.

Reaction ~ Days + (1 + Days | Subject)

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():

c(Reaction, Memory) ~ Days + (1 + Days | Subject)

La partie droite permet de définir les prédicteurs. L’intercept est généralement implicite, de sorte que les deux écritures ci-dessous sont équivalentes.

c(Reaction, Memory) ~ Days + (1 + Days | Subject)
c(Reaction, Memory) ~ 1 + Days + (1 + Days | Subject)

Rappels de syntaxe

Si l’on veut fitter un modèle sans intercept (why not), il faut le spécifier explicitement comme ci-dessous.

c(Reaction, Memory) ~ 0 + Days + (1 + Days | Subject)

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

c(Reaction, Memory) ~ Days + (1 | Subject)
c(Reaction, Memory) ~ Days + (Days | Subject)

Le premier modèle ci-dessus contient seulement un intercept variable, qui varie par Subject. Le deuxième modèle contient également un intercept variable, mais aussi une pente variable pour l’effet de Days.

Rappels de syntaxe

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

c(Reaction, Memory) ~ Days + (1 + Days || Subject)

Les modèles précédents postulaient une fonction de vraisemblance Gaussienne. Ce postulat peut être changé facilement en spécifiant la fonction de vraisemblance souhaitée via l’argument family.

brm(Reaction ~ 1 + Days + (1 + Days | Subject), family = lognormal() )

Mise en pratique - absentéisme expérimental

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

Mise en pratique - absentéisme expérimental

\[ \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.

Mise en pratique - absentéisme expérimental

prior1 <- c(
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 1), class = b),
    prior(cauchy(0, 1), class = sd),
    prior(lkj(2), class = cor)
    )
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"
    )

Mise en pratique - absentéisme expérimental

mod1 %>%
    plot(
        combo = c("dens_overlay", "trace"), pars = c("^b_", "^cor_"), widths = c(1, 1.5),
        theme = theme_bw(base_size = 16, base_family = "Open Sans")
        )