連續(xù)性危險(xiǎn)因素暴露的人群歸因危險(xiǎn)度精確計(jì)算及統(tǒng)計(jì)軟件實(shí)現(xiàn)*

上海市疾病預(yù)防控制中心(200336)

劉曉俠# 王春芳# 汪 晶 虞慧婷 鄭 楊 劉丹妮 楊群娣 施 燕△

人群歸因危險(xiǎn)度(population attributable fraction，PAF)是一種廣泛用于評估人群暴露對健康影響的流行病學(xué)指標(biāo)[1]。其計(jì)算過程是通過對實(shí)際人群中危險(xiǎn)因素暴露的分布與理論最小分布比較，若人群中危險(xiǎn)因素暴露降低到理論最小分布，估計(jì)疾病或死亡降低的比例[2]。從數(shù)據(jù)類型看，人群危險(xiǎn)因素暴露主要分為兩大類：分類變量和連續(xù)性變量。當(dāng)人群危險(xiǎn)因素暴露水平為分類變量時(shí)，人群歸因危險(xiǎn)度計(jì)算過程比較簡單[3]；但人群中危險(xiǎn)因素暴露呈連續(xù)性分布的情況也廣泛存在，如人群中呈正態(tài)(近似正態(tài))分布的收縮壓，目前的研究通常采用離散求和的方式來近似估計(jì)其對疾病的人群歸因危險(xiǎn)度，估算過程較繁瑣復(fù)雜，尚無可供直接實(shí)現(xiàn)此類計(jì)算的統(tǒng)計(jì)軟件或程序模塊，其應(yīng)用的普及受到一定影響。本文將通過理論推導(dǎo)，實(shí)現(xiàn)危險(xiǎn)因素暴露呈正態(tài)分布時(shí)人群歸因危險(xiǎn)度的快速精確計(jì)算，并在Excel和RStudio軟件中以收縮壓(SBP)導(dǎo)致缺血性心臟病(IHD)為例，編制程序快速實(shí)現(xiàn)其人群歸因危險(xiǎn)度計(jì)算，旨在為疾病負(fù)擔(dān)研究中提高相關(guān)指標(biāo)的運(yùn)算效率和通過Excel/RStudio等統(tǒng)計(jì)軟件實(shí)現(xiàn)這一過程操作提供幫助。

人群歸因危險(xiǎn)度原理

對于人群中危險(xiǎn)因素暴露呈連續(xù)性分布造成的疾病負(fù)擔(dān)，在不同的性別和年齡組，通過比較當(dāng)前的暴露分布和理論上的最小風(fēng)險(xiǎn)暴露分布，人群歸因危險(xiǎn)度表示為[4]：

(1)

RR(x)：危險(xiǎn)因素暴露水平為x時(shí)的相對危險(xiǎn)函數(shù)；

p1(x)：危險(xiǎn)因素當(dāng)前暴露密度函數(shù)；

p2(x)：危險(xiǎn)因素反事實(shí)暴露密度(如：理論最小風(fēng)險(xiǎn)暴露分布)函數(shù)；

l：危險(xiǎn)因素在人群中可能的最低暴露值；

h：危險(xiǎn)因素在人群中可能的最高暴露值。

人群歸因危險(xiǎn)度計(jì)算

(2)

(3)

其中Φ(*)為標(biāo)準(zhǔn)正態(tài)分布函數(shù)。

做變換x=σ1t+μ1，則：

因此，

(4)

同理，

(5)

將公式(4)和(5)代入公式(1)：

(6)

在得到公式(6)的結(jié)果后，人群歸因危險(xiǎn)度計(jì)算可以通過統(tǒng)計(jì)軟件快速實(shí)現(xiàn)，極大地方便了其他研究者對此類人群歸因危險(xiǎn)度計(jì)算的實(shí)現(xiàn)。本文結(jié)合Excel和RStudio軟件的操作和程序語句以及常見實(shí)例計(jì)算其人群歸因危險(xiǎn)度，為更好地通過統(tǒng)計(jì)軟件實(shí)現(xiàn)人群歸因危險(xiǎn)度計(jì)算提供幫助。

實(shí)例應(yīng)用

本文的實(shí)例數(shù)據(jù)來源于文獻(xiàn)[6]。研究的目的是結(jié)合文獻(xiàn)[6]中2013年北京市戶籍人口分年齡組、性別的收縮壓(SBP)值和不同年齡收縮壓導(dǎo)致缺血性心臟病(IHD)的RR值以及公式(4)、(5)、(6)，利用Excel和RStudio軟件計(jì)算2013 年北京市戶籍人口分年齡組、性別的收縮壓導(dǎo)致缺血性心臟病的人群歸因危險(xiǎn)度，并簡述其計(jì)算操作步驟。表1為變量重命名、描述及賦值。

