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\:}\]

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.

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\)

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

Le facteur de Bayes est la nouvelle p-valeur

Attention à ne pas interpréter le BF comme un rapport des chances a posteriori (posterior odds)…

Le BF est un facteur qui nous indique de combien notre rapport des chances a priori (prior odds) doit changer, au vu des données. Il ne nous dit pas quelle est l’hypothèse la plus probable, sachant les données (sauf si les prior odds sont de \(1/1\)). Par exemple :

Double sens

Dans le cadre bayésien, le terme de prior peut faire référence soit :

• À la probabilité a priori ou a posteriori d’un modèle (par rapport à un autre modèle), c’est à dire \(\Pr(\mathcal{M}_{i})\). Voir ce blogpost.

• Aux priors qu’on définit sur les paramètres d’un modèle, par exemple \(\beta \sim \mathrm{Normal}(0, 1)\). Voir ce blogpost.

Comparaison de modèles

Deux problèmes récurrents à éviter en modélisation : le sous-apprentissage (underfitting) et sur-apprentissage (overfitting). Comment s’en sortir ?

• Utiliser des priors pour régulariser, pour contraindre le posterior (i.e., accorder moins de poids à la vraisemblance).
• Utiliser des critères d’information (e.g., AIC, WAIC).

Overfitting

ppnames <- c("afarensis", "africanus", "habilis", "boisei",
        "rudolfensis", "ergaster", "sapiens")
brainvolcc <- c(438, 452, 612, 521, 752, 871, 1350)
masskg <- c(37.0, 35.5, 34.5, 41.5, 55.5, 61.0, 53.5)

d <- data.frame(species = ppnames, brain = brainvolcc, mass = masskg)

d %>%
    ggplot(aes(x = mass, y = brain, label = species) ) +
    geom_point() +
    ggrepel::geom_label_repel(hjust = 0, nudge_y = 50, size = 5) +
    xlim(30, 70)

Overfitting

mod1.1 <- lm(brain ~ mass, data = d)
(var(d$brain) - var(residuals(mod1.1) ) ) / var(d$brain)

[1] 0.490158

mod1.2 <- lm(brain ~ mass + I(mass^2), data = d)
(var(d$brain) - var(residuals(mod1.2) ) ) / var(d$brain)

[1] 0.5359967

mod1.3 <- lm(brain ~ mass + I(mass^2) + I(mass^3), data = d)
(var(d$brain) - var(residuals(mod1.3) ) ) / var(d$brain)

[1] 0.6797736

Overfitting

mod1.4 <- lm(brain ~ mass + I(mass^2) + I(mass^3) + I(mass^4), data = d)
(var(d$brain) - var(residuals(mod1.4) ) ) / var(d$brain)

[1] 0.8144339

mod1.5 <- lm(brain ~ mass + I(mass^2) + I(mass^3) + I(mass^4) +
    I(mass^5), data = d)
(var(d$brain) - var(residuals(mod1.5) ) ) / var(d$brain)

[1] 0.988854

mod1.6 <- lm( brain ~ mass + I(mass^2) + I(mass^3) + I(mass^4) +
    I(mass^5) + I(mass^6), data = d)
(var(d$brain) - var(residuals(mod1.6) ) ) / var(d$brain)

[1] 1

Overfitting

Overfitting

Overfitting

Underfitting

\[ \begin{align} v_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha \\ \end{align} \]

mod1.7 <- lm(brain ~ 1, data = d)

Théorie de l’information

Il nous faut une autre mesure des capacités de prédiction pour évaluer nos modèles. Idéalement, on voudrait pouvoir mesurer la “distance” entre notre modèle et le “full model” (i.e., la “réalité”, la nature)… mais on ne sait (encore) pas faire cela.

Incertitude (exemple)

Exemple de prédiction météorologique. Si on imagine que la probabilité qu’il pleuve ou qu’il fasse beau (à Grenoble) est de \(p_{1} = 0.3\) et \(p_{2} = 0.7\).

p <- c(0.3, 0.7) # distribution des probabilités de pluie et soleil à Grenoble
- sum(p * log(p) ) # équivalent à sum(p * log(1 / p) )

[1] 0.6108643

p <- c(0.01, 0.99) # distribution des probabilités de pluie et soleil à Abu Dhabi
- sum(p * log(p) ) # toujours équivalent à sum(p * log(1 / p) )

[1] 0.05600153

Divergence

On a donc un moyen de quantifier l’incertitude. Comment utiliser cette mesure pour quantifier la distance entre notre modèle et la réalité ?

Divergence

\[ D_{KL}(p, q) = \sum_{i} p_{i}\big(\log(p_{i}) - \log(q_{i})\big) = \sum_{i} p_{i} \log\bigg(\frac{p_{i}}{q_{i}}\bigg) \]

Par exemple, supposons que la “véritable” distribution de nos évènements (soleil vs pluie) soit \(p_{1} = 0.3\) et \(p_{2} = 0.7\). Si nous pensons plutôt que ces évènements arrivent avec une probabilité \(q_{1} = 0.25\) et \(q_{2} = 0.75\), quelle quantité d’incertitude avons-nous ajoutée ?

p <- c(0.3, 0.7)
q <- c(0.25, 0.75)

sum(p * log(p / q) )

[1] 0.006401457

# NB : La divergence n'est pas symétrique...
sum(q * log(q / p) ) 

[1] 0.006164264

Divergence

La divergence n’est pas symétrique (ce n’est pas une distance)…

Entropie croisée et divergence

Entropie croisée : \(H(p, q) = \sum_{i} p_{i} \log (q_{i})\)

sum(p * (log(q) ) )

- sum(p * (log(q) - log(p) ) )

[1] 0.006401457

sum(p * log(p / q) )

[1] 0.006401457

Vers la déviance…

OK, mais nous ne connaissons pas la distribution “cible” (la réalité), à quoi cela peut donc nous servir ?

Vers la déviance…

On peut donc utiliser \(\mathbb{E} [\log(q_{i})]\) et \(\mathbb{E} [\log(r_{i})]\) comme estimateurs de la distance relative entre chaque modèle et notre distribution cible. On a seulement besoin de la log-probabilité moyenne des modèles. Comme on ne connaît pas la distribution cible, cela veut dire qu’on ne peut pas interpréter cette quantité en termes absolus mais seulement en termes relatifs. Ce qui nous intéresse c’est \(\mathbb{E} [\log(q_{i})] - \mathbb{E} [\log(r_{i})]\) (i.e., des différences de divergences).

Vers la déviance (exemple d’application)

Exemple d’une distribution cible uniforme \(p\), et 4 modèles \(q1, q2, q3, q4\) (exemple tiré de cet article).

p <- c(1/3, 1/3, 1/3) # distribution cible (distribution uniforme avec p = 1/3)
q1 <- c(0.36, 0.48, 0.16) # modèle q1
q2 <- c(0.2, 0.6, 0.2) # modèle q2
q3 <- c(0.1, 0.8, 0.1) # modèle q3
q4 <- c(0.3, 0.4, 0.3) # modèle q4

Vers la déviance (exemple d’application)

divergence_q1 <- sum(p * log(p / q1) ) # divergence du modèle q1
divergence_q2 <- sum(p * log(p / q2) ) # divergence du modèle q2
divergence_q3 <- sum(p * log(p / q3) ) # divergence du modèle q3
divergence_q4 <- sum(p * log(p / q4) ) # divergence du modèle q4

# vecteur de divergences
divergences <- c(divergence_q1, divergence_q2, divergence_q3, divergence_q4)

# vecteur de negative log-scores
neg_log_scores <- c(-sum(log(q1) ), -sum(log(q2) ), -sum(log(q3) ), -sum(log(q4) ) )

Déviance

Pour approximer la valeur de \(\mathbb{E} [\log(q_{i})]\), on peut utiliser la déviance d’un modèle, qui est une mesure du fit relatif du modèle.

\[D(q) = -2 \sum_{i} \log(q_{i})\]

où \(i\) indice chaque observation et \(q_{i}\) est la vraisemblance de chaque observation.

# on standardise la variable "mass"
d$mass.s <- scale(d$mass)

# on "fit" un modèle de régression linéaire Gaussien
mod1.8 <- lm(formula = brain ~ mass.s, data = d)

# calcul de la déviance
-2 * logLik(mod1.8)

'log Lik.' 94.92499 (df=3)

Déviance

# on récupère les paramètres estimés via MLE (intercept et pente)
alpha <- coef(mod1.8)[1]
beta <- coef(mod1.8)[2]

