基于队列角度的高龄人口死亡率预测

Mortality Projections for the Elderly Population Based on Cohort Perspective

摘要: 传统的死亡率预测模型大多以时期效应为出发点，使用人类死亡率数据库(HMD)中的年龄–时期数据。而队列预测通过捕捉不同的队列效应来提高对死亡率发展动态的捕捉，而这也被传统的年龄–时期角度出发的模型所忽视。跟经典的死亡率预测模型相比，本文提出了一种基于年龄–队列数据的C-STAD死亡率预测方法，并对日本和加拿大的成年男性队列死亡率进行分析，通过两次样本外检验来评估模型预测的准确性。结果表明，使用队列数据构建的死亡率模型的误差远小于使用时期数据的误差，具有良好的预测效果。最后文章对两个国家的死亡率从队列角度和时期角度进行预测。

Abstract: Most of the traditional mortality prediction models take period effects as the starting point and use age-period data from the Human Mortality Database (HMD). Cohort prediction improves the capture of mortality dynamics by capturing different cohort effects, which are also ignored by traditional age-period models. Compared with the classical mortality prediction model, this paper proposes a C-STAD mortality prediction method based on age-cohort data, and analyzes the adult male mortality cohort in Japan and Canada, and evaluates the accuracy of the model prediction through two out-of-sample validations. The results show that the error of the mortality model based on cohort data is far less than that of the period data, and it has a good prediction effect. Finally, the mortality rates of the two countries are predicted from the perspective of cohort and period.

1. 引言

老龄化进程伴随着世界人口的快速增长：几乎每个国家都经历着老年人数量和比例的增长，人类死亡率在上个世纪大幅提高。在过去的一个世纪里，我国人口的预期寿命都从20世纪初的35岁左右增加到2018年统计的76.7岁 [1]，但这一明显的好处也给如医疗保健和养老金供应等机构带来了额外的压力。疾病发生率和死亡率的建模和预测在一般的流行病学与人口研究中，特别是对保险和养老金行业来说，是一个具有根本重要性的问题，因此越来越需要能够准确描述和预测死亡率的创新模型。精算师至少从1924年起就开始对死亡率进行预测，以应对死亡率提高对生活年金和养恤金的不利财政影响(Pollard, 1987) [2]。然而，直到最近30年，人们才提出并使用了更复杂的预测死亡率的统计方法(Booth等，2008) [3]。

绝大多数传统的随机动态死亡率模型都从年龄–时期角度出发，即使有考虑队列因素的模型，如对我国人口死亡率数据模拟效果较好的APC模型，但该模型通过提取HMD中生命表的对角线队列死亡率模式的方式来获得队列死亡率，而生命表数据由调整后的死亡率数据构成，对80岁以后的年龄使用逻辑函数，已失去实际观测数据的真实性。相比之下，从队列角度出发考虑死亡率，使用HMD队列数据，追踪不同时期同一组群体从出生年到死亡年龄的变化，分析结果更为客观。更普遍地说，队列死亡率的分析和预测值得探讨的主要原因有两个：首先，真实出生队列中的存活率不同于假设时期死亡率不变情况下的存活率，这是因为1) 速度效应；2) 队列效应；3) 选择(Borgan等，2019) [4]；第二，队列死亡率的发展是实际观察到的，它们可能不同于在周期生命表中假设的合成队列。对年龄–队列数据的分析也展现了与基于年龄–时期数据研究的不同的关于死亡率发展的视角。另外，Goldstein等(2006) [5] 指出，在死亡率稳步下降的人口中，预期寿命可以被视为队列寿命的滞后指标；Shkolnikov等(2011) [6] 的研究表明，自1840年以来，生于1870年至1920年的最佳实践期女性的预期寿命(即全国人口中观察到的最高预期寿命)增长速度几乎是最佳实践期女性预期寿命增长速度的两倍，最后，Borgan和Keilman (2019) [4] 证明日本和意大利女性与斯堪的纳维亚女性在预期寿命方面的差异，在考虑队列数据时就消失了，甚至出现了逆转。考虑到这些因素，根据年龄–队列数据预测队列死亡率比从年龄时期预测中提取队列模式表现出更合理的趋势。

