Introduction à la modélisation statistique bayésienne

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

Principes de l’analyse bayésienne :

• On dispose d’un ensemble de données à analyser
• On suppose un modèle génératif défini par un ensemble de paramètres
• On dispose d’une connaissance a priori quant à la valeur de ces paramètres

Rappels

Inférence bayésienne : On infère (ou plutôt on déduit) la probabilité que le paramètre ait telle ou telle valeur sachant les données (et le prior) via le théorème de Bayes.

\[ \color{purple}{p(\theta \given y)} = \frac{\color{orangered}{p(y \given \theta)}\ \color{steelblue}{p(\theta)}}{\color{green}{p(y)}} \propto \color{orangered}{p(y \given \theta)}\ \color{steelblue}{p(\theta)} \]

Rappels

Les étapes de l’analyse de données bayésienne :

1. Définir le modèle - Identifier les paramètres du modèle génératif, définir une distribution a priori pour ces paramètres.

Le modèle Beta-Binomial

Pourquoi ce modèle ?

• Le modèle Beta-Binomial couvre de nombreux problèmes de la vie courante :

  • Réussite / échec à un test
  • Présence / absence d’effets secondaires lors du test d’un médicament

• C’est un modèle simple

  • Un seul paramètre
  • Solution analytique

Loi de Bernoulli

S’applique à toutes les situations où le processus de génération des données ne peut résulter qu’en deux issues mutuellement exclusives (e.g., un lancer de pièce). À chaque essai, si on admet que \(\Pr(\text{face}) = \theta\), alors \(\Pr(\text{pile}) = 1 - \theta\).

Schéma de Bernoulli

Si l’on dispose d’une suite de lancers \(\{Y_i\}\) indépendants et identiquement distribués (i.e., chaque lancer a une distribution de Bernoulli de probabilité \(\theta\)), l’ensemble de ces lancers peut être décrit par une distribution binomiale.

Schéma de Bernoulli

Sachant que les essais sont indépendants les uns des autres, la probabilité d’obtenir cette séquence est de \((1 - \theta) \times (1 - \theta) \times (1 - \theta) \times \theta \times \theta\), c’est à dire : \(\theta^{2} (1 - \theta)^{3}\).

Coefficient binomial

Le coefficient binomial nous permet de calculer le nombre de combinaisons possibles résultant en \(y\) succès pour \(n\) lancers de la manière suivante (lu “\(y\) parmi \(n\)” ou “nombre de combinaisons de \(y\) parmi \(n\)”)¹ :

\[ \left(\begin{array}{l} n \\ y \end{array}\right) = C_{y}^{n} = \frac{n !}{y !(n - y) !} \]

# computing the total number of possible configurations in R
choose(n = 3, k = 1)

[1] 3

Loi binomiale

\[ p(y \given \theta) = \Pr(Y = y \given \theta) = \left(\begin{array}{l} n \\ y \end{array}\right) \theta^{y}(1 - \theta)^{n - y} \]

La loi binomiale nous permet de calculer la probabilité d’obtenir \(y\) succès sur \(n\) essais, pour un \(\theta\) donné. Exemple de la distribution binomiale pour une pièce non biaisée (\(\theta = 0.5\)), indiquant la probabilité d’obtenir \(n\) faces sur 10 lancers (en R: dbinom(x = 0:10, size = 10, prob = 0.5)).

Générer des données à partir d’une distribution binomiale

library(tidyverse)
set.seed(666) # for reproducibility

sample(x = c(0, 1), size = 500, prob = c(0.4, 0.6), replace = TRUE) %>% # theta = 0.6
        data.frame() %>%
        mutate(x = seq_along(.), y = cummean(.) ) %>%
        ggplot(aes(x = x, y = y) ) +
        geom_line(lwd = 1) +
        geom_hline(yintercept = 0.6, lty = 3) +
        labs(x = "Nombre de lancers", y = "Proportion de faces") +
        ylim(0, 1)

Définition du modèle (likelihood)

Fonction de vraisemblance (likelihood)

• Nous considérons \(y\) comme étant le nombre de succès
• Nous considérons le nombre d’observations \(n\) comme étant une constante
• Nous considérons \(\theta\) comme étant le paramètre de notre modèle (i.e., la probabilité de succès)

Vraisemblance versus probabilité

On lance à nouveau une pièce de biais \(\theta\) (où \(\theta\) représente la probabilité d’obtenir Face). On lance cette pièce deux fois et on obtient une Face et un Pile.

