Introduction à l’analyse bayésienne

Ce cours présente les concepts de base de l’approche bayésienne pour l’estimation de paramètres, en préparation à l’utilisation des modèles hiérarchiques bayésiens durant les deux prochains cours.

Nous commençons avec une présentation de la théorie de la probabilité conditionnelle et du théorème de Bayes, un sujet qui n’est pas particulier à l’inférence bayésienne, mais qui sera utile pour comprendre la logique de cette approche.

Probabilité conditionnelle

Inférence bayésienne

Exemple de régression bayésienne

Visualiser et vérifier l’ajustement

Supposons qu’un nouveau test vise à dépister une maladie en mesurant la concentration d’une certaine protéine dans le sang. Cette concentration est plus élevée en moyenne chez les personnes atteintes que chez celles non-atteintes, mais ne permet pas de départager parfaitement les deux groupes. Les études cliniques ont permis d’estimer les deux propriétés suivantes de ce test:

La sensibilité du test, soit la probabilité d’obtenir un résultat positif si la maladie est présente, est de 95%.

La spécificité du test, soit la probabilité d’obtenir un résultat négatif si la maladie est absente, est de 99%.

À partir de cette information, peut-on déterminer quelle est la probabilité d’être atteint de la maladie, si on reçoit un résultat positif?

Tel que nous verrons ci-dessous, cette probabilité dépend de la prévalence générale de la maladie dans la population testée. Supposons donc que 0.2% de la population testée soit atteinte.

Avec cette information, nous pouvons diviser la population testée en 4 groupes selon que la maladie soit présente (\(M_1\)) ou non (\(M_0\)) et selon que le test soit positif (\(T_+\)) ou négatif (\(T_-\)). Sur 10 000 personnes recevant le test, nous pouvons calculer qu’en moyenne:

0.2% x 10 000 = 20 seront atteintes.

Sur 20 personnes atteintes, en moyenne 95% x 20 = 19 recevront un résultat positif, et 1 recevra un résultat négatif.

Sur 9980 personnes non-atteintes, en moyenne 1% x 9980 = 99.8 (arrondi à 100) recevront un résultat positif et les autres 9880 recevront un résultat négatif.

Donc, sur les 119 personnes recevant en moyenne un résultat positif, 19 / 119 = environ 16% sont atteintes de la maladie. Parmi celles recevant un résultat négatif, 1 / 9881 = environ 0.01% sont atteintes.

Dans la section suivante, nous reprenons ce calcul en utilisant le concept de probabilité conditionnelle.

Si \(x\) et \(y\) sont deux variables aléatoires, la probabilité conditionnelle \(p(y|x)\) est la probabilité d’une valeur de \(y\) pour une valeur donnée de \(x\) (on dit aussi la probabilité de \(y\) sachant \(x\)).

Dans notre exemple, nous avons deux variables: \(M\) représente la présence ou absence de la maladie et \(T\) représente le résultat positif ou négatif du test. La sensibilité du test est la probabilité de \(T_+\) si on sait que la personne est atteinte, donc \(P(T_+ | M_1)\), égale à 0.95. La spécificité quant à elle représente la probabilité conditionnelle \(P(T_- | M_0)\), égale à 0.99.

La probabilité conjointe d’obtenir à la fois une certaine valeur de \(x\) et une certaine valeur de \(y\), notée \(p(x, y)\), peut être calculée de deux façons: (1) la probabilité d’obtenir \(x\) multipliée par la probabilité d’obtenir \(y\), sachant qu’on a obtenu \(x\); ou (2) la probabilité d’obtenir \(y\) multipliée par la probabilité d’obtenir \(x\), sachant qu’on a obtenu \(y\).

\[p(x, y) = p(x) p(y | x) = p(y) p(x | y)\]

Dans notre exemple, voici les deux façons de calculer la probabilité d’être atteint de la maladie et d’obtenir un test positif:

\[p(M_1, T_+) = p(M_1) p(T_+ | M_1) = p(T_+) p(M_1 | T_+)\]

