曾麗娟

(中山大學(xué)新華學(xué)院，廣東 廣州 510642)

大量的醫(yī)學(xué)研究[1]表明，大氣污染對(duì)人體呼吸系統(tǒng)會(huì)造成較大的損害，空氣污染物的增加在一定程度上使得呼吸系統(tǒng)疾病的患病率上升.因此，研究空氣污染與呼吸疾病患病的關(guān)系，具有重要的意義.

目前，不少的研究工作者對(duì)空氣污染物與呼吸疾病的關(guān)系進(jìn)行了分析[2-3].李寧等[3]采用時(shí)間序列方法建立呼吸系統(tǒng)疾病日門診量的自回歸模型，并分別將污染物SO2, NO2,PM10和PM2.5加入到模型中，得到對(duì)應(yīng)4種污染物的回歸方程.然而，這種方法側(cè)重于對(duì)病人數(shù)的單一建模，而且忽略了空氣污染物的滯后效應(yīng)，具有一定的局限性.相比之下，時(shí)間序列滯后回歸模型同時(shí)考慮了自變量在當(dāng)前時(shí)刻以及滯后時(shí)刻對(duì)因變量的影響，用該模型來分析空氣污染物與呼吸疾病患病人數(shù)的關(guān)系，更能體現(xiàn)空氣污染物對(duì)人體呼吸系統(tǒng)損害的滯后效應(yīng).

因此，本文將建立空氣污染物與患呼吸疾病人數(shù)的時(shí)間序列滯后回歸模型，并運(yùn)用隨機(jī)搜索變量選擇方法(Stochastic Search Variable Selection,SSVS)[4-5]進(jìn)行空氣污染變量及滯后變量的選擇.SSVS的基本思想是通過建立一個(gè)層次貝葉斯正態(tài)混合模型，利用Gibbs抽樣[6]從后驗(yàn)分布中抽樣，借助示性變量來判別變量是否被選入模型，最終得到后驗(yàn)概率最高的子模型.相比于傳統(tǒng)的變量選擇準(zhǔn)則，如AIC準(zhǔn)則，SSVS避免了計(jì)算所有可能子模型的模型信息，具有潛在的優(yōu)越性.

1 滯后回歸模型的統(tǒng)計(jì)結(jié)構(gòu)

(1)

2 滯后回歸模型的貝葉斯推斷

為了將SSVS運(yùn)用于滯后回歸模型的自變量選擇，引入示性變量γ=(γ1′,…,γp′)′,γi=(γi,0,γi,1,…,γi,d)′,i=1,…,p.示性變量γi,j的作用是：當(dāng)所選擇模型包含xi,j項(xiàng)時(shí)，則γi,j=1，反之，則γi,j=0.在下面的推斷中，根據(jù)自變量順序和滯后階數(shù)大小對(duì)示性變量進(jìn)行重新排序，γ=((γ1,γ2,…,γd+1),…,(γ(p-1)(d+1)+1,…,γp(d+1))),其中,γj,j=1,2,…,p(d+1)取值為0或1.下面為模型參數(shù)設(shè)定相應(yīng)的先驗(yàn)信息：

1)系數(shù)β的每個(gè)分量看成來自由不同方差的兩個(gè)正態(tài)分布組成的混合正態(tài)分布：

(2)

3)對(duì)于正態(tài)分布，當(dāng)均值已知方差未知時(shí)，σ的共軛先驗(yàn)分布為逆Gamma分布族，即σ2|γ～I(xiàn)G(vγ/2,vγλγ/2),為了簡(jiǎn)化計(jì)算，超參數(shù)vγ和λγ一般選取υγ=υ和λγ=λ.

根據(jù)模型的似然函數(shù)和參數(shù)的先驗(yàn)分布，可以得到參數(shù)β,σ2和γ的條件后驗(yàn)密度：

π(β|Y,σ2,γ)∝f(Y|β,σ2)π(β|γ)∝

