《1 前言》

1 前言

在SARS流行初期, 即从宏观上研究其发生、发展、演化规律并进行预测分析, 对于防治传染性疾病大规模流行的战略决策、稳定社会秩序具有重要意义。笔者基于生物学领域生长预测模型, 以广东、北京、全国SARS流行数据为例, 初步分析了SARS流行的一般规律, 并进行了预测分析, 供有关方面参考、探讨。

《2 特征参量辨识》

2 特征参量辨识

特征参量主要指此次SARS感染累计人数、爆发峰值时间和终止时间、爆发峰值时刻的总感染人数及当日感染人数。这些都是系统发展演化的特征参量, 对于局势判断极为重要。

广义的Logistic生长模型可以表示为[1],

dΝdt=rΝα[1-(ΝΚ)β]γ(1)

式中, r是增殖系数;α, β, γ是正的实常数;N是状态变量, 为时间的函数;K是饱和值。此模型表现形式一般为N-t空间中呈现不同形式的S型曲线, 对应初始缓慢增长、加速、减速和稳定四个阶段, 反映了一般生命历程, 但是不存在解析解。

一般地, 随系统演化内部驱动机制不同, 曲线拐点位置会不同, 各阶段持续时间、速率也不同。其特征参量中峰值爆发状态变量N (t) inf (在拐点上) 为

Νinf=Κ(1+βγα)-(1/β)(2)

式 (2) 只有在大于Ν*=(1+βγα-1)-(1/β)Κ时才有意义。

最大增长速率为 (dN/dt) max:

(dΝdt)max=(dΝdt)Νinf=rΚα(αα+βγ)α/β(βγα+βγ)γ(3)

爆发峰值时间 (也在拐点上) tinf

tinf=1rΚα-1[(1+βγα)(α-1)/β1-α+(1+βγα)(α-1-β)/β1-α+β+γ(γ+1)2!(1+βγα)(α-1-2β)/β1-α+2β+]-1r[Ν01-α1-α+γΝ01-α+βΚβ(1-α+β)+γ(γ+1)2!Ν1-α+2βΚ2β(1-α+2β)+]α1(4)

针对不同的生命演化系统, 相应于式 (1) 存在多个变种表达[2,3]

首先进行数据预处理:

1) 运用线性内插的方法, 将北京4月19日至4月25日公布的截至每日晚8时的SARS累计感染人数, 转换为截至每日上午10时的等效数据, 以便于与以后所有公布的每日上午10时数据一致。

2) 考虑到每日报告新增病例在随后公布的数据中又有排除, 因此一律将排除病例数从最近的累计报告病例数中减掉, 并由累计病例数推算每日实际新增病例数, 以保证数据准确性和累计数据的非降性。

基于可以利用的并且首先进入稳定期的广东, 以及将进入稳定期的北京、全国每日累计确诊SARS人数序列, 运用优化算法, 设定残差平方和最小为目标函数, 特征参量r, α, β, γ, K为控制变量, 约束条件为K>N (t) , N (t) 为最近一次累计感染人数, 按式 (1) 进行参数识别。

模拟结果均表明, 无论是地区性还是全国性的SARS流行都很好地遵从生长曲线 (表1、图1、图2) 。但是, 由于广东和全国公布的数据为峰值爆发时间以后的序列, 因此由这些短序列 (特别对于广州) 还不能预测到全部的特征参量, 尤其是对应于式 (2) 、式 (3) 和式 (4) 的特征量。对于北京, 采用已公布的所有数据序列 (2003年4月19日—2003年6月12日) , 采用广义的Logistic生长模型进行拟合, 所得特征参量为N (t) inf=953, (dN/dt) max=118, tinf=10.08。结合表1, 即说明最终的感染人数将为2 527人, 峰值爆发时间是4月28日, 峰值爆发时的累计感染人数为953人, 其时每日最大感染人数为118人。

表1 广东、北京、全国SARS广义的Logistic生长模型特征参量辨识

Table 1 The predicted results of characteristic parameters of the generalized Logistic growth model for Guangdong, Beijing and Mainland China 2003 SARS