La probabilité marginale d’une variable \(y\), \(p(y)\), est sa probabilité si on ignore la valeur des autres variables. Si on ne connaît pas directement \(p(y)\), mais qu’on connaît \(p(y, x)\) pour chaque valeur possible d’une autre variable \(x\), alors \(p(y)\) correspond à la somme des probabilités conjointes de \(x\) et \(y\) pour chaque valeur de \(x\). (La marginalisation représente l’action de faire cette somme des différentes possibilités pour les variables autres que \(y\).)

\[p(y) = \sum_x p(y, x) = \sum_x p(y|x) p(x)\]

Si \(x\) est une variable continue, le principe est le même, mais la somme devient une intégrale:

\[p(y) = \int p(y, x) \text{d}x = \int p(y|x) p(x) \text{d}x\]

En revenant à notre exemple, nous pouvons exprimer chaque élément du tableau comme une probabilité conjointe d’une valeur de \(M\) ou de \(T\) (cases intérieures), ou comme une probabilité marginale (dans le cas des sommes de chaque ligne et de chaque colonne).

Nous pouvons ainsi remplir les cases du tableau en utilisant la relation entre probabilités conditionnelles, conjointes et marginales. Par exemple, pour la rangée \(T_+\):

\(p(M_1, T_+) = p(M_1) p(T_+ | M_1) = 0.002 \times 0.95 = 0.0019\)

\(p(M_0, T_+) = p(M_0) p(T_+ | M_0) = 0.998 \times 0.01 = 0.00998\) (\(p(T_+ | M_0)\) est la probabilité d’obtenir un faux positif, donc le complément de la spécificité, égal à 1%.)

\(p(T_+) = p(M_1, T_+) + p(M_0, T_+) = 0.019 + 0.00998 \approx 0.119\)

Voici le tableau complet:

Pour deux variables \(x\) et \(y\), il peut être plus facile de calculer \(p(y|x)\) que \(p(x|y)\), ou vice versa. Dans notre exemple, nous connaissions la probabilité d’avoir un résultat positif si la maladie est présente \(p(T_+ | M_1)\), mais nous cherchions la probabilité que la maladie soit présente si le résultat est positif: \(p(M_1 | T_+)\). Le théorème de Bayes nous indique comment “inverser” la probabilité conditionnelle sans avoir à remplir un tableau comme celui ci-dessus.

Rappelons-nous qu’il y a deux façons de calculer la probabilité conjointe \(p(x, y)\) d’une valeur de \(x\) et d’une valeur de \(y\):

\[p(x, y) = p(x) p(y | x) = p(y) p(x | y)\]

Si on divise les deux parties à droite par \(p(y)\), on obtient le théorème de Bayes:

\[p(x|y) = \frac{p(x) p(y | x)}{p(y)}\]

Celui-ci nous indique qu’on peut calculer la distribution de probabilité de \(x\) conditionnelle à \(y\) si on connaît: (1) la distribution de probabilité de \(y\) conditionnelle à \(x\) et (2) la distribution de probabilité marginale de \(x\). Pour ce qui est du dénominateur \(p(y)\), celui-ci peut être obtenu en faisant la somme (ou l’intégrale) de \(p(x) p(y|x)\) sur l’ensemble des valeurs possibles de \(x\).

Dans notre exemple, l’application du théorème montre que \(p(M_1 | T_+)\) = 16%, tel que déterminé précédemment.

\[p(M_1 | T_+) = \frac{p(T_+ | M_1) p(M_1)}{p(T_+)}\]

–

\[p(M_1 | T_+) = \frac{p(T_+ | M_1) p(M_1)}{p(T_+ | M_1) p(M_1) + p(T_+ | M_0) p(M_0)}\]

\[p(M_1 | T_+) = \frac{0.95 \times 0.002}{0.95 \times 0.002 + 0.01 \times 0.998} = 0.16\]

Notez donc qu’avec un résultat positif, la probabilité d’être atteint est multipliée par 80 (16% vs. 0.2%), mais il demeure plus probable de ne pas être atteint. Autrement dit, si une maladie est assez rare, la plupart des résultats positifs seront des faux positifs. C’est pourquoi certains tests ne sont pas réalisés sans présence d’autres symptômes qui augmenteraient la probabilité d’être atteint avant le test. La question de réaliser un test ou non est une décision éthique parfois difficile, car une absence de détection ou une fausse alerte ont toutes deux des effets néfastes.

