# Base de R INLA
install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)TD2 - Introduction INLA
Ce TP est construit à partir des idées explorées durant la session FinistR2023
Introduction
Présentation de la méthodologie INLA (integrated nested Laplace approximation) en utilisant les packages {INLA}.
Ressources :
IMPORTANT : Utiliser le cluster 10 !!! INLA n’est pas installé sur les autres.
Installation de INLA
(Ici inutile car INLA est déja sur le cluster 10)
La réussite en maths
Chargement de packages
# Chargement des packages
library(tidyverse)
library(INLA)
library(lme4) # pour comparer avec l'approche frequentiste Le jeu de données awards contient le nombre de récompenses (num_awards) obtenues en fonction des résultats en mathématique (math) pour une classe de 200 élèves.
La réponse mesurée est un comptage.
# Chargement des donnees
load("data/awards.RData")
head(awards)- On souhaite voir la note de mathématique est liée au nombre de récompenses obtenues. Proposer une représentation graphique
#| eval: false
#| echo: false
#| label: awards_pres
awards |> ggplot() + aes(x= math, y = num_awards) + geom_point()- Proposer un modèle.
- Ajustement du modèle
#| eval: false
#| label: awards_first_model
# Application de la formule avec un modèle de Poisson
m.pois.lin1 <- inla(num_awards ~ math,
data = awards, family = "poisson")
summary(m.pois.lin1)
#Vous pouvez éventuellement comparer à ce que donneraitjags- Quelle est la distribution a posteriori pour l’effet des mathématiques
#| echo: true
#| label: awards_first_post
m.pois.lin1$marginals.fixed$math |>
ggplot() + aes(x=x, y = y) + geom_line() +
stat_function(fun = dnorm, args = list(mean = mean(m.pois.lin1$summary.fixed[2,1]), sd = m.pois.lin1$summary.fixed[2,2]), col = "red")- Comparer avec les résultats obtenus en utilisant un modèle linéaire généralisé estimé par max de vraisemblance
- Pour prendre en compte le problème classique de surdispersion dans les données de comptages, il peut être intéressant de rajouter un effet aléatoire de la manière suivante. Identifier tous les hyper paramètres et toutes les variables latentes.
\[Y_i\overset{ind}{\sim}\mathcal{P}(\exp(\mu+\alpha\cdot x_i+E_i)) \quad \text{avec}\quad E_i\overset{ind}{\sim}\mathcal{N}(0,\sigma^2).\]
#| echo: true
#| label: awards_ranef
awards <- awards |> mutate(ID = 1:nrow(awards))
m.pois.lin_surd <- inla(num_awards ~ math + f(ID, model = "iid"),
data = awards, family = "poisson")
summary(m.pois.lin_surd)Discuter les distributions a posteriori des effets fixes
on peut regarder la loi a posteriori de la précision sur les effets alétaoires
m.pois.lin_surd$summary.hyperpar
m.pois.lin_surd$marginals.hyperpar$`Precision for ID` |> ggplot() + aes(x=x, y = y) + geom_line()- Si une variable aléatoire V est strictement positive et admet une densité \(f\), quelle est la densité de 1/V ? En déduire la loi a posteriore de la variance de l’effet aléatoire.
m.pois.lin_surd$marginals.hyperpar$`Precision for ID` |> as.data.frame() |> mutate(u = 1/x, fu = x^2 ) |>
ggplot() + aes(x=u, y = fu) + geom_line()- On peut accéder aux effets aléatoires a posteriori, par exemple pour l’élève 1.
#| echo: true
#| label: awards_post_ranef
m.pois.lin_surd$marginals.random$ID$index.1
m.pois.lin_surd$marginals.random$ID$index.1 |> ggplot() + aes(x=x, y = y) + geom_line()- On peut comparer un modèle à effets aléatoires généralisés en fréquentiste. Que constatez vous ?
fit4.glmm<-glmer(num_awards ~ math + (1|id), family = poisson, data = awards)
summary(fit4.glmm)Quels sont les priors utilisés par inla pour les précisions ? Vous pourrez utiliser la commande
inla.doc("iid"). Quelle différence entre prior et postérior ?On souhaite changer le prior en une loi gamma de parametre 1,2 (Attention les priors sont codés dans le format interne d’inla cf Table 5.3 Gómez-Rubio (2020)
m.pois.lin_surd2 <- inla(num_awards ~ math +
f(ID, model = "iid",
hyper = list(prec = list(prior = "loggamma", param = c(3, 1)))),
data = awards, family = "poisson")
summary(m.pois.lin_surd2)- On peut demander à inla de calculer des critères de sélection de modèle type DIC ou WAIC
m.pois.lin_surd2 <- inla(num_awards ~ math +
f(ID, model = "iid",
hyper = list(prec = list(prior = "loggamma", param = c(3, 1)))),
control.compute = list(dic = TRUE, waic = TRUE),
data = awards, family = "poisson")
summary(m.pois.lin_surd2)- Intégrer un effet non linéaire de la note en mathématique
m.pois.nonlin <- inla(num_awards ~ f(math, model = "rw1"),
control.compute = list(dic = TRUE, waic = TRUE),
data = awards, family = "poisson")
summary(m.pois.nonlin)
awards |> mutate(pred_nonlin = m.pois.nonlin$summary.fitted.values[,1],
pred_lin = m.pois.lin1$summary.fitted.values[,1]) |>
pivot_longer(cols = c("pred_nonlin", "pred_lin"), names_to = "method", values_to = "pred") |>
arrange(math) |>
ggplot() + geom_line(aes(x=math, y = pred, col = method)) +
geom_point(aes(x=math, y = num_awards)) + ggtitle("Fit with linear and non linear effect (RW1)")- Il existe un modèle RW2, en vous référrant au chapitre Mixed-effects Models du livre Gómez-Rubio (2020), expliquez la différence entre RW1 et RW2. Dans cet exemple, quelle différence observez vous ?
m.pois.nonlin2 <- inla(num_awards ~ f(math, model = "rw2"),
control.compute = list(dic = TRUE, waic = TRUE),
data = awards, family = "poisson")
summary(m.pois.nonlin2)
awards |> mutate(pred_nonlin = m.pois.nonlin$summary.fitted.values[,1],
pred_nonlin2 = m.pois.nonlin2$summary.fitted.values[,1],
pred_lin = m.pois.lin1$summary.fitted.values[,1]
) |>
pivot_longer(cols = c("pred_nonlin","pred_nonlin2", "pred_lin"), names_to = "method", values_to = "pred") |>
arrange(math) |>
ggplot() + geom_line(aes(x=math, y = pred, col = method)) +
geom_point(aes(x=math, y = num_awards)) + ggtitle("Fit with linear and non linear effect (RW1)")References
Gómez-Rubio, Virgilio. 2020. Bayesian Inference with INLA. Chapman; Hall/CRC.