《表1》


类别

特征参量
残差平方和数据源

r
αβγK*

广东
4.65E-103.26191 491.5992115 614.45781 516198.83122003-04-26~06-12

北京
56.05710.11373.64891.06122 5277 668.02752003-04-19~06-12

全国
486.1013-0.125512.82551.63625 3554 482.88162003-04-26~06-12

*截至2003-06-12广东、北京、全国实际累计感染SRAS人数分别为:1511, 2523, 5328

应当说明的是, 上述模拟结果都应当视为期望值, 并且对于北京和全国SARS而言, 其误差来源主要是由4月26日以前的最初几日统计数字准确度较差造成的。相对来说, 由于笔者没有全部掌握广东前期SARS流行数据, 其计算结果不够理想, 但就饱和值来说, 还是可以的。

《3 过程预测》

3 过程预测

上述广义的Logistic生长模型, 由于不能通过积分给出N (t) 的解析解, 因此尽管它能较精确描述与自然生命历程有关的系统演化特征, 但是不能直接用于过程预测。

作为广义的Logistic生长模型的特例, 历史上先后发展了多种预测模型, 例如早期的指数模型、Verhulst-Pearl模型、Gompertz模型、Von Bertalanffy模型、Richards模型、Weibull模型以及Morgan-Mercer-Flodin模型等。对于一个事件演化的预测, 最终运用或发展何种预测模型, 最好的办法是首先对宏观趋势作估计, 再经过各主要模型模拟, 兼顾选择残差平方和最小者作为过程预测模型, 但是始终不要认为残差平方和最小者一定是最好的预测模型。

《图1》

图1 SARS广义的Logistic生长模型 (北京)

图1 SARS广义的Logistic生长模型 (北京)   

《图2》

图2 SARS广义的Logistic生长模型 (全国) Fig.2 A realization of the generalized Logistic growth model of Mainland China 2003 SARS

图2 SARS广义的Logistic生长模型 (全国) Fig.2 A realization of the generalized Logistic growth model of Mainland China 2003 SARS  

经比较, 采用Gompertz函数[4]进行过程预测效果较好。Gompertz函数是典型的S型曲线, 对应于生命历程同样有4个阶段:初始缓慢增长阶段、加速阶段、减速阶段、稳定阶段。超Gompertz函数或广义Gompertz函数为

dΝdt=limβ0rΝΚβγ(Κβ-Νββ)γ=rΝ[ln(ΚΝ)]γ(5)

γ=1时, 上式即是经典的Gompertz函数, 其解为

Ν(t)=Κ(Ν0Κ)e-rt(6)

经典的Gompertz函数也可以由如下的Richards增长方程, 通过右边除以β, 接着取β→0的极限值得到

dΝdt=rΝ[1-(ΝΚ)β](7)

式 (6) 拐点处的N (t) inftinf值分别为

Νinf=Κe-1(8)tinf=ln[ln(ΚΝ0)]r(9)

有时, 经典的Gompertz函数也可用下列代数式直接给出:

Ν(t)=Κexp(1-aexp(-bt))(10)

式中b为群集系数。

此外, 笔者也分别运用了如下Pearl模型、Von Bertalanffy模型、修正指数模型:

Ν(t)=Κ1+aexp(-bt)(11)Ν(t)=Κ{1-13exp[-a(t-b)]}3(12)Ν(t)=Κ+abt(13)

上述各式特征参量中饱和值K的确定可采用两种方法:

1) 首先由广义的Logistic生长模型, 经优化处理得到 (约束条件:K>N (t) , N (t) 为最近一次累计感染人数) , 再代入上述各式, 再次用无约束优化算法, 求得参量a, b, c等。

2) 直接运用各模型函数, 视所有模型参量K, a, b, c等为变数, 经约束优化求解得到, 其中约束条件为K>N (t) , N (t) 为最近一次累计感染人数。

