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.
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.
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}\]
où \(\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\) où \(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}\) où \(\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 :
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\).
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}\)
Calculer \(\hat{\theta}_{k + 1}\) à l’aide de la formule : \(\hat{\theta}_{k + 1} = (X^{'} W_kX)^{- 1}X^{'} W_k\:\mathbf{z_k}\)
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})\) où \(\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\) où \(\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}\) où \(\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}\) où \(\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é.
| 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.
| 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.
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.