Stata：高维倾向得分法

👇 连享会 · 推文导航 | www.lianxh.cn

• 🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令
• 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享
• 💹 计量：回归分析 | 交乘项-调节  | IV-GMM | 时间序列 | 面板数据  | 空间计量 | Probit-Logit | 分位数回归
• ⛳ 专题：SFA-DEA | 生存分析  | 爬虫 | 机器学习 | 文本分析
• 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching
• 🔨 工具：工具软件 | Markdown | Python-R-Stata
• 🎧 课程：公开课-直播 |  计量专题 | 关于连享会

作者：陈佳慧 (浙江大学)
邮箱：[e-mail]

编者按：本文主要摘译自以下文章，特此致谢！

Ishimaru M. Introduction to High-dimensional Propensity Score Analysis[J]. Annals of Clinical Epidemiology, 2020, 2(4): 85-94. -PDF-Tazare J, Smeeth L, Evans S J W, et al. Implementing high‐dimensional propensity score principles to improve confounder adjustment in UK electronic health records[J]. Pharmacoepidemiology and Drug Safety, 2020, 29(11): 1373-1381. -PDF-

---

目录

• 1. 引言

• 2. 理论背景

• 3. 命令介绍

• 4. Stata 实操

• 5. 相关推文

---

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：

1. 引言

在因果推断中，倾向得分匹配法 (propensity score matching，PSM) 经常被用来解决样本自选择等问题，但在观测样本个体特征参数较多时，从对照组中找到与实验组相似的个体便十分困难。高维倾向得分法 (high-dimensional propensity score approach，hdps) 克服了上述局限性，可以很容易地使用大量数据信息来减少混杂因素(confounders) 引起的偏差。

目前，高维倾向得分法 (hdps) 被广泛应用于流行病学研究中。其相关命令有 SAS、R、以及 Stata 等版本。在这里，我们主要关注高维倾向得分法 (hdps) 的理论内容和 Stata 操作。

2. 理论背景

高维倾向得分法 (hdps) 的关键假设是，可以从大量样本数据中获得不可测度的混杂变量代理信息 (proxy information for important unmeasured confounders)。不同于传统倾向得分法，高维倾向得分的自变量是自动地 (而非人为地) 从样本数据中筛选出来，以最小化混杂偏差。

例如，在评估治疗方式 A (treatment A) 与死亡率的因果关系中，对于患者的日常活动 (activities of daily living，ADL) 这一不可测度的混杂因素，hdps 通过控制与 ADL 相关的其他变量，例如是否中风 (stroke)、是否有康复记录 (rehabilitation)、是否有抗痴呆药物处方 (anti-dementia)  三个可测变量，间接控制了 ADL 带来的干扰偏差。

hdps 的计算，主要分为以下几步：

1. 确定数据维度：将样本数据划分为不同的维度，确定维度数量 (通常设定为 200)，每个维度下的子信息集都有一个特定的编码系统 (code system)；
2. 识别备选协变量：列出每个维度下的变量编码，并按照症状发生率 (prevalence) 对其排序。通常选取发病率最靠前的 100-200 个变量编码被识别为备选协变量；
3. 计算编码出现的次数：对于每个观测样本 (例如一位患者的诊断记录)，计算这个样本下出现的所有备选协变量代码的次数；
4. 估计备选协变量的偏差：

• 首先，根据公式 计算乘性偏差项。其中 和 分别表示变量编码在处理组和对照组出现的频率， 表示未调整的相对发病率；
• 然后，根据 的对数绝对值来比较备选协变量。

1. 选择最终的协变量：将 的对数绝对值降序排列，通常选出前 500 个备选协变量进入到倾向得分估计模型中；
2. 估计倾向得分：接下来，与传统的 PS 方法相同，使用多变量逻辑回归 (multivariate logistic regression) 估计高维倾向得分。逻辑回归中包含由步骤 5 选出的变量、以及预先设定的的变量 (hdps and predefined covariates)。回归得到的 hdps 可以用来估计治疗效应 (treatment effect)。

hdps 方法的局限性：