Vraisemblance versus probabilité

# Représentation graphique de la fonction de vraisemblance de theta pour y = 1 et n = 2

y <- 1 # nombre de faces
n <- 2 # nombre d'essais

data.frame(theta = seq(from = 0, to = 1, length.out = 1e3) ) %>%
  mutate(likelihood = dbinom(x = y, size = n, prob = theta) ) %>%
  ggplot(aes(x = theta, y = likelihood) ) +
  geom_area(color = "orangered", fill = "orangered", alpha = 0.5) +
  labs(x = expression(paste(theta, " - Pr(face)") ), y = "Vraisemblance (likelihood)")

Vraisemblance versus probabilité

Si on calcule l’aire sous la courbe de cette fonction, on obtient :

\[\int_{0}^{1} 2 \theta(1 - \theta) \mathrm{d} \theta = \frac{1}{3}\]

f <- function(theta) {2 * theta * (1 - theta) }
integrate(f = f, lower = 0, upper = 1)

0.3333333 with absolute error < 3.7e-15

Vraisemblance versus probabilité

Notons que la vraisemblance de \(\theta\) pour une donnée particulière est égale à la probabilité de cette donnée pour cette valeur de \(\theta\). Cependant, la distribution de ces vraisemblances (en colonne) n’est pas une distribution de probabilité. Dans l’analyse bayésienne, les données sont considérées comme fixes et la valeur de \(\theta\) est considérée comme une variable aléatoire.

Définition du modèle (prior)

Comment définir un prior dans le cas du lancer de pièce ?

La distribution Beta

\[ \begin{align} \color{steelblue}{p(\theta \given a, b)} \ &\color{steelblue}{= \mathrm{Beta}(\theta \given a, b)} \\ & \color{steelblue}{= \theta^{a - 1}(1 - \theta)^{b - 1} / B(a, b)} \\ & \color{steelblue}{\propto \theta^{a - 1}(1 - \theta)^{b - 1}} \end{align} \]

où \(a\) et \(b\) sont deux paramètres tels que \(a \geq 0\), \(b \geq 0\), et \(B(a, b)\) est une constante de normalisation.

Interprétation des paramètres du prior Beta

• On peut exprimer l’absence de connaissance a priori par \(a = b = 1\) (distribution orange).
• On peut exprimer un prior en faveur d’une absence de biais par \(a = b \geq 2\) (distribution verte).
• On peut exprimer un biais en faveur de Face par \(a > b\) (distribution bleue).
• On peut exprimer un biais en faveur de Pile par \(a < b\) (distribution violette).

Interprétation des paramètres du prior Beta

Le niveau de certitude augmente avec la somme \(\kappa = a + b\).

• Aucune idée sur la provenance de la pièce : \(a = b = 1\) -> prior plat.
• En attendant le début de l’expérience, on a lancé la pièce 10 fois et observé 5 “Face” : \(a = b = 5\) -> prior peu informatif.
• La pièce provient de la banque de France : \(a = b = 50\) -> prior fort.

Interprétation des paramètres du prior Beta

Supposons que l’on dispose d’une estimation de la valeur la plus probable \(\omega\) du paramètre \(\theta\). On peut reparamétriser la distribution Beta en fonction du mode \(\omega\) et du niveau de certitude \(\kappa\) :

\[ \begin{align} a &= \omega(\kappa - 2) + 1 \\ b &= (1 - \omega)(\kappa - 2) + 1 &&\mbox{pour } \kappa > 2 \end{align} \]

Prior conjugué

Formellement, si \(\mathcal{F}\) est une classe de distributions d’échantillonnage \(p(y|\theta)\), et \(\mathcal{P}\) est une classe de distributions a priori pour \(\theta\), alors \(\mathcal{P}\) est conjuguée à \(\mathcal{F}\) si et seulement si :

\[ p(\theta \given y) \in \mathcal{P} \text{ for all } p(\cdot \given \theta) \in \mathcal{F} \text{ and } p(\cdot) \in \mathcal{P} \]

(p.35, Gelman et al., 2013). En d’autres termes, un prior est appelé conjugué si, lorsqu’il est converti en une distribution a posteriori en étant multiplié par la fonction de vraisemblance, il conserve la même forme. Dans notre cas, le prior Beta est un prior conjugué pour la vraisemblance binomiale, car le posterior est également une distribution Beta.

Dérivation analytique de la distribution a posteriori

