基于高斯隨機場下的疾病風險模型

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

（1.燕山大學a.里仁學院基礎教學部；b.理學院統計學系，河北 秦皇島 066004；2.河北科技師范學院 數理系，河北 秦皇島 066004）

1 研究背景及統計框架

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

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

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

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

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

我們定義區域Ai上k階層人口的平均相對風險：

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

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

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

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

2 空間統計模型

2.1 高斯隨機場模型

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

由于：

我們只需計算E[RiRj]：

所以：

但R(x)的分布函數不能得到，空間統計的有關學者通常認為該分布近似對數正態分布，特別是當研究區域Ai相對較小時，近似情況也較好。在這種近似分析下：

服從多元正態分布，其均值、方差分別為：

(見文獻[1])，當區域范圍逐漸減小，我們取極限結果，進一步近似得到E[Si]≈α

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

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

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

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

2.2 算法

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

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

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

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

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

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

對于二次型有配方公式：

其中：

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

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

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

區域對數相對風險Si服從均值向量為α1n協方差矩陣為τ-1H(D)的多元正態分布，則以α、τ和D為條件，Si的滿條件分布：

經計算得到：

其中：

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

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

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

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

3 小結

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

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