1. 并非所有不可测度的混杂因素都可以使用 hdps 来调整，只有与数据样本中变量相关的不可测度混杂因素才可以使用 hdps 方法来控制；
2. hdps 方法可能出现过度调整 (over-adjustment) 的情况，从而造成共线性和估计中的统计无效性。但是有研究证明这种情况出现的概率非常小。

3. 命令介绍

hdps 命令主要用于：

1. 执行 hdps 算法的数据处理和变量选取步骤；
2. 根据图形评估所选协变量的性质。

hdps 命令需要通过 github 外部途径进行安装。因此在安装 hdps 命令之前，要预先安装 github 命令。具体如下：

net install github, from("https://haghish.github.io/github/")
github install johntaz/hdps

// 或者通过 lxhget 命令安装
lxhget hdps.pkg, install replace

hdps 算法包含多个步骤，其命令语法结构如下：

hdps setup(string, varname), [save(string) study(string) patid(varname) EXPosure(varname) OUTcome(varname)]
hdps prevalence, [top(numlist max=1 integer) nofilter] 
hdps recurrence
hdps prioritize, method(string) top(numlist >0 integer) [zerocell] 
hdps graphics varlist(min = 2) [if], type(string) [DIMension(varname) pr(numlist max=1) twoway_options]

• hdps setup()：数据维度，多个数据维度使用 () 隔开，每个 (string, varname) 中包含数据文件名称和数据维度变量名称；[save() study() patid() EXPosure() OUTcome()] 为可选项，用于声明关键协变量，并将初步处理结果存储到指定路径。其中，
• save()：在指定路径下存储处理结果文件；
• study：声明关键解释变量文件名称前缀；
• patid()：声明患者 id 变量；
• exp()：声明患者处理类型变量 (即处理组或控制组)；
• out()：声明患者的治疗结果变量。
• hdps prevalence：计算每个数据维度下的症状发病率 (feature prevalence)。其中，
• top() 指定发病率排名前列的症状个数；
• nofilter 表示不指定具体的症状个数。
• hdps recurrence：根据症状出现与否生成二值的 hdps 协变量 (binary covariates)，并计算每个编码出现的频次；
• hdps prioritize：根据协变量的排序，从备选协变量中选择出最终使用的协变量。其中，
• method() 指定排序方法，bross 表示的 Bross formula 方法，exposure 表示 exposure-based approach；
• top() 指定按照选定方法进行排序的协变量个数；
• zerocell 对 Bross formula 方法的单元格进行 0.1 修正，这种做法通常在较少的结果变量情形中使用。
• hdps graphics：一个独立的命令，根据图形结果评估所选协变量的性质。其中，
• varlist() 指定需要作图的变量，最小变量个数为 2；
• if 为可选的条件语言；
• type() 指定图形评估方法，共有 3 种方法可供选择，boss 选项可以检测 Boss values 的分布，prevalence 选项用来比较不同组别中编码出现的频次，strength 选项用来比较 hdps 协变量与治疗变量和 hdps 协变量与结果变量之间的关联强度；
• DIMension() 指定数据维度变量，该选项仅在 type() 项为 prevalence 或 strength 时才可选；
• 可选项 pr() 指定在图形中使用虚线作出频次比率及其倒数，默认值为 2，该选项仅在 type() 项为 prevalence 时才可选；
• 可选项 twoway_options 表示任何 twoway() 函数的可选项都是允许的。

4. Stata 实操

我们使用 Tazare 等 (2020) 的文章中英国电子健康数据集 (Electronic Health Records，EHRs) 来演示 hdps 算法。代码和数据可通过以下命令获取：

. lxhget hdps.pkg, des

-------------------------------------------------------------------------
package hdps from https://file.lianxh.cn/data/h/hdps
-------------------------------------------------------------------------
INSTALLATION FILES                                 (type net install hdps)
      hdps.ado
      hdps_graphics.ado
      hdps_prevalence.ado
      hdps_prioritize.ado
      hdps_recurrence.ado
      hdps_setup.ado

