實驗用家蠅壽命分布的兩個參數(shù)

實驗用家蠅壽命分布的兩個參數(shù)

家蠅壽命模型建立家蠶是中醫(yī)學和生物實驗中最常用的實驗動物之一。許多新藥物的效果首先通過對家蠶的實驗初步證實。在這些實驗中，有必要處理家蠶的生命數(shù)據(jù)。在文獻中，由于不知道家蠶生命分布的分布類型，所以有兩種常見的方法。其中之一是作為正態(tài)分布來處理的。而采用非參數(shù)方法進行處理。正態(tài)分布家蠶的生命數(shù)據(jù)，并武斷地進行處理。例如，在兩份長度測試中使用t-驗證（見），這是不可取的方法。t-驗證可以在大樣本中確保第一次錯誤的概率，但驗證的效率不一定得到保證。在小樣本中，甚至第一次錯誤的概率，而是以參數(shù)方法來確保每個人都失去的信息。任何一種學習方法都是從所有空白的地方開始的，這尤其是對像家蠶這樣的實驗不太可能少另一方面，從大量研究事實來看，動物、植物和材料的死亡必須基于一定的機制。因此，生命分布的類型相對固定，因此可以使用這種方法。我們直接基于數(shù)據(jù)來尋找家蠅壽命的分布類型.數(shù)據(jù)來源于某種藥物對家蠅壽命影響的實驗數(shù)據(jù),實驗分為對照組,處理組Ⅰ(低劑量用藥組)和處理組Ⅱ(高劑量用藥組).在每個組里將雌雄分別處理,因此共有6個組.本文的是這樣安排的,先用Kaplan-Meier方法給出了生存函數(shù)的非參數(shù)估計及相應的死亡力度函數(shù)的估計,然后據(jù)此導出家蠅的參數(shù)生存模型并對模型中的參數(shù)給出了生物學的解釋,之后討論刪失樣本下模型參數(shù)的估計方法,最后通過圖形直觀地描述了模型對數(shù)據(jù)的擬合狀況,文章的結(jié)尾是一些注釋.生存函數(shù)的描述設(shè)X是表示壽命的隨機變量,則相應的生存函數(shù)定義為Sx(t)=P(X>t).生存函數(shù)的估計采用Kaplan-Meier乘積限估計的方法(見例如).壽命試驗的數(shù)據(jù)記錄為:(X1,δ1),…,(Xn,δn),其中δk=0或1視xk是截尾壽命還是真實壽命而定.為得到估計,將數(shù)據(jù)(X1,δ1),…,(Xn,δn)整理為(t1,n1,d1),…,(tk,nk,dk),其中tj表示死亡時間(假設(shè)t1≤t2≤…≤tk),dj是在ti時刻死亡的個體數(shù)目,nj表示tj時刻面臨風險的個體數(shù)目(包括在tj時刻死亡的個體數(shù)目).下面的六個圖形是我們的六組實驗數(shù)據(jù)的生存函數(shù)的Kaplan-Meier估計(橫軸為時間,單位是天),圖中的“”表示在該處有刪失的數(shù)據(jù).從這些圖中看出,所得六組實驗的生存函數(shù)已經(jīng)具有基本的形狀,曲線的躍度較小,因此,我們可以認為,所取樣本量用于尋找生存函數(shù)的分布類型是合適的.除生存函數(shù)外,死亡力度是生存分析中一個具有基本重要性的概念,其定義為死亡力度在時間點(ti+ti+1)/2處的估計為:這一估計的直觀解釋是:ni個個體在時間區(qū)間(ti,ti+1)內(nèi)平均單位時間內(nèi)死亡的頻率.根據(jù)頻率估計概率的思想,這樣的估計是直觀合理的.下面是對數(shù)死亡力度的圖(橫軸是時間,單位是天),圖中“+”為對數(shù)死亡力度對時間的散點,直線是這些散點對時間擬合的回歸直線.為0時的gomprtz模型從圖2中可以很明顯地看出,死亡力度函數(shù)的對數(shù)可以用線性函數(shù)來擬合,即:根據(jù)知道,生存函數(shù)為其中λ=eb,這是正半軸上的Gompertz分布.當a=0時,理解為α趨于0時的極限,即指數(shù)分布S(t)=e-λx.此外,還應有α≥0,否則會有S(+∞)>0,這是不合理的.相應的密度函數(shù):.同樣,當α=0時,理解為α趨于0時的極限,即指數(shù)分布f(t)=λe-λx.很顯然,在Gompertz模型中,因此,在此模型中,參數(shù)λ是初始的死亡力度,家蠅的性別,羽化時的體重等因素都可能對λ的值發(fā)生影響;而參數(shù)α則是死亡力度在單位時間內(nèi)增加的百分比,也就是說,在Gompertz生存模型中,在單位時間內(nèi),個體的死亡力度以固定的百分比增加,死亡加速因子是常數(shù).我們在這里應該注意到,習慣上,Gompertz模型常寫成位置-刻度模型的形式.但是在討論半軸上的Gompertz模型時,這樣的寫法是不可取的,這會丟掉指數(shù)分布這一重要特例.換句話說,正半軸上的Gompertz模型不再具有位置-刻度的形式.刻失數(shù)據(jù)下gomprtz的二乘估計根據(jù)正半軸上Gompertz分布的兩個參數(shù)的生物學解釋,尋找這兩個參數(shù)的矩估計是不現(xiàn)實的,我們這里敘述刪失數(shù)據(jù)下Gompertz分布的極大似然估計和簡單最小二乘估計.(1)極大似然最大解關(guān)于Gompertz模型的極大似然估計,在文獻中有比較詳細的敘述,正半軸上Gompertz模型的極大似然估計與所述相似.為了后文說明及引用方便,我們將其簡要敘述如下.設(shè)樣本為X1,…,Xr,Xr+1,…,Xn,其中,前r個為觀察到的真實壽命,后n-r個為刪失壽命.這時模型的參數(shù)的對數(shù)似然函數(shù)為:相應的似然方程為:如果則似然方程的解為,=0.否則,從第一個似然方程得代入第二個似然方程,得:該方程的解只能利用數(shù)值方法得到(比如說Newton-Raphson方法).從該方程得到α的極大然似估計為,相應就得到λ的極大似然估計為關(guān)于極大似然估計,可以證明:(1)方程(3)的解存在且唯一;(2)>0的充分必要條件是(2)全實軸gomprtz模型的最小二乘估計正如上面所述,似然方程的求解只有通過計算機的幫助才能實現(xiàn),其次,似然方程的NewtonRaphson方法求解依賴于良好的迭代初值的選取,因此,發(fā)展一種簡單的估計參數(shù)的方法是必要的,這可以得到良好的用來求解似然方程的初值.在全實軸的Gompertz模型中,對生存函數(shù)取兩次對數(shù)就得到線性函數(shù),據(jù)此導出參數(shù)的最小二乘估計.當討論正半軸的Gompertz模型時,情形稍有不同,此時的最小二乘估計充分利用了對數(shù)死亡力度函數(shù)是時間的線性函數(shù)這一特性.根據(jù)最小二乘估計的公式有:其中,這樣,得到參數(shù)的簡單最小二乘估計為:下面的表給出了這兩種方法在六組數(shù)據(jù)下的數(shù)值結(jié)果.從表中可以看出,兩種方法估計的結(jié)果實際上相差不遠.參數(shù)模型估計的擬合情況下面給出的圖形表示了Kaplan-Meier估計與參數(shù)模型估計的擬合情況.在圖中,階梯函數(shù)是Kaplan-Meier估計的結(jié)果而光滑曲線是參數(shù)模型估計的結(jié)果,從圖中可以看出,參數(shù)模型與數(shù)據(jù)的擬合非常良好.mprtz生存模型與蘇氏最優(yōu)模型的關(guān)系從數(shù)據(jù)出發(fā),給出了實驗用家蠅生存函數(shù)的分布類型,正是正半軸上的Gompertz分布.Gompertz分布的重要性是勿容置疑的,比如就有這方面的論述.但一般說來,在生存分析中引入Gompertz分布主要是用來估計Weibull分布的參數(shù),而不是直接用來作為生存模型使用.我們在這里不但給出了正半軸的Gompertz生存模型的實例而且對其中的參數(shù)給出了生物學上的解釋.從這一點來看,正半軸的Gompertz分布與全實軸上的Gompertz分布還是有著些許的不同.此外,正半軸的Gompertz分布與Weibull分布又有著某些相似性,這就是,正半軸的Gompertz生存模型的死亡力度函數(shù)的對數(shù)是年齡的線性函數(shù),而Weibull生存模型的死亡力度函數(shù)的對數(shù)是對數(shù)年齡的線性函數(shù).從這一點出發(fā),假如說我們考慮這樣的生存模型類,即死亡力度的對數(shù)可以表示為α+bg(x)的情形,其中g(shù)(x)是已知函數(shù),α,b是未知參數(shù),則指數(shù)生存模型,正半軸上的Gompertz生存模型以及Weibull生存模型皆屬此類.在藥效實驗中,需要對家蠅施以一定劑量的藥物,以觀察該藥對家蠅壽命的影響,這實際上是個假設(shè)檢驗問題.根據(jù)本文前面對參數(shù)λ和α給出的生物學解釋,λ是生命初始時刻的死亡力度,這是不會被任何處理影響的量,因此所謂藥物效果僅僅是作用在參數(shù)α上,這樣來看,藥物的效果是影響了死亡加速因子α,因而相應的假設(shè)檢驗問題就只需在參數(shù)α上施行就行了.由于篇