新传染病的Logistic模型研究

1、传染病的Logistic模型与统计软件R摘要：从理论上说明Logistic模型可以描述传染病的流行规律，用统计软件R实现对模型中3个参数的最小二乘估计。充分利用Logistic模型曲线的载荷信息，求解其速度函数的一阶和二阶导数，得到了传染病传播过程的渐增期，快增期和缓增期。实证分析表明Logistic 模型比一般文献中的SIR模型等更加实用： 模型中的参数容易确定，易得出定量结果，拟和精度高。 关键词：Logistic模型 ；增长率； 关键期： SARS； 最小二乘法 中图分类号：0241.81 MR.分类号：62M45 文献标识码：A Logistic Model of Infectious

2、 Disease and Statistics with RXU Fu-xia1.2 SHI Dao-ji1 DONG Yong-quan2(1. School of Science,Tanjin University,Tianjin 300072,China2.Math.Department of Tangshan Teacher College,Tangshan 063000,China)Abstract:This paper indicates in theory that it is practicable using logistic model to describe the po

3、p disease law, about 3 parameters of the model can be estimated with R. Make full use of the information of logistic curve, we obtained the gradually growth period, fast growth period, slow growth period by fist and second derivative of velocity funtion of increasing or growth process of logistic cu

4、rve, then solves an example with the model.It is emphasized that logistic model is more practical than SIR or some others because of the formers parameters is determined more easily to get the quantitative results whose fit error is minor.Key word:logistic model; growth factor ; key period; SARS; me

5、thod of the minimum squars 长期以来，建立传染病的数学模型来描述传染病的传播过程、分析受感染人数的变化规律、预报传染病高潮的到来等，一直是各国有关专家和官员关注的课题。一般的传染病数学模型都是以SI（人群分为易感染者Susceptible和已感染者Infective）模型、SIS（人群分为易感染者、已感染者和治愈后重新变为易感染者）模型、SIR（人群分为易感染者、已感染者和治愈后获得免疫力的移出者Removed ）模型为基础进行研究的，这也是经典的研究传染病模型的方法。这些方法的一个共同特点是：模型假设越来越接近实际，所用数学方法越来越深奥，模型的数学表达式越来越复杂

6、冗长，甚至变得不可求解（虽然具有理论意义），忽略了数学模型的实用性和应用的广泛性。例如SIR模型，其理论虽然精确，但模型的解析表达式无法求出，只能通过讨论解的性质来确定其中的参数，这给解决实际问题带来困难。Logistic分布是一个具有较大实用价值的连续型分布，最初起源于生物群体的Logistic增生曲线呈S型，在经济学、政治学、人口统计学和植物种群生态学等方面都有广泛的应用。文8用logistic模型分析了SARS的传播规律，但未从理论上阐明该模型适用于研究一般的传染病的传播规律，而且由于用logistic模型得到的结果比用纯粹的统计回归模型抛物线模型还要差（前者的残差标准误大），这使得人们

7、对用logistic模型描述传染病的传播规律产生怀疑。本文通过进一步的研究发现：用logistic模型拟合的误差大是因为参数估计不准确所致，为此我们应用最新的自由（免费）统计软件R，由原来估计两个参数改为对三个参数的逐步逼近，得到了拟合精度高于文8中抛物线模型的logistic模型。进一步充分利用Logistic模型曲线的载荷信息，通过求解其速度函数的一阶和二阶导数，得到了传染病传播过程的渐增期，快增期和缓增期，对传染病的logistic模型进行了全面深刻的理论分析和实证说明，为预测和控制传染病的蔓延提供了方法。本文的另一个目的在于通过对模型参数的估计，介绍免费的统计软件R的强大统计作图功能，

8、尤其适用于生物医学方面的统计工作，文9是这方面的专著。 1建立模型要预报未来若干天的传染病感染人数，最重要的影响因素是初始时刻的传染病例数和今后若干天的病例增长率（新感染病例率减去治愈率）。记初始时刻的病例数为，天后的病例数为，传染病增长率为(>0)，则预报公式为： <1>将上述离散型的模型<1>式连续化，就得到指数增长模型：记时刻的传染病例数为，当考察一个国家或一个很大的地区时，是较大的整数，可将视为连续可微函数，则到时间内传染病例数的增量为： <2> 于是，满足如下的微分方程： <3> 解之得： <4> 上式表明传染病例数将

