连续性危险因素暴露的人群归因危险度精确计算及统计软件实现*

上海市疾病预防控制中心(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等常用统计软件中简单、灵活地快速实现其计算非常重要，本文的研究结果为相关计算提供了一个便捷、普适的计算工具。