在最早使用队列视角的研究中，Richards等(2006) [7] 采用Currie等(2004) [8] 的二维p样条模型来完成英格兰和威尔士队列的死亡率分析，Eilers等(2007) [9] 完成部分观察队列分年龄死亡人数分布，并假设潜在分布在年龄上是平滑的，Chiou等 [10] 提出了一种利用瑞典死亡率数据预测队列对数风险函数的功能数据分析方法。最近，Oeppen (2008) [11] 和Bergeron Boucher等(2017) [12] 提出了Lee-Carter模型的修正版本，使用主成分分析预测了生命表死亡分布，Pascariu等(2019) [13] 提出用统计矩构建最大熵死亡率模型预测死亡分布，Rizzi等(2020) [14] 提出使用惩罚复合链接模型，Basellini [15] 等在2020年提出使用分段函数构造队列分年龄死亡人数分布。

本文基于Basellini提出的分年龄死亡人数分布的方法，通过一个转换函数将观察数据的分布和“标准”分布进行转换，改变参考“标准”分布的年龄轴来捕捉死亡率动态。变换函数的特征是五个描述寿命变化和寿命变异性的参数，它采用分段函数的形式，称为C-STAD模型。本文首先选取日本和加拿大男性队列数据，取共有的出生年部分，即1887年至1953年，然后建立模型，求出参数，最后对队列死亡率进行分析和预测。另外，由于L-C是完全基于时期效应的模型，二维P样条模型是唯一一种从队列角度来预测队列死亡率的模型，所以本文选择以上两个模型作为对比，说明C-STAD方法的预测效果。除了提供关于人口死亡率经验的非常详细的描述外，C-STAD还提供了关于死亡时平均年龄、中位数和模态年龄等指标(Canudas-Romo等，2010) [16] 以及寿命变化。图1显示了日本男性部分队列的死亡模态年龄分布，从图中可以看到死亡模态年龄发生了显著的右移。因此，分年龄死亡人数分布提供了关于寿命和寿命变化的关键信息，这些信息不能直接从死亡率或生存函数中获得。

Figure 1. Age-at-death distributions for Japanese males

图1. 日本男性分年龄死亡人数分布

2. 预备知识

2.1. L-C模型

L-C模型由Lee等提出，主要思路是将死亡率的变化分解为时间因子 $t$ 和年龄因子 $x$。 $m_{x,t}$ 表示在时间 $t$ 时年龄为 $x$ 岁的人群的中心死亡率， $m_{x,t}$ 满足以下函数关系

$ln{m}_{x,t}={\alpha }_x+{\beta }_x{\kappa }_t+{\epsilon }_{x,t}$ (1)

且为了使参数估计得到唯一解，L-C模型需满足以下条件

${\displaystyle {\sum }_t{\kappa }_t}=0,{\displaystyle {\sum }_x{\beta }_x}=1$ (2)

其中， ${\alpha }_x$ 为年龄因子，反映分年龄死亡率自然对数的平均水平， $k_t$ 为时间因子，称为死亡指数，反映人口死亡随时间变化的速度； ${\beta }_x$ 表示年龄因子对时间因子的敏感度， ${\epsilon }_{x,t}$ 为误差项。

2.2. 基尼系数

基尼系数是社会科学中最常见的多样性或不平等的统计指标 [6]，表达式如下：

$G_0=1-\frac{1}{e\left(0\right){\left[l\left(\text{0}\right)\right]}^{\text{2}}}\cdot {\displaystyle \underset{0}{\overset{\infty }{\int }}{\left[l\left(x\right)\right]}^2\text{d}x}$ (3)

基尼系数最初在经济学中用于定量测定收入分配差异程度(Gini, 1914) [17]，后来成为统计学的常用指标之一，被用来衡量社会经济群体死亡率水平的可变性(Leclerc等，1990) [18]。近年来，该系数已被用来评估死亡率预测(Basellini等2019) [19] 和衡量人口内部和人口之间的寿命不平等(Gigliarano等，2017) [20]。基尼系数在0 (完全相等)和1 (完全不相等)的极限之间变化。对于寿命分布，等于0时表示所有人都在相同年龄死亡，如果所有人都在0岁死亡，而一个人无限长寿死亡，则基尼系数等于1。

2.3. 二维p样条法

P样条法(Eilers等，1996) [21] 是目前GLMs中公认的一种平滑方法，二维P样条法由一维B样条法扩展而来。本文通过使用带有惩罚的二维B-样条构建二元回归的基础，通常被称为P样条。对P样条法可以简单总结为：1) 使用P样条作为回归的基础，2) 通过对回归系数的差罚来修改对数似然 [8]。二维P样条是唯一一种从队列角度来预测队列死亡率的模型。