ANCILLARY FILES                                    (type net get hdps)
      100_hdps_procedure.do
      101_propensity_score_analysis.do
      clinical_dim.dta
      clinical_dim_ever.dta
      cohort.dta
      therapy_dim.dta
      therapy_dim_ever.dta
-------------------------------------------------------------------------

. lxhget hdps.pkg, replace

Step 1：根据患者接触初级医疗服务的程度，在 EHRs 中确定三个数据维度，分别为：clinical，referral 和 prescription。其中，clinical 和 referral 维度的数据使用 ICD-10 编码系统，而 prescription 维度的数据使用 BNF 编码系统。

step 2：运行 hdps setup 命令声明数据维度和关键变量。

. use cohort.dta, replace // 打开了 EHRs 的患者信息数据集
. capt mkdir output       // 创建 output 文件夹

. hdps setup (clinical_dim, icd10 ever) (therapy_dim, bnf), study(example) ///
>   save(./output) patid(patid) exp(trt) out(outcome)

Data dimensions identified (code variable):
   Dimension 1:     clinical_dim (icd10)
   Dimension 2:     therapy_dim (bnf)

Note: 'ever' option specified at least once
Ever dimensions:
   Dimension 1:     clinical_dim_ever (icd10)
   
Output folder:
./output

其中，hdps setup 命令中 (clinical_dim, icd10 ever) (therapy_dim, bnf) 声明两种数据编码系统，前者为 clinical 和 referral 维度的 ICD-10 编码系统，后者为 prescription 维度的 BNF 编码系统。同时将数据结果储存在 output 文件夹下的 example_cohort_info.dta 文件中。该文件包含了已声明的患者 id、患者处理类型和患者的治疗结果变量 (即 patid，trt 和 outcome 变量)。

step 3：运行 hdps prevalence 命令，对于每个数据维度，选择症状出现频次排名前 100 的编码，并将选取结果以及与患者 id 匹配过的编码信息保存到 dta 文件中。

.   hdps prevalence, top(100)

Identifying most prevalent features:
Selecting top 100 from each dimension
   Dimension 1:     Completed: selected 100 features
   Dimension 2:     Completed: selected 100 features

Incorporating 'ever' information:
   Dimension 1:     Completed

Output files:
(1) example_feature_prevalence.dta
(2) example_patient_totals.dta

step 4：运行 hdps recurrence 命令建立备选的 hdps 协变量池，计算每个编码在选定的窗口期中出现的频次，并将匹配到患者 id 的协变量信息保存到 dta 文件中。

. hdps recurrence

Loading data:
Completed
Generating HDPS covariates and assessing feature recurrence:
Progress: 0%...20%...40%...60%...80%...Completed
Number of binary HDPS covariates created:
600

Output file:
(1) example_hdps_covariates.dta

2 个编码系统各包含排名前 100 的编码，每个编码被指定了 3 个指示性的二值变量 (dimension_code_ever 表示症状编码是否出现过，dimension_code_spor 表示症状编码出现次数是否大于中位数，dimension_code_freq 表示症状编码出现次数是否大于 75% 分位数)。因此 hdps 协变量池中的变量个数为 600。

step 5：运行 hdps prioritize 命令使用 Bross formula 方法对池中协变量进行排序，分别选取排名前 50 和前 100 的变量作为最终使用的 hdps 协变量。

. hdps prioritize, method(bross) top(50 100)

Ranking HDPS covariates:
Prioritizing using the Bross formula:
Progress: 0%...20%...40%...60%...80%...Completed

Forming hd-PS cohort(s) based on top ranked covariates:
Selecting: 50, and 100.

Output files:
(1) example_bias_info.dta
(2) example_hdps_covariates_top_50.dta
(3) example_hdps_covariates_top_100.dta

step 6：运行 hdps graphics 命令检测 step 5 确定的 hdps 协变量的性质。

首先，载入偏差信息数据集 example_bias_info.dta，生成数据维度指示变量 dim，并使用 hdps graphics 命令作出排名前 100 个 hdps 协变量的 Boss values 分布图。

