Stata:高维倾向得分法-hdps

发布时间:2021-12-18 阅读 2926

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc, ihelp, rdbalance, gitee, installpkg

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

PDF下载 - 推文合集

作者:陈佳慧 (浙江大学)
邮箱sunny22961@163.com

编者按:本文主要摘译自以下文章,特此致谢!

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. 引言

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

目前,高维倾向得分法 (hdps) 被广泛应用于流行病学研究中。其相关命令有 SASR、以及 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. 估计备选协变量的偏差:
  • 首先,根据公式 BiasM=Pc1(RRcd1)+1Pc0(RRcd1)+1 计算乘性偏差项。其中 Pc1 和 Pc0 分别表示变量编码在处理组和对照组出现的频率,RRcd 表示未调整的相对发病率;
  • 然后,根据 BiasM 的对数绝对值来比较备选协变量。
  1. 选择最终的协变量:将 BiasM 的对数绝对值降序排列,通常选出前 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() 项为 prevalencestrength 时才可选;
    • 可选项 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

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh