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

Introduction

Cinq problèmes, cinq jeux de données. Le but est de comprendre et d’analyser ces données pour répondre à une (ou plusieurs) question(s) théorique(s).

Vous devrez écrire le modèle mathématique, puis fitter ce modèle en utilisant brms.

Ensuite, vous devrez évaluer le modèle, interpréter les résultats, et écrire un paragraphe de résultats (de type article) pour décrire vos analyses et vos conclusions.

Les problèmes sont classés par ordre croissant de difficulté. Vous pouvez travailler individuellement ou par groupe, et des propositions de correction sont disponibles à la suite des énoncés.

Problème n°1

Peut-on prédire la taille d’un individu par la taille de ses parents ?

library(tidyverse)
library(imsb)

d1 <- open_data(parents)
head(d1, 10)

   gender height mother father
1       M   62.5     66     70
2       M   64.6     58     69
3       M   69.1     66     64
4       M   73.9     68     71
5       M   67.1     64     68
6       M   64.4     62     66
7       M   71.1     66     74
8       M   71.0     63     73
9       M   67.4     64     62
10      M   69.3     65     69

Problème n°2

Les données suivantes documentent le naufrage du titanic. La colonne pclass indique la classe dans laquelle chaque passager voyageait (un proxy pour le statut socio-économique), tandis que la colonne parch indique le nombre de parents et enfants à bord.

Peut-on prédire la survie d’un passager grâce à ces informations ?

d2 <- open_data(titanic)
head(d2, 10)

   survival pclass gender age parch
1         0  upper   male  22     0
2         1  lower female  38     0
3         1  upper female  26     0
4         1  lower female  35     0
5         0  upper   male  35     0
6         0  lower   male  54     0
7         0  upper   male   2     1
8         1  upper female  27     2
9         1  upper female   4     1
10        1  lower female  58     0

Problème n°3

Ce jeu de données recense des informations sur le diamètre (colonne diam) de 80 pommes (chaque pomme étant identifiée par la colonne id), poussant sur 10 arbres différents (colonne tree). On a mesuré ce diamètre pendant 6 semaines successives (colonne time).

Que peut-on dire de la pousse de ces pommes, tout en considérant les structures de dépendance existant dans les données (i.e., chaque pomme poussait sur un arbre différent) ?

d3 <- open_data(apples)
head(d3, 10)

   tree apple id time diam
1     1     1  1    1 2.90
2     1     1  1    2 2.90
3     1     1  1    3 2.90
4     1     1  1    4 2.93
5     1     1  1    5 2.94
6     1     1  1    6 2.94
7     1     4  4    1 2.86
8     1     4  4    2 2.90
9     1     4  4    3 2.93
10    1     4  4    4 2.96

Problème n°4

Ces données recensent le nombre de candidatures pour 6 départements (colonne dept) à Berkeley (données disponibles dans le paquet rethinking). La colonne admit indique le nombre de candidatures acceptées et la colonne reject le nombre de candidatures rejetées (la colonne applications est simplement la somme des deux), en fonction du sexe des candidats (applicant.gender).

On veut savoir s’il existe un biais lié au sexe dans l’admission des étudiants à Berkeley.

d4 <- open_data(admission)
head(d4, 10)

   dept applicant.gender admit reject applications
1     A             male   512    313          825
2     A           female    89     19          108
3     B             male   353    207          560
4     B           female    17      8           25
5     C             male   120    205          325
6     C           female   202    391          593
7     D             male   138    279          417
8     D           female   131    244          375
9     E             male    53    138          191
10    E           female    94    299          393

Problème n°5

Le dilemme du tramway (trolley problem) est une expérience de pensée qui permet d’étudier les déterminants des jugements de moralité (i.e., qu’est-ce qui fait qu’on juge une action comme morale, ou pas ?).

Sous une forme générale, ce dilemme consiste à poser la question suivante : si une personne peut effectuer un geste qui bénéficiera à un groupe de personnes A, mais, ce faisant, nuira à une personne B (seule); est-il moral pour la personne d’effectuer ce geste ?

Voir ce lien pour plus d’informations.

Problème n°5

Généralement, on fait lire des scénarios aux participants de l’étude, dans lesquels un individu doit prendre une décision dans une situation similaire à celle décrite à la slide précédente. Par exemple, imaginons que Denis ait le choix entre ne rien faire et laisser un train tuer cinq personnes, ou faire dérailler ce train mais tuer une personne. Ensuite, on demande aux participants de juger de la moralité de l’action choisie par Denis, sur une échelle de 1 à 7.

Des études antérieures ont montré que ces jugements de moralité sont grandement influencés par trois mécanismes de raisonnement inconscients :

Le principe d’action : un préjudice causé par une action est jugé moralement moins acceptable qu’un préjudice causé par omission.
Le principe d’intention : un préjudice causé comme étant le moyen vers un but est jugé moralement moins acceptable qu’un préjudice étant un effet secondaire (non désiré) d’un but.
Le principe de contact : un préjudice causé via contact physique est jugé moralement moins acceptable qu’un préjudice causé sans contact physique.

Problème n°5

Ce jeu de données comprend 12 colonnes et 9930 lignes, pour 331 individus. L’outcome qui nous intéresse est response, un entier pouvant aller de 1 à 7, qui indique à quel point il est permis (moralement) de réaliser l’action décrite dans le scénario correspondant, en fonction de l’âge (age) et genre (male) du participant (id).

On se demande comment les jugements d’acceptabilité sont influencés par les trois principes décrits slide précédente. Ces trois principes correspondent aux trois variables, action, intention, et contact (dummy-coded).

d5 <- open_data(morale)
head(d5)

  response     id age male action intention contact
1        4 96;434  14    0      0         0       1
2        3 96;434  14    0      0         0       1
3        4 96;434  14    0      0         0       1
4        3 96;434  14    0      0         1       1
5        3 96;434  14    0      0         1       1
6        3 96;434  14    0      0         1       1

Propositions de réponses

Réponse possible problème n°1

La taille de la mère a l’air “plus” prédictive de la taille d’un individu, et ce d’autant plus si cet individu est une femme…

d1 %>%
    gather(parent, parent.height, 3:4) %>%
    ggplot(aes(x = parent.height, y = height, colour = parent, fill = parent) ) +
    geom_point(pch = 21, size = 4, color = "white", alpha = 1) +
    stat_smooth(method = "lm", fullrange = TRUE) +
    facet_wrap(~ gender)

Réponse possible problème n°1

On peut fitter plusieurs modèles avec brms::brm(), et les comparer en utilisant le WAIC.

library(brms)

d1$gender <- ifelse(d1$gender == "F", -0.5, 0.5)
d1$mother <- scale(d1$mother) %>% as.numeric
d1$father <- scale(d1$father) %>% as.numeric

p1 <- c(
    prior(normal(70, 10), class = Intercept),
    prior(cauchy(0, 10), class = sigma)
    )

m1 <- brm(
    height ~ 1 + gender,
    prior = p1,
    data = d1
    )

p2 <- c(
    prior(normal(70, 10), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 10), class = sigma)
    )

m2 <- brm(
    height ~ 1 + gender + mother + father,
    prior = p2,
    data = d1
    )

Réponse possible problème n°1

p3 <- c(
    prior(normal(70, 10), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 10), class = sigma)
    )

m3 <- brm(
    height ~ 1 + gender + mother + father + gender:mother,
    prior = p3,
    data = d1
    )

p4 <- c(
    prior(normal(70, 10), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 10), class = sigma)
    )

m4 <- brm(
    height ~ 1 + gender + mother + father + gender:father,
    prior = p4,
    data = d1
    )

Réponse possible problème n°1

m1 <- add_criterion(m1, "waic")
m2 <- add_criterion(m2, "waic")
m3 <- add_criterion(m3, "waic")
m4 <- add_criterion(m4, "waic")

model_comparison_table <- loo_compare(m1, m2, m3, m4, criterion = "waic") %>%
  data.frame %>%
  rownames_to_column(var = "model")

weights <- data.frame(weight = model_weights(m1, m2, m3, m4, weights = "waic") ) %>%
  round(digits = 3) %>%
  rownames_to_column(var = "model")

left_join(model_comparison_table, weights, by = "model")

  model elpd_diff  se_diff  elpd_waic se_elpd_waic   p_waic se_p_waic     waic
1    m3  0.000000 0.000000  -93.24063     5.048439 5.317052  1.361108 186.4813
2    m2 -0.227792 1.424670  -93.46842     5.348427 4.860775  1.418229 186.9368
3    m4 -1.521233 1.646477  -94.76187     5.402762 5.941276  1.744878 189.5237
4    m1 -9.551493 4.784951 -102.79213     4.241509 2.641759  0.647044 205.5843
    se_waic weight
1 10.096877  0.496
2 10.696854  0.395
3 10.805523  0.108
4  8.483018  0.000

Réponse possible problème n°1

summary(m3)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 + gender + mother + father + gender:mother 
   Data: d1 (Number of observations: 40) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept        65.97      0.38    65.21    66.74 1.00     5277     2285
gender            3.57      0.75     2.07     5.09 1.00     5819     2990
mother            1.73      0.39     0.95     2.49 1.00     5693     3249
father            0.59      0.37    -0.13     1.31 1.00     5315     2508
gender:mother    -1.04      0.79    -2.68     0.53 1.00     5321     2607

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     2.35      0.29     1.86     2.98 1.00     4163     2755

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Réponse possible problème n°1

m3 %>%
    plot(
        pars = "^b_",
        combo = c("dens_overlay", "trace"), widths = c(1, 1.5),
        theme = theme_bw(base_size = 14, base_family = "Open Sans")
        )