連續性危險因素暴露的人群歸因危險度精確計算及統計軟件實現*

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

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

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

人群歸因危險度原理

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

(1)

RR(x)：危險因素暴露水平為x時的相對危險函數；

p1(x)：危險因素當前暴露密度函數；

p2(x)：危險因素反事實暴露密度(如：理論最小風險暴露分布)函數；

l：危險因素在人群中可能的最低暴露值；

h：危險因素在人群中可能的最高暴露值。

人群歸因危險度計算

(2)

(3)

其中Φ(*)為標準正態分布函數。

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

因此，

(4)

同理，

(5)

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

(6)

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

實例應用

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

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

1.Excel中操作步驟

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

圖1 Excel中數據集及PAF計算結果

Excel中在L：M列對公式(4)和(5)的計算過程需調用函數NORMDIST(x，mean，standard_dev，cumulative)，EXP(number)和LN(number)。以Item 1為例，在單元格L2中運行命令如下：=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))，輸出結果為2.1266。同理，單元格M2中輸出結果為1.1702。由此，單元格N2中公式(6)的輸出結果為0.4497。

2.RStudio中操作步驟

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

attach(data)

a<-RRper10mmhg^0.1#向量，人群不同性別年齡組的危險因素單位水平的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)，輸出結果四舍五入保留4位小數

dataPAF<-cbind(data，PAF)

View(dataPAF)

輸出結果展示見圖2。

圖2 RStudio中PAF計算結果展示

小 結

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

疾病負擔研究越來越受到研究者的重視，人群歸因危險度是疾病負擔計算比較難的部分，在Excel/RStudio等常用統計軟件中簡單、靈活地快速實現其計算非常重要，本文的研究結果為相關計算提供了一個便捷、普適的計算工具。