人群归因分数（2）：多分类变量、控制变量

---

专注系列化、高质量的R语言教程

推文索引 | 联系小编 | 付费合集

---

本文接上篇继续介绍人群归因分数（PAF）的计算，主要介绍风险因子是多分类变量、含控制变量两种情况。

本篇目录如下：

• 1 多分类变量

• 2 含控制变量

  • 2.1 错误计算方法

  • 2.2 正确计算方法1

  • 2.3 正确计算方法2（条件PAF）

  • 2.4 正确计算方法3

• 3 总结

1 多分类变量

上篇举的例子中，风险因子phys是二分变量。有读者问到风险因子是多分类变量时该如何计算PAF。

示例数据如下：

library(causalPAF)
library(tidyverse)
data <- strokedata %>%
  mutate(phys = as.numeric(phys) - 1,
         ahei3tert = as.numeric(ahei3tert) - 1)

ahei3tert变量表示样本饮食习惯的AHEI得分，原数据经过上述转换后，取值范围为0、1、2，数值越大表示饮食习惯越不健康。

样本分布情况：

table <- data %>%
  group_by(ahei3tert) %>%
  summarise(pop = length(case),
            pop_prop = pop/16623,
            case_num = sum(case),
            case_prop = mean(case))
table
##   ahei3tert   pop pop_prop case_num case_prop
##       <dbl> <int>    <dbl>    <dbl>     <dbl>
## 1         0  5114    0.308     2457     0.480
## 2         1  5824    0.350     2844     0.488
## 3         2  5685    0.342     2934     0.516

多分类情况下，PAF的计算公式仍然如下，只不过分母相加的项变多了：

## 计算相对风险度
(rr <- table$case_prop/table$case_prop[1])
## [1] 1.000000 1.016398 1.074200

## 人群暴露水平分布
(p <- table$pop_prop)
## [1] 0.3076460 0.3503579 0.3419960

## 计算PAF
1 - rr[1]/(sum(rr * p))
## [1] 0.03018201

2 含控制变量

控制变量一般刻画的是样本个体或家庭的特征，如年龄、性别、受教育程度、经济状况等，它们可能会间接影响样本的患病风险，通常不看作是风险因子，而是混杂因素（confounders）。

esex变量表示样本性别，1表示女性，2表示男性。在本例中，将其看作控制变量。

2.1 错误计算方法

以phys、esex为分组变量，展示样本的分布情况：

table <- data %>%
  group_by(esex, phys) %>%
  summarise(pop = length(case),
            pop_prop = pop/16623,
            case_num = sum(case),
            case_prop = mean(case))
table
##    esex  phys   pop pop_prop case_num case_prop
##   <int> <dbl> <int>    <dbl>    <dbl>     <dbl>
## 1     1     0   759   0.0457      344     0.453
## 2     1     1  5985   0.360      2973     0.497
## 3     2     0  1538   0.0925      699     0.454
## 4     2     1  8341   0.502      4219     0.506

将上表第一行的情况看作理想暴露水平（患病人口比例最小），然后计算PAF：

rr <- table$case_prop/table$case_prop[1]
p <- table$pop_prop

1 - rr[1]/sum(rr * p) 
## [1] 0.08512351

但这种计算方法是错误的。