. clear all
. use ./output/example_bias_info.dta, replace 
. gen dimension=substr(code,1,2)
. encode dimension, gen(dim)
. drop dimension 
. hdps graphics abs_log_bias rank if rank<=100, type(bross)
. graph export "./output/bross_formula_distribution.png", width(2000) replace

接着，我们演示 type() 项为 prevalence 时的不同处理组之间协变量出现的频次。

. hdps graphics ce_strength cd_strength if rank<=100, type(strength) dim(dim)     ///
>   legend(order(1 "Clinical" 2 "Therapy") title("Data Dimensions",size(small))   ///
>   cols(3) rows(1) symxsize(*0.4) size(small))                 ///
>   ylabel(0.2 0.5 1 2 4, labsize(medsmall) angle(horizontal))  ///
>   xlabel(0(0.2)1, labsize(medsmall) angle(horizontal))        ///
>   xscale(log) yscale(log) ytitle("Strength of covariate-treatment association") ///
>   xtitle("Strength of covariate-outcome association")
. graph export "./output/covariate_associations.png", width(2000) replace

最后，我们演示 type() 项为 strength 时不同处理组之间的关联强度。

. hdps graphics ce_strength cd_strength if rank<=100, type(strength) dim(dim)     ///
>     legend(order(1 "Clinical" 2 "Therapy") title("Data Dimensions",size(small)) ///
>     cols(3) rows(1) symxsize(*0.4) size(small))                 ///
>     ylabel(0.2 0.5 1 2 4, labsize(medsmall) angle(horizontal))  ///
>     xlabel(0(0.2)1, labsize(medsmall) angle(horizontal))        ///
>     xscale(log) yscale(log) ytitle("Strength of covariate-treatment association") ///
>     xtitle("Strength of covariate-outcome association")

. graph export "./output/covariate_associations.png", width(2000) replace

至此，我们演示了 hdps 算法包的所有命令，需要说明的是，step 1 至 step 5 的子命令是按步骤顺序进行的，不能单独从 step 2 或之后的步骤开始运行。step 6 虽然是一个独立的命令，但也需要 step 5 的运行结果才可以继续。

在 step 6 的图形结果中，我们可以得出经过 hdps 算法后的协变量在不同组患者中是随机分布的，不再受到混杂因素的影响。因此，现在可以进行最后一个阶段——传统的倾向得分分析。在最后的倾向得分分析中，我们分别进行了使用和不使用 hdps 协变量的结果，以形成直观的对比。

其中，不使用 hdps 协变量的倾向得分结果如下：

. use cohort.dta, replace
. mkspline age_c = age, cubic nknots(4) // 对 age 变量进行限制三次样条，节点数量设置为4
. logit trt age_c1 age_c2 age_c3 female ses smoke alc bmicat nsaid_rx cancer hyper 

Logistic regression                                     Number of obs = 10,000
                                                        LR chi2(11)   =   8.04
                                                        Prob > chi2   = 0.7099
Log likelihood = -6750.4373                             Pseudo R2     = 0.0006
------------------------------------------------------------------------------
         trt | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      age_c1 |      0.003      0.009     0.30   0.763       -0.015       0.021
      age_c2 |     -0.018      0.027    -0.66   0.511       -0.070       0.035
      age_c3 |      0.084      0.108     0.77   0.441       -0.129       0.296
      female |     -0.050      0.042    -1.19   0.234       -0.132       0.032
         ses |      0.017      0.038     0.45   0.655       -0.058       0.093
       smoke |      0.018      0.068     0.26   0.798       -0.116       0.151
         alc |      0.017      0.065     0.26   0.798       -0.111       0.144
      bmicat |     -0.031      0.046    -0.67   0.503       -0.121       0.060
    nsaid_rx |      0.073      0.045     1.62   0.104       -0.015       0.160
      cancer |      0.070      0.051     1.38   0.167       -0.029       0.170
       hyper |      0.019      0.042     0.45   0.654       -0.063       0.100
       _cons |      0.271      0.384     0.71   0.480       -0.481       1.023
