 Insea 2025             Statistiques Bayésiennes

TP3: Modèles Bayésiens Hiérarchiques (I)

                 Author: Hicham Janati

Objectifs:

Découvrir la librairie PyMC
Implémenter les premiers modèles bayésiens et faire le diagnostic de convergence
Interpréter les résultats

# Installation si nécessaire
!pip install pymc arviz numpy pandas matplotlib seaborn ipywidgets

On importe les librariries et un crée un générateur aléatoire:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import pymc as pm
import arviz as az

# Configuration pour de meilleurs graphiques
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_context('notebook')

seed = 42
rng = np.random.default_rng(seed)

1. Modèle Poisson-Gamma hiérarchique

Dans le cadre de la modélisation du nombre de sinistres, il n’est pas pratique de considérer un $\lambda_i$ spécifique à chaque individu car les données individuelles contiennent très souvent très peu d’observations. Ici, on souhaite donc regrouper les assurés en utilisant leurs informations individuelles. Nous avons une variable age qui donne l’âge du conducteur en 4 catégories:

age = 0 (< 30 ans)
age = 1 (entre 30 et 50 ans)
age = 2 (entre 50 et 60 ans)
age = 3 (supérieur à 60 ans)

Pour tenir compte des différences entre les catégories, on modélise chaque catégorie par un taux de sinistre spécifique $\textcolor{red}{\lambda_j}$ avec $j\in \{0, 1, 2, 3\}$. Pour prendre en compte leur similarité, les $\textcolor{red}{\lambda_j}$ sont modélisés avec une loi a priori commune $\text{Gamma}(\textcolor{purple}{\alpha}, \textcolor{purple}{\beta})$.

Si des données historiques peuvent être utilisées, alors $\textcolor{purple}{\alpha}$ et $ $ sont choisis (constantes a priori) avec la méthode des moments (comme en TD1)
Sinon, on les modélise comme des variables aléatoires avec une loi apriori $\pi$ assez vague (grande variance, ou uniforme).

Le deuxième cas définit une structure bayésienne à deux niveaux: 1. Le nombre de sinistre $\textcolor{blue}{N}$ dépend de $\textcolor{red}{\lambda_j}$: $\textcolor{blue}{N} | \textcolor{red}{\lambda_j} \sim \mathcal{P}(\textcolor{red}{\lambda_j})$ 2. Le taux de sinistre $\textcolor{red}{\lambda_j}$ dépend de $ $ et $ : | , (, )$ avec $\textcolor{purple}{\alpha}$ et $\textcolor{purple}{\beta}$ et $\textcolor{purple}{\alpha}, \textcolor{purple}{\beta} \sim \pi$.

C’est un modèle hiérarchique. On considère une loi a priori Uniforme(0, 10).