Soit un prior défini par : \(\ \color{steelblue}{p(\theta \given a, b) = \mathrm{Beta}(a, b) = \frac{\theta^{a - 1}(1 - \theta)^{b - 1}}{B(a, b)} \propto \theta^{a - 1}(1 - \theta)^{b - 1}}\)

Un exemple pour digérer

On observe \(y = 7\) réponses correctes sur \(n = 10\) questions. On choisit un prior \(\mathrm{Beta}(1, 1)\), c’est à dire un prior uniforme sur \([0, 1]\). Ce prior équivaut à une connaissance a priori de 0 succès et 0 échecs (i.e., prior plat).

Un exemple pour digérer

Influence du prior sur la distribution postérieure

Cas \(n < a + b, (n = 10, a = 4, b = 16)\).

Influence du prior sur la distribution postérieure

Cas \(n = a + b, (n = 20, a = 4, b = 16)\).

Influence du prior sur la distribution postérieure

Cas \(n > a + b, (n = 40, a = 4, b = 16)\).

Ce qu’il faut retenir

The posterior distribution is always a compromise between the prior distribution and the likelihood function (Kruschke, 2015).

Ce qu’il faut retenir

Plus on a de données, moins le prior a d’influence dans l’estimation de la distribution a posteriori (et réciproquement). Attention : Lorsque le prior accorde une probabilité de 0 à certaines valeurs de \(\theta\), le modèle est incapable d’apprendre (ces valeurs sont alors considérées comme “impossibles”)…

La vraisemblance marginale

\[ \text{Posterior} = \frac{\text{Likelihood} \times \text{Prior}}{\text{Marginal likelihood}} \propto \text{Likelihood} \times \text{Prior} \]

\[ p(\theta \given \text{data}) = \frac{p(\text{data} \given \theta) \times \ p(\theta)}{p(\text{data})} \propto p(\text{data} \given \theta) \times p(\theta) \]

La vraisemblance marginale

Petit problème : \(p(\text{data})\) s’obtient en calculant la somme (pour des variables discrètes) ou l’intégrale (pour des variables continues) de la densité conjointe \(p(\text{data}, \theta)\) sur toutes les valeurs possibles de \(\theta\). Cela se complique lorsque le modèle comprend plusieurs paramètres continus…

La vraisemblance marginale

Trois méthodes pour résoudre (contourner) ce problème :

• Solution analytique \(\longrightarrow\) Utilisation d’un prior conjugué (e.g., le modèle Beta-Binomial).

• Solution discrétisée \(\longrightarrow\) Calcul de la solution sur un ensemble fini de points (grid method).

• Solution approchée \(\longrightarrow\) On échantillonne “intelligemment” l’espace conjoint des paramètres (e.g., méthodes MCMC, Cours n°05).

Distributions discrètes

Distributions continues

Problème : Cette solution est très contraignante. Idéalement, le modèle (likelihood + prior) devrait être défini à partir de l’interprétation que l’on peut faire des paramètres de ces distributions, et non pour faciliter les calculs…

La distribution postérieure, grid method

La distribution postérieure, grid method

La distribution postérieure, grid method

La distribution postérieure, grid method

La distribution postérieure, grid method

La distribution postérieure, grid method

Problème du nombre de paramètres… En affinant la grille on augmente le temps de calcul :

• 3 paramètres avec une grille de \(10^3\) noeuds = une grille de \(10^9\) points de calcul
• 10 paramètres avec une grille de \(10^3\) noeuds = une grille de \(10^{30}\) points de calcul

Échantillonner la distribution postérieure

Pour échantillonner (de manière intelligente) une distribution postérieure, on peut utiliser différentes implémentations des méthodes MCMC (e.g., Metropolis-Hastings, Hamiltonian Monte Carlo) que l’on discutera au Cours n°05. En attendant, on va travailler avec des échantillons de la distribution postérieure i) pour s’habituer en préparation aux méthodes MCMC et ii) car c’est plus simple de calculer une moyenne ou un intervalle de crédibilité sur des échantillons plutôt qu’en calculant des intégrales.

p_grid <- seq(from = 0, to = 1, length.out = 1000) # creates a grid
prior <- rep(1, 1000) # uniform prior
likelihood <- dbinom(x = 12, size = 20, prob = p_grid) # computes likelihood
posterior <- (likelihood * prior) / sum(likelihood * prior) # computes posterior
samples <- sample(x = p_grid, size = 1e3, prob = posterior, replace = TRUE) # sampling
hist(samples, main = "", xlab = expression(theta) ) # histogram

La distribution postérieure, résumé