De la même façon, on pourrait calculer \(p(M_1 | T_-)\) = 0.01%. Donc avec un résultat négatif, la probabilité d’être atteint est divisée par 20 par rapport à la prévalence générale de la maladie.

En réarrangeant les termes du théorème de Bayes:

\[p(M_1 | T_+) = \frac{p(T_+ | M_1)}{p(T_+)} p(M_1)\]

on voit qu’il s’agit d’une méthode pour réviser une probabilité initiale \(p(M_1)\) en fonction d’une nouvelle information donnée par le résultat du test, pour obtenir une probabilité \(p(M_1 | T_+)\). Cette idée est à la base de l’inférence bayésienne.

Dans la section précédente, nous avons vu comment le théorème de Bayes permet de calculer la probabilité d’être atteint d’une maladie en fonction de sa prévalence générale \(p(M_1)\) et du résultat d’un test de dépistage, \(T_+\) ou \(T_-\):

\[p(M_1 | T_+) = \frac{p(T_+ | M_1)}{p(T_+)} p(M_1)\] \[p(M_1 | T_-) = \frac{p(T_- | M_1)}{p(T_-)} p(M_1)\]

Considérons maintenant \(M_1\) comme une hypothèse selon laquelle un patient donné est atteint de la maladie. La probabilité a priori de \(M_1\) (avant le test) est égale à \(p(M_1)\). Après le test, la probabilité a posteriori de \(M_1\) est \(p(M_1 | T_+)\) ou \(p(M_1 | T_-)\), selon le résultat.

Selon l’interprétation fréquentiste, les probabilités représent la fréquence d’événements après un grand nombre de répétitions d’une observation ou d’une expérience. Dans ce cas-ci, nous pouvons assigner une probabilité à \(M_1\) car le patient provient d’une population et la maladie a une certaine fréquence dans cette population.

L’interprétation fréquentiste est à la base de la plupart des cours d’introduction aux statistiques, car elle permet notamment de définir des tests d’hypothèse et des intervalles de confiance. Dans cette approche, on peut associer une probabilité aux statistiques basées sur les données, comme la moyenne d’un échantillon \(\bar{x}\), mais pas aux paramètres d’un modèle comme la moyenne de la population \(\mu\). Quand on définit un intervalle de confiance à 95% autour de \(\bar{x}\), ce n’est pas cet intervalle particulier qui a une probabilité de 95% de contenir \(\mu\) (après l’échantillonnage, l’intervalle et \(\mu\) sont tous les deux fixes), mais c’est 95% des échantillons possibles de \(x\) qui produiraient un intervalle contenant la valeur de \(\mu\).

Selon l’interprétation bayésienne, les probabilités représentent notre incertitude sur la valeur d’une quantité. On peut donc parler d’une distribution de probabilité même pour une valeur présumée fixe, ex.: un paramètre d’un modèle.

Historiquement, les débats entre les deux approches ont souvent été acrimonieux. Aujourd’hui, les mêmes statisticiens peuvent employer l’approche fréquentiste ou l’approche bayésienne selon la nature du problème. Cependant, il faut s’assurer de toujours interpréter les résultats en fonction de l’approche utilisée. Il faut se rappeler, par exemple, qu’un intervalle de confiance fréquentiste ne représente pas une distribution de probabilité du paramètre, ou qu’une vérification de l’ajustement d’un modèle bayésien n’est pas équivalente à un test d’hypothèse nulle.

Supposons que nous avons une série d’observations d’une variable \(y\), que nous représentons par un modèle incluant un paramètre ajustable \(\theta\). Dans l’approche bayésienne, nous assignons une distribution de probabilité a priori à \(\theta\), \(p(\theta)\), représentant l’incertitude sur la valeur du paramètre avant d’avoir observé les données. La probabilité des observations \(y\), conditionnelle à une valeur de \(\theta\) donnée, est donnée par la fonction de vraisemblance \(p(y|\theta)\).

À partir de cette information, nous pouvons utiliser le théorème de Bayes pour déduire \(p(\theta | y)\), soit la distribution a posteriori de \(\theta\) après avoir observé \(y\).