# calcul de la log-vraisemblance
ll <- sum(dnorm(
    x = d$brain,
    mean = alpha + beta * d$mass.s,
    sd = sigma(mod1.8),
    log = TRUE
    ) )

# calcul de la déviance
(-2) * ll

[1] 95.2803

Log-pointwise-predictive density

Dans le monde fréquentiste, on multiplie le log-score par \(-2\) car la différence de deux déviances suit une loi de \(\chi^{2}\), ce qui est utile pour tester l’hypothèse nulle. Mais sans besoin de tester l’hypothèse nulle (avec une forme prédéfinie), on peut très bien travailler directement avec le log-score \(S(q) = \sum_{i} \log(q_{i})\), qu’on traite comme une estimation de \(\mathbb{E} [\log(q_{i})]\).

In-sample and out-of-sample

La déviance a le même problème que le \(R^{2}\), lorsqu’elle est calculée sur l’échantillon observé. Dans ce cas, on l’appelle déviance in-sample.

In sample and out of sample deviance

\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i}, 1) \\ \mu_{i} &= (0.15) x_{1, i} - (0.4) x_{2, i} \\ \end{align} \]

On a réalisé ce processus 10.000 fois pour cinq modèles de régression linéaire de complexité croissante. Les points bleus représentent la déviance calculée sur l’échantillon d’apprentissage et les points noirs la déviance calculée sur l’échantillon de test.

Régularisation

Une autre manière de lutter contre le sur-apprentissage (overfitting) est d’utiliser des priors “sceptiques” qui vont venir ralentir l’apprentissage réalisé sur les données (i.e., accorder plus de poids au prior).

Régularisation

Comment décider de la précision du prior ? Est-ce que le prior est “assez” régularisateur ou pas ? On peut diviser le jeu de données en deux parties (training et testing) afin de choisir le prior qui produit la déviance out-of-sample la plus faible. On appelle cette stratégie la validation croisée (cross-validation).

Critères d’information

On mesure ici la différence entre la déviance in-sample (en bleu) et la déviance out-of-sample (en noir). On remarque que la déviance out-of-sample est presque exactement égale à la déviance in-sample, plus deux fois le nombre de paramètres du modèle…

Akaike information criterion

L’AIC fournit une approximation de la déviance out-of-sample :

\[\text{AIC} = D_{\text{train}} + 2p = -2 \text{lppd} + 2p \approx D_{\text{test}}\]

où \(p\) est le nombre de paramètres libres (i.e., à estimer) dans le modèle. L’AIC donne donc une approximation des capacités de prédiction (out-of-sample) du modèle.¹

Deviance information criterion

Une condition d’application de l’AIC est que les priors soient plats ou dépassés par la vraisemblance (e.g., lorsqu’on a beaucoup de données). Le DIC est un indice qui ne requiert pas cette condition, en s’accommodant de priors informatifs.

Widely applicable information criterion

Une condition d’application de l’AIC et du DIC est que la distribution a posteriori soit une distribution gaussienne multivariée. Le WAIC relâche cette condition et il est souvent plus précis que le DIC.

Widely applicable information criterion

library(brms)
data(cars)

priors <- c(
  prior(normal(0, 100), class = Intercept),
  prior(normal(0, 10), class = b),
  prior(exponential(0.1), class = sigma)
  )

mod1 <- brm(
  formula = dist ~ 1 + speed,
  prior = priors,
  data = cars
  )

Widely applicable information criterion

\[\text{lppd}(y, \Theta) = \sum_{i} \log \frac{1}{S} \sum_{s} p(y_{i} \given \Theta_{s})\]

# pointwise log-likelihood (S samples * N observations)
ll <- log_lik(mod1) %>% data.frame()

(lppd <-
  ll %>% 
  pivot_longer(
    cols = everything(),
    names_to = "i",
    values_to = "loglikelihood"
    ) %>%
  # log-likelihood to likelihood
  mutate(likelihood = exp(loglikelihood) ) %>%
  # pour chaque observation
  group_by(i) %>%
  # logarithme de la vraisemblance moyenne
  summarise(log_mean_likelihood = log(mean(likelihood) ) ) %>%
  # on prend la somme de ces valeurs
  summarise(lppd = sum(log_mean_likelihood) ) %>%
  ungroup() %>%
  pull(lppd) )

[1] -206.6047

Widely applicable information criterion

