楊 洋 ，宋向東 ，陸 瑤 ，劉麗靜

（1.燕山大學(xué)a.里仁學(xué)院基礎(chǔ)教學(xué)部；b.理學(xué)院統(tǒng)計學(xué)系，河北 秦皇島 066004；2.河北科技師范學(xué)院 數(shù)理系，河北 秦皇島 066004）

1 研究背景及統(tǒng)計框架

設(shè)A表示研究區(qū)域；Ai,i=1,…,n表示區(qū)域A上的一個分割；Yi和Ei分別表示各區(qū)域疾病發(fā)生或死亡事件的觀測值和期望值。對于稀有的和非傳染性疾病，傳統(tǒng)框架中通常假設(shè)：

Yi～Poisson(RiEi) (1)其中Ri表示區(qū)域Ai的相對風(fēng)險。區(qū)域Ai上相對風(fēng)險Ri的極大似然估計稱為標準死亡率(standardized mortality ratio簡稱為SMR)，SMRi=Yi/Ei。SMR對于區(qū)域總體相對風(fēng)險的估計很有參考價值，但是：

因此對于Ei較小的情況來說，相對風(fēng)險Ri的極大似然估計的方差很大，也就是說，對Ei較小的地區(qū)，相對風(fēng)險Ri的極大似然估計SMRi將有很高的不準確性，它可能會使面積大但人口小的區(qū)域的估計值偏高，從而掩蓋了真實區(qū)域的風(fēng)險格局。

我們提出了一種新的統(tǒng)計框架，并將其用于模擬相對風(fēng)險。假設(shè)面臨風(fēng)險的k階層人口數(shù)量服從一個泊松過程，λk(x)表示其強度，則發(fā)病人數(shù)也服從一個泊松過程，其強度為λk(x)×pk(x)，其中pk(x)表示面臨風(fēng)險的k階層人口在x水平處的發(fā)病率。我們再假設(shè)pk(x)=pk×Rk(x)，其中Rk(x)表示k階層人口在x水平處的相對風(fēng)險，pk表示k階層人口的標準參考發(fā)病率。

在通常的應(yīng)用中，可通過一些不太充分的數(shù)據(jù)獲取區(qū)域人口密度。通過身份證上相關(guān)信息，如性別、住址，我們可以構(gòu)造一個分段階層化的人口密度函數(shù)fik(x)，用以表示區(qū)域Ai上k階層人口的空間分布密度。于是對于每個Ai，有λk(x)=Nik×fik(x)，其中Nik表示區(qū)域Ai上k階層人口數(shù)量。通過Nik有Yik～binomial(Nik,pik)，其中對k求和得到：

我們定義區(qū)域Ai上k階層人口的平均相對風(fēng)險：

那么 Yi～Poisson(RiEi)，其中 Ri=ΣkwikRik，wik=Nikpk/Ei，Ei=ΣkNikpk。區(qū)域Ai上的相對風(fēng)險Ri是Rik的加權(quán)平均，權(quán)重wik為區(qū)域Ai上k階層人口所占比例。

由于數(shù)據(jù)的不充分，我們不可能估計得到k階層人口相對風(fēng)險，因此我們假設(shè)對所有的x，k，都有Rk(x)=R(x)，即我們將區(qū)域總體風(fēng)險視為一個連續(xù)的風(fēng)險曲面，且有：

其中fi(x)=Σkwikfik(x)是階層化人口密度的加權(quán)平均，權(quán)重wik為區(qū)域Ai上k階層人口所占比例。

我們將上面的統(tǒng)計框架與之前介紹的傳統(tǒng)疾病地圖進行比較，后者說明式(1)中的Ri代表區(qū)域Ai上的每個人共有的相對風(fēng)險，這要求空間相對風(fēng)險R(x)在整個Ai范圍之內(nèi)不存在空間差異，且與fik(x)的形式無關(guān)。而將區(qū)域總體風(fēng)險視為一個連續(xù)的風(fēng)險曲面有很多好處，尤其是它形象有清晰地描述了區(qū)域的相對風(fēng)險，準確地模擬R(x)能夠得到Ri和Rj間的協(xié)方差，i≠j，而且還可以進一步的分析R(x)，例如在研究點的附近再建立一個風(fēng)險曲面模型，從而研究相對風(fēng)險的近似情況。

2 空間統(tǒng)計模型

2.1 高斯隨機場模型

我們假設(shè)空間相對風(fēng)險R（x）是一個連續(xù)的隨機場，則利用，R(x)區(qū)域相對風(fēng)險Ri可由前面的式(4)得到。特別的，我們假設(shè)S(x)=logR(x)來自一個平穩(wěn)高斯隨機場(Gaussian Random Field,簡稱 GRF)，其均值為 α，臨界方差為 σ2，相關(guān)函數(shù)為γ(x,y)，則區(qū)域相對風(fēng)險Ri的均值和協(xié)方差可得：

由于：

我們只需計算E[RiRj]：

所以：

但R(x)的分布函數(shù)不能得到，空間統(tǒng)計的有關(guān)學(xué)者通常認為該分布近似對數(shù)正態(tài)分布，特別是當(dāng)研究區(qū)域Ai相對較小時，近似情況也較好。在這種近似分析下：

