灰色模型在城乡居民心脏病死亡趋势中的建模研究

Modeling Study on Grey Model in the Trend of Heart Disease Mortality among Urban and Rural Residents

摘要: 为对城乡居民未来5年的心脏病死亡率作出预测，本文通过分析2008年至2021年城乡心脏病死亡率的发展趋势，结合Matlab软件，使用MDNGM(1, 1)模型对城乡居民心脏病死亡率建模，并分别与GM(1, 1)模型、DNGM(1, 1)模型的预测结果作对比。结果表明，MDNGM(1, 1)模型对城乡居民心脏病死亡率的预测精度分别为99.797%和99.157%。因此，MDNGM(1, 1)模型可以很好地拟合与预测城乡心脏病死亡率，为预防和控制心脏病提供了科学的预前理论依据。

Abstract: In order to forecast the heart disease mortality rates for urban and rural residents over the next five years, by analyzing trends in heart disease mortality rates in urban and rural areas from 2008 to 2011 and employing Matlab, the MDNGM(1, 1) model was utilized to establish a model for heart disease mortality in this paper. Subsequently, the predictions were compared with those generated by the GM(1, 1) and DNGM(1, 1) models. The results indicated that the MDNGM(1, 1) model achieved a prediction accuracy of 99.797% and 99.157% for urban and rural residents, respectively. Consequently, the MDNGM(1, 1) model effectively models and predicts the urban and rural heart disease mortality rates, offering a scientific advance theoretical basis for the prevention and control of heart disease.

1. 引言

世界心脏病联盟(WHF)发布的报告显示，心血管病长期以来一直是全球居民死亡的主要原因。2021年，有2050万人因心血管病去世，约占全球死亡总人数的三分之一；《中国心血管健康与疾病报告2022》指出，中国心血管病患病率处于持续上升阶段，推算2022年患病人数约为3.3亿[1]。2019年中国死因监测数据显示，心脏病死亡率为160.26/10万，位于死亡顺位第2位，仅次于恶性肿瘤。近几年心脏病死亡率趋势研究成为热门话题，如韦琴等人通过Joinpoint回归模型描述分析了2004~2019年中国城乡人群心脏病死亡率的时间变化趋势[2]；丁璐璐等人通过启东市居民死因登记监测系统收集了1990~2019年启东市居民心脏病死亡资料，并使用SAS 9.2软件采用时间序列分析中的ARIMA模型进行了趋势预测[3]；赵明日等人通过利用2019全球疾病负担研究中1990~2019年中国居民归因于高盐饮食的IHD死亡和DALY数据，应用年龄–时期–队列(APC)模型探讨了年龄–时期–队列效应[4]等。而通过梳理相关文献发现，灰色预测模型适宜处理历史数据少、波动小的预测问题，在趋势预测中受到了众多学者的青睐，但灰色模型预测精度相对较低。因此，本文基于2008年~2021年城乡心脏病死亡率，使用MDNGM(1, 1)模型进行建模，并对心脏病的发展趋势进行预测，旨在为未来制定心脏病防治策略提供科学的理论依据。

2. 方法

2.1. GM(1, 1)模型

灰色预测模型在建模过程中有着所需数据量较少，精度较高，同时预测模型的样本分布不需要有规律性的优点[5]。随着灰色理论的不断发展，GM(1, 1)模型的使用越来越广泛。建立灰色GM(1, 1)模型的步骤如下：

设原始时间序列数据为$X^{\left(0\right)}=\left\{X^{\left(0\right)}\left(i\right),i=1,2,\cdots ,n\right\}$ ，原始时间序列多表现为无规律性和波动性。为减少原始序列的波动性，生成一次累加序列：

$X^{\left(1\right)}=\left\{X^{\left(1\right)}\left(k\right),k=1,2,\cdots ,n\right\}$ ,

其中，$X^{\left(1\right)}\left(k\right)={\displaystyle \sum _{i=1}^k{X}^{\left(0\right)}}\left(i\right)$ ，$\left(k=1,2,\cdots ,n\right)$ 。

建立一阶微分方程：

$\frac{\text{d}{X}^{\left(1\right)}}{\text{d}t}+a{X}^{\left(1\right)}=b$ .

生成$X^{\left(1\right)}$ 的紧邻均值序列：

$Z^{\left(1\right)}=\left(z^{\left(1\right)}\left(2\right),z^{\left(1\right)}\left(3\right),\cdots ,z^{\left(1\right)}\left(n\right)\right)$ ,

其中，$z^{\left(1\right)}\left(k\right)=\frac{1}{2}\left[X^{\left(1\right)}\left(k\right)+X^{\left(1\right)}\left(k-1\right)\right],k=2,3,\cdots ,n$ 。

