上QQ阅读APP看书,第一时间看更新
第二节 传染病疫情预测技术
预测分为定性预测及定量预测。
定性预测指由熟悉业务知识及经验丰富的工作人员,根据熟悉的历史资料和直观材料,通过知识及经验去分析判断事物的未来发展趋势,给出预测的判断。定性预测虽然方法简单,不需要借助复杂的数学公式,但是会带有一定的主观性。
量预测是指通过分析历史数据,结合数学及统计学知识,模拟历史数据的走势规律,建立合适的数学模型,用来推测数据的未来发展变化。常用的定量预测方法有回归分析法、Box-Jenkins模型法、灰色动态模型法、动力学模型法、人工神经网络(artificial neural networks,ANN)等。其中较流行,常用于传染病发病率预测的模型方法有:①Box-Jenkins模型方法:对于传染病的预测而言,该模型可以很好地控制发病的长期趋势、季节性、周期性等因素的影响,该方法最受学者的青睐。②灰色动态模型法:不需要太多数据、运算方便、原理简单以及可检验,所以自提出以来,应用范围已拓展到众多科学领域,并取得了显著成果。③动力学模型法:通过建立模型,拟合模型基本性态和数值,发现疾病流行传播的影响因素,分析其对疾病发展趋势的影响,为疾病今后发展提供预测预警。与传统的统计方法相比较,传染病动力学模型是从疾病传播机制中反映发病过程,一旦与实际流行病学调查结果相结合,可更加真实反映疾病流行状况。④人工神经网络方法,该方法能考虑到传染病发病率数据受多种因素影响,不规则、混沌等非线性特征而得到广泛的应用。
通常情况下,定性预测与定量预测配合使用预测效果更好。
近年来,传染病网络直报系统的建立,使得传染病相关数据资料库得到了进一步完善,有了较全面的传染病数据资料,通过对数据的统计分析,来挖掘与传染病相关的信息,对传染病的预防控制将会给予很大帮助。由传染病网络直报系统的数据对传染病进行预测分析,是对疾病未来发生、发展和流行趋势的研究,一个好的分析和预测是传染病预防与控制的重要环节,给出准确的预测,才能制订相应的预防措施。传染病的预测可以帮助及早发现疾病流行发展的趋势,增加预防控制工作的预见性,对传染病的治疗、预防及制定卫生决策都有很重要的现实意义。
在国外,一些学者在传染病建模预测方面做了许多研究,如,Korthals等拟合季节性自回归移动平均(SARIMA)模型,以时间序列揭示了肺结核发病的季节趋势。我国使用传染病预测方法对传染病分析预测起步较晚,20世纪80年代以后才逐步发展起来,后来慢慢成为疾病监测工作不可缺少的一部分。
此处重点介绍以下四种模型:
一、回归分析
回归分析是在分析自变量和因变量之间相关关系的基础上,把影响因变量的各种因素(自变量)找出来,建立变量之间的回归方程式(回归模型),然后利用样本数据估计模型参数,并对模型进行误差检验,如果模型确定,就可以用模型对因素的变化值进行预测。回归预测是回归方程的重要应用之一,当因变量为计量变量时,根据自变量的多少,可分为简单线性回归(simple linear regression)和多重线性回归(multiple linear regression),当因变量为分类变量时,如疾病发生与否(发生、不发生),则可采用logistic回归族。回归分析预测在医学统计中应用广泛且技术较成熟,能够综合考虑多种因素的共同作用,但模型需要的数据量较大,并且对样本的分布有较高的要求。这里重点介绍常用的三种回归分析预测方法:简单线性回归、多重线性回归和logistic回归。
(一)简单线性回归
1.简单线性回归的概念回归分析中,最简单的模型仅包含一个应变量( Y)和一个自变量( X),若两变量间的变化呈线性趋势,则选用直线回归方程来描述其变化规律,称之为简单线性回归,又称直线回归。总体直线回归方程式为:
Y= α+ βX+ε
其中ε为随机误差, β称为总体回归系数,即直线的斜率,含义为当 X每增加或减小1个单位时,因变量Y平均改变 β个单位。α称为截距,为 X=0时, Y的估计值。
样本资料的直线回归方程一般形式为:
为 y的估计值, a和 b分别为样本截距和回归系数。线性回归分析要求数据满足线性、独立性、正态性等方差的前提假设。一般要求因变量y是来自正态总体的随机变量,自变量x可以是正态随机变量,也可以是精确测量和严密控制的值。
2.直线回归分析预测的一般步骤
(1)绘制散点图:首先通过将 n个观察单位的变量( x, y)在直角坐标系中绘制散点图,判断x与y是否呈直线趋势。若呈线性趋势,则可拟合直线回归方程。
(2)求回归方程的回归系数 b和截距 a:方程 中的 a和 b是两个待定常数,根据样本实测( x, y)计算 a和 b的过程就是求回归方程的过程。利用最小二乘法原理,当各实测点到回归直线的纵向距离的平方和( )最小,方程能较好地反映各点的分布规律,求出 b和 a。
具体计算公式如下:
式中 l XY为 X、 Y的离均差积和, l XX为 X的离均差平方和;
(3)对回归方程进行假设检验:回归方程检验的目的是检验求得的回归方程在总体中是否成立,即是否样本代表的总体也有直线回归关系。我们知道即使 X、 Y的总体回归系数 β为零,由于抽样误差的原因,其样本回归系数 b也不一定为零,因此,须作 β是否为零的假设检验,方法有以下两种:
方差分析:其基本思想是将因变量 Y的总变异 SS 总分解为 SS 回归和 SS 剩余,然后利用 F检验来判断回归方程是否成立。
t检验:其基本思想是利用样本回归系数 b与总体均数回归系数进行比较来判断回归方程是否成立。
(4)利用回归方程进行预测预报:把预报因子(即自变量 x)代入回归方程对预报量(即因变量y)进行估计,即可得到个体y值的容许区间。
实际上,在传染病预测预警中的很多场合下,需要考察因变量如何依赖多因素变化而变化的规律,再通过这种规律对传染病的趋势进行预测预警。如手足口病的流行受到儿童EV71疫苗接种率、是否为流行季节、儿童人口密度等多种因素的影响,此时,所需要的统计分析方法称为“多重回归分析”,当因变量为计量变量时,称为“多重线性回归”。
(二)多重线性回归
1.多重线性回归的概念
当模型包含一个因变量与多个自变量的线性回归分析,称为多重线性回归,它是简单线性回归分析的推广。以 Y为因变量, X 1. X 2.…、 X i分别代表 i个自变量。对于总体而言,多重线性回归模型的表达式为:
其中 β 0为总体截距, β 1、 β 2、…、 β i分别为各个自变量所对应的总体偏回归系数, ε为随机误差,常假定其服从正态分布。偏回归系数 β i( i = 1、2、…、i)表示在其他自变量固定不变的情况下, X i每改变1个单位时,因变量 Y平均改变的单位数。对于样本而言,多重线性回归模型为:
其中 表示 Y的估计值, b 0、 b 1、 b 2、…、 b i为样本的截距和偏回归系数,它们是相应总体参数的估计值。多重线性回归分析同样要求数据满足线性、独立、正态等方差的前提假设,且要求样本量足够大,至少是自变量个数的10~20倍。
2.多重线性回归预测分析的步骤
(1)根据预测目标,确定自变量和因变量:
明确预测的具体目标,也就确定了因变量。如预测具体目标是各县区手足口病的发病率,那么手足口病发病率Y就是因变量。通过查阅资料,结合专业理论和经验,寻找与预测目标的可能相关影响因素,即自变量。
(2)判断变量间是否存在多重共线性:
多重共线性是指线性回归模型中的解释变量之间由于存在精确相关关系或高度相关关系而使模型估计失真或难以估计准确,如气候因素中的温度、湿度、降雨量等之间可能有显著的相关关系,即可能会存在多重共线性。严重的多重共线性会大大影响模型的预测结果。诊断协变量间的共线性是回归分析的重要环节,常用指标有简单相关系数、容忍度、方差膨胀因子法和条件指数等,一般来说,简单相关系数大于0.8、容忍度小于0.1、方差膨胀因子大于30、特征根趋于0、条件指数大于100,提示存在的共线性。
多重共线性问题可以在建模过程中借助软件来处理,如统计模型处理时可以通过逐步回归,在R语言中,用step()命令,选择最小AIC信息统计量的模型,可以排除引起共线性的变量,是解决多重共线性的比较常用方法。另一种方法多重共线性的解决方法是岭回归,R语言中的MASS包提供了lm.ridge()函数用于实现岭回归。
(3)筛选自变量,建立回归预测模型:
淘汰那些对因变量影响无统计学意义的自变量,从中选出主要的影响因素,使拟合的多重线性回归模型精简且具有更高的价值,即变量的筛选和建模过程。
变量筛选方法主要分为:以回归模型拟合优劣准则判断的全局择优法和基于统计检验准则的局部择优法。局部择优法筛选变量时常用的方法有后退法、前进法和逐步法。变量的筛选可以通过R语言、SAS软件等进行,为了选出最优模型,通常需要用到一些评判指标,如Cp、AIC、BIC以及Adjusted R2等,Adjusted R2越接近1说明模型拟合得越好,其他三个指标则是越小越好,在实际使用中,AIC和BIC应用更广,当样本量非常大时,AIC倾向于选择更多的变量,BIC倾向于选择出正确的那个子模型。
(4)检验回归预测模型,计算预测误差:
回归预测模型是否可用于实际预测,取决于对回归预测模型的检验和对预测误差的计算。回归方程只有通过各种检验,且预测误差较小,才能将回归方程作为预测模型进行预测。
(5)计算并确定预测值:
利用回归预测模型计算预测值,并对预测值进行综合分析,确定最后的预测值。
(三)logistic回归
前面介绍的简单线性回归和多重线性回归因变量均为计量资料,当因变量为分类变量时,如疾病发生与否(发生、不发生)、疾病严重程度(轻、中、重),则采用logistic回归,根据因变量的多少,可分为二分类logistic回归和多分类logistic回归。logistic回归由于因变量Y不满足多重线性回归模型的应用条件,因此,须对其先进行logit变换,其模型表达式为:
其中 P为发生概率,取值在[0,1]之内,发生概率 P与未发生概率1- P之比为优势(odds), β 1、 β 2、…、 β i、称为Logistic偏回归系数,当X i为二分类变量或无序分类变量时,偏回归系数 β i表示当其他自变量固定不变时,变量X i的一个类别和另一个类别或参照水平比值比的对数值。当X i为定量变量或有序多分类变量,偏回归系数 β i表示当其他自变量固定不变时,变量X i每改变一个单位,因变量y发生与不发生概率之比的对数值。
logistic回归模型预测分析的步骤与多重线性回归预测分析类似,首先要根据预测目标,确定自变量和因变量,进行自变量筛选,然后通过原始数据估计方程中的各个偏回归系数 β i,建立logistic回归模型,并对模型及各偏回归系数做假设检验,AIC和BIC准则也可用于推断logistic回归模型的拟合优度,对于同一组数据,进行变量筛选时,AIC和BIC值越小,表示模型的拟合效果越好。
(四)分析法实例
下面以我们曾经开展过的急性脑炎、脑膜炎病因预测为例,介绍回归预测分析的方法,具体见参考文献3。该案例利用广西2007—2012年急性脑炎、脑膜炎综合征监测的数据,根据541例急性脑炎、脑膜炎综合征监测病例的相关症状、体征和血液、脑脊液检测结果,对引起急性脑炎、脑膜炎的可能病因进行预测。具体做法如下:
1.首先确定预测的具体目标,也就是因变量。引起脑炎、脑膜炎的原因很多,常见的有细菌性、病毒性、真菌性等,临床上往往需要根据病因采取有针对性的治疗措施,如细菌性脑膜炎须尽早使用抗生素治疗,真菌性脑炎、脑膜炎抗真菌治疗对改善患者的预后非常有效,病毒性脑炎无特异性治疗办法,以对症支持治疗为主。由于病原菌的培养需要较长的时间周期,而且阳性率不高,大部分医院也不开展乙脑的病毒特异性抗体检测,如果临床用药依赖病原学检测结果,在实际工作中缺少可操作性。如果能找到一种简易的方法判断病例是由哪种原因引起,对临床的治疗非常有帮助,如通过患者的年龄,发热、呕吐等症状,脑膜刺激征等体征,血液和脑脊液等常规检测结果,来预测急性脑炎、脑膜炎的病因。因此,因变量为急性脑炎、脑膜炎病因,我们根据病原检测结果,分为乙脑、其他病毒性脑炎、细菌或真菌性脑炎、阴性组共4个组,可以看出,这样的因变量为无序多分类变量,且有多个自变量,因此,采用多项logistic回归(multinomial logistic regression)方法。
2.采用单因素分析的方法分析各种类型脑炎脑膜炎各个自变量指标之间的差异,涉及的自变量包括发热、是否高热(体温≥39℃)、头痛、腹泻、呕吐、喷射性呕吐、精神萎靡、嗜睡、意识障碍、抽搐、颈项强直、角弓反张、脑膜刺激征、皮肤瘀点瘀斑、血液白细胞计数、中性粒细胞计数、脑脊液蛋白含量、脑脊液白细胞计数、葡萄糖、氯化物,详见表4-1,将 P≤0.05的自变量将纳入多项logistic回归分析。
表4-1 各种病原急性脑炎、脑膜炎病例的临床特征及血液、脑脊液检测结果
*中位数(四分位数间距), P值为Kruskal-Wallis秩和检验。
3.由于患者的血液和脑脊液检测结果,如血液白细胞计数、中性粒细胞计数、脑脊液蛋白含量、脑脊液白细胞计数、葡萄糖、氯化物等指标为非独立变量,可能存在多重共线性,因此,先采用因子分析的方法对这些变量进行分析(表4-2),最后将因子得分纳入多项logistic回归分析。
表4-2 541例急性脑炎/脑膜炎病例血液和脑脊液生化检测结果探索性因子分析
*因子载荷值,正交旋转方法为方差极大化(varimax),因子1和因子2的累积方差贡献率为32.8%。
4.建立多项logistic回归模型。将症状、体征和血液、脑脊液检测指标的因子得分作为自变量,以乙脑、其他病毒性脑炎、细菌或真菌性脑炎、阴性组作为因变量,采用后退法筛选有模型变量,先将全部自变量纳入回归方程,然后采用AIC法筛选出最优模型,选择AIC值最小的模型。
最后纳入回归模型的变量见表4-3。结果得出年龄、意识改变、血液白细胞和中性粒细胞因子得分、脑脊液因子得分可以用于急性脑炎脑膜炎综合征病因预测。
表4-3 多项logistic回归模型预测急性脑炎/脑膜炎(阴性组 n=176为参照组)
如果把病原检测阴性组与乙脑、其他病毒性脑炎合并为一组,另一组为细菌/真菌性脑膜炎,因变量变成了一个二分类变量,最终只有年龄、意识障碍和脑脊液因子得分进入最终的logistic模型(表4-4)。即当患者年龄≥18岁,且出现意识障碍、脑脊液炎性改变时,更多可能为细菌或真菌感染。
表4-4 细菌/真菌性脑炎、脑膜炎logistic回归模型(病原检测阴性组与乙脑、其他病毒性脑炎合并为一组,并作为参照组)
续表
对模型的诊断效果(预测效果)进行评价。采用诊断特征曲线(ROC曲线)下面积对模型的诊断效果(预测效果)进行评价,本案例中ROC曲线下面积为72.12%,显示模型的预测效果尚可但并非最优。
二、Box-Jenkins模型
所谓时间序列是按时间顺序所排列的数据集合。时间序列中常含有确定因素和不规则因素,确定因素指时间序列中含有长期趋势、季节变化等,不规则因素则主要指由随机因素引起的偶然变化。时间序列分析建模,就是根据观察指标数据的变化规律来建立适合的数学模型,从而推测指标未来的变化趋势。
Box-Jenkins模型是一种精确度较高的短期预测时间序列模型。它将预测对象随时间变化形成的序列看作是一个随机序列,并呈现一定的规律性,可以用数学模型近似描述。在医学卫生领域中,传染病的发病会受到许多因素的影响,而且影响因素之间又存在着错综复杂的联系,很难运用结构式的因果模型加以解释;而数据之间的这种相互依存关系恰恰是研究对象最重要和最有用的特性,这时根据其自身的变动规律建立时间序列的动态模型则是一种行之有效的方法。Box-Jenkins模型是具有代表性的时间序列分析和预测方法。可以揭示研究对象与其他对象随着时间的发展变化其数量关系的变化规律,此模型应用较其他模型更为广泛,也能够取得很好的预测效果。
(一)Box-Jenkins模型方法
20世纪60年代,美国的Box及英国的Jenkins给出了他们的研究成果:时间序列分析、预测及控制方法,该种方法被称为Box-Jenkins建模方法,也称为B-J模型或ARMA(auto regression moving average)模型。基本的ARMA模型有三种:自回归模型(autoregressive model),简称AR模型;移动平均模型(moving average model),简称MA模型;自回归移动平均模型简称ARMA模型。
ARMA(p,q)模型的表达式为:
φ p为滞后项的系数,模型自回归的阶数为p,θ q为滞后的随机误差项的系数,模型移动平均的阶数为q;{ε t}为零均值白噪声序列,条件E(X sε t)=0,对任意的s<t说明当期的随机干扰与过去的序列值无关。
ARMA模型是一种比较成熟的模型,适用于短期预测,建模的数据要求是平稳时间序列,如果在解决实际问题的预测中,所用数据非平稳,就不能直接用RMA模型建模。这种情况下需要对数据做差分处理,再用ARIMA(p,d,q)模型建模。d为把非平稳的时间序列转化成平稳时间序列时对原时间序列进行差分的次数。ARIMA(p,d,q)的数学表达式为:
X t表示t时刻的序列,B是后移算子 =1-B,p,d,q分别表示自回归阶数、差分阶数和移动平均阶数。
SARIMA模型是ARIMA模型的推广,对于包含有季节性和趋势性的非平稳序列,如果可以通过逐期差分和季节差分使序列平稳化,就需要运用更复杂的模型描述,这就是SARIMA(p,d,q)(P,D,Q) s模型。该模型的一般形式为:
Φ( B s) φ( B)(1- B) d(1- B s) D X t=Θ( B s) θ( B)ε t,
φ( B)=1- φ 1 B- φ 2 B 2-…- φ p B p,
θ( B)=1- θ 1B-θ 2B 2-…- θ q B q,
Φ( B s)=1- φ s,1B s-φ s,2 B 2s-…- φ s , P B Ps,
Θ(B s)=1- θ s,1B s-θ s,2 B 2s-…- θ s , Q B Qs,
其中X t是一个t时刻非稳定的时间序列,ε t是一个白噪声,d是一般差分次数,D是季节差分次数,s代表周期,当随机事件的发展变化随时间表现出季节性时,若以月为单位,则s=12。p为一般模型的自回归阶数,P为季节模型的自回归阶数,q为一般模型的移动平均阶数,Q为季节模型的移动平均阶数,B是后移算子(BX t=X t-1)。
SARIMA模型的建模步骤为:
第一步:序列平稳化
将数据分为两部分,一部分用来建立模型,另一部分用来检验模型预测效果。建模前要检验数据的平稳性,如果不平稳,可通过差分将非平稳的时间序列转化为平稳的时间序列,可用增广的迪基-富勒(ADF)检验法检验序列是否平稳。
第二步:模型识别
通过观察平稳序列的自相关函数(autocorrelation function,ACF)和偏相关函数(partial autocorrelation function,PACF)图,确定SARIMA(p,d,q)(P,D,Q) s模型中可能的p、q、P及Q值,从而确定可能的模型。
第三步:参数估计与模型诊断
参数估计采用非条件最小二乘,选取各参数有统计意义的预测模型,对其进行残差检验,残差的Box-Jenkins Q检验的pro值大于0.05,则残差序列通过白噪声检验,模型是适应的。通过赤池信息准则(akaike information criterion,AIC)及施瓦茨准则(schwaz criterion,SC)确定拟合效果最好的模型(AIC与SC越小,模型拟合优度越高)
第四步:预测应用
用所建模型进行预测,并根据预测值与实际值的平均绝对误差(mean absolute error,MAE)、平均绝对百分比误差(mean absolute percentage error,MAPE)及均方根误差(root mean square error,RMSE)来衡量模型的预测精度,三个指标值越小说明模型预测精度越高。
其中AIC准则为:
AIC越小模型拟合效果越好。
SC准则为:
SC也是越小模型拟合效果越好。
MAE的计算公式为:
MAPE的计算公式为
RMSE的计算公式为:
其中 为预测值,X t是实际观测值,n是观测值的个数。
(二)Box-Jenkins模型实例分析
某地2004年1月至2014年6月肺结核(tuberculosis,TB)月发病率见图4-1。
图4-1 某地2004年1月至2014年6月肺结核(TB)月发病率
月发病率数据分为两部分,2004年1月至2013年12月的数据用来建立模型,2014年1—6月数据用来检验模型预测效果。对建模部分数据做平稳性ADF检验,结果 P>0.5,显示数据非平稳,于是对建模部分数据做一阶逐期差分,则建模中d=1,对差分后的数据做自相关ACF及偏相关PACF图发现数据还存在12个月的季节周期,继续做周期为12的一阶季节差分,则建模中D=1,两次差分后的数据再做ADF检验,结果如表4-5。
表4-5 建模部分数据一阶普通差分一阶季节差分后数据的单位跟ADF检验
由ADF检验 P值0.000 1远小于0.05,可以看出,两次差分后数据已平稳,在此基础上再做自相关及偏相关函数图,如图4-2。
图4-2 自相关及偏相关函数图
由自相关及偏相关图分析时间序列的相关性,为了考虑更全面,建立了九个可能的模型ARIMA(1,1,1)(1,1,1) 12、ARIMA(1,1,2)(1,1,1) 12、ARIMA(1,1,3)(1,1,1) 12、ARIMA(2,1,1)(1,1,1) 12、ARIMA(2,1,2)(1,1,1) 12、ARIMA(2,1,3)(1,1,1) 12、ARIMA(3,1,1)(1,1,1) 12、ARIMA(3,1,2)(1,1,1) 12、ARIMA(3,1,3)(1,1,1) 12.其中,通过残差序列白噪声检验的有6个模型,见表4-6。表中给出了相应的AIC和SC的值,AIC及SC越小模型拟合效果越好,比较以下6个模型,ARIMA(1,1,2)(1,1,1) 12有最小的AIC和次小的SC值,因此最终选择该模型去拟合建模数据。
表4-6 6种ARIMA模型及对应的AIC和SC值
续表
拟合图,如图4-3(由于做了两次差分,损失了部分数据,所以拟合部分数据只有2006年3月至2013年12月的拟合数据):
图4-3 拟合图
由图可看出该模型拟合效果较好,用该模型对2014年1—6月发病率数据做预测,计算得RMSE是2.58,MAE是2.14,MAPE是9.51,这三个指标都较小,表明SARIMA模型可用于该地肺结核发病率预测,预测效果较好(以上数据分析及绘图采用MATLAB2015b及Eviews7.2软件)。
任何一种预测模型都是对实际情况的抽象,不免存在局限性及不完备性,B-J模型方法一样也有缺点,由于卫生资源、社会条件等在不断变化,该模型方法长期预测性能会下降。
三、灰色系统理论
灰色系统理论是由我国学者邓聚龙教授于1982年创立。所谓灰色性是指一个系统的层次和结构关系模糊,具有随机的动态变化性,指标数据也是不全的或不确定的。具有灰色性的系统被称为灰色系统。
灰色系统理论重点研究“小样本,贫信息”的不确定性问题及“外延明确,内涵不明确”的对象,该类问题是概率统计与模糊数学中难以解决的。不同于经典的统计学,灰色预测模型在实验观测数据及其分布方面没有特殊的要求和限制,也不需要大量的历史资料数据,甚至4个数据即可建模预测,较好地解决了一些常见的预测方法难以解决的问题。
由于灰色模型拥有不需要太多数据、运算方便、原理简单以及可检验等许多优点,所以自提出以来,应用范围已拓展到众多科学领域,并取得了显著成果。近年来,灰色系统理论在传染病预测分析中也有一定应用,并取得了较好的预测效果。
灰色系统理论在预防医学中的应用集中于几个传统的灰色预测模型、GM(1,1)、GM(1,N)等。
(一)灰色预测模型
灰色数列预测模型的相关理论和建模方法:在建模前,先对原始数据进行建模可行性检验和生成处理,通过累加生成或相关生成强化序列规律性,弱化数据随机性,使系统灰色量白化并呈现一定的规律性,然后建立动态微分方程模型并做出预报。
1.灰色系统理论的主要建模预测步骤
数据检验与预处理、模型建立、模型检验几个步骤。
(1)数据检验与预处理:
首先,为了验证模型的适用性与可行性,需要对已知数据列做必要的检验。准光滑性检验、准指数律检验和级比检验是常用的三种检验方法。
设参考数据为 x (0)=( x (0)(1), x (0)(2),…, x (0)( n)), x (0)的一阶累加生成序列为 x (1)=( x (1)(1), x (1)(2),…, x (1)( n))。
准光滑性检验: ;
准指数规律检验: ;
级比检验 。
如果所有的级比σ( k)都落在可容覆盖( )内,则数列 x (0)可以作为模型GM(1,1)的数据进行灰色预测。否则,对原始数列 x (0)应用数据变换处理,通过数据变换使其落入可容覆盖( )内。常用的变换方法为平移变换,选取合适的常数 c,对 x (0)作平移变换,即令 y (0)( k)= x (0)( k)+ c, k=1,2…, n,使数列 y (0)=( y (0)(1), y (0)(2), …, y (0)( n))的级比 。
三种检验关系为:
因此这三种检验方法是相通的。
(2)模型建立:
基于累加生成处理后的数据,构建灰色预测模型如下:
其中, M为微分方程的阶, N为方程中变量个数,以上方程也被称为灰色动力学模型,通常记作GM(M,N)。模型参数可由最小二乘法估计。常用的灰色预测模型有:GM(1,1)、GM(2,1)、GM(0,N)、GM(1,N)。
(3)模型检验:
对生成的模型可信度检验,常用的方法有残差检验,级比偏差检验,灰色关联度检验和后验差检验。
1)残差检验:
记残差序列为 E=( e(1), e(2),…, e( n)),
其中 。
计算相对误差序列: Rel=( ε 1, ε 2,…, ε n),其中 。
当 ε k<0.2时,可认为模型达到一般要求;当 ε k<0.1时,则认为模型达到较高要求。
2)级比偏差值检验:
首先根据参考数据 x (0)( k-1), x (0)( k)计算出级比σ( k),然后由发展系数 a求出相应的级比偏差
若 λ( k)<0.2,则可认为模型达到一般要求;若 λ( k)<0.1,则认为模型达到较高要求。
3)灰色关联度检验:
设 x (0)=( x (0)(1), x (0)(2),…, x (0)( n))为系统参考序列, =( (1), (2),…, ( n))为 x (0)的模拟序列,定义 x (0)与 的灰色绝对关联度为 ,
4)后验差检验:
设初始数据序列 x (0)及残差序列 E的方差分别为 和 ,计算后验差比值
和小误差概率 ,
其中 m为满足 的例数。对于给定的 P 0>0,当 P> P 0时,称模型为小概率合格模型。可以参照下表4-7判断模型预测精度。
表4-7 模型精度评价标准
2.外推预测
若模型预测精度较好,直接应用预测模型外推预测。
灰色GM(1,1)模型是灰色预测模型GM(M,N)中应用最为广泛的一个模型,模型构建简单,对样本量要求低,4个数据即可建模。GM(1,1)模型通过构建一个一阶微分方程描述系统行为特征,其应用过程虽含有大量复杂的矩阵代数,但却可在Matlab环境下通过编写简单代码轻易实现。
模型包括三个基本组成部分:累加生成、逆累加生成、灰色建模。主要预测步骤如下:
(1)步骤一:
为弱化原始序列随机性,增强其规律性,对原始序列进行一阶累加生成处理,生成离散光滑序列。
设x (0)为非负原始序列:x (0)=(x (0)(1),x (0)(2),…,x (0)(n)),其一阶累加生成序列为:x (1)=(x (1)(1),x (1)(2),…,x (1)(n)),其中 。
(2)步骤二:
基于生成数据建立GM(1,1)预测模型,计算模型参数(常用最小二乘法拟合)。
累加生成处理后的原始数列,弱化了系统中“坏”数据即随机因素的干扰,使处理后的数据呈现指数增长规律。根据累加生成序列x (1)建立如下一阶线性微分方程模型
记模型参数列 为 = [a,b] T,记一阶累加生成序列x (1)的紧邻均值生成序列z (1)=(z (1)(2),z (1)(3),…,z (1)(n)),其中 ,
且
其中,a称为灰色模型的发展系数,反映了原始数列x (0)及其一阶累加生成序列x (1)的发展趋势;b称为模型的协调系数,反映系统数据间的变化关系。由最小二乘法计算模型参数
其中, 。
(3)步骤三:
根据建立的GM(1,1)模型得到预测表达式。
由以上结果,可得到微分方程 的解为 。
所以其时间响应式为 ,
其中c为未知积分常数。常用的确定c的方法有:
1) (假定x (1)(1)=x (0)(1));
2) (假定 x (1)(1)= x (0)( n));
3) (最小二乘法)。
经典的GM(1,1)采用 ,即假定x (1)(1)=x (0)(1),则 ,
其离散形式为 。
其还原模型为 。
当GM(1,1)模型的发展系数| a|≥2时,该模型无意义,即当| a|<2时,模型才有意义。模型预测效果随着 a取值的不同也有所不同。要得到较高拟合精度,往往要求 。
(二)GM(1,1)模型实例分析
基于2014年至2018年5年的某院某病住院人数数据,见表4-8,建立GM(1,1)模型。
表4-8 2014年至2018年各年某病住院人数
x (0)=( x (0)(1), x (0)(2),…, x (0)( n))=(309,238,227,242,243), x (0)的一阶累加生成序列为:
x (1)=( x (1)(1), x (1)(2),…, x (1)( n))=(309,547,774,1 016,1 259),
采用这五年的某病住院人数数据,建立的GM(1,1)模型的参数 a=-0.013, b=227.54,由于| a| < 2,所以模型精度较高。采用该模型拟合2014年至2018年某病患者数据情况,见表4-9。
表4-9 2014年至2018年某病患者数据模型拟合结果
由此表可看出拟合残差小,平均相对误差为1.9,后验差比值C=0.17,说明模型拟合精度高,建立得好,模型能够较好地预测未来几年该院某病的发病人数。
采用该模型预测2017年至2021年帕金森患者人数情况,如表4-10。
表4-10 2017年至2021年帕金森患者模型预测人数
由表预测结果可以看出,未来五年,该院某病住院患者逐年小幅增加,住院人数的年预测,能够帮助该院合理配置医生及药品物资情况提供一定的科学参考。
类似于许多常用的预测方法,灰色系统预测方法主要从数据上获取信息探讨疾病的统计规律。它缺乏详细地讨论疾病发生、发展的各种影响因素,因此须综合各方面因素考虑预测结论;当采用小样本数据做预测时,有时由于数据的偏倚也会导致预测结果的偏离;随着时间的推移,模型的灰度越来越大,拟合误差也会随之增大,模型在实际应用时,不宜做长期预测,要不断加入新信息,剔除老信息,改进模型。
四、传染病动力学
传染病动力学是依据疾病的传播机制,依据影响疾病传播发展的因素,建立反映疾病的传播特征、发展过程和流行趋势的数学模型,该模型通过研究数学模型的基本形态、数值模拟,是理论流行病学的研究方法之一。近年来,随着传染病动力学研究的深入发展,越来越多的模型被应用于各种传染病。除去性病、艾滋病等一些特殊的传染病研究,大部分模型都是适用于各类传染病。根据模型选入的因素,探讨不同因素对于疾病传播、发展流行的影响,得到控制疾病传播的关键因素,为传染性疾病流行趋势的预测提供数量依据和理论支持。
随着传染病动力学模型在实践中的深入探索,在忽略种群人口变动的前提下,可用于描述一些病程较短的疾病。
依据疾病是否存在潜伏性,传染病动力学模型可分为SI模型、SIR模型、SIRS模型、SIS模型、SEIR模型和SEIRS模型。
SI模型适用于无疾病潜伏期,且患病后难以治愈的传染性疾病;SIS模型则适用于无疾病潜伏期,但患病后可以治愈的传染性疾病;SIR模型适用于无疾病潜伏期,但传染病患者康复后可获得终身免疫性的传染性疾病。SIRS模型则适用于无疾病潜伏期,但传染病患者康复后只能获得暂时免疫力的传染性疾病,传染病感染者在前期愈后获得暂时性免疫力,一段时间后会有部分痊愈者失去免疫力,而再次变成易感者而染病。SEIR模型则适用于存在疾病潜伏期,且染病者痊愈后可获得永久性免疫力的传染性疾病。SEIRS模型通常存在疾病潜伏期,且染病者愈后仅获得暂时性免疫力的传染性疾病,康复者在经过一段时间后会失去免疫力而再次因感染患病。
若考虑到出生率、死亡率和人口流动等种群动力学变化,考虑是否存在先天性免疫力以及疾病是否存在垂直传染的事宜,动力学模型还可添加出生率、死亡率或人口流动等这些种群动力学因素,建立更加复杂的传染病动力学模型,以期全面细致揭示疾病的传播机制,描述疾病流行过程,为预测传染性疾病流行趋势提供进一步的理论依据和可靠基础。
1926年,为了研究1965—1966年伦敦黑死病及1906年孟买瘟疫的流行规律,Kermack与McKendrick构建了SIR仓室模型。在1932年又提出了SIS仓室模型。并且在仓室模型的基础上,提出“阈值理论”,应用该理论判断疾病是否会在某地区流行。
(一)Kermack-McKendrick的SIR仓室模型
长期以来,传染病动力学模型是以“仓室”模型为主。该模型于1927年由Kermack与McKendrick首次提出,到目前为止,“仓室”模型的理论依然不断发展并广泛应用。
就某地区某种传染病而言,所谓的SIR仓室模型就是依据该地区全部人群分为三类,每一类代表一仓室,即共有三个仓室。
其中,易感者类(susceptibles),数量记为 S( t),表示 t时刻,具有易感性且未染病者的人数;传染者(infectives)类,数量记为 I( t),表示 t时刻具有传染力且已被感染为患者的人数;和移出者类(removed),数量记为 R( t),表示 t时刻已从染病者类移出的人数。
假定该地区总人口为 N( t),在 t时刻,则有 N( t)= S( t)+ I( t)+ R( t)。Kermack-McKendrick的SIR仓室模型的建立,是依据以下三个假设为前提条件:
考虑研究环境的开放性。即不考虑人口出生率、死亡率和人口流动等动力学因素。忽略出生和死亡对疾病产生的影响,假定环境中的种群总数固定不变,是唯一常数: N( t)= k或 S( t)+ I( t)+ R( t)= K( 2)。
每一个患者与易感者接触,且具有传染力。假定传染病患者与易感者接触后就一定能被传染。假定在 t时刻的单位时间里,该环境内易感者总数 S( t)与一个传染病患者能使易感者患病的数量呈正比,称为感染率,也就是比例系数为 β,故, t时刻单位时间内所有传染病患者使易感者染病的总数,即新发病的人数为 βS( t) I( t)。
t时刻,单位时间内移出者与患者数量呈正比。单位时间内移出者在所有染病者中所占的比例,称其为移出率系数,为了使用方便本研究也将其简称为移出率,通用 γ表示。若从传染病染病者中移出部分仅为康复者,没有死亡等其他个体,移出率 γ又可称为恢复率系数或恢复率。在 t时刻单位时间内,传染病患者数量与移出者数量呈正比,移出率比例系数为 γ, t时刻单位时间内移出者的数量为 γI( t)。
基于上述三个基本假定条件,绘制传染病传播框架图,如图4-4所示,描述传染病易感者从染病到移出的状态变化过程。
图4-4 传染病易感者染病到移出的变化过程(愈后具有免疫力)
建立每一仓室人口变化率的平衡方程式,得到如下模型:
上述这一模型一般用于通过病毒为介质传播的疾病,疾病患者康复后,机体会获得免疫力,如流感、麻疹和水痘等。
(二)Kermack-McKendrick的SIS仓室模型
若传播介质为细菌,且疾病患者痊愈后不具有免疫力,可再次被感染的疾病,如淋病、脑炎等,则适用于以下SIS模型。SIS模型是于1932年,由Kermack与McKendrick提出,适用于描述细菌性染病者,且痊愈后不具有免疫力的疾病传播过程,模型的传播框架图,见图4-5所示。
图4-5 传染病易感者染病到移出的变化过程(愈后不具有免疫力)
假定传染病染病者在愈后不具有免疫力,会再次成为易感者,若该模型的假设与SIR模型相同,SIS模型则为:
假设环境中的总人数为 N,且传染性疾病是通过人与人之间接触传播,传染性疾病染病者在单位时间内与其他人接触的次数为 n,接触频率为( n/ N)。假设每次染病者通过接触传播疾病的概率为 β 0, β 0被称为有效接触概率。接触者若为易感者,则会造成疾病的传播流行。
基本再生数是指在疾病发病初期,种群中存在易感者,一个染病者在其患病期间内能够感染易感者的人数,称为基本再生数,通常用 R 0表示,是区分传染性疾病是否流行的阈值。若 R 0<1,表示一个染病者在其患病期间内能够传染的总人数小于1人,说明该疾病会慢慢消亡;若 R 0>1,则表示该疾病将会持续存在,并成为一种地方病。
不同模型有着不同的适用范围和优势,其预测精度往往也是不同的。如果简单地将预测误差较大的一些方法舍弃掉,将会丢失一些有用的信息。若将这些方法进行组合,可增大信息量,能更好地进行预测。组合预测将各种预测效果进行总体性综合考虑,比单个预测模型更系统、更全面。组合预测在国外称为Combination forecasting,Combined forecasting或Combination of forecasts等,在国内也称为综合预测、结合预测或复合预测。组合的主要目的就是当单一的预测方法无法预测或达不到预期效果时,较大限度地综合利用各种方法所提供的信息,尽可能地提高预测精度。流行的组合模型有:ARIMA和神经网络模型组合、灰色模型与神经网络组合等。
除了以上的常用预测模型以外,还有很多其他的模型,例如:logistic回归流行潜势模型,马尔可夫数学模型等。各个预测方法有其特点,不同传播途径的传染病有着不同的流行特征,同一种传染病由于受自然环境和社会因素的影响,在不同地区的流行趋势也不尽一样,对一组随机变量的动态过程同时进行观测,并将其作为整体加以研究,根据资料数据的特点选用适合的方法,并用多种方法进行预测,能寻求到一个基于这些单一模型而又博采众长的模型,或者尝试着将多种方法综合起来进行分析,势必能更系统、更全面地反映动态现象的内部规律性和未来趋势。但是,模型是实际原型的模型,而不是原型本身,决定事物未来发展的因素众多。有时由于所掌握的资料不够准确和完整或所用方法不够恰当,会影响预测的准确性。利用各种数学统计模型探讨某个实际问题时,关键是找到符合实际的好模型。随着对模型研究的深入,会有更多的预测经验上升到理论并模型化,会有更多的预测方法渗透到对传染病未来流行趋势的预测中,为更好地防制传染病和制定各项卫生管理决策服务。
传染病的发生发展是一种复杂的现象,对其进行准确可靠的预测,需要不断纳入新方法、新知识,不断提高传染病模型的模拟精度,使流行强度与传播区域预测更加准确。
(董柏青 谢艺红 马金凤 朱金辉 郑彦玲)