一种种植病流行的时间动态模型

一种种植病流行的时间动态模型

近年来，生存分析理论和应用的研究相对活跃。许多基于生命分析的统计模型（如生命表方法和cox模型）也用于生物研究和应用。作者运用生存分析理论中生存函数的推导方法,导出了一个描述植病流行在时间动态上的新模型。该模型与以往若干植病流行模型的主要区别在于:一是其解析形式不同;二是其推导基于生存分析理论和植病流行的机理分析,具有一定的生物学意义。因此,该模型可用于描述一些植物病害引致病害流行的规律分析。1子菌素质中病原菌的累积与消除的关系植物病害是病原物和寄主植物通过寄主作用构成的生物系统。植病的流行,无论是轻微的还是严重的,爆发的还是迟延的,都是一个随着时间和空间而发展变化的动态过程。因而,植病流行模型大多以时间t为自变量,其流行规律主要指病情指数x(t)与时间t的一个关系式(确定性或非确定性的)。在一定条件下,病原物的传播体在着落于寄主植物体表后可引致植物发病,发病的过程实际上也是菌量的积累过程。但是,附着萌发的孢子不一定全部都能侵入,即使侵入了也不一定全部都能定殖而引致发病。也就是说,孢子也有死亡和侵染未成功的情形。这一点可以理解为在发病过程中,有一部分菌量被消除了(降雨等因素也可造成部分菌量的消除)。基于这一分析,下面来推导植物受到病害侵染的生存函数S(t)。设植株个体从0时刻开始受到侵染,到t时刻发病。菌量在其某一时间间隔(t′,t′+dt′)内所积累下来的量为δdt′,并按速率υ消除,至t时刻其残留量为δexp{-υ(t-t′)}dt′所以,在(0,t)时间间隔内,总的菌量积累为∫t0δexp{-υ(t-t′)}dt′=δυ[1-exp{-υt}]再设在t时刻的侵染率λ(t)与此刻植物所积累存活的菌量成正比,则有λ(t)=k⋅δυ[1-exp{-υt}]=βυ[1-exp{-υt}]‚其中β=kδ,故得生存函数为S(t)=exp{-∫t0λ(u)du}=exp{-βυ[t-1υ(1-exp{-υt})]}‚从而得累积发病率,即病情指数为x(t)=1-S(t)=1-exp{-βυ[t-1υ(1-exp{-υt})]},(1)其密度函数为dxdt=βυ{1-exp{-υt}exp{-βυ[t-1υ(1-exp{-υt})]}.(2)式(1)即为利用生存分析理论和菌量的积累与消除机制所导出的一个新的植病流行模型,式(2)即为病情指数x(t)所满足的微分方程。显然,病害的增长速度dx/dt与菌量的积累δ和其消除率υ有关,这一点在其它植病流行模型中是没有得到反映的。此外,该模型从其解析形式上也是以往的一些植病流行模型所不能替代的,从而为研究人员在选择合适的植病流行模型时提供了余地。2模型的参数估计和应用2.1性化法最小二乘估计先把式(1)等价地写成ln11-x=-βυ2+βυt+βυ2exp{-υt},记y=ln11-x,则对一组试验数据(t1,x1),(t2,x2),…,(tn,xn)可得到非线性回归模型yi=-βυ2+βυti+βυ2exp{-υti}+εi,(3)其中,误差εi满足E(εi)=0,Var(εi)=σ,E(εiεj)=0,i≠j,i,j=1,2,…,n。考虑到模型(3)无法线性化,而且一些常规的迭代法(如Gauss-Newton法和Levenberg-Marquart法等)又强烈地依赖于参数初值的选取。因此,本文先采用部分线性化方法求得待估参数的一个初始估计,然后以此构造参数搜索区域,再借助于正交表安排一个较小的参数代表点集网格,在此网格上计算残差平方和,取有最小残差平方和的那个估计再次作为新的参数初始值。如此继续下去,直至达到网格精度,求出待估参数的最小二乘估计。该法的主要步骤如下:(1)选取参数初始值θ0=(β0,υ0)T。由于在植病流行中有0≤t≤1,则取模型(3)的线性近似回归模型yi=-βυ2+βυti+εi,i=1,2,⋯,n(4)中θ=(β,υ)T的LSE作为θ0=(β0,υ0)T。若记式(4)的LSE为ˆA,ˆB,则有{ˆB=n∑i=1(ti-ˉt)(xi-ˉx)n∑i=1(ti-ˉt)2,ˉt=1nn∑i=1ti,ˉx=1nn∑i=1xi‚ˆA=ˉx-ˆBˉt‚即有{βυ=ˆB‚-βυ2=ˆA‚解得{υ0=-ˆBˆA‚β0=-ˆB2ˆA‚从而得到θ的初始值为θ0=(β0,υ0)Τ=(-ˆB2/ˆA,-ˆB/ˆA)Τ。(2)构造参数搜索区域Dh(θ0)。这里构造为Dh(θ0)=[β0/h,β0h]×[υ0/h,υ0h],其中“×”为集合的乘积。(3)在一定精度下,把每个参数分成t个水平,以此把搜索区域Dh(θ0)分成网格状,则有t2个网格点组成点集ˉDh(θ0)⊂Dh(θ0).然后,借助正交表从ˉDh(θ0)中选取有代表性的点集ˉD′h(θ0),在ˉD′h(θ0)上计算残差平方和S(θ),其中S(θ)=n∑i=1[yi-(-βυ2+βυti+βυ2exp{-υti})]2.若有S(θ′)=minS(θ)θ∈ˉD′h(θ0)且θ′达到网格精度,则停止搜索,有ˆθLSE=θ′;否则以θ′为新的初始值,在θ′附近加密网格进行第二轮计算,直至达到精度为止。此法相对于Guass-Newton等迭代法而言,对初值变化不敏感,故不必为获得优良的初值而作过多的考虑。另外,该法只需计算函数值而不需计算微商,具有计算简单等优点。2.2以0为中心的初步参数估计以合肥地区1994年水稻纹枯病流行的一组数据(0.42,0.1858),(0.54,0.2900),(0.60,0.3241),(0.67,0.4058),(0.72,0.4369)为例,按前述方法进行参数估计并建立植病流行模型。首先把数据组(ti,xi)转化成(ti,ln11-xi),即为(0.42,0.2056),(0.54,0.3425),(0.60,0.3917),(0.67,0.5205),(0.72,0.5743)。然后计算出估计值ˆA=-0.3267,ˆB=1.2434,于是得到参数θ的初始值θ0=(4.7323,3.8060)T。取h=1.8,得到参数搜索区域为Dh(θ0)=[4.7323/h,4.7323h]×[3.8060/h,3.8060h]=[2.6291,8.5181]×[2.1144,6.8510].再取t=7,把区间[β0/h,β0],[β0,β0/h],[υ0/h,υ0],[υ0,υ0h]各三等份,这样每个参数就在以θ0为中心分成了7个水平,得到下列参数水平表。使用正交表Ft2+12(t2)=F25(72)选取设计ˉD′h(θ0),其中F25(7)为F25(72)=(11131517222426⋯⋯71737577),然后在ˉD′h(θ0)上计算残差平方和S(θ),经计算,当ˆθ=(4.7323,3.8060)Τ时,有S(ˆθ)=minS(θ),θ∈ˉDh′(θ0)代入(1)中,得到该病害流行模型为x(t)=1-exp{-1.2434[t-0.2627(1-exp{-3.8060t})]}.该模型的理论值与实测值的拟合情况如下:根据表2计算得残差平