3 Données agrégées de mortalité
Cette note traite de l’agrégation des données de mortalité et de l’utilisation de celles-ci 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 données individuelles mentionnés dans la précédente note ou qu’elles travaillent directement sur des 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 motivé. 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 et leurs extensions établi. Enfin, les données agrégées de la Human Mortality Database fourniront un exemple de données agrégées utile et en libre 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 d’échelles de temps ou 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.
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{aligned}\ell(\boldsymbol{\theta}) = & \underset{i = 1}{\overset{n}{\sum}} \left\lbrace \delta_i\ln\mu_{x_i + t_i, y_i + t_i, s_i}(\boldsymbol{\theta}) - \int_{u = 0}^{t_i} \mu_{x_i + u, y_i + u, s_i}(\boldsymbol{\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}(\boldsymbol{\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}(\boldsymbol{\theta})\text{d}u \right\rbrace \notag\\ = & \underset{x, y, s}{\sum} \left\lbrace \ln\mu_{x, y, s}(\boldsymbol{\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}(\boldsymbol{\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}(\boldsymbol{\theta}) - \mu_{x,y,s}(\boldsymbol{\theta}) e^c_{x,y,s}\right\rbrace\end{aligned} \tag{3.1}\]
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{aligned} 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{aligned}\]
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.
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 présent.
Dans l’Équation 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’Équation 3.1 se réécrit :
\[\ell(\boldsymbol{\theta}) = \underset{\chi}{\sum} \left\lbrace d_\chi\ln\mu_\chi(\boldsymbol{\theta}) - \mu_\chi(\boldsymbol{\theta}) e^c_\chi\right\rbrace \tag{3.2}\]
où \(\chi\) désigne ici toutes les combinaisons possibles des variables discrétisées, 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(\boldsymbol{\theta}) = \underset{\chi}{\sum} \left\lbrace d_\chi\ln\mu_\chi(\boldsymbol{\theta}) - \mu_\chi(\boldsymbol{\theta}) e^c_\chi\right\rbrace\]
dont la somme porte maintenant sur les combinaisons possibles des variables explicatives discrétisées 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 de passer aux notations algébriques. Nous désignerons dans ce qui suit les matrices par des lettres majuscules et les vecteurs par des caractères gras.
L’équation Équation 3.2 se réécrit ainsi :
\[\ell(\boldsymbol{\theta}) = \ln\boldsymbol{\mu}(\boldsymbol{\theta})^T\mathbf{d} - \boldsymbol{\mu}(\boldsymbol{\theta})^T\mathbf{e^c} \tag{3.3}\]
en notant à l’aide de \(A^T\) la transposée d’une matrice \(A\) (qui peut être comme ici un vecteur). On se place ici dans le cas où le logarithme de la force de mortalité s’exprime comme une fonction linéaire des paramètres du modèle \(\ln\mu(\boldsymbol{\theta}) = X\boldsymbol{\theta}\) où \(X\) est appelé matrice de modèle.
Les modèles linéaires constituent une classe de modèle beaucoup plus large que les ajustements affines avec lesquelles ils ne doivent pas être confondus. Les polynômes, les splines - sommes de polynômes par morceaux - ainsi que les fonctions échelons, qui permettent de modéliser des discontinuités, font partie des modèles linéaires.
La dérivée de la fonction de vraisemblance s’écrit :
\[\frac{\partial \ell}{\partial\boldsymbol{\theta}} = X^T\left(\mathbf{d} - \mathbf{e^{X\boldsymbol{\theta}}} * \mathbf{e^c} \right) \tag{3.4}\]
Dans le modèle log-linéaire \(\mu(\boldsymbol{\theta}) = e^{X\boldsymbol{\theta}}\) et la force de mortalité ne peut ainsi prendre que des valeurs positives, en cohérence 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\boldsymbol{\theta}^2} = - X^T W_\boldsymbol{\theta} X\]
en notant \(W_\boldsymbol{\theta} = \text{Diag}({\mathbf{e^{X\boldsymbol{\theta}}} * \mathbf{e^c}})\) la matrice diagonale dont la diagonale correspond au vecteur \(\mathbf{e^{X\boldsymbol{\theta}}} * \mathbf{e^c}\) (\(*\) désigne ici le produit de Hadamard ou produit élément par élément de deux vecteurs ou matrices).
3.2.2 Intervalles de confiance asymptotiques
Les propriétés de l’estimateur du maximum de vraisemblance font que la matrice de variance-covariance de l’estimateur \(\hat{\boldsymbol{\theta}}\) peut être asymptotiquement estimée par :
\[\widehat{\text{Var}}(\hat{\boldsymbol{\theta}}) \simeq (X^{'} W_{\hat{\boldsymbol{\theta}}} X)^{- 1}\]
Il est possible de montrer qu’asymptotiquement, lorsque le nombre d’observation devient très grand, la distribution de \(\hat{\boldsymbol{\theta}}\) tend vers \(\mathcal{N}(\boldsymbol{\theta}_0, \mathcal{I}^{- 1})\) où \(\boldsymbol{\theta}_0\) désigne la vraie valeur (inconnue) du paramètre \(\boldsymbol{\theta}_0\) et
\[\mathcal{I} = - \mathbb{E}\left(\left.\frac{\partial^2 \ell}{\partial\boldsymbol{\theta}^2}\right|_{\boldsymbol{\theta} = \boldsymbol{\theta}_0}\right)_{} = X^{'} \mathbb{E}(W_{\boldsymbol{\theta}_0}) X.\]
La quantité \(\mathcal{I}\) est appelée matrice d’information de Fisher. Celle-ci est en général inconnue, car dépendant du paramètre \(\boldsymbol{\theta}_0\) inconnu, mais peut en pratique être approchée par la matrice d’information empirique \(\hat{\mathcal{I}} = X^{'} W_{\hat{\boldsymbol{\theta}}} X\) qui se calcule quant à elle directement à partir des données. C’est donc celle-ci qui sera utilisée pour obtenir un estimateur de la variance de \(\hat{\boldsymbol{\theta}}\).
Par bilinéarité de la variance \(\text{Var}(X\hat{\boldsymbol{\theta}}) = X\text{Var}(\hat{\boldsymbol{\theta}})X^{'}\), un résultat connu également sous le nom de loi de propagation des erreurs. S’en déduit un estimateur de la variance associée à \(\ln \hat{\boldsymbol{\mu}} = X\hat{\boldsymbol{\theta}}\) :
\[\widehat{\text{Var}}(X\hat{\boldsymbol{\theta}}) = X(X^{'} W_{\hat{\boldsymbol{\theta}}} X)^{- 1}X^{'}\]
Cela permet d’obtenir l’intervalle de confiance (asymptotique) suivant :
\[\ln\boldsymbol{\mu} \in \left[X\hat{\boldsymbol{\theta}} \pm \Phi^{- 1}\left(\frac{1 - \alpha}{2}\right) \sqrt{\text{diag}(X(X^{'} W_{\hat{\boldsymbol{\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.3 Estimateur des taux bruts
Un cas particulier intéressant est celui du modèle saturé qui inclut autant de paramètres que d’observations. La matrice de modèle \(X\) correspondante est la matrice identité et le vecteur de paramètres est alors directement \(\boldsymbol{\theta} = \ln\mu\).
En reprenant l’Équation 3.4, on peut alors calculer explicitement \(\boldsymbol{\mu}\) :
\[0 = \frac{\partial \ell}{\partial\boldsymbol{\theta}} = X^{'}\left(\mathbf{d} - \mathbf{e^{X\boldsymbol{\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é.
Comme \(\boldsymbol{\mu}*\mathbf{e^c} = \mathbf{d}\), la matrice de variance-covariance associée se simplifie dans ce cas en \(\widehat{\text{Var}}(\ln\hat{\boldsymbol{\mu}}) = \text{Diag}(\mathbf{1} / \mathbf{d})\) 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.
L’estimateur des taux bruts possède un certain nombre de propriétés intéressantes :
Il s’agit de l’estimateur le plus précis parmi les estimateurs se basant sur les données discrétisées. Il convient de noter que cette précision s’entend par rapport aux données d’origine et traduit en grande partie du surapprentissage.
Il s’agit de l’estimateur possédant le plus grand nombre de paramètres indépendants parmi les estimateurs se basant sur les données discrétisées.
Il s’agit de l’estimateur possédant les intervalles de confiance les plus larges.
L’estimateur des taux bruts permet également de passer du cadre des modèles de durées à celui de la régression pondérée. En effet, il a été montré que le logarithme des taux bruts suit asymptotiquement une loi normale centrée sur le logarithme de la force de mortalité :
\[\ln\frac{\mathbf{d}}{\mathbf{e^c}} \sim \mathcal{N}(\ln\boldsymbol{\mu}, \text{Diag}(\mathbf{1} / \mathbf{d}))\] Cela suggère ainsi d’utiliser \(y = \ln\mathbf{d} - \ln\mathbf{e^c}\) comme vecteur d’observation et \(w = d\) comme vecteur de poids afin de construire un estimateur et des intervalles de confiance pour \(\ln\boldsymbol{\mu}\) puis pour \(\boldsymbol{\mu}\) en passant à l’exponentielle. Il convient toutefois de noter que cette approximation normale n’est valable qu’asymptotiquement. Pour un échantillon de taille fini, elle induit un biais qui peut être très conséquent. En effet, \(y\) et \(w\) sont des fonctions croissantes du nombre de décès \(d\) et l’estimateur construit sur la base de l’approximation normale aura ainsi tendance à donner plus de poids aux observations présentant un nombre de décès élevés et donc à surestimer le taux de mortalité. On évitera donc de recourir à cette approximation à moins de travailler sur des populations présentant un volume de données très conséquent.
Supposer que le taux brut, et non son logarithme, suit une loi normale, est fréquent dans les études actuarielles et conduit parfois à des taux estimés négatifs, ce qui peut aisément être évité.
3.2.4 Résolution numérique : l’algorithme de Newton
L’estimateur du maximum de vraisemblance \(\hat{\boldsymbol{\theta}}\) de \(\boldsymbol{\theta}\) s’obtient en annulant la dérivée de la fonction de vraisemblance. Si l’équation É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{\boldsymbol{\theta}}_k)_{k \ge 0}\) convergeant vers \(\hat{\boldsymbol{\theta}}\).
Ceux-ci sont définis récursivement par :
\[\begin{aligned} \hat{\boldsymbol{\theta}}_{k + 1} &= \hat{\boldsymbol{\theta}}_k - \left(\left.\frac{\partial^2 \ell}{\partial\boldsymbol{\theta}^2}\right|_{\boldsymbol{\theta} = \hat{\boldsymbol{\theta}}_k}\right)^{- 1} \left.\frac{\partial \ell}{\partial\boldsymbol{\theta}}\right|_{\boldsymbol{\theta} = \hat{\boldsymbol{\theta}}_k} \\ &= \hat{\boldsymbol{\theta}}_k + (X^{'}W_kX)^{- 1} X^{'}\left(\mathbf{d} - \mathbf{e^{X\hat{\boldsymbol{\theta}}_k}} * \mathbf{e^c} \right) \\ &= (X^{'} W_kX)^{- 1}X^{'} W_k \mathbf{z_k} \end{aligned}\]
en notant \(W_k = \text{Diag}(\exp(X\boldsymbol{\theta}_k) * \mathbf{e_c})\) et \(z_k = X\hat{\boldsymbol{\theta}}_k + \mathbf{1} - \mathbf{d} / (\exp(X\boldsymbol{\theta}_k) * \mathbf{e_c})\).
L’algorithme de Newton peut être utilisé pour trouver le zéro d’une fonction \(\mathbf{f} : \mathbb{R}^p\rightarrow \mathbb{R}^n\) suffisamment régulière et une fois dérivable. Il s’agit d’un algorithme récursif qui à partir d’un vecteur \(\mathbf{x}_0\) construit une suite \((\mathbf{x}_k)_{k \ge 0}\) convergeant vers le vecteur \(\mathbf{x}\) vérifiant \(\mathbf{f}(\mathbf{x}) = \mathbf{0}\).
Le passage de \(\mathbf{x}_k\) à \(\mathbf{x}_{k + 1}\) s’appuie sur un développement de Taylor à l’ordre 1 au voisinage de \(\mathbf{x}_k\) :
\[ \small \mathbf{0} = \mathbf{f}(\mathbf{x}) \simeq \mathbf{f}(\mathbf{x}_k) + J(\mathbf{x}_k)(\mathbf{x} - \mathbf{x}_k)\quad\text{où}\quad J(\mathbf{x}_k) = \left(\frac{\partial f_i}{\partial x_j}\right)_{1\le i\le n, 1\le j\le p} \text{ est la matrice Jacobienne de } \mathbf{f}\text{ évaluée en } \mathbf{x}_k. \]
Cela permet d’obtenir une expression récursive des termes de la suite : \(\mathbf{x}_{k + 1} = \mathbf{x}_k - J^{- 1}(\mathbf{x}_k) \mathbf{f}(\mathbf{x}_k)\).
L’initialisation de cet algorithme peut être réalisée de différente manières, mais l’une des plus utilisées consiste à remplacer \(X\hat{\boldsymbol{\theta}}_0\) par \(\ln\mathbf{d} - \ln\mathbf{e^c}\) qui correspond au logarithme du taux brut de mortalité introduit à la section précédente. Il n’est pas nécessaire ici de connaître la valeur de \(\hat{\boldsymbol{\theta}}_0\) qui mène à ce résultat et cette valeur n’existe d’ailleurs pas sauf cas particulier. 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}\). Ainsi, avec ce choix d’initialisation, utiliser de l’algorithme de Newton revient à partir de l’approximation normale associée à l’estimateur des taux bruts puis à raffiner itérativement le choix des vecteurs d’observation et de poids à chaque étape.
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{\boldsymbol{\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{\boldsymbol{\theta}}_1 = (X^{'} W_0X)^{- 1}X^{'} W_0\: \mathbf{z}_0\).
Calculer \(W_k = W_{\hat{\boldsymbol{\theta}}_k}\) et \(\mathbf{z_k} = X\hat{\boldsymbol{\theta}}_k + \frac{\mathbf{d}}{\mathbf{e^{X\hat{\boldsymbol{\theta}}_k}} *\:\mathbf{e^c}} - \mathbf{1}\)
Calculer \(\hat{\boldsymbol{\theta}}_{k + 1}\) à l’aide de la formule : \(\hat{\boldsymbol{\theta}}_{k + 1} = (X^{'} W_kX)^{- 1}X^{'} W_k\:\mathbf{z_k}\)
Si \(\ell(\hat{\boldsymbol{\theta}}_{k + 1}) - \ell(\hat{\boldsymbol{\theta}}_k) < \epsilon\) avec par exemple \(\epsilon = 10^{- 12}(\mathbf{d}^T\mathbf{1})\), arrêter l’algorithme. Sinon incrémenter \(k\) et retourner à l’étape 2. Rappelons que \(\ell(\boldsymbol{\theta}) = \ln\boldsymbol{\mu}(\boldsymbol{\theta})^{T}\mathbf{d} - \boldsymbol{\mu}(\boldsymbol{\theta})^{T}\mathbf{e^c}\).
La convergence de cet algorithme très performant s’obtient en général au bout de quelques itérations seulement.
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 (GLMs 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}(\boldsymbol{\theta})\) et d’autre part au nombre total d’années d’observations \(\mathbf{e^c}\) qui est ici considéré comme déterministe. Formellement, on écrit : \(\mathbf{D} | \mathbf{e^c} \sim \mathcal{P}\left(\boldsymbol{\mu}(\boldsymbol{\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 apparaît naturellement lorsqu’on modélise un nombre d’évènements observés.
La vraisemblance (conditionnelle) associée à ce modèle s’écrit :
\[\begin{aligned} \ell(\boldsymbol{\theta} | \mathbf{e^c}) = & \ln \left\lbrace \underset{\chi}{\prod} \mathbb{P}(D_\chi = d_\chi | e^c_\chi) \right\rbrace = \ln \left\lbrace \underset{\chi}{\prod} e^{- \mu_\chi(\boldsymbol{\theta}) e^c_\chi} \frac{(\mu_\chi(\boldsymbol{\theta}) e^c_\chi)^{d_\chi}}{d_\chi!}\right\rbrace \\ = & \underset{\chi}{\sum} \left\lbrace - \mu_\chi(\boldsymbol{\theta}) e^c_\chi + d_\chi \ln\mu_\chi(\boldsymbol{\theta}) + d_\chi \ln e^c_\chi + \ln(d_\chi!)\right\rbrace \end{aligned} \tag{3.5}\]
Cette somme comprend 4 termes dont les 2 derniers ne dépendent pas de \(\boldsymbol{\theta}\). Les 2 premiers termes coïncident quant à eux avec ceux de l’Équation 3.2. Cela signifie notamment que les estimateurs de \(\boldsymbol{\theta}\) obtenus en maximisant les expressions de l’Équation 3.2 et l’Équation 3.5 sont identiques.
Ce nouveau modèle appartient à la famille des modèles linéaires généralisés (GLMs). 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 conditionnellement à l’exposition au risque \(\mathbf{e^c}\).
La partie modélisée : dans le cadre des modèles linéaires, elle s’écrira en notation matricielle \(X\boldsymbol{\theta}\) où \(\boldsymbol{\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} | \mathbf{e^c}) = \ln(\boldsymbol{\mu}(\boldsymbol{\theta}) \mathbf{e^c}) = \ln\boldsymbol{\mu}(\boldsymbol{\theta}) + \ln \mathbf{e^c}\) où \(\ln\boldsymbol{\mu}(\boldsymbol{\theta}) = X\boldsymbol{\theta}\) correspondra à la partie modélisée. Par ailleurs dans le cadre de la loi de Poisson \(\text{var}(\mathbf{D} | \mathbf{e^c}) = \mathbb{E}(\mathbf{D} | \mathbf{e^c}) = \boldsymbol{\mu}(\boldsymbol{\theta})\mathbf{e^c}\).
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 de valeur fixe et unitaire.
L’équivalence, à une constante près des vraisemblances de l’Équation 3.3 et de l’Équation 3.5, signifie que tous les résultats issus de la théorie des GLMs sont immédiatement transposables au modèle de durée défini par l’Équation 3.3. Cela comprend notamment la méthode d’estimation, les intervalles de confiance, la validation des résidus du modèle ou encore les tests d’hypothèse.
Notons toutefois que les hypothèses sous-jacentes des deux méthodes sont très différentes. Le GLM Poisson repose sur une hypothèse très forte sur la distribution du nombre de décès quand le modèle de durée de l’Équation 3.3 suppose uniquement :
l’indépendance entre les décès des individus
le fait que ces individus sont identiques à l’exception des variables explicatives prises en compte dans le modèle
la constance par morceaux de la force de mortalité.
Si l’on suppose que le nombre de décès observés \(\mathbf{D} | \mathbf{e^c}\) suit une loi de Poisson de paramètre \(\boldsymbol{\mu}(\boldsymbol{\theta}) \mathbf{e^c}\) conditionnellement à \(\mathbf{e^c}\) où \(\boldsymbol{\mu}(\boldsymbol{\theta})\) est la force de mortalité et \(\mathbf{e^c}\) l’exposition centrale au risque, on retombe sur l’estimateur de l’Équation 3.3. Les méthodes issues du cadre des GLMs sont ainsi directement applicables à ce modèle de durée.
3.3 Données de population
Cette ultime section introduit à travers les données de la Human Mortality Database un exemple de donné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 | 2021 | 331 806 167 | 3 464 231 |
| Russia | 1959 | 2014 | 143 829 248 | 1 878 039 |
| Japan | 1947 | 2021 | 122 923 519 | 1 439 856 |
| Germany | 1956 | 2020 | 83 165 834 | 985 572 |
| U.K. | 1922 | 2020 | 66 938 304 | 689 629 |
| France | 1816 | 2020 | 65 368 119 | 654 599 |
| Italy | 1872 | 2020 | 59 438 188 | 742 842 |
| Republic of Korea | 2003 | 2020 | 51 356 413 | 304 948 |
| Spain | 1908 | 2020 | 47 364 869 | 492 447 |
| Ukraine | 1959 | 2013 | 45 306 382 | 662 368 |
| Poland | 1958 | 2019 | 38 400 047 | 409 709 |
| Canada | 1921 | 2020 | 37 941 481 | 307 205 |
| Australia | 1921 | 2020 | 25 586 764 | 162 567 |
| Taiwan | 1970 | 2019 | 23 596 830 | 175 546 |
| Chile | 1992 | 2020 | 17 859 399 | 125 841 |
| Netherlands | 1850 | 2020 | 17 441 678 | 168 678 |
| Belgium | 1841 | 2021 | 11 568 653 | 112 331 |
| Greece | 1981 | 2019 | 10 722 972 | 124 965 |
| Czechia | 1950 | 2021 | 10 528 215 | 139 891 |
| Sweden | 1751 | 2022 | 10 489 651 | 94 737 |
| Portugal | 1940 | 2021 | 10 332 833 | 124 802 |
| Hungary | 1950 | 2020 | 9 751 911 | 141 002 |
| Belarus | 1959 | 2018 | 9 484 490 | 120 053 |
| Austria | 1947 | 2019 | 8 879 702 | 83 386 |
| Switzerland | 1876 | 2021 | 8 703 647 | 71 192 |
| Israel | 1983 | 2016 | 8 546 145 | 44 203 |
| Hong Kong | 1986 | 2020 | 7 492 088 | 50 666 |
| Bulgaria | 1947 | 2021 | 6 879 290 | 148 995 |
| Denmark | 1835 | 2022 | 5 903 405 | 59 435 |
| Finland | 1878 | 2022 | 5 556 551 | 63 041 |
| Norway | 1846 | 2022 | 5 458 444 | 45 774 |
| Slovakia | 1950 | 2019 | 5 454 025 | 53 234 |
| New Zealand | 1948 | 2021 | 5 113 692 | 34 908 |
| Ireland | 1950 | 2020 | 4 982 855 | 32 387 |
| Croatia | 2001 | 2020 | 4 048 049 | 57 023 |
| Lithuania | 1959 | 2020 | 2 795 648 | 43 547 |
| Slovenia | 1983 | 2019 | 2 088 355 | 20 588 |
| Latvia | 1959 | 2019 | 1 913 855 | 27 719 |
| Estonia | 1959 | 2019 | 1 326 978 | 15 401 |
| Luxembourg | 1960 | 2021 | 640 012 | 4 489 |
| Iceland | 1838 | 2021 | 372 550 | 2 333 |
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 |
|---|---|---|---|---|---|
| Russia | 1984 | 60 | Femme | 934 590 | 12 169 |
| U.S.A. | 1990 | 80 | Femme | 598 330 | 31 825 |
| U.K. | 1939 | 8 | Homme | 342 063 | 486 |
| Belgium | 1892 | 21 | Femme | 53 458 | 354 |
| Russia | 2002 | 47 | Homme | 1 114 120 | 22 316 |
| Netherlands | 1883 | 1 | Homme | 57 651 | 3 396 |
| Austria | 1949 | 65 | Femme | 35 650 | 869 |
| Spain | 1943 | 9 | Homme | 278 393 | 629 |
| France | 1939 | 26 | Homme | 305 908 | 1 330 |
| U.S.A. | 2012 | 11 | Femme | 2 030 056 | 200 |
| U.S.A. | 1939 | 28 | Femme | 1 103 152 | 3 052 |
| U.S.A. | 1934 | 16 | Homme | 1 180 467 | 2 698 |
| Italy | 1896 | 20 | Femme | 256 214 | 1 857 |
| Germany | 2020 | 63 | Femme | 559 950 | 3 778 |
| France | 1842 | 51 | Femme | 178 298 | 2 990 |
| France | 1955 | 1 | Homme | 394 736 | 1 691 |
| Canada | 1974 | 9 | Homme | 229 096 | 102 |
| Denmark | 1884 | 41 | Femme | 11 497 | 111 |
| Germany | 1974 | 38 | Femme | 573 055 | 758 |
| U.K. | 1981 | 80 | Femme | 168 753 | 12 515 |
3.3.2 Données HMD pour la France
Les données HMD pour la France couvrent la période allant de 1816 à 2020. 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.
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 équivalent au précédent mais qui relève du cadre théorique des GLM. Les avantages associés sont un algorithme d’estimation performant, et la disponibilité d’intervalles de confiance, de méthodes de validation ou de tests d’hypothèses associés.