------------------------------------------------------------------------------

. // 对不同组别估计倾向得分
. predict pscore // 保存倾向得分
. hist pscore, by(trt) name(investigator, replace)
. graph export "./output/propensity_score_dist_investigator.png", width(2000) replace

. // 根据组别绘制倾向得分的概率密度图，并将其保存至 figs 文件夹中
. gen wts = 1/ps if trt == 1
. replace wts = 1/(1-ps) if trt == 0 // 根据组别倾向得分生成权重
. logistic outcome i.trt [pw=wts], robust // 根据处理结果估计倾向得分

Logistic regression                                     Number of obs = 10,000
                                                        Wald chi2(1)  =   0.12
                                                        Prob > chi2   = 0.7271
Log pseudolikelihood = -13468.506                       Pseudo R2     = 0.0000
------------------------------------------------------------------------------
             |               Robust
     outcome | Odds ratio   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         trt |
 Study Drug  |      1.015      0.042     0.35   0.727        0.935       1.101
       _cons |      1.483      0.047    12.30   0.000        1.393       1.579
------------------------------------------------------------------------------

从图形中可以看出，处理组 (Study Drug) 患者的倾向得分分布与控制组 (Comparator) 不是对称的，因而不具有可比性。

使用 hdps 协变量的倾向得分结果如下：

. tempname john      // 创建一个临时内存区域 "john" 来储存临时创建的文件
. postfile `john' str4(num_vars) float(es ll ul) using "./output/hdpsResults", replace 
. // 指定临时文件的存储位置和名称
. foreach v in 50 100 {
  2.     use "cohort.dta", replace
  3.     merge 1:1 patid using "./output/example_hdps_covariates_top_`v'.dta", assert(match) nogen
  4.     mkspline age_c = age , cubic nknots(4)
  5.     set matsize 200
  6.     logit trt age_c1 age_c2 age_c3 female ses smoke alc bmicat nsaid_rx cancer hyper d1* d2*
  7.     // 包含了 hdps 协变量、针对组别的倾向得分估计
.        predict pscore
  8.     hist pscore, by(trt) name(top_`v', replace)
  9.     graph export "./output//propensity_score_dist_top_`v'.png", width(2000) replace
 10.     gen wts = 1/ps if trt == 1
 11.     replace wts = 1/(1-ps) if trt == 0
 12.     logistic  outcome i.trt [pw=wts], robust
 13.     // 包含了 hdps 协变量的倾向、针对处理效应的倾向得分估计
.        mat b =r(table)
 14.     mat list b
 15.     local es = b[1,2]
 16.     local ll = b[5,2]
 17.     local ul = b[6,2]
 18.     post `john' ("`v'") (`es') (`ll') (`ul') 
 19. }   // 循环语句分别针对排名前 50 和前 100 的协变量估计两个组别的倾向得分
. postclose `john'   // 结束追加数据

可以看出，与不包括 hdps 协变量的倾向得分估计相比，该图形中控制组和对照组的倾向得分分布更对称，因此更具有可比性。

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
  lianxh psm, m
安装最新版 lianxh 命令：
  ssc install lianxh, replace

• 专题：倍分法DID
• 面板PSM+DID如何做匹配？
• 专题：断点回归RDD
• 当PSM遇上RDD：rddsga命令详解
• 专题：PSM-Matching
• Stata：psestimate-倾向得分匹配(PSM)中协变量的筛选
• 伍德里奇先生的问题：PSM-分析中的配对——小蝌蚪找妈妈
• Stata：psestimate-倾向得分匹配(PSM)中匹配变量的筛选
• Stata+PSM：倾向得分匹配分析简介
• Stata-从匹配到回归：精确匹配、模糊匹配和PSM
• Stata：PSM-倾向得分匹配分析的误区

New！ Stata 搜索神器：lianxh 和 songbl  GIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
  . ssc install lianxh
  . ssc install songbl
👉  使用：
  . lianxh DID 倍分法
  . songbl all

🍏 关于我们

• 连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。
• 直通车： 👉【百度一下：连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。