Chapitre 3 Mise en oeuvre pratique

3.1 Introduction

Dans ce chapitre, nous découvrirons brms un outil très pratique pour faire de la statistique bayésienne sans trop d’efforts. Dans la version enrichie en ligne à https://oliviergimenez.github.io/statistique-bayes/logiciels.html, vous trouverez aussi une introduction à NIMBLE. NIMBLE et brms sont deux packages R qui implémentent pour vous les algorithmes MCMC. Concrètement, il vous suffit de spécifier une vraisemblance et des priors pour que le théorème de Bayes s’applique automatiquement. Grâce à une syntaxe proche de celle de R, les deux packages rendent cette étape relativement simple, même pour des modèles complexes.

3.2 La syntaxe du package`brms`

brms signifie Bayesian Regression Models using Stan. Ce package permet de formuler et d’estimer des modèles de régression (voir la section suivante et les Chapitres 5 et 6) de manière intuitive grâce à une syntaxe proche de celle du package lme4 (la référence R pour les modèles mixtes), tout en s’appuyant sur Stan, un logiciel de référence en statistique bayésienne. Le package est en constant développement, rendez-vous à https://paul-buerkner.github.io/brms/. Vous pouvez obtenir de l’aide via https://discourse.mc-stan.org/.

Pour utiliser brms, on commence par préparer les données :

dat <- data.frame(y = 19, n = 57)

Sans oublier de charger brms :

library(brms)

La vraisemblance est binomiale dans notre exemple fil rouge. Dans brms, on peut exprimer cela simplement :

bayes.brms <- brm(
  y | trials(n) ~ 1, # le nombre de succès est fonction d'un intercept
  family = binomial("logit"), # famille binomiale avec fonction de lien logit
  data = dat, # données utilisées
  chains = 3, # nombre de chaînes MCMC
  iter = 2000, # nombre total d'itérations par chaîne
  warmup = 300, # nombre d'itérations burn-in
  thin = 1 # pas de sous-échantillonnage (chaque itération est conservée)
)

La syntaxe est relativement simple mais nécessite quelques explications. L’argument y | trials(n) ~ 1 permet de spécifier un modèle dans lequel on a $y$ succès parmi $n$ essais, et on estime une ordonnée à l’origine ou intercept seulement, le 1 après ~. Pourquoi un intercept ici ? Pourquoi pas directement la survie $\theta$ ? Parce qu’on utilise family = binomial("logit") à la ligne d’après pour spécifier à brms que la variable réponse suit une loi binomiale. Autrement dit on a un modèle linéaire généralisé (voir le Chapitre 6) avec $\text{logit}(\theta) = \beta$ et on estime $\beta$ l’intercept. Les arguments iter = 2000, warmup = 300 et chains = 3 stipulent à brms d’utiliser 300 itérations pour l’adaptation (burn-in), et les 1700 suivantes pour l’inférence, avec 3 chaînes.

Jetons un coup d’oeil aux résultats :

summary(bayes.brms)
#>  Family: binomial 
#>   Links: mu = logit 
#> Formula: y | trials(n) ~ 1 
#>    Data: dat (Number of observations: 1) 
#>   Draws: 3 chains, each with iter = 2000; warmup = 300; thin = 1;
#>          total post-warmup draws = 5100
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept    -0.69      0.28    -1.24    -0.15 1.00     1588     2307
#> 
#> 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).

Cette commande affiche un tableau récapitulatif des estimations postérieures pour chaque paramètre du modèle. On y trouve :

Estimate est la moyenne a posteriori.
Est.Error est l’écart-type de la distribution a posteriori.
l-95% CI et u-95% CI sont les bornes de l’intervalle de crédibilité à 95%.
Le diagnostic de convergence Rhat.
Bulk_ESS est la taille effective d’échantillon (Tail_ESS est une autre mesure de la taille effective d’échantillon qu’on n’utilisera pas ici).

La moyenne a posteriori vaut -0.69 bien loin de la proportion de ragondins qui ont survécu à l’hiver ($19/57 \approx 0.33$). Comme toujours dans R et l’implémentation des modèles linéaires généralisés (voir Chapitre 6), les estimations des paramètres sont données sur l’échelle de la fonction de lien. Ici l’intercept estimé est exprimé sur l’échelle logit. Pour le convertir en probabilité de survie (entre 0 et 1), on extrait d’abord les valeurs générées dans la distribution a posteriori de l’intercept $\beta$ avec la fonction brms::as_draws_matrix() :