3. 数据

由于中国人口死亡率数据统计体系起步较晚，至今未存在有效的队列死亡率数据，故本文选取HMD中的日本和加拿大的男性队列死亡率数据进行分析和说明。为了更方便的比较拟合效果，我们选取日本加拿大队列数据的共有的出生年部分，即1887年至1953年。另外，本文将年龄从110岁扩展到120岁，以便观察分布的后期趋势。以日本数据为例，图2显示了原死亡率数据、扩展的死亡率数据和经过平滑处理的死亡率数据。

Figure 2. Mortality for Japan males

图2. 日本死亡率

本文数据的处理和作图均由R语言完成。

4. C-STAD模型

C-STAD模型，即队列分段转换分年龄死亡人数模型，首先采用队列死亡率数据，对于模态年龄前后分别使用不同的函数进行计算，其次通过转换函数将标准分布和观察分布进行转换，本节具体介绍C-STAD模型的结构，以及模型的参数估计与预测方法。

4.1. C-STAD模型

假设有分布 $f\left(x\right)$， $g\left(x\right)$， $x\ge 60$，令 $g\left(x\right)$ 为观察分布，是实际值形成的分布， $f\left(x\right)$ 为标准分布， $t\left(x;\theta \right)$ 为关于年龄 $x$ 和参数向量 $\theta$ 的转换函数，令

$g\left(x\right)=f\left(t\left(x;\theta \right)\right)$ (4)

其中 $\theta =\left(s,b_L,c_L,d_L,b_U\right)$,定义转换函数 $t\left(x;\theta \right)$ 如下：