由最小二乘法可得${\left(a,b\right)}^{\text{T}}={\left(B^{\text{T}}B\right)}^{-1}{B}^{\text{T}}Y$ ，其中

$B=\left[\begin{array}{cc}-z^{\left(1\right)}\left(2\right)& 1\\ -z^{\left(1\right)}\left(3\right)& 1\\ ⋮& ⋮\\ -z^{\left(1\right)}\left(n\right)& 1\end{array}\right],Y=\left[\begin{array}{c}{X}^{\left(0\right)}\left(2\right)\\ X^{\left(0\right)}\left(3\right)\\ ⋮\\ X^{\left(0\right)}\left(n\right)\end{array}\right]$ .

通过求解后可得出微分方程[6]，即

${\widehat{X}}^{\left(1\right)}\left(k\right)=\left(X^{\left(0\right)}\left(1\right)-\frac{b}{a}\right){\text{e}}^{-a\left(k-1\right)}+\frac{b}{a}$.

累减还原可得原始序列的预测值：

${\widehat{X}}^{\left(0\right)}\left(k\right)={\widehat{X}}^{\left(1\right)}\left(k\right)-{\widehat{X}}^{\left(1\right)}\left(k-1\right)$ .

2.2. DNGM(1, 1)模型

DNGM(1, 1)模型也称为基于白化微分方程参数直接估计法的灰色预测模型，该模型直接从微分方程的时间响应函数出发，通过求解参数$a$ 、$b$ 和$c$ 的估计值，并代入时间响应函数，通过累减还原得到近似非齐次指数序列的表达式，从而得到预测结果。一般而言，从差分方程到微分方程的过程中，模型参数估计值$\widehat{\alpha }$ 、$\widehat{\beta }$ 和$\widehat{\gamma }$ 存在跳跃性误差。然而，DNGM(1, 1)模型避免了此类误差，其得到的模拟值不存在模型误差，只存在计算误差[7]-[12]。具体建模过程如下：

设原始时间序列数据为$X^{\left(0\right)}=\left\{X^{\left(0\right)}\left(i\right),i=1,2,\cdots ,n\right\}$ ，$X^{\left(1\right)}$ 为序列$X^{\left(1\right)}$ 的一次累加生成(1-AGO)序列，即$X^{\left(1\right)}=\left\{X^{\left(1\right)}\left(t\right),t=1,2,\cdots ,n\right\}$ 。其中，$X^{\left(1\right)}\left(t\right)={\displaystyle \sum _{i=1}^t{X}^{\left(0\right)}}\left(t\right)$ ，$\left(t=1,2,\cdots ,n\right)$ 。

建立DNGM(1, 1)模型的白化微分方程：

$\frac{\text{d}{X}^{\left(1\right)}\left(t\right)}{\text{d}t}+a{X}^{\left(1\right)}\left(t\right)=bt+c$ ,

其中$a$ 、$b$ 、$c$ 为待定参数，对白化微分方程进行求解运算，可得：

$X^{\left(1\right)}\left(t\right)=\left(X^{\left(1\right)}\left(1\right)-\frac{b}{a}-\frac{c}{a}+\frac{b}{a^2}\right){\text{e}}^{-a\left(t-1\right)}+\frac{b}{a}t+\frac{c}{a}-\frac{b}{a^2}$ .

进行累减还原后，可得最终响应式：

${\widehat{X}}^{\left(0\right)}\left(t\right)=X^{\left(1\right)}\left(t\right)-X^{\left(1\right)}\left(t-1\right)=\left(1-{\text{e}}^a\right)\left(X^{\left(1\right)}\left(1\right)-\frac{b}{a}-\frac{c}{a}+\frac{b}{a^2}\right){\text{e}}^{-a\left(t-1\right)}+\frac{b}{a}$.

由上式可得：

$X^{\left(1\right)}\left(t+1\right)={\text{e}}^{-a}{X}^{\left(1\right)}\left(t\right)+\frac{b}{a}\left(1-{\text{e}}^{-a}\right)t+\left(1-{\text{e}}^{-a}\right)\left(\frac{c}{a}-\frac{b}{a^2}\right)+\frac{b}{a}$ .

令$\alpha ={\text{e}}^{-a}$ ，$\beta =\frac{b}{a}\left(1-{\text{e}}^{-a}\right)$ ，可求出误差平方和为：

$S={\displaystyle \sum _{t=1}^{n-1}{\left[X^{\left(1\right)}\left(t+1\right)-\widehat{\alpha }{X}^{\left(1\right)}\left(t\right)-\widehat{\beta }t-\widehat{\gamma }\right]}^2}$ .