draws_fit <- as_draws_matrix(bayes.brms)

Puis on applique la fonction logistique réciproque plogis() sur chacune de ces valeurs pour obtenir tout un tas de valeurs simulées dans la distribution a posteriori de la survie $\theta$ :

beta <- draws_fit[,'Intercept'] # sélectionne la colonne intercept
theta <- plogis(beta)  # conversion logit -> [0,1]

On obtient ainsi une estimation directe de la moyenne a posteriori de la probabilité de survie, accompagnée de son intervalle de crédibilité à 95% :

mean(theta)
#> [1] 0.3365172
quantile(theta, probas = c(2.5,97.5)/100)
#>        0%       25%       50%       75%      100% 
#> 0.1282870 0.2936775 0.3338043 0.3776052 0.5752788

Ou plus directement avec la fonction posterior::summarise_draws() :

summarise_draws(theta)
#> # A tibble: 1 × 10
#>   variable   mean median     sd    mad    q5   q95  rhat ess_bulk ess_tail
#>   <chr>     <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
#> 1 Intercept 0.337  0.334 0.0616 0.0621 0.240 0.442  1.00    1588.    2307.

3.3 Visualisation

Pour visualiser la distribution a posteriori de la probabilité de survie, il suffit d’utiliser (Figure 3.1) :

draws_fit %>%
  ggplot(aes(x = theta)) +
  geom_histogram(color = "white", fill = "steelblue", bins = 30) +
  labs(x = "Probabilité de survie", y = "Fréquence")

$Histogramme de la distribution a posteriori de la probabilité de survie ($\theta$).$

Figure 3.1: Histogramme de la distribution a posteriori de la probabilité de survie ($\theta$).

Dans brms, on peut évaluer la convergence des chaînes MCMC (Figure 3.2) :

plot(bayes.brms)

Histogramme de la distribution a posteriori et trace plot de la probabilité de survie sur l’échelle logit (b). Dans l’histogramme, l’axe des abscisses représente les valeurs possibles de l’intercept (échelle logit) et l’axe des ordonnées la fréquence des valeurs simulées. Dans le trace plot, l’axe des abscisses correspond au numéro de l’itération MCMC et l’axe des ordonnées aux valeurs simulées de l’intercept (échelle logit).

Figure 3.2: Histogramme de la distribution a posteriori et trace plot de la probabilité de survie sur l’échelle logit (b). Dans l’histogramme, l’axe des abscisses représente les valeurs possibles de l’intercept (échelle logit) et l’axe des ordonnées la fréquence des valeurs simulées. Dans le trace plot, l’axe des abscisses correspond au numéro de l’itération MCMC et l’axe des ordonnées aux valeurs simulées de l’intercept (échelle logit).

Ce graphique affiche les trace plots (à droite) ainsi que les densités a posteriori (à gauche).

Au passage, pour déterminer de la longueur de la période de pré-chauffage ou burn-in, il suffit de faire tourner brms avec warmup = 0 pour quelques centaines ou milliers d’itérations et d’examiner la trace du paramètre pour décider du nombre d’itérations à utiliser pour atteindre la convergence.

Un avantage majeur des méthodes MCMC est qu’elles permettent d’obtenir la distribution a posteriori de n’importe quelle fonction des paramètres en appliquant cette fonction aux valeurs tirées dans les distributions a posteriori de ces paramètres. A noter qu’ici on estime l’intercept $\beta$ et on a donc déjà utilisé cette idée pour obtenir la distribution a posteriori de la probabilité de survie en appliquant la fonction logit réciproque. Comme autre exemple, disons que j’aimerais calculer l’espérance de vie des ragondins, celle-ci étant donnée par $\lambda = -1/\log(\theta)$ :

beta <- draws_fit[,'Intercept'] # sélectionne la colonne intercept
theta <- plogis(beta)  # conversion logit -> [0,1]
lambda <- -1 / log(theta) # transforme survie en espérance de vie
summarize_draws(lambda) # résumé des tirages : moyenne, médiane, intervalles
#> # A tibble: 1 × 10
#>   variable   mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
#>   <chr>     <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
#> 1 Intercept 0.931  0.911 0.162 0.154 0.700  1.22  1.00    1588.    2307.