\[p(\theta | y) = \frac{p(y | \theta) p(\theta)}{p(y)}\]

Le dénominateur \(p(y)\) est obtenu en faisant la somme (ou l’intégrale) de \(p(y | \theta) p(\theta)\) pour l’ensemble des valeurs possibles de \(\theta\). La plupart du temps, cette quantité ne peut pas être calculée exactement, mais elle est approximée par les méthodes de Monte-Carlo que nous verrons au prochain cours.

Voici un exemple simple visant à illustrer le principe de l’inférence bayésienne.

Supposons que dix lancers d’une pièce de monnaie produisent la série de valeurs suivantes (0 = pile, 1 = face): 0,0,0,0,0,1,1,1,0,0. Nous cherchons à estimer \(p\), la probabilité d’obtenir “face” pour cette pièce.

Si \(y\) est le nombre de “face” obtenus sur \(n\) lancers, alors la fonction de vraisemblance est donnée par la distribution binomiale: \(y \sim \text{Bin}(n, p)\).

Les différents panneaux du graphique ci-dessous montrent la distribution a posteriori de \(p\) (ligne pleine) après chaque lancer dans la séquence. Dans ce cas, la distribution a priori indiquée par la ligne pointillée était très diffuse, accordant une probabilité égale à chaque valeur possible de \(p\). Après 10 lancers, le maximum de la distribution a posteriori est égal à la proportion de “face” dans les données (0.3), qui est aussi le maximum de vraisemblance.

Le graphique ci-dessous montre la même inférence, mais avec une distribution a priori beaucoup plus concentrée autour de 0.5. Cette distribution représente l’idée qu’une pièce a beaucoup plus de chance d’être équilibrée (50% face) et que les déviations autour de cette valeur sont généralement mineures. Dans ce cas, la distribution a posteriori se déplace en suivant les données, mais reste beaucoup plus près de celle a priori.

Finalement, le graphique ci-dessous compare la distribution a posteriori (ligne orange) en fonction de chaque distribution a priori (ligne pointillée). La vraisemblance (identique dans les deux cas) est indiquée par une ligne pleine et a été normalisée pour être comparée aux distributions.

Pour la distribution a priori diffuse (à gauche), la distribution a posteriori est exactement proportionnelle à la vraisemblance. Lorsque la distibution a priori est plus concentrée que la vraisemblance (à droite), la distribution a posteriori se situe entre les deux, mais plus près de la distribution a priori.

Dans certains cas, les connaissances antérieures peuvent nous donner une idée assez spécifique de la distribution a priori à utiliser. Par exemple, la distribution a posteriori obtenue par une étude peut servir de distribution a priori pour une étude subséquente.

Le plus souvent toutefois, nous pouvons utiliser une distribution assez diffuse pour pénaliser les valeurs très implausibles des paramètres, sans trop contraindre l’analyse. En anglais, le terme weakly informative prior décrit ce type de choix.

Une des critiques courantes de l’inférence bayésienne est que l’assignation d’une distribution a priori ajoute un biais à l’analyse. Cependant, si le choix de cette distribution est justifiée par le besoin de pénaliser des valeurs trop extrêmes des paramètres, le rôle de la distribution a priori n’est pas si différent de celui d’un effet aléatoire qui resserre les moyennes de groupes vers la moyenne générale, ou du paramètre de lissage dans un modèle additif qui pénalise les courbes trop complexes. Toutes ces méthodes sont des exemples de régularisation, c’est-à-dire l’imposition de contraintes permettant de contrôler le risque de surajustement dans un modèle complexe, sans avoir à fixer complètement certains paramètres et certains effets.

On peut considérer le choix de la distribution a priori comme une supposition qui s’ajoute aux autres suppositions du modèle, comme le choix de la distribution représentant les observations, le choix des prédicteurs et interactions à inclure ou non dans le modèle, etc. Ultimement, c’est le modèle au complet qui doit être validé en fonction de sa capacité à reproduire les caractéristiques des observations, incluant des observations autres que celles utilisées pour réaliser son ajustement.

