Introduction à la modélisation statistique bayésienne

Un cours en R et Stan avec 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\:}\]

Null Hypothesis Significance Testing (NHST)

On s’intéresse aux différences de taille entre hommes et femmes. On va mesurer 100 femmes et 100 hommes.

set.seed(19) # pour reproduire les résultats
men <- rnorm(100, 175, 10) # 100 tailles d'hommes
women <- rnorm(100, 170, 10) # 100 tailles de femmes
t.test(men, women) # test de student pour différence de moyenne

    Welch Two Sample t-test

data:  men and women
t = 2.4969, df = 197.98, p-value = 0.01335
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.7658541 6.5209105
sample estimates:
mean of x mean of y 
 175.1258  171.4825

Null Hypothesis Significance Testing (NHST)

On va simuler des t-valeurs issues de données générées sous l’hypothèse d’une absence de différence entre hommes et femmes.

nsims <- 1e4 # nombre de simulations
t <- rep(NA, nsims) # initialisation d'un vecteur vide

for (i in 1:nsims) {
    
    men2 <- rnorm(100, 170, 10) # 100 tailles d'hommes
    women2 <- rnorm(100, 170, 10) # 100 tailles de femmes (de même distribution)
    t[i] <- t.test(men2, women2)$statistic # on conserve la t-valeur
    
}
# une autre manière de réaliser la même opération, sans boucle for
t <- replicate(nsims, t.test(rnorm(100, 170, 10), rnorm(100, 170, 10) )$statistic)

Null Hypothesis Significance Testing (NHST)

data.frame(t = t) %>%
    ggplot(aes(x = t) ) +
    geom_histogram() +
    labs(x = "Valeur de t", y = "Nombre d'échantillons")

Null Hypothesis Significance Testing (NHST)

data.frame(x = c(-5, 5) ) %>%
    ggplot(aes(x = x) ) +
    stat_function(fun = dt, args = list(df = t.test(men, women)$parameter), size = 1.5) +
    labs(x = "Valeur de t", y = "Densité de probabilité")

Null Hypothesis Significance Testing (NHST)

alpha <- 0.05 # seuil de significativité
abs(qt(p = alpha / 2, df = t.test(men, women)$parameter) ) # valeur de t critique
[1] 1.972019

Null Hypothesis Significance Testing (NHST)

tobs <- t.test(men, women)$statistic # valeur de t observée
tobs %>% as.numeric
[1] 2.496871

P-valeur

Une p-valeur est une aire sous la courbe (une intégrale) sous la distribution de statistiques de test sous l’hypothèse nulle (i.e., étant admis que l’hypothèse nulle est vraie). La p-valeur indique la probabilité d’observer la valeur de la statistique de test (e.g., une valeur de t) observée, ou une valeur plus extrême, sous l’hypothèse nulle.

\[p[\mathbf{t}(\mathbf{x}^{\text{rep}}; \mathcal{H}_{0}) \geq t(x)]\]

t.test(men, women)$p.value
[1] 0.01334509
tvalue <- abs(t.test(men, women)$statistic)
df <- t.test(men, women)$parameter
2 * integrate(dt, tvalue, Inf, df = df)$value
[1] 0.01334509
2 * (1 - pt(abs(t.test(men, women)$statistic), t.test(men, women)$parameter) )
         t 
0.01334509

Intervalles de confiance

Les intervalles de confiance peuvent être interprétés comme des régions de significativité (ou des régions de compatibilité avec l’hypothèse nulle). Par conséquent, un intervalle de confiance contient la même information qu’une p-valeur et s’interprète de manière similaire.

On ne peut pas dire qu’un intervalle de confiance a une probabilité de 95% de contenir la vraie valeur (i.e., la valeur dans la population) de \(\theta\) (cf. the inverse fallacy), contrairement à l’intervalle de crédibilité bayésien.

Un intervalle de confiance à 95% représente une degré de “recouvrement” (coverage). Le “95%” fait référence à une propriété fréquentiste (i.e., sur le long-terme) de la procédure, mais ne fait pas référence au paramètre \(\theta\). Autrement dit, sur le long-terme, 95% des intervalles de confiance à 95% que l’on pourrait calculer (dans une réplication exacte de notre expérience) contiendraient la valeur du paramètre dans la population (i.e., la “vraie” valeur de \(\theta\)). Cependant, nous ne pouvons pas dire qu’un intervalle de confiance en particulier a une probabilité de 95% de contenir la “vraie” valeur de \(\theta\)… soit ce dernier contient la vraie valeur de \(\theta\), soit il ne la contient pas.

Facteur de Bayes

On compare deux modèles :

  • \(\mathcal{H}_{0}: \mu_{1} = \mu_{2} \rightarrow \delta = 0\)
  • \(\mathcal{H}_{1}: \mu_{1} \neq \mu_{2} \rightarrow \delta \neq 0\)

\[ \underbrace{\dfrac{p(\mathcal{H}_{0} \given D)}{p(\mathcal{H}_{1} \given D)}}_{\text{posterior odds}} = \underbrace{\dfrac{p(D \given \mathcal{H}_{0})}{p(D \given \mathcal{H}_{1})}}_{\text{Bayes factor}} \times \underbrace{\dfrac{p(\mathcal{H}_{0})}{p(\mathcal{H}_{1})}}_{\text{prior odds}} \]

\[\text{evidence}\ = p(D \given \mathcal{H}) = \int p(\theta \given \mathcal{H}) p(D \given \theta, \mathcal{H}) \text{d}\theta\]

L’évidence en faveur d’un modèle correspond à la vraisemblance marginale d’un modèle (le dénominateur du théorème de Bayes), c’est à dire à la vraisemblance moyennée sur le prior… Ce qui fait du facteur de Bayes un ratio de vraisemblances, pondéré par (ou moyenné sur) le prior.

Facteur de Bayes, exemple d’application

On lance une pièce 100 fois et on essaye d’estimer la probabilité \(\theta\) (le biais de la pièce) d’obtenir Face. On compare deux modèles qui diffèrent par leur prior sur \(\theta\).

\[ \begin{align} \mathcal{M_{1}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(6, 10) \\ \end{align} \]

\[ \begin{align} \mathcal{M_{2}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(20, 12) \\ \end{align} \]

Facteur de Bayes, exemple d’application


\[ \begin{align} \mathcal{M_{1}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(6, 10) \\ \end{align} \]

\[ \begin{align} \mathcal{M_{2}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(20, 12) \\ \end{align} \]


\[ \text{BF}_{12} = \dfrac{p(D | \mathcal{M_{1}} )}{p(D | \mathcal{M_{2}} )} = \dfrac{\int p(\theta | \mathcal{M_{1}} ) p(D | \theta, \mathcal{M_{1}} ) \text{d}\theta}{\int p(\theta | \mathcal{M_{2}}) p(D | \theta, \mathcal{M_{2}} ) \text{d}\theta} = \dfrac{\int \mathrm{Binomial}(n, \theta)\mathrm{Beta}(6, 10)\text{d}\theta}{\int \mathrm{Binomial}(n, \theta)\mathrm{Beta}(20, 12) \text{d}\theta} \]

Facteur de Bayes, exemple d’application