注意, 所有优化算法中, 目标函数均为残差平方和最小。这两种方法所得过程模拟结果, 将因参量不同而不同, 但一般差别不会太大, 可以接受。k值也可以因为采取了不同防治策略而为时变[5]。亦即采用不同的外部干预形式, 事件演化的道路将极为不同。

上述各模型预测结果见表 2、图 3 (K不固定) , 结果显示经典Gompertz模型预测效果最好, 并且预测曲线趋于稳定由快到慢的顺序依次为:Pearl模型>Von Bertalanffy模型>Gompertz模型>修正指数模型。

表2 广东、北京、全国SARS几个生长模型模拟结果

Table 2 The predicted results of Guangdong, Beijing and Mainland China 2003 SARS using growth models

《表2》


生长模型
广东
(2003.4.26—6.12)
北京
(2003.4.26—6.12)
全国
(2003.4.26—6.12)

Pearl模型

a
0.13718.09823.4952

b
0.10740.19170.1544

K
1 529.10392 523.00005 404.1002


Von Bertalanffy模型

a
0.12470.13380.1276

b
-16.78656.7905-2.3612

K
1 516.31832 523.00005 385.6542


Gompertz模型

a
0.12772.94220.82524

b
0.10990.13880.12768

K
1 514.97142 537.03795 442.85


修正指数模型*

a
-180.1327-1844.3256-3192.1477

b
0.85180.89780.8946

K
1 516.47402 560.39345 404.8864

下面以经过预处理的北京SARS全部观察数据 (2003-04-19~2003-06-12) 为例, 介绍经典Gompertz模型过程预测结果。

首先运用广义的Logistic生长模型调优得到的饱和值K=2 527, 直接辨识经典Gompertz函数式 (6) 中的参量得r=0.1458, 其中N0=105, 此时残差平方和为36 404.1042, 并且如以四舍五入计, 6月11日事件将结束。另外的特征参量为:N (t) inf=929.61, tinf=7.94, 后者对应于爆发峰值时间4月26日。

《图3》

图3 SARS感染累计人数N (t) 实际值与各模型预测值比较

图3 SARS感染累计人数N (t) 实际值与各模型预测值比较  

Fig.3 Predicted population size of Beijing 2003 SARS using different models and optimization projects

如果对式 (10) 中的所有参量进行辨识, 如表2所示, 则得K=2 537, a=2.9422, b=0.1388, 此时残差平方和为22 529.6951, 并且如以四舍五入计, 6月13日事件将结束。

经典Gompertz函数模拟结果见表3, 可以看出与实际数据较吻合。

《4 结论》

4 结论

1) SARS感染累计的人数演化过程总体上可以划分为初始缓慢增长、加速、减速和稳定四个阶段, 遵从广义生长模型。

2) 就实际效果而言, 为了能运用具有准确性、一致性的数据进行较高精度预测, 卫生突发事件的统计数字及时、全面、准确地公布极其重要。

3) 根据最新发展情况, 适时更新预报非常重要 (如可以运用Bayesian Method以及耦合时变参数[5]) 。预测结果的准确性有赖于决不放松警惕, 不出现“反弹”。对于复杂的长期发展过程, 需要运用分段函数进行预测。

4) 生长模型预测中, 原始数据的预处理很重要;预测的峰值爆发点、饱和值、趋稳时间与实际事件过程的吻合程度, 可作为判断生长模型预测结果满意与否的标志。

5) 对于具有四阶段的生长过程, 为了做出较好的预测, 所需要的数据列的长度及所要求数据列覆盖哪些生长阶段仍需要做进一步研究。深入理解事件发展过程的物理机制, 是发展合适预测模式的关键, 如广义的流行病模型 (SIR) 是有前景的[5]。进一步工作应以疑似、确诊、排除 (包括排除、治愈、死亡) 三者协同演化为主线建立合适的模型。

  

表3 经典Gompertz函数模拟结果Table 3 The predicted results of Beijing 2003 SARS using Gompertz function  

《图4》

表3 经典Gompertz函数模拟结果Table 3 The predicted results of Beijing 2003 SARS using Gompertz function