π(σ2|Y,β,γ)∝f(Y|β,σ2)π(σ2|γ)∝

π(γj|β,σ2,γ(-j),Y)∝f(Y|β,σ2)π(β|γ(-j),γj)π(σ2|γ(-j),γj)π(γ(-j),γj)，令

aj=f(Y|β,σ2)π(β|γ(-j),γj=1)π(σ2|γ(-j),γj=1)π(γ(-j),γj=1)，

bj=f(Y|β,σ2)π(β|γ(-j),γj=0)π(σ2|γ(-j),γj=0)π(γ(-j),γj=0)，

3 滯后回歸模型的Gibbs抽樣

SSVS利用Gibbs抽樣從γ的后驗(yàn)分布中產(chǎn)生樣本{γl,l=1,2,…,L}，并篩選出MCMC估計(jì)中γj接近1的項(xiàng)，作為模型選擇的變量.具體抽樣步驟如下：

1)設(shè)定初始狀態(tài)，給出起始點(diǎn)φ(0)=(β(0),(σ2)(0),γ(0))；

上述過程完成了第一輪抽樣，第二輪以φ(1)=(β(1),(σ2)(1),γ(1))作為初始值，重復(fù)上述過程.在足夠次數(shù)的迭代后，樣本序列收斂到一個(gè)獨(dú)立于初始值的平穩(wěn)分布.對(duì)于得到的長(zhǎng)度為k的Gibbs序列ψ(k)=(φ(1),φ(2),…,φ(k)),舍棄前m次的預(yù)迭代，則馬爾科夫鏈的實(shí)現(xiàn)值為ψ(m+1,…,k)=(φ(m+1),φ(m+2),…,φ(k))，取其平均值作為參數(shù)的估計(jì)值.

4 香港空氣污染與呼吸疾病數(shù)據(jù)分析

4.1 數(shù)據(jù)的選取

4.2 貝葉斯計(jì)算

根據(jù)香港空氣污染與呼吸道疾病相關(guān)文獻(xiàn)[7],SO2,TSP,NO,NO2是造成患呼吸疾病的主要因素，故xi,t-j,i=1,…,4,j=0,1,2，對(duì)應(yīng)的γ先驗(yàn)選取為p(γ=1)=U(0.75,0.95),其余的γ先驗(yàn)選取為p(γ=1)=0.5.τj取模型系數(shù)β最小二乘估計(jì)的標(biāo)準(zhǔn)差的均值，(σ2)0取方差的最小二乘估計(jì)值，超參數(shù)cj=10,j=1,…,21.R取單位陣I.SSVS選擇模型變量的Gibbs抽樣過程分成兩個(gè)步驟，在每個(gè)步驟的抽樣中均造兩條馬爾可夫鏈，燃燒期為10 000次，抽樣總次數(shù)為20 000次.

2)在變量縮減模型下進(jìn)行隨機(jī)變量搜索，對(duì)1)中第一階段隨機(jī)搜索篩選后的自變量進(jìn)行重新排序：x1,t-2,x2,t-0,x2,t-1,x2,t-2,x3,t-0,x4,t-0,x4,t-2,x7,t-0,對(duì)應(yīng)的γ編號(hào)為γ1,…,γ8.超參數(shù)τj≈0.03,j=1,2,…,8,(σ2)0≈190,σ2|γ～I(xiàn)G(2,95).表1給出了幾種模型變量選擇結(jié)果,可見，具有最大后驗(yàn)概率的滯后回歸模型包含的自變量有x7,t-0,x4,t-2,x2,t-0,x4,t-0，即濕度、NO2及其二階滯后值、TSP對(duì)呼吸疾病影響顯著.

表1 幾種具有較大后驗(yàn)概率的模型變量選擇結(jié)果

