人群归因分数
人群归因分数（Population Attributable Fraction, PAF）是指总人群发病率中归因于暴露的部分。其计算过程是通过对实际人群中危险因素暴露的分布与理论最小分布比较，若人群中危险因素暴露降低到理论最小分布，估计疾病或死亡降低的比例。当危险因素是连续型变量时，可以通过以下公式计算人群归因分数：
其中，RR(x)是危险因素暴露水平为x时的相对危险函数；P1(x)是危险因素当前暴露密度函数；P2(x)是危险因素暴露水平为理论最小值时的密度函数；m为危险因素的最高暴露值。
人群归因分数的计算
示例数据使用的数据集如下，涉及的变量有：
value：PM2.5的年平均暴露值
status：是否患病，0代表否、1代表是
time：观察到结局发生的时间（年），或者随访结束的时间
2.1数据导入
pacman::p_load(survival) #安装并加载包setwd("~") #设置工作路径PAFdata <- read.csv("PAFdata.csv") #导入数据
2.2构建模型，提取相关值
res.cox <- coxph(Surv(time, status) ~ value, data = PAFdata) #构建模型HR1 <- summary(res.cox)$conf.int[1,1] #提取模型的HR值，以HR值代表RR值进行PAF计算HRlower1 <- summary(res.cox)$conf.int[1,3] #提取模型的HR值置信区间下限值HRupper1 <- summary(res.cox)$conf.int[1,4] #提取模型的HR值置信区间上限值miu1_1 <- mean(PAFdata$value) #计算value的均值sigma1_1 <- sd(PAFdata$value) #计算value的标准差lowerlimit1 <- min(PAFdata$value) #计算value的最小值higherlimit1 <- max(PAFdata$value) #计算value的最大值
2.3拟合危险因素平均暴露水平为理论最小值的数据集
set.seed(1234) #设置随机种子，确保结果的可重复性data2$expouse.F = rnorm(n = 22298, mean = 15, sd = 1) #构建的数据集与真实数据集的样本量均为22298，均值15微克/立方米为中国PM2.5年平均浓度限值一级标准，标准差为1miu2_1 <- mean(data2$expouse.F) #提取expouse.F的均值sigma2_1 <- sd(data2$expouse.F) #提取expouse.F的标准差
2.4计算真实数据的暴露分布
F1 <- function(HR1, miu1_1, sigma1_1, miu2_1, lowerlimit1,higherlimit1){return(pnorm ((miu2_1 - miu1_1)/sigma1_1) - pnorm((lowerlimit1 - miu1_1)/sigma1_1) + HR1^(miu1_1 - miu2_1)*exp((sigma1_1 * log (HR1))^ 2 *0.5)*(pnorm((higherlimit1 - miu1_1)/sigma1_1 - sigma1_1 * log (HR1))-pnorm((miu2_1 - miu1_1)/sigma1_1 - sigma1_1* log(HR1))))}
2.5计算拟合数据的暴露分布
F2 <- function (HR1, miu2_1, sigma2_1, lowerlimit1,higherlimit1){return(pnorm (0)- pnorm ((lowerlimit1 - miu2_1)/sigma2_1)+ exp((sigma2_1* log(HR1))^2 * 0.5)* (pnorm((higherlimit1 - miu2_1)/sigma2_1 - sigma2_1 * log (HR1))-pnorm(- sigma2_1* log(HR1))))}
2.6计算PAF值
PAF1 <- round ((F1 (HR1,miu1_1,sigma1_1,miu2_1,lowerlimit1,higherlimit1 ) - F2 (HR1,miu2_1,sigma2_1,lowerlimit1,higherlimit1))/F1 (HR1,miu1_1,sigma1_1,miu2_1,lowerlimit1,higherlimit1),4) #4表示保留四位小数
2.7计算PAF的上下限区间
#计算PAF的置信区间上限值F1 <- function(HRlower1, miu1_1, sigma1_1, miu2_1, lowerlimit1,higherlimit1){return(pnorm ((miu2_1 - miu1_1)/sigma1_1) - pnorm((lowerlimit1 - miu1_1)/sigma1_1) + HRlower1^(miu1_1 - miu2_1)*exp((sigma1_1 * log (HRlower1))^ 2 *0.5)*(pnorm((higherlimit1 - miu1_1)/sigma1_1 - sigma1_1 * log (HRlower1))-pnorm((miu2_1 - miu1_1)/sigma1_1 - sigma1_1* log(HRlower1))))}F2 <- function (HRlower1, miu2_1, sigma2_1, lowerlimit1,higherlimit1){return(pnorm (0)- pnorm ((lowerlimit1 - miu2_1)/sigma2_1)+ exp((sigma2_1* log(HRlower1))^2 * 0.5)* (pnorm((higherlimit1 - miu2_1)/sigma2_1 - sigma2_1 * log (HRlower1))-pnorm(- sigma2_1* log(HRlower1))))}PAFlower1 <- round ((F1 (HRlower1,miu1_1,sigma1_1,miu2_1,lowerlimit1,higherlimit1 ) - F2 (HRlower1,miu2_1,sigma2_1,lowerlimit1,higherlimit1))/F1 (HRlower1,miu1_1,sigma1_1,miu2_1,lowerlimit1,higherlimit1),4)#计算PAF的置信区间上限值F1 <- function(HRupper1, miu1_1, sigma1_1, miu2_1, lowerlimit1,higherlimit1){return(pnorm ((miu2_1 - miu1_1)/sigma1_1) - pnorm((lowerlimit1 - miu1_1)/sigma1_1) + HRupper1^(miu1_1 - miu2_1)*exp((sigma1_1 * log (HRupper1))^ 2 *0.5)*(pnorm((higherlimit1 - miu1_1)/sigma1_1 - sigma1_1 * log (HRupper1))-pnorm((miu2_1 - miu1_1)/sigma1_1 - sigma1_1* log(HRupper1))))}F2 <- function (HRupper1, miu2_1, sigma2_1, lowerlimit1,higherlimit1){return(pnorm (0)- pnorm ((lowerlimit1 - miu2_1)/sigma2_1)+ exp((sigma2_1* log(HRupper1))^2 * 0.5)* (pnorm((higherlimit1 - miu2_1)/sigma2_1 - sigma2_1 * log (HRupper1))-pnorm(- sigma2_1* log(HRupper1))))}PAFupper1 <- round ((F1 (HRupper1,miu1_1,sigma1_1,miu2_1,lowerlimit1,higherlimit1 ) - F2 (HRupper1,miu2_1,sigma2_1,lowerlimit1,higherlimit1))/F1 (HRupper1,miu1_1,sigma1_1,miu2_1,lowerlimit1,higherlimit1),4)
PAF置信区间的上下限计算只是把公式中的HR值替换为HRlower1和HRupper1
2.8查看PAF值及其置信区间
PAF1PAFlower1PAFupper1
计算结果表明，38.7%（95%CI:32.03%，44.83%）的病例发生可归因于PM2.5，将PM2.5控制到年平均浓度一级标准，可减少38.7%的病例发生。
参考文献：