服從多元正態(tài)分布，其均值、方差分別為：

(見文獻[1])，當(dāng)區(qū)域范圍逐漸減小，我們?nèi)O限結(jié)果，進一步近似得到E[Si]≈α

其中γ(Ai，Aj)是分別在區(qū)域Ai和Aj上隨機選擇兩地點的協(xié)方差的均值，并假定對數(shù)相對風(fēng)險曲面的區(qū)域均值Si的分布近似于

的分布，以上的近似處理更方便計算。

S(x)的相關(guān)結(jié)構(gòu)決定了Ri的相關(guān)構(gòu)成。我們假設(shè)GRF是均勻各向同性的，那么 γ(x，y)=ρ(d)，其中 d=||x-y||，Wackernagel在文獻[2]中定義該函數(shù)為：

通常要求假設(shè)的風(fēng)險曲面既要和子區(qū)域的空間相關(guān)性相符合，又要和總體區(qū)域的風(fēng)險相符合。因此上面的高次函數(shù)只依賴一個參數(shù)，并且我們只考慮了兩區(qū)域的距離，而沒有考慮區(qū)域的寬度廣度。小區(qū)域的范圍也限制了我們分析相互關(guān)系的信息，因此在分析光滑地圖時，小范圍的特點不會被夸大描述。

2.2 算法

我們首先討論之前的論述。為提高參數(shù)的可解釋性和計算效果，我們用來描述 σ2，其中=median{Var[Si]，i=1，…，n}，則其中 c=median{γ(Ai，Aj)，i=1，…，n}。因為的后驗分布比σ2對參數(shù)D的依賴更小，因此增強了Markov鏈的收斂性。 定義H（D）為以 γ(Ai，Aj)/c為元的 n×n階矩陣，因此cov(Si，Sj)=H(D)ij。假設(shè) α 的先驗為正態(tài)分布，τ的先驗為gamma分布，后驗密度為：

但它不能用分析方法化簡，因此我們轉(zhuǎn)而使用Markov chain Monte Carlo(MCMC)方法。范圍參數(shù)D離散的先驗會彈性丟失一些，但允許預(yù)先對其進行大量的計算，包括矩陣的建立和矩陣的逆。式(7),每帶入一次，需計算一次值D，用到了數(shù)值積分的計算。

MCMC所需要的滿條件分布在下面給出。

應(yīng)用多元正態(tài)分布的性質(zhì)和H(D)對D依賴性的減弱，Si|S-i的分布是正態(tài)的且：

其中 Hi是矩陣 H 去掉了第 i行 i列后的（n-1）×（n-1）階矩陣，hi是矩陣H第i列同時又去掉了第i個元生成的，由這些先驗條件分布可得每一個Si的滿條件分布。

假設(shè)均值水平參數(shù)α的先驗為正態(tài)分布N(mα，vα)，則α的條件分布 p(α|S)∝：

對于二次型有配方公式：

其中：

假設(shè) τ～Ga(a，b)，則 τ的條件分布：

從而得出條件分布為Ga(a'，b')，其中：

令πj=p(D=Dj)，j=1，…，k表示D的k種可能取值下的先驗概率。那么D對S、α和τ的條件分布是

區(qū)域?qū)?shù)相對風(fēng)險Si服從均值向量為α1n協(xié)方差矩陣為τ-1H(D)的多元正態(tài)分布，則以α、τ和D為條件，Si的滿條件分布：

經(jīng)計算得到：

其中：

雖然從Si值的角度來說，我們已計算出估計值，但我們需要進一步分析對數(shù)相對風(fēng)險曲面S(x)。

利用來自后驗樣本的每一組參數(shù)集和S(x)的條件分布，就可以生成一個m維向量S(x)在x處的后驗樣本，因此我們只需得到S(x)的條件分布即可。S(x)的條件分布是多元正態(tài)分布，為得到條件分布的期望向量和協(xié)方差陣，我們需要計算S(x)和Sj的協(xié)方差。該值可利用2.1中的近似計算得到：

其中γ(x，Aj)表示x處的S(x)與區(qū)域 Aj上的隨機點間的自相關(guān)均值。于是：

其中 K是一個由 Kij=γ(Ai，xj)/c構(gòu)成的n×m階矩陣，G是由Gij=γ(xi，xj)/c構(gòu)成的m×m階矩陣。由此，來自S(x)后驗分布的樣本以及對應(yīng)的R(x)即可得。

3 小結(jié)

本文中，我們提出了一種用于空間疾病地圖中模擬空間差異的方法，將空間疾病風(fēng)險模型R(x)建立在高斯隨機場下，并得出了空間疾病風(fēng)險R(x)的計算方法，對建立我國疾病風(fēng)險具有一定的實踐意義。

[1]J.P.Chiles,P.Delfiner.Geostatistics:Modeling Spatial Uncertainty[M].New York:Wiley,2003.

[2]H.Wackrnagel. Multivariate Geostatistics:an Introduction with Applications[M].New York:Springer,1995.