表1 變量重命名、描述及賦值

1.Excel中操作步驟

根據(jù)表1，將不同的性別、年齡組相關(guān)參數(shù)存儲至Excel中，數(shù)據(jù)集格式見圖1中A：K列。

圖1 Excel中數(shù)據(jù)集及PAF計(jì)算結(jié)果

Excel中在L：M列對公式(4)和(5)的計(jì)算過程需調(diào)用函數(shù)NORMDIST(x，mean，standard_dev，cumulative)，EXP(number)和LN(number)。以Item 1為例，在單元格L2中運(yùn)行命令如下：=NORMDIST((G2-E2)/F2，0，1，TRUE)-NORMDIST((I2-E2)/F2，0，1，TRUE)+K2^(E2-G2)*EXP((F2*LN(K2))^2*0.5)*(NORMDIST((J2-E2)/F2-F2*LN(K2)，0，1，TRUE)-NORMDIST((G2-E2)/F2-F2*LN(K2)，0，1，TRUE))，輸出結(jié)果為2.1266。同理，單元格M2中輸出結(jié)果為1.1702。由此，單元格N2中公式(6)的輸出結(jié)果為0.4497。

2.RStudio中操作步驟

對數(shù)據(jù)集以圖1中 A：K列形式存儲的入Excel文件data.xlsx，RStudio可直接將其導(dǎo)入。點(diǎn)擊菜單File-Import Dataset-From Excel-Browse-選擇data.xlsx文件，點(diǎn)擊open-Import，導(dǎo)入原始數(shù)據(jù)集data.xlsx。導(dǎo)入成功后，運(yùn)行命令如下：

attach(data)

a<-RRper10mmhg^0.1#向量，人群不同性別年齡組的危險(xiǎn)因素單位水平的RR值

F1<-function(a，miu1，sigma1，miu2，lowerlimit，higherlimit){

return(pnorm((miu2-miu1)/sigma1)-pnorm((lowerlimit-miu1)/sigma1)+a^(miu1-miu2)*exp((sigma1*log(a))^2*0.5)*(pnorm((higherlimit-miu1)/sigma1-sigma1*log(a))-pnorm((miu2-miu1)/sigma1-sigma1*log(a))))

}#編輯公式(4)

F2<-function(a，miu2，sigma2，lowerlimit，higherlimit){

return(pnorm(0)-pnorm((lowerlimit-miu2)/sigma2)+exp((sigma2*log(a))^2*0.5)*(pnorm((higherlimit-miu2)/sigma2-sigma2*log(a))-pnorm(-sigma2*log(a))))

}#編輯公式(5)

PAF<-round((F1(a，miu1，sigma1，miu2，lowerlimit，higherlimit)-F2(a，miu2，sigma2，lowerlimit，higherlimit))/F1(a，miu1，sigma1，miu2，lowerlimit，higherlimit)，4)

}#編輯公式(6)，輸出結(jié)果四舍五入保留4位小數(shù)

dataPAF<-cbind(data，PAF)

View(dataPAF)

輸出結(jié)果展示見圖2。

圖2 RStudio中PAF計(jì)算結(jié)果展示

小 結(jié)

人群歸因危險(xiǎn)度是定量描述危險(xiǎn)因素對人群致病效應(yīng)大小的統(tǒng)計(jì)指標(biāo)，受到流行病學(xué)家和統(tǒng)計(jì)學(xué)家的高度重視[7]。對危險(xiǎn)因素連續(xù)性暴露引起的疾病負(fù)擔(dān)研究，如常見的危險(xiǎn)因素暴露呈正態(tài)(近似正態(tài))分布時(shí)，本文通過積分計(jì)算，實(shí)現(xiàn)了與理論最小風(fēng)險(xiǎn)暴露分布平均值相關(guān)的連續(xù)性危險(xiǎn)因素的人群歸因危險(xiǎn)度的精確計(jì)算，同時(shí)將其計(jì)算過程通過常見的Excel和RStudio統(tǒng)計(jì)軟件結(jié)合實(shí)例快速實(shí)現(xiàn)，以方便更多的研究人員使用。

疾病負(fù)擔(dān)研究越來越受到研究者的重視，人群歸因危險(xiǎn)度是疾病負(fù)擔(dān)計(jì)算比較難的部分，在Excel/RStudio等常用統(tǒng)計(jì)軟件中簡單、靈活地快速實(shí)現(xiàn)其計(jì)算非常重要，本文的研究結(jié)果為相關(guān)計(jì)算提供了一個(gè)便捷、普適的計(jì)算工具。