'data.frame': 504 obs. of 8 variables:
$ actor : int 1 1 1 1 1 1 1 1 1 1 ...
$ recipient : int NA NA NA NA NA NA NA NA NA NA ...
$ condition : int 0 0 0 0 0 0 0 0 0 0 ...
$ block : int 1 1 1 1 1 1 2 2 2 2 ...
$ trial : int 2 4 6 8 10 12 14 16 18 20 ...
$ prosoc_left : int 0 0 1 0 1 1 1 1 0 0 ...
$ chose_prosoc: int 1 0 0 1 1 1 0 0 1 1 ...
$ pulled_left : int 0 1 0 0 1 1 0 0 0 0 ...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
Introduction
\[ \begin{align} \color{orangered}{y_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)}\\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \beta_{1} \times x_{1i} + \beta_{2} \times x_{2i}} \\ \end{align} \]
Le modèle linéaire Gaussien qu’on a vu aux Cours n°03 et n°04 est caractérisé par un ensemble de postulats, entre autres choses :
- Les résidus sont distribués selon une loi Normale.
- La variance de cette distribution Normale est constante (postulat d’homogénéité de la variance).
- Les prédicteurs agissent sur la moyenne de cette distribution.
- La moyenne suit un modèle linéaire ou multi-linéaire.
Introduction
Cette modélisation (le choix d’une distribution Normale) induit plusieurs contraintes, par exemple :
- Les données observées sont définies sur un espace continu.
- Cet espace n’est pas borné.
Comment modéliser certaines données qui ne respectent pas ces contraintes ? Par exemple, la proportion de bonnes réponses à un test (bornée entre 0 et 1), un temps de réponse (restreint à des valeurs positives et souvent distribué de manière approximativement log-normale), un nombre d’avalanches…
Introduction
Nous avons déjà rencontré un modèle différent : le modèle Beta-Binomial (cf. Cours n°02).
\[ \begin{align} \color{orangered}{y_{i}} \ &\color{orangered}{\sim \mathrm{Binomial}(n, p)} \\ \color{steelblue}{p} \ &\color{steelblue}{\sim \mathrm{Beta}(a, b)} \\ \end{align} \]
- Les données observées sont binaires (e.g., 0 vs. 1) ou le résultat d’une somme d’observations binaires (e.g., une proportion).
- La probabilité de succès (obtenir 1) a priori se caractérise par une distribution Beta.
- La probabilité de succès (obtenir 1) ne dépend d’aucun prédicteur.
Introduction
Cette modélisation induit deux contraintes :
- Les données observées sont définies sur un espace discret.
- Cet espace est borné.
Comment pourrait-on ajouter des prédicteurs à ce modèle ?
Modèle linéaire généralisé
\[ \begin{align} \color{orangered}{y_{i}} \ &\color{orangered}{\sim \mathrm{Binomial}(n, p_{i})} \\ \color{black}{f(p_{i})} \ &\color{black}{= \alpha + \beta \times x_{i}} \\ \end{align} \]
Objectifs :
- Rendre compte de données discrètes (e.g., échec/succès) générées par un processus unique.
- Introduire des prédicteurs dans le modèle.
Deux changements par rapport au modèle Gaussien :
- L’utilisation d’une distribution de probabilité Binomiale.
- Le modèle linéaire ne sert plus à décrire directement un des paramètres de la distribution, mais une fonction de ce paramètre (on peut aussi considérer que le modèle Gaussien était formulé avec une fonction identité).
Fonction de lien
Les fonctions de lien ont pour tâche de mettre en correspondance l’espace d’un modèle linéaire (non borné) avec l’espace d’un paramètre potentiellement borné comme une probabilité, définie sur l’intervalle \([0, 1]\).
Fonction de lien
Les fonctions de lien ont pour tâche de mettre en correspondance l’espace d’un modèle linéaire (non borné) avec l’espace d’un paramètre potentiellement borné comme une probabilité, définie sur l’intervalle \([0, 1]\).
Régression logistique
La fonction logit du GLM binomial (on parle de “log-odds”) :
\[\text{logit}(p_{i}) = \log \left(\frac{p_{i}}{1 - p_{i}} \right)\]
La cote d’un évènement (“odds” en anglais) est le ratio entre la probabilité que l’évènement se produise et la probabilité qu’il ne se produise pas. Le logarithme de cette cote est prédit par un modèle linéaire.
\[ \log \left(\frac{p_{i}}{1 - p_{i}} \right) = \alpha + \beta \times x_{i} \]
Pour retrouver la probabilité d’un évènement, il faut utiliser la fonction de lien inverse, la fonction logistique (ou logit-inverse) :
\[ p_{i} = \frac{\exp(\alpha + \beta \times x_{i})}{1 + \exp(\alpha + \beta \times x_{i})} \]
Complications induites par la fonction de lien
Ces fonctions de lien posent des problèmes d’interprétation : Un changement d’une unité sur un prédicteur n’a plus un effet constant sur la probabilité mais la modifie plus ou moins en fonction de son éloignement à l’origine. Quand \(x = 0\), une augmentation d’une demi-unité (i.e., \(\Delta x = 0.5\)) se traduit par une augmentation de la probabilité de \(0.25\). Puis, chaque augmentation d’une demi-unité se traduit par une augmentation de \(p\) de plus en plus petite…
Complications induites par la fonction de lien
Deuxième complication : cette fonction de lien “force” chaque prédicteur à interagir avec lui même et à interagir avec TOUS les autres prédicteurs, même si les interactions ne sont pas explicites…
Dans un modèle Gaussien, le taux de changement de \(y\) en fonction de \(x\) est donné par \(\partial(\alpha + \beta x)~/~\partial x = \beta\) et ne dépend pas de \(x\) (i.e., \(\beta\) est constant).
Dans un GLM binomial (avec une fonction de lien logit), la probabilité d’un évènement est donnée par la fonction logistique :
\[p_{i} = \frac{\exp(\alpha + \beta \times x_{i})}{1 + \exp(\alpha + \beta \times x_{i})}\]
Et le taux de changement de \(p\) en fonction du prédicteur \(x\) est donné par :
\[ \frac{\partial p}{\partial x} = \frac{\beta}{2(1 + \cosh(\alpha + \beta \times x))} \]
On voit que la variation sur \(p\) due au prédicteur \(x\) est fonction du prédicteur \(x\), et dépend également de la valeur de \(\alpha\)… !
Exemple de régression logistique : La prosocialité chez le chimpanzé
Régression logistique
- pulled_left : 1 lorsque le chimpanzé pousse le levier gauche, 0 sinon.
- prosoc_left : 1 lorsque le levier gauche est associé à l’option prosociale, 0 sinon.
- condition : 1 lorsqu’un partenaire est présent, 0 sinon.
Régression logistique
Le problème
On cherche à savoir si la présence d’un singe partenaire incite le chimpanzé à appuyer sur le levier prosocial, c’est à dire l’option qui donne de la nourriture aux deux individus. Autrement dit, est-ce qu’il existe une interaction entre l’effet de la latéralité et l’effet de la présence d’un autre chimpanzé sur la probabilité d’actionner le levier gauche.
Les variables
Observations (
pulled_left) : Ce sont des variables de Bernoulli. Elles prennent comme valeur 0/1.Prédicteur (
prosoc_left) : Est-ce que les deux plats sont sur la gauche ou sur la droite ?Prédicteur (
condition) : Est-ce qu’un partenaire est présent ?
Régression logistique
\[ \begin{align} \color{orangered}{L_{i}} \ &\color{orangered}{\sim \mathrm{Binomial}(1, p_{i})} \\ \text{(équivalent à)} \quad \color{orangered}{L_{i}} \ &\color{orangered}{\sim \mathrm{Bernoulli}(p_{i})} \\ \color{black}{\text{logit}(p_{i})} \ &\color{black}{= \alpha} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(0, \omega)} \\ \end{align} \]
Modèle mathématique sans prédicteur. Comment choisir une valeur pour \(\omega\)… ?
Prior predictive check
On écrit le modèle précédent avec brms et on échantillonne à partir du prior afin de vérifier que les prédictions du modèle (sur la base du prior seul) correspondent à nos attentes.
Prior predictive check
# récupère les échantillons (sur la base) du prior
prior_draws(x = mod1.1) %>%
# applique la fonction de lien inverse
mutate(p = brms::inv_logit_scaled(Intercept) ) %>%
ggplot(aes(x = p) ) +
geom_density(fill = "steelblue", adjust = 0.1) +
labs(x = "Probabilité a priori de tirer le levier gauche", y = "Densité de probabilité")