LinkMath #3 Données agrégées de mortalité

Cette note traite de l’agrégation des données de mortalité et de l’utilisation des données agrégées dans les modèles statistiques. Elle s’adresse à un public maîtrisant certains concepts de statistiques et théorie des probabilités (fonction de répartition, densité de probabilité) ainsi que les notations présentées dans la note LinkMath #1 et les résultats de la note LinkMath #2. Cette note sera particulièrement pertinente pour les personnes intéressées par les aspects de modélisation, qu’elles aient au préalable réalisé les traitements sur des données individuelles mentionnés dans la note précédente ou qu’elles travaillent directement sur la base de données de portefeuille ou de population déjà agrégées.

Le passage des données individuelles aux données agrégées, qui se fait en supposant une force de mortalité constante par morceaux sera d’abord justifié. La notion d’exposition centrale au risque, qui apparaît naturellement lorsqu’on procède à cette agrégation, sera en même temps introduite. L’estimateur du maximum vraisemblance de la force de mortalité sera ensuite présenté et le lien sera avec le cadre théorique des modèles linéaires généralisés établi. Enfin, les données agrégées de la Human Mortality Database fourniront un exemple à la fois utile et facile d’accès.

3.1 Agrégation des données

Agréger les données de mortalité signifie regrouper l’expérience associée à des individus de caractéristiques proches. Cela a pour but principal de réduire la taille de la base de données ce qui permettra d’optimiser l’estimation de modèles statistiques. L’agrégation des données produira également des indicateurs agrégés offrant une meilleure visualisation des données.

3.1.1 Représentation âge - temps calendaire

Pour illustrer la méthode d’agrégation des données dans un cadre concret, les caractéristiques de l’individu à la date de début d’observation seront supposées consister en un triplet \(\chi_i = (x_i, y_i, s_i)\) correspondant respectivement à l’âge de début d’observation, au temps calendaire de début d’observation et au sexe de l’individu. L’approche présentée est toutefois généralisable à n’importe quelle nombre de variables explicatives.

Le diagramme de Lexis permet de représenter les données de mortalité dans un repère utilisant le temps calendaire en abscisse et l’âge en ordonnée. Une trajectoire individuelle correspond ainsi à une diagonale dans ce repère. En effet chaque individu vieillit d’un an lorsqu’une année calendaire s’écoule. La Figure 3.1 présente dans ce nouveau repère les trajectoires individuelles utilisées comme exemple dans la note LinkMath #2. Cette représentation servira à illustrer le calcul de l’exposition au risque introduit dans la suite de cette section.

Trajectoires individuelles représentées dans un repère âge - temps calendaire

Figure 3.1: Trajectoires individuelles représentées dans un repère âge - temps calendaire

3.1.2 Vraisemblance associée aux données agrégées