Cas analytique :

a <- b <- 1 # paramètres du prior Beta
n <- 9 # nombre d'observations
y <- 6 # nombre de succès
p_grid <- seq(from = 0, to = 1, length.out = 1000)
posterior <- dbeta(p_grid, y + a, n - y + b) # plot(posterior)

p_grid <- seq(from = 0, to = 1, length.out = 1000)
prior <- rep(1, 1000) # uniform prior
likelihood <- dbinom(x = y, size = n, prob = p_grid)
posterior <- (likelihood * prior) / sum(likelihood * prior) # plot(posterior)

samples <- sample(x = p_grid, size = 1e4, prob = posterior, replace = TRUE) # hist(samples)

La distribution postérieure, résumé

Méthode analytique

• La distribution postérieure est décrite explicitement
• Le modèle est fortement contraint

Utiliser les échantillons pour résumer la distribution postérieure

Estimation de la tendance centrale : À partir d’un ensemble d’échantillons d’une distribution postérieure, on peut calculer la moyenne, le mode, et la médiane. Par exemple pour un prior uniforme, 10 lancers, et 3 Faces.

mode_posterior <- find_mode(samples) # en bleu
mean_posterior <- mean(samples) # en orange
median_posterior <- median(samples) # en vert

Utiliser les échantillons pour résumer la distribution postérieure

Quelle est la probabilité que le biais de la pièce \(\theta\) soit supérieur à 0.5 ?

sum(samples > 0.5) / length(samples) # équivalent à mean(samples > 0.5)

[1] 0.112

sum(samples > 0.2 & samples < 0.4) / length(samples)

[1] 0.5482

Highest density interval (HDI)

Highest density interval (HDI) :

• Le HDI indique les valeurs du paramètre qui sont les plus probables (sachant les données et le prior)
• Plus le HDI est étroit et plus le degré de certitude est élevé
• La largeur du HDI diminue avec l’augmentation du nombre de mesures

Highest density interval (HDI)

Highest density interval (HDI)

library(imsb)

set.seed(666)
p_grid <- seq(from = 0, to = 1, length.out = 1e3)
pTheta <- dbeta(p_grid, 3, 10)
massVec <- pTheta / sum(pTheta)
samples <- sample(x = p_grid, size = 1e4, replace = TRUE, prob = pTheta)

posterior_plot(samples = samples, credmass = 0.89)

Region of practical equivalence (ROPE)

Cette procédure permet d’accepter ou de rejeter une valeur nulle (null value). La région d’équivalence pratique ou region of practical equivalence (ROPE) définit un intervalle de valeurs qu’on considère comment étant “équivalentes” à la valeur nulle. La figure ci-dessous résume les décisions possibles issues de cette procédure (Kruschke, 2018).

Region of practical equivalence (ROPE)

La valeur du paramètre (e.g., \(\theta = 0.5\)) est rejetée si le HDI est entièrement hors de la ROPE. La valeur du paramètre (e.g., \(\theta = 0.5\)) est acceptée si le HDI est entièrement dans la ROPE. Si le HDI et la ROPE se chevauchent, on ne peut pas conclure…

posterior_plot(samples = samples, rope = c(0.49, 0.51) ) +
    labs(x = expression(theta) )

Comparaison de modèles

On lance une pièce 200 fois et on obtient 115 “Faces”. Est-ce que la pièce est biaisée ? Nous construisons deux modèles et essayons de savoir lequel rend le mieux compte des données.

Comparaison de modèles

Le facteur de Bayes (Bayes factor) fait le rapport des vraisemblances (marginales) des deux modèles.

\[\frac{p(\mathcal{M}_{0} \given \text{data})}{p(\mathcal{M}_{1} \given \text{data})} = \color{green}{\frac{p(\text{data} \given \mathcal{M}_{0})}{p(\text{data} \given \mathcal{M}_{1})}} \times \frac{p(\mathcal{M}_{0})}{p(\mathcal{M}_{1})}\]

Model checking

Les deux rôles de la fonction de vraisemblance :

• C’est une fonction de \(\theta\) pour le calcul de la distribution postérieure : \(\mathcal{L}(\theta \given y, n)\)
• Lorsque \(\theta\) est connu / fixé, c’est une distribution de probabilité : \(p(y \given \theta, n) \propto \theta^y(1 - \theta)^{(n - y)}\)

samples <- rbinom(n = 1e4, size = 10, prob = 0.6)

Model checking