9、按指数规律无限增长（），这不符合实际情况，因为随着传染病例数的增加，政府、医疗部门和居民等都不会无动于衷，任其泛滥，而是采取各种有效措施来预防和控制传染病的迅速蔓延，传染病增长率会随着传染病人数的增加而逐渐减少。将传染病增长率表示为传染病例数的函数，按照上面的分析应是的减函数，一个最简单的假设是是的线性函数： <5> 这里表示初始时刻的增长率，显然，对于任意的，增长率. 为了确定系数的意义，引入最大传染病例数，当时增长率应为零，即，代入<5>式得到，所以传染病增长率可表示为： <6> 其中,可以根据传染病统计数据确定，因子（）体现了对传染病例数增长的阻滞作用

10、。在<6>式的假设下， <3>式变为： <7>非线性微分方程<7>式可用分离变量法求解得： <8> <8>式即为所求的传染病Logistic数学模型。2模型分析 (1) 上述模型应用的关键是3个参数和的确定：表示最大传染病例数，它应该比传染病的实际统计数据中的最大值大一点，但具体取多大合适，只能凭经验或试验的方法估计；表示初始时刻的传染病例数，由于传染病的数据是在疾病爆发后统计出来的，初始数据不容易精确得到，当然第3个参数初始传染病增长率也就不可能得到准确的估计，文8采用了将取定为统计数据的第一个值，对其它两个参数用最小二

11、乘法估计，所以误差较大（还有一个原因是我们已经将离散的问题连续化，所以取定了一个离散值作为初始值，自然会造成误差的增大）。本文考虑到3个参数之间的关系，对它们实行同步最小二乘估计，减小了估计的误差，用统计软件R可以实现这个过程，祥见第3部分的讨论。(2) Logistic曲线的特点是开始增长缓慢，而在以后某一范围内迅速增长，达到某限度后，增长又缓慢下来，曲线略呈拉长的S型 。下面对方程<8>式进行解析，充分了解上述Logistic方程所载荷的信息，以期对传染病的传播规律进行深入研究。()求方程的一阶导数得到传染过程的速度函数： <9>再求方程<9>生长速度函

12、数的一阶导数，并令其等于0： <10>解之得： <11> 可见当时，传染病传播速度最快，为其高峰期。() 求方程<9>生长速度函数的二阶导数，并令其等于0： <12>求解得：， <13> 这是速度函数的两个拐点，加上高峰点，这样Logistic曲线有三个关键点，其横坐标分别是：， ， 这三个点对应者传染病传染过程的三个关键时期：始盛期，高峰期，盛末期，也可用速度函数的两个拐点将Logistic曲线的生长过程分为渐增期，快增期和缓增期通过以上分析就可以比较全面地了解和把握传染病传播的规律和过程，并给于科学的解释。3 实证分析以2003年

13、全国大学生数学建模竞赛A题为例，对附件2北京市疫情数据，建立形如<8>式的数学模型，并确定其中的参数，和，在这个模型中要求，令：第天的累计病例数 ：第天的死亡病例累计数：第天的治愈出院累计数：第天的实际病例数，（4月20日为起始日，5月15日实际病例数达到最大，为终止日，见附表）则有： <14> 净增长率： <15>下面用R语言实现对模型中3个参数的最小二乘估计，并对模型的拟合结果进行分析：(1) 将时间序号（附表中第2列，取名time）和实际病例数（附表中第6列，取名bingli）以文本文档chuan保存在D盘上，然后输入命令：> read.tabl