Voici les données (extraites et filtrées à partir de https://www.kaggle.com/datasets/saisatish09/insuranceclaimsdata?select=dataCar.csv)

df = pd.read_csv("http://hichamjanati.github.io/data/claims_age.csv", index_col=0)
print(df.shape)
df.head()

Nous avons donc les données de 10205 assurés. On peut commencer par voir la taille de chaque groupe:

df.age.value_counts()

Question 1: On remarque que la catégorie d’âge 1 (30-40 ans) est la plus grande avec 4610 assurés. Celle des > 60 ans est la plus petite avec 1400 assurés. À quoi peut-on s’attendre concernant la qualité de l’estimation de chaque $\textcolor{red}{\lambda_j}$ ?

On regarde la distribution du nombre de sinistre déclarés par assuré:

df.numclaims.value_counts()

Question 2: On remarque que plus de 90% des assurés ne déclarent jamais de sinistres. Seulement 50/10250 ont déclaré 2 ou 3 sinistres. Quel est l’ordre de grandeur (ou fourchette de valeurs) des $\textcolor{red}{\lambda_j}$ auquel on peut s’attendre ?

Question 3: On note les nombres de sinistres de chaque groupe d’âge $j$ par $\textcolor{blue}{N}_1^j, \dots, \textcolor{blue}{N}_{n_j}^j | \textcolor{red}{\lambda_j} \sim \mathcal{P}(\textcolor{red}{\lambda_j})$. Ainsi, d’après la distribution ci-dessus des catégories d’âge: $n_0 = 2377$, $n_1 = 4610$ etc… Représentez le graphe probabiliste de ce modèle et déterminez la formule de la loi a posteriori jointe en fonction des lois de l’énoncé.

Question 4: Implémentez le modèle hiérarchique avec pymc et faites le diagnostic MCMC

Question 5: Calculez les bonnes probabilités de type $\mathbb P(\textcolor{red}{\lambda_j} < \textcolor{red}{\lambda_k})$ pour déterminer si certains groupes d’âge ont des risques différents ou non. Commenter

2. Bayesian Poisson regression

En plus de l’âge du conducteur, nous avons également la catégorie d’âge du véhicule (veh_age) et la valeur du véhicule en ‘$’ (veh_value) (divisée par 10000). Une structure bayésienne hiérarchique n’est plus possible (à moins de diviser la variable veh_value en catégories et créer tous les croisements de catégories age x veh_age x veh_value possibles chacune avec son taux $\textcolor{red}{\lambda_j}$). Une meilleur approche est de considérer une regression linéaire où on prédit le taux de sinistre avec une combinaison linéaire des variables: $\textcolor{red}{\lambda} = \textcolor{blue}{\beta_0} + \textcolor{blue}{\beta_1} \text{age} + \textcolor{blue}{\beta_2} \text{veh\_age} + \textcolor{blue}{\beta_3} \text{veh\_value}$. Or $\textcolor{red}{\lambda} > 0$, ce qui n’est pas respecté ici. On utilise un modèle linéaire généralisé où c’est $\log(\textcolor{red}{\lambda})$ qui est expliquée:

\[ N | \textcolor{red}{\lambda} \sim \mathcal P(\textcolor{red}{\lambda})\] \[\log(\textcolor{red}{\lambda}) = \textcolor{blue}{\beta_0} + \textcolor{blue}{\beta_1} \text{age} + \textcolor{blue}{\beta_2} \text{veh\_age} + \textcolor{blue}{\beta_3} \text{veh\_value}\]

Avec une loi a priori $\textcolor{blue}{\beta_0}, \textcolor{blue}{\beta_1}, \textcolor{blue}{\beta_2}, \textcolor{blue}{\beta_3} \sim \mathcal{N}(0, 1)$.

Question 6: On suppose que $\textcolor{blue}{\beta_1} = 0.2$. Si toutes les variables sauf l’âge ne changent pas, quel est l’effet de passer à une catégorie d’âge supérieur (càd que l’âge passe de 0 à 1, ou 1 à 2 ou 2 à 3) sur $\textcolor{red}{\lambda}$ ? Répondez à la question en terme de pourcentage de changement.

df = pd.read_csv("http://hichamjanati.github.io/data/claims_reg.csv", index_col=0)
df.head()

Question 7: Complétez le modèle bayésien ci-dessous et faites le diagonostic de convergence. En pratique, on définit $\textcolor{red}{\lambda}$ comme l’exp de la combinaison linéaire.

coords = dict(var_name=["intercept", "age", "veh_age", "veh_value"]) # dictionnaire qui sert à nommer les variables
with pm.Model(coords=coords) as reg_model:
    # beta = pm.Normal("beta", mu=0, sigma=1, dims="var_name") # vecteur des betas de taille 4 nommé selon "var_name" de coords
    # TO DO
    #
    #
    # lambda_ =
    lambda_ = pm.Gamma("lambda", 0.1, 0.1)
    N = pm.Poisson("N", mu=lambda_, observed=df["numclaims"])
    trace = pm.sample()

Question 8: Interprétez les valeurs et HDI obtenus pour chaque coefficient de regression $\beta_i$.

az.summary(trace)

Question 9: On souhaite vérifier que le modèle fit bien les données. Pour cela on peut utiliser les échantillons MCMC ($\beta$) pour générer des données $N_i$ (`pm.sample_posterior_predictive`) et comparer la vraisemblance avec la distribution des données générées. Qu’en pensez-vous ?

with reg_model:
    pm.sample_posterior_predictive(trace, extend_inferencedata=True) # ce paramètre = True fait qu'on modifie l'objet `trace` en rajoutant les samples de la predictive posterior
az.plot_ppc(trace)

Question 10 Now time to break it ! Pour bien cerner ce qui explique le bon fit des données de la question précédente, réduisez la complexité du modèle (le simplifier en enlevant des variables) jusqu’à ce que la fonction prédictive s’éloigne des données. Que peut-on en déduire ?

Question 11: Vus les résultats obtenus, il se peut que certains coefficients soient biaisés par le choix restrictif de l’a priori Gaussien avec variance égale à 1. On considère à présent $\sigma$ comme une variable aléatoire avec un a priori gaussienne positive (tronquée, `pm.HalfNormal`) avec un hyperparamètre $\sigma = 1$. Implémentez ce modèle. les résultats ont-ils changé considérablement ? Interpréter.

Question 12: On peut évaluer la qualité de ces modèles (et choisir le meilleur) avec le critère bayésien LOO (Leave-one-out).

LOO consiste à évaluer la log-vraisemblance de prédiction d’un échantillon $i$ après l’avoir enlevé des données, autrement dit, si on note $N_1, \dots, N_n$ les observations, alors $N_{-i}$ représente toutes les données sauf $N_i$, on note donc $N_{-i} = \{N_1, \dots, N_{i-1}, N_{i+1}, \dots N_n\}$. La fonction de prédiction (en log-probabilité) pour des données nouvelles est dit “expected log predictive density (ELPD)”: \[ ELPD_i = \log p(N_i | N_{-i}) \] Le critère LOO évalue la log-vraisemblance pour tous les échantillons (qui sont enlevés et prédits avec le reste à tour de rôle): \[ ELPD = \sum_{i=1}^n ELPD_i = \sum_{i=1}^n \log p(N_i | N_{-i}) \] Avec la loi des probabilités totale: \[ p(N_i | N_{-i}) = \int p(N_i | \beta, N_{-i}) p(\beta | N_{-i}) \mathrm d\beta = \int p(N_i | \beta) p(\beta | N_{-i}) \mathrm d\beta \]
Cette intégrale est approchée par Importance sampling (IS) avec la loi a posteriori full $p(\beta | N_1, \dots N_n)$ et une approximation des poids IS avec la loi de Pareto généralisée qui doit avoir des moments finies sinon IS est instable. Sa variance est finie si son paramètre (scale) $k < 0.5$. Sa moyenne est finie si $k < 1$. arviz nous donne l’estimation du ELPD ainsi que la qualité de l’estimation (k pour chaque $i$). Il faut d’abord calculer les loglikelihoods avec pm:

with reg_model:
    pm.compute_log_likelihood(trace)
az.plot_loo(trace)

Pour comparer tous les modèles vus dans ce notebook, on peut utiliser la fonction az.compare qui prend en argument un dictionnaire avec des noms des modèles en keys et les objets models en valeurs. Complétez ce code avec vos modèles et analysez, arviz trie les modèles par défaut du meilleur au pire:

models_dict = {"reg_model": reg_model, ...}
az.compare(models_dict)

TP3: Modèles Bayésiens Hiérarchiques (I)

Objectifs:

1. Modèle Poisson-Gamma hiérarchique

Question 1: On remarque que la catégorie d’âge 1 (30-40 ans) est la plus grande avec 4610 assurés. Celle des > 60 ans est la plus petite avec 1400 assurés. À quoi peut-on s’attendre concernant la qualité de l’estimation de chaque \(\textcolor{red}{\lambda_j}\) ?

Question 2: On remarque que plus de 90% des assurés ne déclarent jamais de sinistres. Seulement 50/10250 ont déclaré 2 ou 3 sinistres. Quel est l’ordre de grandeur (ou fourchette de valeurs) des \(\textcolor{red}{\lambda_j}\) auquel on peut s’attendre ?

Question 4: Implémentez le modèle hiérarchique avec pymc et faites le diagnostic MCMC

Question 5: Calculez les bonnes probabilités de type \(\mathbb P(\textcolor{red}{\lambda_j} < \textcolor{red}{\lambda_k})\) pour déterminer si certains groupes d’âge ont des risques différents ou non. Commenter

2. Bayesian Poisson regression

Question 7: Complétez le modèle bayésien ci-dessous et faites le diagonostic de convergence. En pratique, on définit \(\textcolor{red}{\lambda}\) comme l’exp de la combinaison linéaire.

Question 8: Interprétez les valeurs et HDI obtenus pour chaque coefficient de regression \(\beta_i\).

Question 10 Now time to break it ! Pour bien cerner ce qui explique le bon fit des données de la question précédente, réduisez la complexité du modèle (le simplifier en enlevant des variables) jusqu’à ce que la fonction prédictive s’éloigne des données. Que peut-on en déduire ?

Question 12: On peut évaluer la qualité de ces modèles (et choisir le meilleur) avec le critère bayésien LOO (Leave-one-out).