Afin de procéder à l’agrégation des données, la force de mortalité sera supposée constante par morceaux entre deux âges entiers et deux années calendaires entières. Plus formellement, on écrira \(\mu_{x + \delta,y + \epsilon,s} = \mu_{x,y,s}\) pour tout couple \((x, y) \in \mathbb{N}\) et tout couple \((\delta, \epsilon) \in [0,1[\).

Sous ces conditions la vraisemblance se réécrit : \[\begin{align} \ell(\theta) = & \underset{i = 1}{\overset{n}{\sum}} \left\lbrace \delta_i\ln\mu_{x_i + t_i, y_i + t_i, s_i}(\theta) - \int_{u = 0}^{t_i} \mu_{x_i + u, y_i + u, s_i}(\theta)\text{d}u \right\rbrace \notag\\ = & \underset{i = 1}{\overset{n}{\sum}} \left\lbrace \left[\underset{x, y, s}{\sum} 1_{\left\lbrace x \le x_i + t_i < x + 1, y \le y_i + t_i < y + 1, s_i = s\right\rbrace}\right]\delta_i\ln\mu_{x_i + t_i, y_i + t_i, s_i}(\theta) \right. \\ & - \left. \int_{u = 0}^{t_i}\left[\underset{x,y,s}{\sum} 1_{\left\lbrace x \le x_i + u < x + 1, y \le y_i + u < y + 1, s_i = s\right\rbrace}\right]\mu_{x_i + u, y_i + u, s_i}(\theta)\text{d}u \right\rbrace \notag\\ = & \underset{x, y, s}{\sum} \left\lbrace \ln\mu_{x, y, s}(\theta) \underset{i = 1}{\overset{n}{\sum}} 1_{\left\lbrace x \le x_i + t_i < x + 1, y \le y_i + t_i < y + 1, s_i = s\right\rbrace} \delta_i \right. \\ & - \left. \mu_{x, y, s}(\theta) \underset{i = 1}{\overset{n}{\sum}} \int_{u = 0}^{t_i} 1_{\left\lbrace x \le x_i + u < x + 1, y \le y_i + u < y + 1, s_i = s\right\rbrace} \text{d}u \right\rbrace \notag\\ = & \underset{x,y,s}{\sum} \left\lbrace d_{x,y,s}\ln\mu_{x,y,s}(\theta) - \mu_{x,y,s}(\theta) e^c_{x,y,s}\right\rbrace \tag{3.1} \end{align}\]

en notant, pour des individus de sexe \(s\), d’âge compris entre \(x\) et \(x + 1\) et pour une date comprise entre \(y\) et \(y + 1\) :

  • \(d_{x,y,s} = \underset{i = 1}{\overset{n}{\sum}} 1_{\left\lbrace x \le x_i + t_i < x + 1, y \le y_i + t_i < y + 1, s_i = s\right\rbrace} \times \delta_i\) le nombre de décès observés

  • \(e^c_{x,y,s} = \underset{i = 1}{\overset{n}{\sum}} \int_{u = 0}^{t_i}1_{\left\lbrace x \le x_i + u < x + 1, y \le y_i + u < y + 1, s_i = s\right\rbrace}\ \text{d}u\) le nombre total d’années d’observations appelé également exposition centrale au risque.

Pour mener à bien ce calcul, deux sommes d’indicatrices portant sur un découpage de l’espace des variables explicative sont tout d’abord introduites dans l’expression de la vraisemblance. Ces deux sommes ont pour valeur 1 car une et une seule des indicatrices est non-nulle simultanément. La seconde étape du calcul est l’interversion de la somme sur les individus et de la somme sur les variables explicatives.

3.1.3 Notion d’exposition centrale au risque

L’exposition centrale au risque associée à une combinaison de variables explicatives correspond au nombre total d’années d’observations d’individus présentant ces variables explicatives. L’intégrale intervenant dans son expression de l’exposition centrale au risque peut être calculée explicitement. En effet :

\[\begin{align*} e^c_{x,y,s} & = \underset{i = 1}{\overset{n}{\sum}} \int_{u = 0}^{t_i}1_{\left\lbrace x \le x_i + u < x + 1, y \le y_i + u < y + 1, s_i = s\right\rbrace}\ \text{d}u \\ & = \underset{i = 1}{\overset{n}{\sum}} 1_{\left\lbrace s_i = s\right\rbrace} \left[\min(t_i, x + 1 - x_i, y + 1 - y_i) - \max(0, x - x_i, y - y_i)\right]^{+}. \end{align*}\]

Pour justifier ce résultat, remarquons que l’intégrande est non-nul si \(x \le x_i + u < x + 1\) et \(y \le y_i + u < y + 1\). En retranchant \(x_i\) (resp. \(y_i\)) à la première (resp. la deuxième) inégalité, on obtient deux nouvelles inégalités sur \(u\). La valeur de l’intégrale correspond par conséquent à la longueur de l’intersection des intervalles \([0,t_i]\) (donné par les bornes de l’intégrale), \([x - x_i, x + 1 - x_i]\) et \([y - y_i, y + 1 - y_i]\). On retrouve ainsi l’expression \(\left[\min(t_i, x + 1 - x_i, y + 1 - y_i) - \max(0, x - x_i, y - y_i)\right]^{+}\) proposée, où \(a^+\) vaut \(a\) si \(a\) est positif et 0 sinon et permet de traiter le cas où l’intersection est vide.

Comme illustré sur la Figure 3.2, la contribution de chaque individu à l’exposition centrale au risque pour la cellule \([x,x + 1] \times [y,y + 1]\) est proportionnelle à la longueur de l’intersection, représentée en bleu, de la trajectoire de l’individu avec la frontière de la cellule (le coefficient de proportionnalité vaut \(\frac{1}{\sqrt{2}}\)). Par ailleurs seuls les décès survenus à l’intérieur de la cellule, en rouge sur la Figure, sont comptabilisés.

Schéma explicatif du calcul de l'exposition au risque

Figure 3.2: Schéma explicatif du calcul de l’exposition au risque

On parle d’exposition centrale pour la distinguer de l’exposition initiale parfois utilisée pour désigner le nombre d’individus présents au début d’une année calendaire. L’exposition centrale ne dépend pas de la cause de sortie étudiée. Ainsi, si l’on voulait réaliser une étude des rachats au sein du portefeuille, l’exposition centrale serait la même que dans le cas des décès.

Dans l’Equation (3.1), le nombre de termes ne dépend plus du nombre d’individus mais uniquement du nombre de combinaisons possibles des variables explicatives. L’agrégation permet ainsi de condenser toute l’information disponible dans les vecteurs \(d_{x,y,s}\) et \(e^c_{x,y,s}\) de taille fixe.

3.1.4 Cas général

Afin d’illustrer plus intuitivement la notion d’exposition centrale au risque, les résultats précédents ont été présentés en prenant en compte 3 variables explicatives dont 2 échelles de temps. Dans un cadre plus général où la vraisemblance dépend d’un nombre quelconque de variables explicatives, l’Equation (3.1) se réécrit :

\[\begin{equation} \ell(\theta) = \underset{\chi}{\sum} \left\lbrace d_\chi\ln\mu_\chi(\theta) - \mu_\chi(\theta) e^c_\chi\right\rbrace \tag{3.2} \end{equation}\]

\(\chi\) désigne ici les combinaisons possibles de ces variables, en faisant l’hypothèse de force de mortalité constante par morceaux comme précédemment pour chaque échelle de temps.

L’hypothèse de force de mortalité constante par morceaux a été présentée ici en prenant des intervalles d’un an pour l’âge comme pour l’année calendaire, ce qui consiste à effectuer un découpage en carrés unitaires du plan âge \(\times\) année calendaire.

Il est en fait possible de construire des intervalles de n’importe quelle longueur - par exemple 1 jour, 1 semaine, 1 mois, 3 mois, 5 ans - voire de combiner des longueurs d’intervalle différentes. Ainsi dans la construction de lois de mortalité d’expérience pour le risque dépendance, où la durée passée dans l’état de dépendance est une variable explicative essentielle pour expliquer la mortalité dans cet état, il est commun de choisir des intervalles mensuels la première année puis annuels les années suivantes.

Plus généralement, le choix des intervalles doit permettre de concilier les deux objectifs adverses suivants :

  • Disposer d’assez d’observations pour chaque combinaison de modalités

  • Conserver un niveau de mortalité homogène sur chaque intervalle afin que l’hypothèse de force de mortalité constante reste pertinente.

Sous l’hypothèse que la force de mortalité est constante par morceaux sur des intervalles d’un an sur l’âge, l’année calendaire, ou d’autres variables explicatives, la vraisemblance associée à l’observation des individus prend une forme agrégée

\[\ell(\theta) = \underset{\chi}{\sum} \left\lbrace d_\chi\ln\mu_\chi(\theta) - \mu_\chi(\theta) e^c_\chi\right\rbrace\]

dont la somme porte maintenant sur les combinaisons possibles des variables explicatives et non plus sur les individus.

Les vecteurs \(d_\chi\) et \(e^c_\chi\) désignent respectivement le nombre de décès observé et le nombre d’années d’observations sommées sur tous les individus pour chaque combinaison \(\chi\) des variables explicatives.

3.2 Résolution des équations de vraisemblance

3.2.1 Dérivées dans le cas de modèles log-linéaires

L’estimateur du maximum de vraisemblance s’obtient en annulant la dérivée de la fonction de vraisemblance. A ce stade, il est plus pratique d’utiliser des notations algébriques. L’équation (3.2) se réécrit ainsi :

\[\ell(\theta) = \ln\boldsymbol{\mu}(\theta)^{'}\mathbf{d} - \boldsymbol{\mu}(\theta)^{'}\mathbf{e^c}\] en notant à l’aide de \(^{'}\) la transposée d’une matrice (ou ici d’un vecteur) et en utilisant des caractères en gras pour indiquer lorsqu’on parle d’un vecteur (sauf pour \(\theta\)).

La dérivée de la fonction de vraisemblance s’écrit :

\[\begin{equation} \frac{\partial \ell}{\partial\theta} = \left(\frac{\partial \boldsymbol{\mu}}{\partial\theta}\right)^{'}\left( \frac{\mathbf{d}}{\boldsymbol{\mu}(\theta)} - \mathbf{e^c} \right). \tag{3.3} \end{equation}\]

Dans le cas où la force de mortalité suit un modèle log-linéaire, \(\ln\mu(\theta) = X\theta\)\(X\) est appelé matrice de modèle. L’Equation (3.3) se réécrit alors :

\[\begin{equation} \frac{\partial \ell}{\partial\theta} = X^{'}\left(\mathbf{d} - \mathbf{e^{X\theta}} * \mathbf{e^c} \right) \tag{3.4} \end{equation}\]

Dans le modèle log-linéaire \(\mu(\theta) = e^{X\theta}\) et la force de mortalité ne peut ainsi prendre que des valeurs positives, ce qui cohérent avec la définition de cette quantité.

La dérivée seconde de la fonction de log-vraisemblance est quant à elle :

\[\frac{\partial^2 \ell}{\partial\theta^2} = - X^{'} W_\theta X\]

en notant \(W_\theta = \text{diag}({\mathbf{e^{X\theta}} * \mathbf{e^c}})\) la matrice diagonale dont les éléments sont les coordonnées du vecteur \(\mathbf{e^{X\theta}} * \mathbf{e^c}\) (\(*\) désigne le produit de Hadamard ou produit élément par élément de deux vecteurs ou matrices).

3.2.2 Résolution numérique : l’algorithme de Newton

L’estimateur du maximum de vraisemblance \(\hat{\theta}\) de \(\theta\) s’obtient en annulant la dérivée de la fonction de vraisemblance. Si l’équation (3.4) n’admet pas de solution explicite dans le cas général, il est néanmoins possible d’utiliser l’algorithme de Newton, appliqué à la dérivée de la fonction de vraisemblance, afin d’obtenir une suite d’estimateurs \((\hat{\theta}_k)_{k \ge 0}\) convergeant vers \(\hat{\theta}\).

Ceux-ci sont définis par :

\[\begin{align*} \hat{\theta}_{k + 1} &= \hat{\theta}_k + \left(\left.\frac{\partial^2 \ell}{\partial\theta^2}\right|_{\theta = \hat{\theta}_k}\right)^{- 1} \left.\frac{\partial \ell}{\partial\theta}\right|_{\theta = \hat{\theta}_k} \\ &= \hat{\theta}_k + (X^{'}W_kX)^{- 1} X^{'}\left(\mathbf{d} - \mathbf{e^{X\hat{\theta}_k}} * \mathbf{e^c} \right) \\ &= (X^{'} W_kX)^{- 1}X^{'} W_k \mathbf{z_k} \end{align*}\]

en notant \(W_k = W_{\hat{\theta}_k}\) et \(\mathbf{z_k} = X\hat{\theta}_k + \frac{\mathbf{d}}{\mathbf{\exp({X\hat{\theta}_k)}} \:*\:\mathbf{e^c}} - \mathbf{1}\).

L’algorithme de Newton peut être utilisé pour trouver le zéro d’une fonction \(f\) suffisamment régulière et dérivable. Il s’agit d’un algorithme récursif qui en partant d’un vecteur de paramètre \(\mathbf{x}_0\) construit une suite \(\mathbf{x}_k\) convergeant vers la valeur \(\mathbf{x}\) telle que \(f(\mathbf{x}) = 0\).

Le passage de \(\mathbf{x}_k\) à \(\mathbf{x}_{k + 1}\) se fait en approchant \(f\) par sa tangente en \(\mathbf{x}_k\) : \(f(\mathbf{x}) \simeq f(\mathbf{x}_k) + f^{'}(\mathbf{x}_k)(\mathbf{x} - \mathbf{x}_k)\)

En utilisant \(f(\mathbf{x}) = 0\), dans cette relation, nous obtenons l’expression du prochain terme de la suite : \(\mathbf{x}_{k + 1} = \mathbf{x}_k + f^{'-1}(\mathbf{x}_k)f(\mathbf{x}_k)\).

L’initialisation de cet algorithme peut être réalisée directement en remplaçant \(X\hat{\theta}_0\) par \(\ln\mathbf{d} - \ln\mathbf{e^c}\). Il n’est pas nécessaire de connaître la valeur de \(\hat{\theta}_0\) qui mène à ce résultat ni même qu’une telle valeur existe. Cela revient à considérer \(W_0\) dont les éléments diagonaux sont les décès observés \(\mathbf{d}\) et \(\mathbf{z}_0 = \ln\mathbf{d} - \ln\mathbf{e^c}\).

Dans le cas où \(\mathbf{d}\)\(\mathbf{e^c}\) présentent des valeurs nulles, \(\ln\mathbf{d}\) et \(\ln\mathbf{e^c}\) ne sont pas définis et pourront être remplacés par une valeur négative par exemple \(- 20\) dans l’expression de \(\mathbf{z}_0\).

L’estimateur du maximum de vraisemblance \(\hat{\theta}\) peut être obtenu numériquement de la manière suivante :

  1. Initialisation : Poser \(W_0 = \text{diag}(\mathbf{d})\), \(\mathbf{z}_0 = \ln\mathbf{d} - \ln\mathbf{e^c}\). Fixer \(k = 1\) et calculer \(\hat{\theta}_1 = (X^{'} W_0X)^{- 1}X^{'} W_0\: \mathbf{z}_0\).

  2. Calculer \(W_k = W_{\hat{\theta}_k}\) et \(\mathbf{z_k} = X\hat{\theta}_k + \frac{\mathbf{d}}{\mathbf{e^{X\hat{\theta}_k}} *\:\mathbf{e^c}} - \mathbf{1}\)

  3. Calculer \(\hat{\theta}_{k + 1}\) à l’aide de la formule : \(\hat{\theta}_{k + 1} = (X^{'} W_kX)^{- 1}X^{'} W_k\:\mathbf{z_k}\)

  4. Si \(\|\hat{\theta}_{k + 1} - \hat{\theta}_k\|^2 < \epsilon\) avec par exemple \(\epsilon = 10^{- 6}\), arrêter l’algorithme. Sinon incrémenter \(k\) et retourner à l’étape 2.

3.2.3 Intervalles de confiance asymptotiques

Un estimateur \(\widehat{\text{var}}(\hat{\theta})\) de la variance de l’estimateur \(\hat{\theta}\) obtenu à la convergence de l’algorithme précédent est donné par :

\[\widehat{\text{var}}(\hat{\theta}) \simeq X(X^{'} W_{\hat{\theta}} X)^{- 1}\]

Il est possible de montrer qu’asymptotiquement, lorsque le nombre d’observation devient très grand, la distribution de \(\hat{\theta}\) tend vers \(\mathcal{N}(\theta_0, \mathcal{I}^{- 1})\)\(\theta_0\) désigne la vraie valeur (inconnue) du paramètre \(\theta_0\) et

\[\mathcal{I} = - \mathbb{E}\left(\left.\frac{\partial^2 \ell}{\partial\theta^2}\right|_{\theta = \theta_0}\right)_{} = X^{'} \mathbb{E}(W_{\theta_0}) X.\] La quantité \(\mathcal{I}\) est appelée matrice d’information de Fisher. Celle-ci n’est en général pas connue mais peut en pratique être approchée par la matrice d’information empirique \(\hat{\mathcal{I}} = X^{'} W_{\hat{\theta}} X\). C’est donc celle-ci qui sera utilisée pour obtenir un estimateur de la variance de \(\hat{\theta}\).

Par bilinéarité de la variance \(\text{var}(X\hat{\theta}) = X\text{var}(\hat{\theta})X^{'}\). S’en déduit un estimateur de la variance associée à \(\ln \hat{\boldsymbol{\mu}} = X\hat{\theta}\) :

\[\widehat{\text{var}}(X\hat{\theta}) = X(X^{'} W_{\hat{\theta}} X)^{- 1}X^{'}\]

Cela permet d’obtenir l’intervalle de confiance (asymptotique) suivant :

\[\ln\boldsymbol{\mu} \in \left[X\hat{\theta} \pm \Phi^{- 1}\left(\frac{1 - \alpha}{2}\right) \sqrt{\text{diag}(X(X^{'} W_{\hat{\theta}} X)^{- 1}X^{'})}\right]\] avec un niveau de confiance \(\alpha\), où \(\Phi\) désigne la fonction de répartition de la loi normale standard.

Un intervalle de confiance sur \(\boldsymbol{\mu}\) peut être obtenu simplement en prenant l’exponentielle des deux bornes de l’intervalle de confiance sur \(\ln\boldsymbol{\mu}\).

3.2.4 Estimateur des taux bruts

Un cas particulier intéressant est celui du modèle saturé qui inclut autant de paramètres que d’observations. Le vecteur de paramètres est alors directement \(\theta = \ln\mu\) et la matrice de modèle \(X\) associée est la matrice identité.

En reprenant l’Equation (3.4), on peut alors calculer explicitement \(\boldsymbol{\mu}\) :

\[0 = \frac{\partial \ell}{\partial\theta} = X^{'}\left(\mathbf{d} - \mathbf{e^{X\theta}} * \mathbf{e^c} \right) = \left(\mathbf{d} - \boldsymbol{\mu} * \mathbf{e^c} \right) \quad \text{et} \quad \boldsymbol{\mu} = \frac{\mathbf{d}}{\mathbf{\:e^c}}.\]

L’estimateur de la force de mortalité associé \(\hat{\boldsymbol{\mu}}\) est dans ce cas appelé taux brut de mortalité.

L’estimateur de la variance précédent se simplifie ici en \(\widehat{\text{var}}(\ln\hat{\boldsymbol{\mu}}) = W_{\hat{\theta}}^{- 1}\) où les éléments diagonaux de \(W_{\hat{\theta}}^{- 1}\) sont les coordonnées du vecteur \(\hat{\boldsymbol{\mu}} * \mathbf{e^c}\) ce qui permet d’obtenir l’intervalle de confiance asymptotique suivant :

\[\ln\boldsymbol{\mu} \in \left[\ln\frac{\mathbf{d}}{\mathbf{e^c}} \pm \Phi^{- 1}\left(\frac{1 - \alpha}{2}\right) \frac{1}{\sqrt{\mathbf{d}}}\right]\] avec un niveau de confiance \(\alpha\), où \(\Phi\) désigne là encore la fonction de répartition de la loi normale standard.

3.2.5 Lien avec le GLM Poisson

Il est possible d’opérer le lien entre le modèle précédent et le cadre des modèles linéaires généralisés (GLM en anglais). Pour ce faire, on suppose que \(\mathbf{D}\), variable aléatoire dont \(\mathbf{d}\) est une réalisation, suit une loi de Poisson et que le paramètre de cette loi, qui représente le nombre moyen de décès observés, est proportionnel d’une part à la force de mortalité \(\boldsymbol{\mu}(\theta)\) et d’autre part au nombre total d’années d’observations \(\mathbf{e^c}\). Formellement, on écrit : \(\mathbf{D} \sim \mathcal{P}\left(\boldsymbol{\mu}(\theta) \mathbf{e^c}\right)\).

La loi de Poisson est une loi discrète ne prenant que des valeurs entières positives. Elle ne dépend que d’un unique paramètre, qui représente à la fois la moyenne et la variance de la loi. Cette loi est un choix naturel pour modéliser un nombre d’évènements observés.

La vraisemblance associée à ce modèle s’écrit :

\[\begin{align} \small \ell^{'}(\theta) = & \ln \left\lbrace \underset{\chi}{\prod} \mathbb{P}(D_\chi = d_\chi) \right\rbrace = \ln \left\lbrace \underset{\chi}{\prod} e^{- \mu_\chi(\theta) e^c_\chi} \frac{(\mu_\chi(\theta) e^c_\chi)^{d_\chi}}{d_\chi!}\right\rbrace \\ = & \underset{\chi}{\sum} \left\lbrace - \mu_\chi(\theta) e^c_\chi + d_\chi \ln\mu_\chi(\theta) + d_\chi \ln e^c_\chi + \ln(d_\chi!)\right\rbrace \tag{3.5} \end{align}\]

Cette somme comprend 4 termes dont les 2 derniers ne dépendent pas de \(\theta\) et donc n’interviennent pas dans la maximisation de la vraisemblance. Les 2 premiers termes coïncident quant à eux avec ceux de l’Equation (3.2). Cela signifie que les estimateurs de \(\theta\) obtenus en maximisant les expressions (3.2) et (3.5) sont identiques.

L’hypothèse d’une distribution de Poisson permet de se ramener au cadre des modèles linéaires généralisés (GLM). Pour rappel, un GLM est défini par trois éléments :

  • La variable observée : ici le nombre de décès \(\mathbf{D}\), supposé suivre une loi de Poisson.

  • La partie modélisée : dans le cadre des modèles linéaires, elle s’écrira en notation matricielle \(X\theta\)\(\theta\) est le vecteur des paramètres à estimer et \(X\) est la matrice de modèle qui décrit le rôle des paramètres dans le modèle.

  • La fonction de lien : elle permet de faire le lien entre l’espérance de la variable observée et la partie modélisée. Ici la fonction \(\ln\), fonction de lien canonique pour la loi de Poisson, sera naturellement utilisée.

On peut alors écrire \(\ln\mathbb{E}(\mathbf{D}) = \ln(\boldsymbol{\mu}(\theta) \mathbf{e^c}) = \ln\boldsymbol{\mu}(\theta) + \ln \mathbf{e^c}\)\(\ln\boldsymbol{\mu}(\theta) = X\theta\) correspondra à la partie modélisée. Par ailleurs dans le cadre de la loi de Poisson \(\text{var}(\mathbf{D}) = \mathbb{E}(\mathbf{D})\).

Le terme \(\ln \mathbf{e^c}\), appelé offset, peut être vu comme un décalage appliqué au vecteur d’observation ou comme un vecteur de paramètres supplémentaires du modèle dont la valeur serait déjà connue.

Le lien entre ce modèle et le précédent permet d’interpréter la force de mortalité \(\boldsymbol{\mu}\) en tant que coefficient de proportionnalité entre le nombre de décès prédit \(\mathbb{E}(\mathbf{D})\) et l’exposition centrale au risque observée \(\mathbf{e^c}\).

Si l’on suppose que le nombre de décès observés \(\mathbf{D}\) suit une loi de Poisson de paramètre \(\boldsymbol{\mu}(\theta) \mathbf{e^c}\)\(\boldsymbol{\mu}(\theta)\) est la force de mortalité et \(\mathbf{e^c}\) l’exposition centrale au risque, on retombe sur le même estimateur que pour l’Equation (3.2).

3.3 Données de population

Cette ultime section introduit à travers les données de la Human Mortality Database un exemple de donnnées agrégées.

3.3.1 Aperçu des données HMD

En France, les données de mortalité sont récoltées et traitées par l’INSEE. Au niveau international, une initiative portée par l’Institut Max Planck pour la recherche démographique et l’Université de Berkeley a permis de créer la Human Mortality Database, une base de données internationale rassemblant les données de plus de 40 pays. Celles-ci, disponibles directement dans un format agrégé, contiennent le nombre de décès observés ainsi que l’exposition centrale au risque pour chaque combinaison des variables suivantes :

  • Pays : la base regroupe les données de plus de 40 pays. Les pays occidentaux, les plus susceptibles de disposer de données historiques de qualité, y sont sur-représentés.

  • Age : avec des tranches d’âge d’un an, entre 0 et 110 ans, les plus de 110 ans étant regroupés dans une même tranche d’âge 110+.

  • Année calendaire par tranche d’un an. Les pays disposent d’un historique plus ou moins long. Ainsi les données pour la Suède remontent à 1751 et pour la France à 1816. A l’opposé, les données pour la Croatie et la Corée ne sont disponibles qu’à partir de 2001 et 2003 respectivement. D’autre part, pour la plupart des pays, il s’écoule entre 6 mois à 3 ans après la fin de la période considérée avant que les données les plus récentes ne soient ajoutées à la base.

  • Sexe : homme ou femme.

La Table 3.1 recense les pays pour lesquels les données sont disponibles, du plus peuplé au moins peuplé.

Table 3.1: Pays dont les données sont disponibles dans la base HMD
Pays Début Fin Population dernière année Décès dernière année
U.S.A. 1933 2019 327 610 213 2 854 838
Russia 1959 2014 143 829 248 1 878 039
Japan 1947 2019 123 834 559 1 381 093
Germany 1956 2017 82 654 471 932 263
U.K. 1922 2018 66 431 622 616 014
France 1816 2018 64 775 395 596 552
Italy 1872 2018 59 873 497 629 345
Republic of Korea 2003 2018 51 309 885 298 820
Spain 1908 2018 46 793 921 425 153
Ukraine 1959 2013 45 306 382 662 368
Poland 1958 2019 38 400 047 409 709
Canada 1921 2019 37 523 508 284 082
Australia 1921 2018 24 998 942 158 494
Taiwan 1970 2019 23 596 830 175 546
Chile 1992 2017 17 553 425 106 388
Netherlands 1850 2019 17 345 093 151 885
Belgium 1841 2020 11 459 028 126 896
Greece 1981 2017 10 756 216 124 501
Czechia 1950 2019 10 672 699 112 362
Sweden 1751 2020 10 355 368 98 124
Portugal 1940 2020 10 298 854 123 358
Hungary 1950 2017 9 787 062 131 674
Belarus 1959 2018 9 484 490 120 053
Austria 1947 2019 8 879 702 83 386
Israel 1983 2016 8 546 145 44 203
Switzerland 1876 2018 8 514 109 67 088
Hong Kong 1986 2019 7 481 456 48 957
Bulgaria 1947 2017 7 075 664 109 791
Denmark 1835 2020 5 831 529 54 645
Finland 1878 2020 5 529 564 55 348
Slovakia 1950 2017 5 439 192 53 914
Norway 1846 2020 5 380 194 40 611
Ireland 1950 2017 4 808 379 30 484
New Zealand 1948 2013 4 450 370 29 528
Croatia 2001 2019 4 066 904 51 794
Lithuania 1959 2019 2 794 758 38 281
Slovenia 1983 2017 2 066 240 20 509
Latvia 1959 2019 1 913 855 27 719
Estonia 1959 2019 1 326 978 15 401
Luxembourg 1960 2019 619 929 4 283
Iceland 1838 2018 352 683 2 254

La Table 3.2 présente quant à lui un extrait de la base. Les colonnes E et D correspondent respectivement à l’exposition centrale et au nombre de décès observés pour chaque combinaison des variables explicatives.

Table 3.2: Extrait des données HMD
Pays Annee Age Sexe E D
U.S.A. 2005 68 Homme 910 654 20 146
Germany 1987 49 Femme 584 214 1 624
U.K. 1998 54 Homme 346 730 2 283
Bulgaria 2006 55 Homme 54 568 843
U.S.A. 2008 69 Femme 1 094 622 16 774
Austria 2000 28 Femme 58 781 21
Belgium 1850 27 Femme 36 351 347
U.K. 1976 57 Homme 283 977 3 991
Ukraine 2003 35 Homme 310 758 2 141
U.S.A. 1965 5 Femme 2 012 717 935
U.S.A. 1965 48 Homme 1 083 273 8 848
U.S.A. 1959 29 Femme 1 161 240 1 140
Poland 2013 41 Homme 261 405 757
Germany 1984 42 Femme 554 425 907
Canada 1981 8 Homme 183 805 55
U.K. 1946 0 Femme 399 972 17 146
France 1838 38 Femme 234 199 2 668
Austria 1948 76 Homme 11 607 1 004
Germany 1979 27 Homme 564 982 790
Poland 1994 52 Femme 175 414 772

3.3.2 Données HMD pour la France

Les données HMD pour la France couvrent la période allant de 1816 à 2018. Les vecteurs de décès et d’exposition centrale par âge, année calendaire et sexe associés à ces données sont représentés sur la Figure 3.3, sous forme de pyramide des âges. L’exposition centrale s’interprète ici comme un nombre moyen de personnes en vie au cours de l’année pour chaque âge et chaque sexe. La représentation sous forme de pyramide est donc naturelle. Les décès sont par contre beaucoup moins fréquemment présentés sous ce format.

Pour l’exposition centrale deux creux qui correspondent à une chute de la natalité durant les deux guerres mondiales peuvent être observés. Le second creux est de surcroît accentué par le contraste avec la période du baby-boom allant de 1947 à 1973 qui voit la natalité exploser.

Décès observés et exposition centrale par âge et par sexe pour la France

Figure 3.3: Décès observés et exposition centrale par âge et par sexe pour la France

Pour certaines données agrégées de population ou issues de portefeuilles d’assurance, il peut arriver que l’exposition centrale au risque \(e^c_{x,y,s}\) ne soit pas disponible mais qu’on dispose à la place de l’effectif par âge et sexe dans le portefeuille au début de chaque année \(E_{x,y,s}\). Dans ce cas, l’approximation \(e^c_{x,y,s} \simeq \frac{E_{x,y,s} + E_{x,y + 1,s}}{2}\) peut être utilisée pour se ramener au cadre présenté dans cette note.

L’agrégation des données individuelles issues de portefeuille d’assurance peut se faire en supposant que la force de mortalité est constante par morceaux entre deux âges entiers et au cours d’une année calendaire. Cela permet de réécrire la vraisemblance associée au modèle en passant d’une somme sur les individus à une somme sur les combinaisons possibles des variables explicatives. L’expression de la vraisemblance ainsi obtenue fait intervenir le vecteur du nombre de décès observé et celui de l’exposition centrale au risque. Cette dernière quantité correspond au nombre d’années d’observation totalisées par l’ensemble des individus du portefeuille.

En faisant l’hypothèse que le nombre de décès suit une loi de Poisson dont le paramètre est proportionnel à la force de mortalité et à l’exposition au risque, on obtient un second estimateur de la vraisemblance qui donne le même résultat que le précédent mais relève du cadre des GLM. Les avantages associés sont un algorithme d’estimation performant, et la disponibilité d’intervalles de confiance ou de tests d’hypothèses associés.

Votre interlocuteur R&D :

Guillaume Biessy

Actuaire certifié IA

Docteur en Mathématiques Appliquées de l’Université Paris-Saclay

Professeur associé à temps partiel à Sorbonne Université

guillaume.biessy@linkpact.fr