$t\left(x;\theta \right)=\{\begin{array}{l}{M}^f+b_L\tilde{x}+c_L{\tilde{x}}^2+d_L{\tilde{x}}^3{{\displaystyle }}_{}如{果}_{}x\le M^g\\ M^f+b_U\tilde{x}{{\displaystyle }}_{}{{\displaystyle }}_{}{{\displaystyle }}_{}{{\displaystyle }}_{}{{\displaystyle }}_{}{{\displaystyle {}_{{}_{}}}}_{}如{果}_{}x\le M^g\end{array}$ (5)

其中 $\tilde{x}=x-s-M^f$，下标L和U是指分段变换(即 $M^g$ )的上下两部分，s表示的是在 $g\left(x\right)$ 和 $f\left(x\right)$ 之间模态年龄的差异， $b_L$ 和 $b_U$ 表示两种分布的模态年龄前后变异性的变化，而 $c_L$ 和 $d_L$ 代表了对于 $M^g$ 以下年龄时 $g\left(x\right)$ 和 $f\left(x\right)$ 在左尾不对称性和重量方面的差异。图3提供了C-STAD模型的图解。图中黑色线条表示 $b_L=b_U=1$， $c_L=d_L=0$，转换函数为 $t\left(x;\theta \right)=x$，即没有转换函数。再将之平移s个单位得到蓝色线条，此时转换函数为 $t\left(x;\theta \right)=x-s$。最后在 $b_L$ 、 $c_L$ 、 $d_L$ 和 $b_U$ 的作用下得到橙色线条，此时转换函数为一般式。

Figure 3. A schematic overview of the C-STAD model

图3. 分段转化示意图

4.2. 标准分布

通过使用二维p样条平滑法对1887~1953年的队列死亡率进行拟合，得到每个队列的分年龄死亡人数分布。具体来说，我们假设给定年龄 $x$ 和队列 $c$ 的观察死亡数 $d_{x,c}$，随机变量 $D_{x,c}$ 遵循泊松分布(Brillinger, 1986) [22]：

$D_{x,c}~P\left(e_{x,c},{\mu }_{x,c}\right)$ (6)

$e_{x,c}$ 表示暴露人数， ${\mu }_{x,c}$ 表示死亡力，使用 $b$ 样条基底的张量积对年龄和队列进行平滑处理，并将暴露量作为偏移量。为了解释不完整的队列，并在数据中保持矩形结构，我们令回归权重

$w=\left({\omega }_{x,c}\right)=\{\begin{array}{l}1,{{\displaystyle }}^{}{{\displaystyle }}^{}如果d_{x,c}和e_{x,c}为观察到的数据\\ 0,{{\displaystyle }}^{}{{\displaystyle }}^{}其他\end{array}$ (7)

再通过贝叶斯信息准则最小化方法选择年龄和队列间的平滑参数，进而通过取均值推导出完整和部分的分布。

最后，我们将标准分布 $f\left(x\right)$ 表示为年龄 $x$ 和系数 ${\beta }_f$ 的等间距 $b$ 样条基 $B\left(x\right)$ 的线性组合：

$f\left(x\right)=\exp \left[B\left(x\right){\beta }_f\right]$ (8)

在最后一步中，我们选择了大量带有较小惩罚项的 $b$ 样条，在从原始数据中去除不必要的随机波动能够在标准分布中保留所有重要的特征。

4.3. 参数的估计与预测

根据估计的标准分布，我们可以得出1887~1953年各队列的C-STAD参数 $\theta =\left(s,b_L,c_L,d_L,b_U\right)$。如第3.2节所述，根据每个队列的可用数据，我们使用三种不同的方法来估计θ；因此，我们将队列分为三组，第一组数据是完全观察到的，第二组大部分数据为观察数据，第三组只有少数队列数据为观察数据。对于前两组首先使用一维p样条捕获队列的死亡率波动，计算每个队列的密度函数和死亡模态年龄估计参数 $\hat{s}$ ( ${\hat{s}}_c=M_c-M_f$ )，然后使用最大似然函数估计其余参数，即：

$\ln L\left(\alpha |d_{x,c}\text{,}{e}_{x,c}\text{,}{w}_{x,c}\text{,}{\hat{s}}_c\text{,}{\beta }_f\right)$ (9)

对于第三组数据，使用 $c_1$ 和 $c_2$ 中估计的 $\hat{s}$ 和 ${\hat{b}}_U$ 来预测第三组的 $\hat{s}$ 和 ${\hat{b}}_U$，其余参数遵循前述方式。

5. C-STAD模型预测性能检验

首先我们使用文章所选数据对日本和加拿大的队列死亡率得到标准分布，然后进行参数估计，最后我们对得到的分布进行样本外验证，测试其拟合效果。

5.1. 标准分布

由日本和加拿大数据可以得到观察数据的分年龄死亡人数分布(数据已经经过平滑处理)，为了计算标准分布，我们将每个队列的分布进行对齐(Ramsay等，2005) [23]，这是一种经常用于功能性数据分析的技术，其目的是对齐所分析曲线的重要特征，对于分年龄死亡人数分布，本文按死亡模态年龄进行对齐。对齐通过分布的平动变换实现，保留了分布的除了模态死亡年龄外的所有特征(图4和图5(左))。对齐分布后求均值可得标准分布(图4和图5(右))。

5.2. 参数估计

估计的参数如图6所示，可以看到对于日本， $\hat{s}$ 和 ${\hat{b}}_U$ 一直处于增大趋势，而 $b_L$ 到后期逐渐趋向于零， $c_L$ 和 $d_L$ 逐渐趋向于零，加拿大也是相似的趋势，说明在拟合后期，转换函数的模态年龄前的函数趋向于系数为1的一次函数，但有别于模态年龄之后的一次函数。

Figure 4. Age-at-death distributions for Japan males

图4. 日本分年龄死亡人数分布

Figure 5. Age-at-death distributions for Canada males

图5. 加拿大分年龄死亡人数分布

5.3. 样本外测试

本节通过对日本和加拿大成年男性数据进行两组预测性样本外验证来评估C-STAD模型的准确性。具体来说，假设收集数据的最后一年为2015 − h，其中h = 5年和10年，然后将C-STAD模型拟合到完全观测的队列 = [1887, …, 1905 − h]，我们预测h年的死亡率。另外，寿命(前者)和寿命不平等(后者)的衡量都有助于评估死亡率预测的准确性(Bohk Ewald等，2017) [24]，故我们选择将60岁时的预期寿命(e60)和60岁时的基尼系数(G60)与样本外观察值进行比较。本文选用均方根误差(RMSE)、平均绝对误差(MAE)和平均绝对百分比误差(MAPE)作为衡量误差的标准 [25]。从表1~4可以看到，不管对于日本还是加拿大，C-STAD对于预期寿命的预测误差和基尼系数总要小于二维P样条和LC模型，e60和G60的各误差值在两组样本外验证中都很低，并且不会随着预测范围的增加而显著增加。

Figure 6. Parameter estimation

图6. 参数估计

Table 1. e60 error value for each model of Japanese males

表1. 日本各模型e60误差值

Table 2. G60 error value for each model of Japan males

表2. 日本各模型G60误差值

Table 3. e60 error value for each model of Canada males

表3. 加拿大各模型e60误差值

Table 4. G60 error value for each model of Canada females

表4. 加拿大各模型G60误差值

5.4. 样本外测试

图7和图8显示了观测值与C-STAD (80%的点置信区间)和二维P样条对三组数据的拟合和补充值之间的比较，从图7和图8(左)可以看到，两个国家的剩余寿命一直处于上升趋势，这种趋势与人们生活水平的提高和医疗条件的改善有直接联系， $c_3$ 的置信区间相对于 $c_2$ 明显变宽，是由于 $c_2$ 的数据大部分来自观测值，而 $c_3$ 只有一小部分为观测值，其余为预测值。就基尼系数而言(图7和图8(右))，日本在1919年到1930年处于停滞状态，死亡率在这段时期内没有得到明显改善，引起两个国家的寿命不平等。但总体来看，两个国家整体是递减状态，说明死亡率逐渐趋向于稳定的正常死亡状态。

下面用时期和队列两种角度对部分群体的死亡率进行预测。首先图9和图10左侧从时期角度对年龄为60、70、80、90岁的死亡进行拟合和预测，死亡率随着年份的后移逐渐下降，说明人们的寿命逐渐延长，预期寿命不断增加，其次两图右侧为以队列为主的60岁以上年龄的死亡率发展，表现出上升状态，符合实际情况，即同一批高龄群体随着时间的推移死亡率逐渐增大，剩余寿命逐渐减小，直至为0。在1887年，死亡率相对较高，后逐渐下降，到1853年时，死亡率明显改善，且上升曲线平坦。

Figure 7. Residual life expectancy (left) and Gini coefficient (right) of 60- year-olds in Japan

图7. 日本60岁群体剩余寿命(左)和基尼系数(右)

Figure 8. Residual life expectancy (left) and Gini coefficient (right) at age 60 in Canada

图8. 加拿大60岁群体剩余寿命(左)和基尼系数(右)

总体而言，从图中可以看出两个国家的死亡率都在发生明显的改善，这也意味着对应的长寿风险随之产生。结合实际，随着人们生活水平的提高和医疗条件的改善，寿命后延已经成为必然趋势，死亡率的降低导致了长寿风险的增加，对各国家保险行业生命表的定制造成一定的困难。如果只从经典的时期效应考虑死亡率的发展，容易使长寿风险的度量和管理变得越来越不可控。而本模型很好的预测了死亡率趋势，大大降低了预测误差，对长寿风险的把控更为精准。

Figure 9. Cohorts mortality rates for selected ages in Japan

图9. 日本部分队列死亡率

Figure 10. Cohorts mortality rates for selected ages in Canada

图10. 加拿大部分队列死亡率

6. 结论

人口老龄化的问题日益加剧，给保险公司和养老金行业带来了巨大的挑战，死亡率直接影响着保险产品的定价和风险，因此对死亡率愈加精准的分析和预测能在很大程度上降低风险，减小损失。本文从对队列角度使用分段函数建立模型，分两组不同范围的时间段进行样本外验证，并与经典的Lee-Carter模型和二维P样条法进行预测效果的比较，全面分析了日本和加拿大高龄人口的死亡率发展趋势，经典的时期效应模型在三个模型中得到最大的预测误差，两个队列效应的模型相比，C-STAD模型表现出明显的预测优势。另外，本文使用剩余寿命和基尼系数两个指标，从剩余寿命和寿命的发展状态对日本和加拿大的死亡率进行说明，使死亡率的变化规律和异常状态可视化。最后本文从时期角度和队列角度分析和预测了死亡率的发展。

由于各国人口的发展规律不同，经典的死亡率模型并没有考虑死亡率数据自身的特点，给出的预测也缺乏有力支持。相比之下，C-STAD模型能捕捉数据内部特点和规律，对每一个年龄的群体给出对应的死亡分布，使得在此基础上对各群体未来死亡率的预测更有依据，由于可以对特定群体进行分析预测，故得到的死亡率准确度更高、更具有现实意义。但不足的是，使用本模型需要大量的队列数据，故对于部分队列数据缺失或队列数据较少的国家并不适用。

基金项目

本文由国家自然科学基金项目(12061066)和甘肃省自然基金项目(20JR5RA528)共同支持。