La deuxième partie du calcul du WAIC est le nombre de paramètres effectif, \(p_{\text{WAIC}}\). On définit \(\text{var}_{\theta}\log[p(y_{i} | \theta)]\) comme la variance de la log-vraisemblance pour chaque observation \(i\) de l’échantillon d’entraînement.

\[p_{\text{WAIC}} = \sum_{i} \text{var}_{\theta}\log[p(y_{i} \given \theta)]\]

(pwaic <-
  ll %>% 
  pivot_longer(
    everything(),
    names_to = "i",
    values_to = "loglikelihood"
    ) %>%
  # pour chaque observation
  group_by(i) %>%
  # variance de la log-vraisemblance
  summarise(var_loglikelihood = var(loglikelihood) ) %>%
  # somme de ces variances
  summarise(pwaic = sum(var_loglikelihood) ) %>%
  ungroup() %>%
  pull(pwaic) )

[1] 3.342853

Widely applicable information criterion

Ensuite, le WAIC est défini par :

\[ \text{WAIC}(y, \Theta) = - 2 \big(\text{lppd} - \underbrace{\sum_{i} \text{var}_{\theta}\log[p(y_{i} \given \theta) ]}_{\text{penalty term}}\big) \]

(WAIC <- -2 * (lppd - pwaic) )

[1] 419.8952

Widely applicable information criterion

Le WAIC est également un estimateur de la déviance out-of-sample. La fonction brms::waic() permet de le calculer directement.

waic(mod1)

Computed from 4000 by 50 log-likelihood matrix

          Estimate   SE
elpd_waic   -209.9  6.5
p_waic         3.3  1.2
waic         419.9 13.0

2 (4.0%) p_waic estimates greater than 0.4. We recommend trying loo instead. 

Critères d’information et régularisation

Le DIC et le WAIC peuvent être conceptualisés (au même titre que l’AIC) comme des approximations de la déviance out-of-sample. On remarque que le WAIC produit des approximations plus précises que le DIC, et que l’utilisation de priors régularisateurs permet de réduire la déviance out-of-sample.

Et ensuite ?

Sélection de modèle : On choisit le meilleur modèle un utilisant un des outils présentés et on base nos conclusions sur les paramètres estimés par ce meilleur modèle.

Comparaison de modèles

On essaye de prédire les kg par gramme de lait (kcal.per.g) avec les prédicteurs neocortex et le logarithme de mass. Nous allons ensuite fitter 4 modèles qui correspondent aux 4 combinaisons possibles de prédicteurs et les comparer en utilisant le WAIC.

library(imsb)
d <- open_data(milk)
d <- milk[complete.cases(milk), ] # removing NAs
d$neocortex <- d$neocortex.perc / 100 # rescaling explanatory variable
head(d)

              clade            species kcal.per.g perc.fat perc.protein
1     Strepsirrhine     Eulemur fulvus       0.49    16.60        15.42
6  New World Monkey Alouatta seniculus       0.47    21.22        23.58
7  New World Monkey         A palliata       0.56    29.66        23.46
8  New World Monkey       Cebus apella       0.89    53.41        15.80
10 New World Monkey         S sciureus       0.92    50.58        22.33
11 New World Monkey   Cebuella pygmaea       0.80    41.35        20.85
   perc.lactose mass neocortex.perc neocortex
1         67.98 1.95          55.16    0.5516
6         55.20 5.25          64.54    0.6454
7         46.88 5.37          64.54    0.6454
8         30.79 2.51          67.64    0.6764
10        27.09 0.68          68.85    0.6885
11        37.80 0.12          58.85    0.5885

Comparaison de modèles

mod2.1 <- brm(
  formula = kcal.per.g ~ 1,
  family = gaussian,
  data = d,
  prior = c(
    prior(normal(0, 100), class = Intercept),
    prior(exponential(0.01), class = sigma)
    ),
  iter = 2000, warmup = 1000,
  backend = "cmdstanr" # on peut changer le backend de brms
  )

