基于高斯隨機(jī)場(chǎng)下的疾病風(fēng)險(xiǎn)模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 高斯隨機(jī)場(chǎng)模型

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

由于：

我們只需計(jì)算E[RiRj]：

所以：

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

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

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

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

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

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

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

2.2 算法

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

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

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

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

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

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

對(duì)于二次型有配方公式：

其中：

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

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

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

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

經(jīng)計(jì)算得到：

其中：

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

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

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

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

3 小結(jié)

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

[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.