Introduction à la modélisation statistique bayésienne

Un cours en R, Stan, et 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

On considère un modèle de régression linéaire gaussien avec un prédicteur continu. Ce modèle contient trois paramètres à estimer : l’intercept \(\alpha\), la pente \(\beta\), et l’écart-type des “résidus” \(\sigma\).

\[ \begin{aligned} \color{orangered}{y_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \beta x_{i}} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(100, 10)} \\ \color{steelblue}{\beta} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 10)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{Exponential}(0.01)} \\ \end{aligned} \]

Rappels

Ce modèle s’implémente simplement via brms::brm().

library(brms)

priors <- c(
  prior(normal(100, 10), class = Intercept),
  prior(normal(0, 10), class = b),
  prior(exponential(0.01), class = sigma)
  )

model <- brm(
  formula = y ~ 1 + x,
  family = gaussian(),
  prior = priors,
  data = df
  )

Régression multiple

On va étendre le modèle précédent en ajoutant plusieurs prédicteurs, continus et/ou catégoriels. Pourquoi ?

  • “Contrôle” des facteurs de confusion (e.g., spurious correlations, simpson’s paradox). Un facteur de confusion est un facteur (une variable aléatoire) qui “perturbe” l’association entre deux variables d’intérêts.

  • Multiples causes : un phénomène peut émerger sous l’influence de multiples causes.

  • Interactions : l’influence d’un prédicteur sur la variable observée peut dépendre de la valeur d’un autre prédicteur.

Associations fortuites

library(tidyverse)
library(imsb)

df1 <- open_data(waffle) # import des données dans une dataframe
str(df1) # affiche la structure des données
'data.frame':   50 obs. of  13 variables:
 $ Location         : chr  "Alabama" "Alaska" "Arizona" "Arkansas" ...
 $ Loc              : chr  "AL" "AK" "AZ" "AR" ...
 $ Population       : num  4.78 0.71 6.33 2.92 37.25 ...
 $ MedianAgeMarriage: num  25.3 25.2 25.8 24.3 26.8 25.7 27.6 26.6 29.7 26.4 ...
 $ Marriage         : num  20.2 26 20.3 26.4 19.1 23.5 17.1 23.1 17.7 17 ...
 $ Marriage.SE      : num  1.27 2.93 0.98 1.7 0.39 1.24 1.06 2.89 2.53 0.58 ...
 $ Divorce          : num  12.7 12.5 10.8 13.5 8 11.6 6.7 8.9 6.3 8.5 ...
 $ Divorce.SE       : num  0.79 2.05 0.74 1.22 0.24 0.94 0.77 1.39 1.89 0.32 ...
 $ WaffleHouses     : int  128 0 18 41 0 11 0 3 0 133 ...
 $ South            : int  1 0 0 1 0 0 0 0 0 1 ...
 $ Slaves1860       : int  435080 0 0 111115 0 0 0 1798 0 61745 ...
 $ Population1860   : int  964201 0 0 435450 379994 34277 460147 112216 75080 140424 ...
 $ PropSlaves1860   : num  0.45 0 0 0.26 0 0 0 0.016 0 0.44 ...

Associations fortuites

On observe un lien positif entre le nombre de “waffle houses” et le taux de divorce…

df1 %>%
  ggplot(aes(x = WaffleHouses, y = Divorce) ) +
  geom_text(aes(label = Loc) ) +
  geom_smooth(method = "lm", color = "black", se = TRUE) +
  labs(x = "Nombre de 'Waffle Houses' (par million d'habitants)", y = "Taux de divorce")

Associations fortuites

On laisse de côté les Waffle Houses. On observe un lien positif entre le taux de mariage et le taux de divorce… mais est-ce qu’on peut vraiment dire que le mariage “cause” le divorce ?

df1$Divorce.s <- scale(x = df1$Divorce, center = TRUE, scale = TRUE)
df1$Marriage.s <- scale(x = df1$Marriage, center = TRUE, scale = TRUE)

df1 %>%
  ggplot(aes(x = Marriage, y = Divorce) ) +
  geom_point(pch = 21, color = "white", fill = "black", size = 5, alpha = 0.8) +
  geom_smooth(method = "lm", color = "black", se = TRUE) +
  labs(x = "Taux de mariage", y = "Taux de divorce")

Associations fortuites

On observe l’association inverse entre le taux de divorce et l’âge médian de mariage.

df1$MedianAgeMarriage.s <- scale(x = df1$MedianAgeMarriage, center = TRUE, scale = TRUE)

df1 %>%
  ggplot(aes(x = MedianAgeMarriage, y = Divorce) ) +
  geom_point(pch = 21, color = "white", fill = "black", size = 5, alpha = 0.8) +
  geom_smooth(method = "lm", color = "black", se = TRUE) +
  labs(x = "Age médian de mariage", y = "Taux de divorce")

Influence du taux de mariage sur le taux de divorce

\[ \begin{aligned} \color{orangered}{D_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \beta_{R} R_{i}} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 10)} \\ \color{steelblue}{\beta_{R}} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 1)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{Exponential}(1)} \\ \end{aligned} \]

priors <- c(
  prior(normal(0, 10), class = Intercept),
  prior(normal(0, 1), class = b),
  prior(exponential(1), class = sigma)
  )

mod1 <- brm(
  formula = Divorce.s ~ 1 + Marriage.s,
  family = gaussian(),
  prior = priors,
  # for prior predictive checking
  sample_prior = TRUE,
  data = df1
  )

Représentation visuelle des priors

Prédictions a priori (pour \(\mu\))

# getting the samples from the prior distribution
prior <- prior_samples(mod1)

# displaying the first six samples
head(prior)
   Intercept           b     sigma
1 -20.080658 -0.46479347 1.6214921
2 -20.123171 -0.46424381 0.3944802
3   5.961326 -0.88470400 0.1494592
4   6.894092 -0.76043929 0.7728310
5  10.770446 -0.28969324 1.1654583
6  -6.228275 -0.03852849 1.2831928

Prédictions a priori (pour \(\mu\))

prior %>% 
  sample_n(size = 1e2) %>% 
  rownames_to_column("draw") %>% 
  expand(nesting(draw, Intercept, b), a = c(-2, 2) ) %>%
  mutate(d = Intercept + b * a) %>% 
  ggplot(aes(x = a, y = d) ) +
  geom_line(aes(group = draw), color = "steelblue", size = 0.5, alpha = 0.5) +
  labs(x = "Taux de mariage (standardisé)", y = "Taux de divorce (standardisé)")

Prédictions a priori (pour \(D_{i}\))

pp_check(object = mod1, prefix = "ppd", ndraws = 1e2) + labs(x = "Taux de divorce", y = "Densité")