mod2.2 <- brm(
  formula = kcal.per.g ~ 1 + neocortex,
  family = gaussian,
  data = d,
  prior = c(
    prior(normal(0, 100), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(exponential(0.01), class = sigma)
    ),
  iter = 2000, warmup = 1000,
  backend = "cmdstanr"
  )

Comparaison de modèles

mod2.3 <- brm(
  formula = kcal.per.g ~ 1 + log(mass),
  family = gaussian,
  data = d,
  prior = c(
    prior(normal(0, 100), class = Intercept),
    prior(exponential(0.01), class = sigma)
    ),
  iter = 2000, warmup = 1000,
  backend = "cmdstanr"
  )

mod2.4 <- brm(
  formula = kcal.per.g ~ 1 + neocortex + log(mass),
  family = gaussian,
  data = d,
  prior = c(
    prior(normal(0, 100), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(exponential(0.01), class = sigma)
    ),
  iter = 2000, warmup = 1000,
  backend = "cmdstanr"
  )

Comparaison de modèles

On peut utiliser la méthode update() qui permet de fitter plus rapidement un nouveau modèle qui ressemble à un modèle déjà existant.

mod2.3 <- update(
  object = mod2.2,
  newdata = d,
  formula = kcal.per.g ~ 1 + log(mass)
  )

mod2.4 <- update(
  object = mod2.3,
  newdata = d,
  formula = kcal.per.g ~ 1 + neocortex + log(mass)
  )

Comparaison de modèles

# calcul du WAIC et ajout du WAIC à chaque modèle

mod2.1 <- add_criterion(mod2.1, "waic")
mod2.2 <- add_criterion(mod2.2, "waic")
mod2.3 <- add_criterion(mod2.3, "waic")
mod2.4 <- add_criterion(mod2.4, "waic")

# comparaison des WAIC de chaque modèle

w <- loo_compare(mod2.1, mod2.2, mod2.3, mod2.4, criterion = "waic")
print(w, simplify = FALSE)

       elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic  se_waic
mod2.4   0.0       0.0     8.4       2.6          3.1    0.8     -16.8   5.1  
mod2.3  -3.9       1.7     4.5       2.1          2.0    0.4      -9.0   4.2  
mod2.1  -4.0       2.4     4.4       1.9          1.3    0.3      -8.8   3.7  
mod2.2  -4.8       2.5     3.6       1.6          1.9    0.3      -7.2   3.3  

Akaike’s weights

Le poids d’un modèle est une estimation de la probabilité que ce modèle fera les meilleures prédictions possibles sur un nouveau jeu de données, conditionnellement au set de modèles considéré.

\[w_{i} = \frac{\exp \left(-\frac{1}{2}\text{dWAIC}_{i}\right)}{\sum_{j=1}^{m} \exp \left(-\frac{1}{2}\text{dWAIC}_{j}\right)}\]

model_weights(mod2.1, mod2.2, mod2.3, mod2.4, weights = "waic") %>% round(digits = 3)

mod2.1 mod2.2 mod2.3 mod2.4 
 0.017  0.008  0.019  0.955 

Moyennage de modèles (model averaging)

Pourquoi ne conserver uniquement le premier modèle et oublier les autres ? Une autre stratégie consisterait à pondérer les prédictions des modèles par leurs poids respectifs. C’est ce qu’on appelle le moyennage de modèles (model averaging).

• Calculer le WAIC de chaque modèle.
• Calculer le poids de chaque modèle.
• Simuler des données à partir de chaque modèle.
• Combiner ces valeurs simulées dans un ensemble de prédictions pondérées par le poids du modèle.

Moyennage de modèles (model averaging)

On peut utiliser la fonction brms::pp_average() qui pondère les prédictions de chaque modèle par leur poids.

# grille de valeurs pour lesquelles on va générer des prédictions
new_data <- data.frame(
  neocortex = seq(from = 0.5, to = 0.8, length.out = 30),
  mass = 4.5
  )

# prédictions du modèle mod2.4
f <- fitted(mod2.4, newdata = new_data) %>%
  as.data.frame() %>%
  bind_cols(new_data)

# prédictions moyennées sur les 4 modèles
averaged_predictions <- pp_average(
  mod2.1, mod2.2, mod2.3, mod2.4,
  weights = "waic",
  method  = "fitted",
  newdata = new_data
  ) %>%
  as.data.frame() %>%
  bind_cols(new_data)

Moyennage de modèles (model averaging)

Voici les prédictions de tous les modèles considérés, pondérés par leur poids respectif. Comme le modèle mod2.4 concentrait quasiment tout le poids, il fait sens que cette prédiction moyennée soit similaire aux prédictions du modèle mod2.4.

R-squared

bayes_R2(mod2.4) %>% round(digits = 3)

   Estimate Est.Error  Q2.5 Q97.5
R2    0.496      0.13 0.167 0.667

posterior_plot(samples = bayes_R2(mod2.4, summary = FALSE)[, 1])

R-squared

posterior_plot(
    samples = bayes_R2(mod2.4, summary = FALSE)[, 1] -
        bayes_R2(mod2.3, summary = FALSE)[, 1],
    compval = 0
    ) + labs(x = "Différence de R2")

Conclusions

Se méfier des interprétations intuitives (et souvent erronées) de la p-valeur, des intervalles de confiance, ou du facteur de Bayes.

Travaux pratiques

# import des données "howell"
d <- open_data(howell) %>% mutate(age = scale(age) )

# on définit une graine (afin de pouvoir reproduire les résultats)
set.seed(666)

# on échantillonne des lignes du jeu de données
i <- sample(1:nrow(d), size = nrow(d) / 2)

# on définit l'échantillon d'entraînement
d1 <- d[i, ]

# on définit l'échantillon de test
d2 <- d[-i, ]

Nous avons maintenant deux dataframes, de 272 lignes chacune. On va utiliser d1 pour fitter nos modèles et d2 pour les évaluer.

Travaux pratiques

Soit \(h_{i}\) les valeurs de taille et \(x_{i}\) les valeurs centrées d’âge, sur la ligne \(i\). Construisez les modèles suivants avec d1, en utilisant brms::brm() et des priors faiblement régularisateurs.

Travaux pratiques

Solution - Question 1

mod3.1 <- brm(
  formula = height ~ 1 + age,
  family = gaussian(),
  data = d1,
  prior = c(
    prior(normal(0, 100), class = Intercept),
    prior(exponential(0.01), class = sigma)
    ),
  backend = "cmdstanr"
  )

mod3.2 <- update(
  mod3.1,
  newdata = d1,
  formula = height ~ 1 + age + I(age^2)
  )

mod3.3 <- update(
  mod3.1,
  newdata = d1,
  formula = height ~ 1 + age + I(age^2) + I(age^3)
  )

Solution - Question 1

mod3.4 <- update(
  mod3.1,
  newdata = d1,
  formula = height ~ 1 + age + I(age^2) + I(age^3) + I(age^4)
  )

mod3.5 <- update(
  mod3.1,
  newdata = d1,
  formula = height ~ 1 + age + I(age^2) + I(age^3) + I(age^4) + I(age^5)
  )

mod3.6 <- update(
  mod3.1,
  newdata = d1,
  formula = height ~ 1 + age + I(age^2) + I(age^3) + I(age^4) + I(age^5) + I(age^6)
  )

Solution - Question 1

# calcul du WAIC et ajout du WAIC à chaque modèle

mod3.1 <- add_criterion(mod3.1, "waic")
mod3.2 <- add_criterion(mod3.2, "waic")
mod3.3 <- add_criterion(mod3.3, "waic")
mod3.4 <- add_criterion(mod3.4, "waic")
mod3.5 <- add_criterion(mod3.5, "waic")
mod3.6 <- add_criterion(mod3.6, "waic")

# comparaison des WAIC de chaque modèle

mod_comp <- loo_compare(mod3.1, mod3.2, mod3.3, mod3.4, mod3.5, mod3.6, criterion = "waic")
print(mod_comp, digits = 2, simplify = FALSE)

       elpd_diff se_diff  elpd_waic se_elpd_waic p_waic   se_p_waic waic    
mod3.4     0.00      0.00  -954.40     12.41         5.74     0.79   1908.79
mod3.5    -0.70      0.47  -955.10     12.48         6.46     0.87   1910.20
mod3.6    -1.29      1.03  -955.69     12.27         7.44     0.93   1911.37
mod3.3   -21.90      6.37  -976.30     11.83         6.11     1.31   1952.60
mod3.2  -132.62     13.63 -1087.02     11.29         5.51     1.28   2174.03
mod3.1  -259.20     15.32 -1213.60     10.91         3.43     0.42   2427.20
       se_waic 
mod3.4    24.83
mod3.5    24.95
mod3.6    24.53
mod3.3    23.66
mod3.2    22.58
mod3.1    21.82

Solution - Question 2

# on crée un vecteur de valeurs possibles pour "age"
age_seq <- data.frame(age = seq(from = -2, to = 3, length.out = 1e2) )

# on récupère les prédictions du modèle pour ces valeurs
mu <- data.frame(fitted(mod3.1, newdata = age_seq) ) %>% bind_cols(age_seq)

# on récupère les prédictions du modèle pour ces valeurs
pred_age <- data.frame(predict(mod3.1, newdata = age_seq) ) %>% bind_cols(age_seq)

# on affiche les dix premières prédictions
head(pred_age, 10)

   Estimate Est.Error     Q2.5    Q97.5       age
1  100.4499  21.13775 58.51867 142.4313 -2.000000
2  101.0368  20.87114 60.20049 141.1053 -1.949495
3  101.8842  20.93198 60.09004 141.7364 -1.898990
4  103.0255  21.18667 60.41278 144.7774 -1.848485
5  104.5247  21.22701 63.20084 145.5188 -1.797980
6  104.9453  20.88182 65.70421 146.6455 -1.747475
7  105.6927  20.74311 64.78789 145.7145 -1.696970
8  106.8630  20.86651 65.64000 148.4519 -1.646465
9  107.6850  20.79897 66.01354 148.0521 -1.595960
10 108.9927  20.77924 67.99611 149.7180 -1.545455

Solution - Question 2

d1 %>%
  ggplot(aes(x = age, y = height) ) +
  geom_ribbon(
    data = mu, aes(x = age, ymin = Q2.5, ymax = Q97.5),
    alpha = 0.8, inherit.aes = FALSE
    ) +
  geom_smooth(
    data = pred_age, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
    stat = "identity", color = "black", alpha = 0.5, size = 1
    ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  labs(x = "Age", y = "Height")

Solution - Question 2

Solution - Question 3

# prédictions moyennées sur les 4 modèles
averaged_predictions_mu <- pp_average(
  mod3.1, mod3.2, mod3.3, mod3.4, mod3.5, mod3.6,
  weights = "waic",
  method  = "fitted",
  newdata = age_seq
  ) %>%
  as.data.frame() %>%
  bind_cols(age_seq)

# prédictions moyennées sur les 4 modèles
averaged_predictions_age <- pp_average(
  mod3.1, mod3.2, mod3.3, mod3.4, mod3.5, mod3.6,
  weights = "waic",
  method  = "predict",
  newdata = age_seq
  ) %>%
  as.data.frame() %>%
  bind_cols(age_seq)

Solution - Question 3

d1 %>%
  ggplot(aes(x = age, y = height) ) +
  geom_ribbon(
    data = averaged_predictions_mu, aes(x = age, ymin = Q2.5, ymax = Q97.5),
    alpha = 0.8, inherit.aes = FALSE
    ) +
  geom_smooth(
    data = averaged_predictions_age, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
    stat = "identity", color = "black", alpha = 0.5, size = 1
    ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  labs(x = "Age", y = "Height", title = "Model-averaged predictions")

Solution - Question 4

# calcul de la log-vraisemblance (log-likelihood) du modèle mod3.1
# chaque ligne est une itération et chaque colonne une observation
log_lik_mod3.1 <- log_lik(mod3.1)

# NB : la déviance possède également une distribution dans le monde bayésien...
dev.mod3.1 <- mean(-2 * rowSums(log_lik_mod3.1) )

# calcul de la log-vraisemblance (log-likelihood) du modèle mod3.2
dev.mod3.2 <- mean(-2 * rowSums(log_lik(mod3.2) ) )

# calcul de la log-vraisemblance (log-likelihood) du modèle mod3.3
dev.mod3.3 <- mean(-2 * rowSums(log_lik(mod3.3) ) )

# calcul de la log-vraisemblance (log-likelihood) du modèle mod3.4
dev.mod3.4 <- mean(-2 * rowSums(log_lik(mod3.4) ) )

# calcul de la log-vraisemblance (log-likelihood) du modèle mod3.5
dev.mod3.5 <- mean(-2 * rowSums(log_lik(mod3.5) ) )

# calcul de la log-vraisemblance (log-likelihood) du modèle mod3.6
dev.mod3.6 <- mean(-2 * rowSums(log_lik(mod3.6) ) )

Solution - Question 4

deviances <- c(dev.mod3.1, dev.mod3.2, dev.mod3.3, dev.mod3.4, dev.mod3.5, dev.mod3.6)
comparison <- mod_comp %>% data.frame %>% select(waic) %>% rownames_to_column()
waics <- comparison %>% arrange(rowname) %>% pull(waic)