Comme le maximum de vraisemblance, l’inférence bayésienne a l’avantage d’être applicable à n’importe quel type de modèle génératif, c’est-à-dire un modèle qui décrit mathématiquement comment les observations sont générées à partir des paramètres. Dans une approche bayésienne, les paramètres de divers types de modèles, incluant tous ceux vus dans ce cours (modèles linéaires généralisés, à effets mixtes, additifs, avec dépendance temporelle ou spatiale) peuvent être estimés avec les mêmes algorithmes; nous discuterons davantage de ces algorithmes au prochain cours. Ces algorithmes produisent la distribution a posteriori conjointe des paramètres du modèle, à partir de laquelle nous pouvons facilement obtenir la distribution de n’importe quelle quantité dérivée du modèle: combinaison des paramètres, prédiction, etc. À titre de comparaison, le maximum de vraisemblance ne nous indique que la valeur de chaque paramètre maximisant la vraisemblance et l’obtention d’intervalles de confiance sur des valeurs dérivées de plusieurs paramètres peut être très laborieuse. En contrepartie, les méthodes bayésiennes produisent plus d’informations, mais demandent beaucoup plus de ressources de calcul.

Du point de vue de la conception du modèle, l’approche bayésienne est aussi plus demandante, car il faut spécifier une distribution a priori pour chaque paramètre ajustable. Comme nous avons mentionné ci-dessus, ces distributions a priori sont néanmoins utiles pour stabiliser l’estimation de modèles complexes; lorsque les données sont limitées, contraindre la valeur des effets peut être une solution préférable à simplement éliminer ces effets du modèle.

Pour illustrer l’approche bayésienne, nous estimerons ici la relation entre le nombre d’espèces de plantes et la superficie de 30 îles de l’archipel des Galapagos, à partir du jeu de données galapagos.csv déjà vu plus tôt dans ce cours.

Dans ce cas, nous choisissons comme distribution a priori pour \(\beta_1\) une distribution exponentielle avec un paramètre de 4: \(\beta_1 \sim \text{Exp}(4)\). La distribution exponentielle a un maximum à 0 et décroît de façon exponentielle; le taux de cette décroissance est donnée par le paramètre ajustable. Ici, avec un paramètre de 4 nous avons environ 2% de probabilité que \(\beta_1 > 1\). Ainsi, les valeurs supérieures à 1 sont jugées improbables mais pas impossibles.

Quant à l’ordonnée à l’origine \(\beta_0\), il s’agit du logarithme du nombre d’espèces moyen lorsque \(\log A = 0\), autrement dit lorsque la superficie est de 1 km\(^2\).

Note: Le package brms que nous utiliserons transforme les prédicteurs d’une régression pour les centrer sur leur valeur moyenne. Cela n’affecte pas les résultats produits pour la régression, mais cela nécessite que nous spécifions une distribution a priori non pas pour \(\beta_0\) tel que défini ici, mais pour une ordonnée à l’origine qui représenterait la valeur moyenne de la réponse si les prédicteurs étaient à leur valeur moyenne. Autrement dit, nous cherchons une distribution a priori du log du nombre d’espèces pour une île ayant une log-superficie moyenne.

Supposons que les valeurs plausibles pour le nombre d’espèces de l’île moyenne sont entre 1 et 1000. En prenant le log de ces valeurs, nous obtenons \(\log(1) = 0\) et \(\log(1000) = 6.91\). Dans ce cas, nous choisissons une distribution a priori normale, avec moyenne de 3 et écart-type de 2: \(\beta_0 \sim \text{N}(3, 2)\). Celle-ci place environ 95% de la probabilité entre -1 et 7.

Pour voir si nos distributions a priori couvrent bien les scénarios plausibles pour notre modèle, il est utile de simuler les prédictions réalisées par ces distributions. Dans le code R ci-dessous, nous créons d’abord un jeu de données sim_df contenant 100 valeurs de \(\beta_0\) et \(\beta_1\) à partir des distributions a priori, avec un index i identifiant la simulation. Nous créons ensuite une grille associant à chaque numéro de simulation une série géométrique de valeurs de la superficie (entre 0.01 km\(^2\) et 10 000 km\(^2\)) pour réaliser les prédictions. Finalement, nous combinons les deux jeux de données et nous calculons la valeur \(\lambda\) prédite pour chaque superficie pour chacune des 100 simulations.