L’espérance de vie est d’un an approximativement. On peut également visualiser la distribution a posteriori de l’espérance de vie (Figure 3.3) :

lambda %>%
  as_tibble() %>%
  ggplot() +
  geom_histogram(aes(x = Intercept), color = "white") +
  labs(x = "Espérance de vie")

Histogramme de la distribution a posteriori de l'espérance de vie. L'axe des abscisses représente les différentes valeurs possibles de l'espérance de vie. L'axe vertical indique le nombre de tirages simulés (Count) pour chaque valeur.

Figure 3.3: Histogramme de la distribution a posteriori de l’espérance de vie. L’axe des abscisses représente les différentes valeurs possibles de l’espérance de vie. L’axe vertical indique le nombre de tirages simulés (Count) pour chaque valeur.

3.4 Les priors

Il y a tout un tas de paramètres qui sont fixés par défaut dans brms, il est important d’en être conscient. Ça concerne les priors en particulier. Dans brms, les priors par défaut sont souvent non informatifs ou faiblement informatifs, mais il est toujours bon de les examiner explicitement. La commande suivante permet d’afficher le résumé des priors utilisés dans un modèle déjà ajusté :

prior_summary(bayes.brms)
#> Intercept ~ student_t(3, 0, 2.5)

Le package brms utilise comme a priori faiblement informatif une loi de Student à 3 degrés de liberté, centrée en 0, avec un écart-type de 2.5. Les 3 degrés de liberté donnent une distribution avec des queues plus épaisses qu’une normale, ce qui donne une certaine robustesse aux valeurs extrêmes. Le centre en 0 traduit une absence d’a priori fort sur la valeur de l’intercept. La largeur 2.5 autorise une variation raisonnablement large de l’intercept sans être complètement non-informatif.

Dans certains cas, il est pertinent de définir soi-même un prior, par exemple pour refléter des connaissances issues de la littérature ou pour contraindre d’avantage l’estimation (prior informatif). Ici, on propose un prior normal centré en 0 avec un écart-type de 1.5 sur l’intercept, on en reparlera au Chapitre 4 :

nlprior <- prior(normal(0, 1.5), class = "Intercept")

On peut ensuite l’utiliser dans la spécification du modèle :

bayes.brms <- brm(y | trials(n) ~ 1,
                  family = binomial("logit"),
                  data = dat,
                  prior = nlprior, # nos propres priors
                  chains = 3,
                  iter = 2000,
                  warmup = 300,
                  thin = 1)

Vous pouvez vérifier que les résultats sont proches de ceux obtenus avec le prior par défaut :

summary(bayes.brms)
#>  Family: binomial 
#>   Links: mu = logit 
#> Formula: y | trials(n) ~ 1 
#>    Data: dat (Number of observations: 1) 
#>   Draws: 3 chains, each with iter = 2000; warmup = 300; thin = 1;
#>          total post-warmup draws = 5100
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept    -0.69      0.28    -1.28    -0.15 1.00     2053     2049
#> 
#> 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).

3.5 Application interactive associée au chapitre

Ce chapitre est complété par une application interactive (shinyApp) basée sur brms, permettant d’explorer l’ajustement d’un modèle bayésien simple, les sorties numériques et les diagnostics de convergence, sans écrire de code supplémentaire. L’application est fournie dans le script 03-brms.R et peut être lancée depuis R avec :

library(shiny)
runApp("03-brms.R")

3.6 En résumé

brms permet de tirer parti des méthodes MCMC sans avoir à écrire le modèle soi-même (la vraisemblance en particulier).
Sa syntaxe est simple et proche de celle de lme4, ce qui le rend particulièrement adapté pour les modèles linéaires généralisés (mixtes ou pas ; voir Chapitre 6).
En contrepartie, brms repose sur des composants préprogrammés (familles de modèles, etc.), et il est important d’être attentif aux choix faits par défaut, notamment en ce qui concerne les distributions a priori.
Ce chapitre propose ainsi une première approche concrète de l’implémentation des modèles bayésiens, avant d’aborder des modèles plus riches, comme les modèles mixtes.

2 Les méthodes MCMC

4 Les distributions a priori