Un cours en R, Stan, et 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\:}\]
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} \]
Ce modèle s’implémente simplement via brms::brm()
.
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.
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 ...
On observe un lien positif entre le nombre de “waffle houses” et le taux de divorce…
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")
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")
\[ \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} \]
# 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
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é)")