14、e("d:/chuan.txt",s=c("time", "bingli")->ex> attach(ex)> bingli.nls<-nls(binglixm/(1+(xm/x0-1)*exp(-r0*time),data=ex,start=c(xm=1991,x0=288,r0=0.2987),trace=T)上述命令中，估计的初始值取的是实际病例数中的最大值；取的是观察为起始日的实际病例数；取自文8的最后结论。事实上，只要我们给3个参数赋予合理的初始值，都可以收敛到最优估计，只不过收敛速度

15、有所不同。运行结果如下：1183751 : 1991.0000 288.0000 0.2987 349977.3 : 1974.5882029 338.2398554 0.1762658 14059.92 : 2051.0846865 298.7330217 0.2051997 13053.82 : 2066.0106971 302.3097350 0.2032687 13053.56 : 2066.0849606 302.2509151 0.2033012 13053.56 : 2066.0865266 302.2527073 0.2033004结果表明：经过5次逼近，得到3个参数的最小二乘

16、估计为，由此得到北京疫情数据的Logistic拟合模型如下： <16>(2) 将时间序号代入上式，可以求出传染病例数的拟合值（理论值），命令如下：> fitted(bingli.nls)拟合结果如下，将其填入附表中最后一列。 1 358.5671 422.8583 495.3330 575.8761 663.9798 758.6999 858.6574 8 962.0930 1066.9780 1171.1672 1272.5718 1369.3220 1459.8950 1543.190415 1618.5489 1685.7238 1744.8173 1796.1998

17、1840.4272 1878.1651 1910.126822 1937.0260 1959.5445 1978.3120 1993.8954 2006.7951(3) 回归的详细结果如下：> summary(bingli.nls)Formula: bingli xm/(1 + (xm/x0 - 1) * exp(-r0 * time)Parameters: Estimate Std. Error t value Pr(>|t|) xm 2.066e+03 1.519e+01 136.03 <2e-16 *x0 3.023e+02 9.107e+00 33.19 <2e

18、-16 *r0 2.033e-01 4.454e-03 45.64 <2e-16 *-Signif. codes: 0 '*' 0.001 '*' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 23.82 on 23 degrees of freedom上述结果表明：3个参数的拟合的误差均在0.001以下，拟合残差的标准误为23.82(4) 回归效果检验画残差的（Q-Q）诊断图> qqnorm(resid(bingli.nls)画出实际病例数散

19、点图并添加由<16>式确定的logistic理论曲线：> plot(time,bingli,lty=1,xlab="时间序号",ylab="病例数")> matlines(time,fitted(bingli.nls)残差的（Q-Q）图显示残差近似服从正态分布，logistic模型的理论曲线与实际数据图则显示模型的拟合误差小。结果表明通过回归检验，拟合效果显著。(5) SARS传播过程的三个关键点将， 代入式<11><13>得SARS传播过程的三个关键点为：传播的始盛期为：，即 4月22日传播的高峰期为：，

20、即 4月29日传播的盛末期为：即5月5日进一步可以得到SARS 的渐增期为4月20日22日，快增期为4月22日5月5日，缓增期为5月5日以后4结论用Logistic模型可以描述传染病的流行趋势，其理论严谨符合实际而且简单明了，用免费的统计软件R可以准确方便地实现模型中参数的最小二乘估计及回归效果的检验。进一步通过对Logistic模型曲线方程的解析，求出传染病传播过程的三个关键点，这样就可以比较全面地分析和把握传染病传播的规律和过程。实证分析表明模型中的参数比一般文献中的SIR模型等容易确定，易得出定量结果，拟合误差小。文章强调了数学模型的实用性和应用的广泛性，例证了一种应用数学的思想：在复杂

21、性和精确性中取折中，在一定的精确度要求下，所用的数学方法越简单越好。附表 （r0时的logistic阻滞增长模型）日期时间序号 已确诊病例累计死亡累计治愈出院累计实际病例数净增长率理论病例数4月20日13391833288358.56714月21日248225434140.4375422.85834月22日358828465140.2416495.33304月23日469335556030.1732575.87614月24日577439646710.1128663.97984月25日687742737620.1356 758.69994月26日798848768640.1339858.6574

22、4月27日8111456789800.1343962.09304月28日91199597810620.08371066.97804月29日101347668311980.12811171.16724月30日111440759012750.06431272.57185月1日1215538210013710.07531369.32205月2 日1316369110914360.04741459.89505月3日1417419611515300.06551543.19045月4日15180310011815850.03591618.54895月5日16189710312116730.05551685.72385月6日17196010713417190.02751744.81735月7日18204911014117980.04601796.19985月8日19213611215218720.04121840.42725月9日20217711416818950.01231878.16515月10日21222711617519360.02161910.12685月11日2222651201