根据最小二乘法原理构造关于$\alpha$ 、$\beta$ 和$\gamma$ 的非齐次方程组，使得误差平方和$S$ 的值最小[13]。

$\{\begin{array}{l}\alpha {\displaystyle \sum _{t=1}^{n-1}{\left(X^{\left(1\right)}\left(t\right)\right)}^2}+\beta {\displaystyle \sum _{t=1}^{n-1}t}{X}^{\left(1\right)}\left(t\right)+\gamma {\displaystyle \sum _{t=1}^{n-1}{X}^{\left(1\right)}}\left(t\right)={\displaystyle \sum _{t=1}^{n-1}{X}^{\left(1\right)}}\left(t\right)X^{\left(1\right)}\left(t+1\right)\hfill \\ \alpha {\displaystyle \sum _{t=1}^{n-1}t}{X}^{\left(1\right)}\left(t\right)+\beta {\displaystyle \sum _{t=1}^{n-1}{t}^2}+\gamma {\displaystyle \sum _{t=1}^{n-1}t}={\displaystyle \sum _{t=1}^{n-1}t}{X}^{\left(1\right)}\left(t+1\right)\hfill \\ \alpha {\displaystyle \sum _{t=1}^{n-1}{X}^{\left(1\right)}}\left(t\right)+\beta {\displaystyle \sum _{t=1}^{n-1}t}+\gamma \left(n-1\right)={\displaystyle \sum _{t=1}^{n-1}{X}^{\left(1\right)}}\left(t+1\right)\hfill \end{array}$ .

根据克拉默法则解该方程组，得到参数$a$ 、$b$ 和$c$ 的估计值为：

$\widehat{a}=-\ln \alpha$ , $\widehat{b}=\frac{a\beta }{1-\alpha }$ , $\widehat{c}=\frac{a\gamma -b}{1-\alpha }+\frac{b}{a}$ .

将参数估计值代入累减还原后的响应式，可得到原始序列的预测值。

2.3. MDNGM(1, 1)模型

马尔可夫链是随机变量$X_1$ ，$X_2$ ，$X_3$ …的一个数列。这些变量的范围，即他们所有可能取值的集合，被称为“状态空间”，而$X_n$ 的值则是在时间$n$ 的状态。如果$X_{n+1}$ 对于过去状态的条件概率分布仅是$X_n$ 的一个函数，则：

$P\left(X_{n+1}=x|X_0,X_1,X_2,\cdots ,X_n\right)=P\left(X_{n+1}=x|X_n\right)$ .

这里$x$ 为过程中的某个状态。上面这个恒等式可以被看作是马尔可夫性质[14]。

在使用DNGM(1, 1)模型对原始序列进行预测时，往往会有偏差。因此，需要用马尔可夫模型进行修正。马尔可夫模型的修正步骤如下[15] [16]：

① 划分状态区间：假设系统有$m$ 种状态$S_1,S_2,\cdots ,S_m$ ；

② 计算状态转移概率：$p_{ij}$ 表示系统从状态$i$ 转移到状态$j$ 的转移概率，

$p_{ij}=\frac{M_{ij}}{M_i}$ ,

在上式中，$M_i$ 表示状态$S_i$ 发生的总次数，$M_{ij}$ 表示状态$S_i$ 转移到状态$S_j$ 发生的次数；

③ 构建状态转移矩阵：状态转移概率矩阵$P$ 可表示为：

$P=\left[\begin{array}{cccc}{p}_{11}& p_{12}& \cdots & p_{1n}\\ p_{21}& p_{22}& \cdots & p_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ p_{n1}& p_{n2}& \cdots & p_{nn}\end{array}\right]$ ;

④ DNGM(1, 1)模型预测值修正：

${\widehat{Y}}^{\left(0\right)}\left(k\right)=\frac{1}{2}\left(m_i+n_i\right){\widehat{X}}^{\left(0\right)}\left(k\right)$ ,

上式中，$m_i$ 表示区间$i$ 的上限，$n_i$ 表示区间$i$ 的下限，${\widehat{X}}^{\left(0\right)}\left(k\right)$ 表示DNGM(1, 1)模型的预测值，${\widehat{Y}}^{\left(0\right)}\left(k\right)$ 表示经马尔可夫修正过后的DNGM(1, 1)模型预测值。

3. 数据建模

3.1. 数据来源

根据国家统计局公开发表的《中国统计年鉴》，整理得2008年至2021年的城乡居民心脏病死亡率数据。

3.2. 精度检验

本文模型检验的指标如下：

• 残差：

$\epsilon =X^{\left(0\right)}\left(k\right)-{\widehat{X}}^{\left(0\right)}\left(k\right)$

• 相对误差：

$\xi =\left|\frac{\epsilon }{X^{\left(0\right)}\left(k\right)}\right|$

• 平均相对误差：

$\overline{\xi }=\frac{1}{n-1}{\displaystyle \sum _{k=2}^n\left|\frac{\epsilon }{X^{\left(0\right)}\left(k\right)}\right|}$

• 相对值：

$\eta =\frac{X^{\left(0\right)}\left(k\right)}{{\widehat{X}}^{\left(0\right)}\left(k\right)}$

3.3. 预测模型对比

根据DNGM(1, 1)模型的预测结果，可对其相对值的跨度范围进行区间划分。不同的区间对应表示不同的转移状态，在区间划分的基础上进行状态转移矩阵概率计算。通过转移状态矩阵中的最优状态，从而对未来的变化趋势做出估算。区间划分过多会导致数据过于冗杂，划分太少又会影响预测结果的精度，一般以3~4个区间为宜[17]。因此，将相对值区间跨度进行四等均分。各年份所处状态如表1、表2所示。

Table 1. List of city status by year

表1. 城市各年份状态列表

年份	2008	2009	2010	2011	2012	2013	2014
状态	3	2	2	4	3	4	4
年份	2015	2016	2017	2018	2019	2020	2021
状态	3	2	2	3	1	1	3

Table 2. List of rural status by year

表2. 农村各年份状态列表

年份	2008	2009	2010	2011	2012	2013	2014
状态	2	3	1	3	1	4	3
年份	2015	2016	2017	2018	2019	2020	2021
状态	2	2	2	2	2	2	3

根据各年份的状态转移列表进行马尔可夫修正，修正后的模型预测值及相对误差分别如表3、表4所示。

Table 3. Urban relative error comparison table

表3. 城市相对误差对照表

年份	真实值 (1/10万)	GM(1, 1)模型		DNGM(1, 1)模型		MDNGM(1, 1)模型
		预测值 (1/10万)	|相对误差|	预测值 (1/10万)	|相对误差|	预测值 (1/10万)	|相对误差|
2008	121.00	121.0000	0.0000	121.0000	0.0000	121.0000	0.0000
2009	128.82	124.7497	0.0316	129.8289	0.0078	128.7254	0.0007
2010	129.19	127.1682	0.0156	130.4383	0.0097	129.3296	0.0011
2011	132.04	129.6336	0.0182	131.2080	0.0063	131.9296	0.0008
2012	131.64	132.1468	0.0038	132.1801	0.0041	131.9818	0.0026
2013	133.84	134.7087	0.0065	133.4078	0.0032	134.1415	0.0023
2014	136.21	137.3203	0.0082	134.9585	0.0092	135.7008	0.0037
2015	136.61	139.9825	0.0247	136.9170	0.0022	136.7116	0.0007
2016	138.70	142.6963	0.0288	139.3905	0.0050	138.2057	0.0036
2017	141.61	145.4627	0.0272	142.5146	0.0064	141.3032	0.0022
2018	146.34	148.2828	0.0133	146.4604	0.0008	146.2407	0.0007
2019	148.51	151.1575	0.0178	151.4439	0.0198	149.0965	0.0039
2020	155.86	154.0880	0.0114	157.7380	0.0120	155.2930	0.0036
2021	165.37	157.0753	0.0502	165.6875	0.0019	165.4390	0.0004
$\overline{\xi }$			0.019792		0.006806		0.002030

Table 4. Rural relative error comparison table

表4. 农村相对误差对照表

年份	真实值 (1/10万)	GM(1, 1)模型		DNGM(1, 1)模型		MDNGM(1, 1)模型
		预测值 (1/10万)	|相对误差|	预测值 (1/10万)	|相对误差|	预测值 (1/10万)	|相对误差|
2008	87.10	87.1000	0.0000	87.1000	0.0000	87.1000	0.0000
2009	112.89	112.7994	0.0008	110.2756	0.0232	112.6879	0.0018
2010	111.34	117.4660	0.0550	116.0815	0.0426	110.5531	0.0071
2011	123.69	122.3256	0.0110	121.9057	0.0144	124.5724	0.0071
2012	119.50	127.3863	0.0660	127.7485	0.0690	121.6645	0.0181
2013	143.52	132.6563	0.0757	133.6098	0.0691	141.1755	0.0163
2014	143.72	138.1443	0.0388	139.4897	0.0294	142.5410	0.0082
2015	144.79	143.8594	0.0064	145.3882	0.0041	143.5163	0.0088
2016	151.18	149.8109	0.0091	151.3055	0.0008	149.3574	0.0121
2017	154.40	156.0086	0.0104	157.2416	0.0184	153.2429	0.0075

续表

2018	162.12	162.4628	0.0021	163.1964	0.0066	161.0952	0.0063
2019	164.66	169.1839	0.0275	166.1702	0.0092	164.0308	0.0038
2020	171.36	176.1831	0.0281	175.1629	0.0222	172.9077	0.0090
2021	188.58	183.4719	0.0271	185.1747	0.0181	189.2254	0.0034
$\overline{\xi }$			0.027543		0.023365		0.008430

经过对比表3、表4可知，MDNGM(1, 1)模型对城乡居民心脏病死亡率的预测精度更高，可将精度分别提升至99.797%和99.157%。为更直观地显示出不同模型的拟合效果，将城乡心脏病死亡率的真实值与以上三种模型的预测值对比作图。分别如图1、图2所示，可以看出，MDNGM(1, 1)模型的拟合结果更加准确，充分说明了MDNGM(1, 1)模型的优越性。

Figure 1. Fitted diagram of urban heart disease mortality rate

图1. 城市心脏病死亡率拟合图

Figure 2. Fitted diagram of rural heart disease mortality rate

图2. 农村心脏病死亡率拟合图

4. 城乡居民心脏病死亡率预测

根据DNGM(1, 1)模型预测值和状态转移矩阵可得到状态区间的预测值，取其区间中点进行马尔可夫修正，对未来5年城乡居民死亡率进行预测，可得到最终DNGM(1, 1)-马尔可夫预测值，结果如表5、表6所示。

Table 5. Prediction of heart disease mortality rate among urban residents

表5. 城市居民心脏病死亡率预测

年份	DNGM(1, 1) 模型预测值	预测区间	MDNGM(1, 1) 预测值
2022	175.7277	[173.6190, 174.8491]	174.2340
2023	188.4085	[186.1476, 187.4665]	186.8070
2024	204.4244	[201.9713, 203.4023]	202.6868
2025	224.6524	[221.9566, 223.5291]	222.7429
2026	250.2004	[247.1980, 248.9494]	248.0737

Table 6. Prediction of heart disease mortality rate among rural residents

表6. 农村居民心脏病死亡率预测

年份	DNGM(1, 1) 模型预测值	预测区间	MDNGM(1, 1) 预测值
2022	189.2055	[176.9071, 183.4820]	180.1946
2023	193.2554	[194.1250, 200.8407]	197.4829
2024	199.3246	[186.3685, 193.2950]	189.8318
2025	205.4130	[206.3374, 214.4361]	210.3867
2026	211.5207	[197.7719, 205.1222]	201.4470

5. 结论与建议

近年来，由于生活方式的改变、环境污染的加剧以及人们心理压力的增加等多种因素的影响，导致心脏病死亡率持续呈上升趋势。特别是在发展中国家，心脏病死亡率更是高居不下。自2008年以来，城乡心脏病死亡率虽有小范围的波动，但总体呈现出上升趋势。城市心脏病死亡率从2008年的121/10万升至2021年的165/10万，涨幅超36%，农村心脏病死亡率从2008年的87.1/10万升至2021年的188.58/10万，涨幅高达116.51%。城乡居民心脏病死亡率的差异，可能是由于经济条件和医疗水平不平衡造成的。表5和表6的预测结果表明，未来五年城乡居民心脏病死亡率将持续上升。因此，关注心脏病的预防和控制，不仅是每个人自身健康的保障，也是全社会共同关心的重要议题。为了遏制心脏病死亡率不断攀升的趋势，政府、医疗机构和个人都应积极采取有效措施，促进心脏健康，减少心脏病的发病和死亡率，以保障人民的健康和生命安全。首先个人应该保持适当的体重、良好的饮食习惯、定期运动、戒烟限酒等，这些都是预防心脏病的重要因素[18]。并定期进行心脏健康检查，包括测量血压、血糖、血脂等指标，及时发现潜在的心脏问题[19]。而对于已经确诊的心脏病患者，应按照医生的建议规范用药，以控制病情发展。同时，也要积极面对压力，保持良好的心理状态。政府应加强对心脏病预防知识的宣传和普及，提高公众对心脏病的认识和防范意识，通过健康教育有效降低心脏病发病率和死亡率。医疗机构应提高医疗服务水平，建立更加完善的心脏病防治体系，加强对心脏病患者的救治和管理，提高救治率和康复率。并开展相关的科学研究，加深对心脏病病因和治疗机制的了解，进一步提高心脏病的预防和治疗水平。