此外，示性變量的抽樣結(jié)果中,γ8,γ7,γ2,γ6,γ3,γ4,γ1,γ5分別為1,0.885 2,0.691 6,0.615,0.502 7,0.482,0.405 4,0.389 1，同樣說明濕度、NO2及其二階滯后值、TSP對(duì)呼吸疾病的影響較為顯著.

可見，SSVS根據(jù)γ值可以一次性準(zhǔn)確地找出模型，表2給出了最優(yōu)模型的貝葉斯估計(jì)結(jié)果.

表2 最優(yōu)滯后回歸模型參數(shù)的MCMC估計(jì)結(jié)果

圖1和圖2為對(duì)呼吸疾病人數(shù)影響最顯著的污染物NO2和TSP的后驗(yàn)密度估計(jì).

圖1 二氧化氮后驗(yàn)密度

圖2 懸浮粒子后驗(yàn)密度

因此，相關(guān)空氣污染與患呼吸疾病人數(shù)的滯后回歸模型表示為：

yt=0.072 74x2,t-0+0.063 55x4,t-0+0.116 9x4,t-2+0.508 6x7,t-0.

(3)

模型表明，濕度是影響患呼吸疾病人數(shù)最顯著的因素，NO2的二階滯后值對(duì)呼吸疾病的影響比當(dāng)前值顯著，TSP也會(huì)影響呼吸疾病患病人數(shù).各種污染物對(duì)呼吸疾病的影響呈現(xiàn)正相關(guān)性，即污染物的增加均會(huì)在一定程度上使得患呼吸疾病人數(shù)上升.模型的變量選擇結(jié)果中，SO2并沒被篩選到模型中，這是因?yàn)樵谙愀劭諝馕廴疚镏?，SO2的濃度在污染水平下[7]，模型體現(xiàn)了香港空氣污染的真實(shí)情況.

表3 兩個(gè)模型的估計(jì)精度對(duì)比

5 結(jié) 論

本文將隨機(jī)搜索變量選擇方法運(yùn)用于空氣污染與呼吸疾病的時(shí)間序列滯后回歸模型的變量選擇中,借助WinBUGS軟件，構(gòu)造基于Gibbs抽樣的數(shù)值計(jì)算,對(duì)模型參數(shù)進(jìn)行估計(jì).SSVS較好地結(jié)合各種空氣污染物以及呼吸疾病人數(shù)的數(shù)據(jù)特征進(jìn)行變量選擇，得出濕度對(duì)香港患呼吸疾病人數(shù)影響最顯著，NO2和TSP也是兩種影響比較顯著的空氣污染影響因素，此外，NO2的二階滯后值也對(duì)呼吸疾病人數(shù)有一定影響.

[1] Zanobetti A, Schwartz J. The effect of fine and coarse particulate air pollution on mortality: a national analysis[J]. Environmental Health Perspectives,2009,117(6):898-903.

[2] 劉迎春，龔潔，楊念念，等．武漢市大氣污染與居民呼吸系統(tǒng)疾病死亡關(guān)系的病例交叉研究[J]．環(huán)境與健康雜志，2012，29(3)：241-244．

[3] 李寧，彭曉武，張本延，等．大氣污染與呼吸系統(tǒng)疾病日門診量的時(shí)間序列分析[J]．環(huán)境與健康雜志，2009，26(12)：1077-1080．

[4] Chen M H,Shao Q M,Ibrahim J G.Monte carlo methods in bayesian computation[M].New York:Springer,2000:267-304

[5] Kwon D, Landi M T, Vannucci M,etal. An efficient stochastic search for bayesian variable selection with High-Dimensional correlated predictors[J]. Computational Statistics & Data Analysis,2011,55(10):2807-2818.

[6] George E I,McCulloch R E. Variable selection via gibbs sampling[J]. Journal of American Statistical Association,1993,88(423):881-889.

[7] Brajer V, Mead R W, Xiao F.Valuing the health impacts of airpollution in Hong Kong[J]. Journal of Asian Economics,2006,17(1):85-102.