Dans un modèle bayésien, il existe deux sources d’incertitude lorsqu’on génère des prédictions :

• Incertitude liée au processus d’échantillonnage
  -> On tire une donnée issue d’une distribution Binomiale
• Incertitude sur la valeur de \(\theta\) elle-même
  -> L’incertitude quant à la valeur de \(\theta\) est représentée par une distribution de probabilité (postérieure)

posterior <- rbeta(n = 1e4, shape1 = 16, shape2 = 10)
samples <- rbinom(n = 1e4, size = 10, prob = posterior)

Prior and posterior predictive checking

Posterior predictive checking

Travaux pratiques

Un analyste qui travaille dans une fabrique de célèbres petits pains suédois a lu un livre qui soulevait une épineuse question… Pourquoi la tartine tombe toujours du côté du beurre ? À défaut de proposer une réponse plausible, il se propose de vérifier cette assertion. La première expérience qu’il réalise consiste à faire tomber une tartine beurrée de la hauteur d’une table. Les résultats obtenus sont disponibles dans le jeu de données tartine1 du paquet imsb.

Récupérer les données

Première tâche : Récupérer les donées.

# importer les données
data <- open_data(tartine1)

# description sommaire des données
str(data)

'data.frame':   500 obs. of  2 variables:
 $ trial: int  1 2 3 4 5 6 7 8 9 10 ...
 $ side : int  1 1 0 1 0 0 1 1 1 0 ...

Questions

Proposition de solution - Question 1

La tartine n’ayant que deux faces, le résultat s’apparente à un tirage sur une loi binomiale de paramètre \(\theta\) inconnu. Quelle est la distribution postérieure du paramètre \(\theta\) ?

# nombre d'essais
nbTrial <- length(data$trial)

# nombre de "succès" (i.e., la tartine tombe du côté du beurre)
nbSuccess <- sum(data$side)

# taile de la grille
grid_size <- 1e3

# génère la grille
p_grid <- seq(from = 0, to = 1, length.out = grid_size)

# prior uniforme
prior <- rep(1, grid_size)

# calcul de la vraisemblance
likelihood <- dbinom(x = nbSuccess, size = nbTrial, prob = p_grid)

# calcul du posterior
posterior <- likelihood * prior / sum(likelihood * prior)

Proposition de solution - Question 2

Calculer le HDI à 95% de la distribution postérieure et en donner une représentation graphique.

samples <- sample(x = p_grid, prob = posterior, size = 1e3, replace = TRUE)
posterior_plot(samples = samples, credmass = 0.95) + labs(x = expression(theta) )

Proposition de solution - Question 3

Peut-on rejeter l’hypothèse nulle selon laquelle \(\theta = 0.5\) ? Non, car le HDI recouvre partiellement la ROPE…

posterior_plot(
  samples = samples, credmass = 0.95,
  compval = 0.5, rope = c(0.49, 0.51)
  ) + labs(x = expression(theta) )

Proposition de solution - Question 4

À ce stade, on ne peut pas conclure. L’analyste décide de relancer une série d’observations afin d’affiner ses résultats.

data2 <- open_data(tartine2)
str(data2)

'data.frame':   100 obs. of  2 variables:
 $ trial: int  1 2 3 4 5 6 7 8 9 10 ...
 $ side : int  0 0 1 0 0 1 1 1 0 0 ...

nbTrial2 <- length(data2$trial) # nombre d'essais
nbSucces2 <- sum(data2$side) # nombre de succès

Proposition de solution - Question 4

On utilise le posterior précédent comme prior de ce nouveau modèle.

mode1 <- find_mode(samples)
prior2 <- dbeta(p_grid, mode1 * (nbTrial - 2) + 1, (1 - mode1) * (nbTrial - 2) + 1)

data.frame(x = p_grid, y = prior2) %>%
  ggplot(aes(x = x, y = y) ) +
  geom_area(alpha = 0.8, fill = "steelblue") +
  geom_line(size = 0.8) +
  labs(x = expression(theta), y = "Densité de probabilité")

Proposition de solution - Question 4 (suite)

likelihood2 <- dbinom(x = nbSucces2, size = nbTrial2, prob = p_grid)
posterior2 <- likelihood2 * prior2 / sum(likelihood2 * prior2)
samples2 <- sample(p_grid, prob = posterior2, size = 1e4, replace = TRUE)

posterior_plot(
  samples = samples2, credmass = 0.95,
  compval = 0.5, rope = c(0.49, 0.51)
  ) + labs(x = expression(theta) )