周晓华：如何用数学语言来描述进行因果分析？

导语

我们的生活中充满因果推断问题，比如打疫苗的效果怎样，一种治疗方案结果如何。但仅仅依靠数据，往往只能观察到两个变量之间的相关性，要如何更进一步，确定它们的因果性呢？如何用数学语言将这些实际问题转化为因果推断问题，进行分析并指导实际决策？

为了更多地推动因果科学学科的发展，聚集国内外因果科学的一线科研工作者，共同讨论因果科学的最新进展，北京大学讲席教授、北京大学公共卫生学院生物统计系系主任、北京大学北京国际数学研究中心生物统计和信息研究事室主任周晓华等发起了泛太平洋因果推断大会。2021泛太平洋因果推断大会将于2021年9月11日-12日全程在线举办，详情见文末，欢迎感兴趣的朋友报名参与！

研究领域：因果科学，因果推断

周晓华、邓宇昊 | 作者

邓一雪 | 编辑

1. 因果推断的数学基础

一种新药在特定人群中是否有效？犯罪率改变是否由政策引起？从公立学校换到私立学校能否提高学习成绩？这些都是因果推断的问题。事实上，如果仅仅依靠数据，我们无法回答上述问题。因为很多时候，数据生成机制对观测者是未知的，我们只能观察到相关性。因果性在相关性的基础上更进一步：相关性可以由观察到的变量联合分布来描述，而对于因果性的分析需要知道数据的生成机制，或者条件改变时变量分布的变化。

如何用数学语言来描述进行因果分析呢？Rubin因果模型包含四个部分：(U,K,Y,Z)，其中，U是个体，K是一组可能的处理集合，Y是从K×U到实数的映射，Z是从U到K的映射，表示个体实际接受的处理。个体u在处理t下的潜在结果记作Y_t(u)。对于个体u来说，处理t相对于处理c的作用是Y_t(u)-Y_c(u)；对于整个人群来说，平均的因果作用就是E[Y_t(u)-Y_c(u)]。因果推断的基本问题是，这两种潜在结果Y_t(u)和Y_c(u)是无法同时观测的，因此我们必须依赖于额外的假设或条件才能识别出因果作用。

统计学家Fisher曾提出随机化的概念。在一个随机化试验中，每个个体都被随机地分配到两种处理方案中的一种，那么平均因果作用就可以识别了。在随机化试验中，潜在结果可以转化为观测结果，即E[Y_t(u)]=E[Y(u)|Z=t]，因此平均因果作用可以用观测数据来估计。

但实际上，随机化条件有可能被破坏。例如，

(1) 非依从性：患者并没有按照原来安排的治疗方案进行治疗；

(2) 缺失数据：没有观察到结局（结局有定义）；

(3) 死亡截断：在收集到结局之前患者死亡（结局无定义）。

2. 非标准条件下的因果推断之非依从性

用Z_i表示第i个个体被随机分配的处理方案，D_i(Z_i)表示个体i是否依从于分配（实际接受的处理），Y_i(Z_i,D_i(Z_i))表示潜在结果。在处理效应稳定假设（SUTVA）——个体之间不会互相干扰——下，Z对Y的个体因果作用为Y_i(1,D_i(1))- Y_i(0,D_i(0))，Z对Y的平均因果作用（意向治疗作用）为E[Y_i(1,D_i(1))- Y_i(0,D_i(0))]。如果我们关心的是D对Y的因果作用，也就是实际接受的处理对结局的因果作用，我们还需要更多的假设。

首先，我们把人群分为四层：依从组c（D_i(z)=z）、永远服药组a（D_i(z)=1）、永不服药组n（D_i(z)=0）、绝不依从组d（D_i(z)=1-z）。只有c组和d组包含了两种处理方案的结局数据，因此只能在这两组中定义因果作用，以c组为例，考虑依从组平均因果作用CACE= E[Y_i(1,1)- Y_i(0,0)]。观察到的数据为：分配方案Z_i、实际接受的处理D_i=D_i(Z_i)、结局Y_i=Y_i(Z_i)。需做如下假设：

(1) 可忽略性假设，即两种潜在结果和分配方案独立；

(2) 单调性假设，即D_i(1)≥D_i(0)，不存在d组；

(3) a组和d组的排他性约束假设，即a组和d组的两种潜在结果相等。

下面我们用一个实际例子来说明带非依从性的因果推断[1]。为了研究流感疫苗的作用，实施一个激励试验。这里，分配的方案是鼓励打疫苗、不鼓励打疫苗，但被鼓励打疫苗的人仍有可能不打疫苗，未被鼓励打疫苗的人也有可能打疫苗。传统的意向治疗分析是有问题的，因为这种做法得到的是鼓励打疫苗的作用，而不是实际打疫苗的作用。用工具变量把人群分层后，估计依从组的因果作用，才能代表打疫苗的真实作用。

周晓华和他的同事解决了随机临床试验中存在非依从性和不可忽略的结局缺失时研究参数的可识别性问题，提出了针对非依从性的贝叶斯分析方法，证明了在不同类型的完全不可忽略缺失数据下（即缺失机制依赖于结局），满足一定条件时，感兴趣的因果参数是可识别的，同时推导出了参数的最大似然估计和矩估计，并分析了它们在有限样本中的性质。当结局存在不可忽略缺失时，或者对于聚类激励试验，研究了多重填补方法[2-5]。

利用高维协变量和观测数据，在不存在强可忽略性假设的情况下，周晓华和合作者提出了异质局部治疗效应的新估计和推断方法[6]。针对两阶段广义线性模型，给出了非凸目标函数下的Lasso估计，并提出了一种协变量特异治疗效果置信区间的构造方法，这种方法同时纠正了由于两个阶段的高维估计而产生的偏差。这项研究成果即将发表在 JRSSB 上。

3. 非标准条件下的因果推断之死亡截断

和非依从性类似，死亡截断也会破坏经典的因果分析假设。在临床试验中，一些个体可能在收集到结局之前就发生死亡，这一现象被称为死亡截断。需要特别强调的是，死亡截断与缺失数据是两个完全不同的问题：前者的结局没有定义，而后者的结局有定义、只不过是未被观察到罢了。

用Z_i表示第i个个体被随机分配的处理方案（假设个体依从于分配方案），S_i(Z_i)表示个体i的潜在存活状态（1表示存活，0表示死亡），Y_i(Z_i)表示潜在结果（如果S_i(Z_i)=1），用X_i表示协变量。仍然利用主层分析的方法，把人群分为四层，用G表示：永远存活组LL（S_i(z)=1）、永远死亡组DD（S_i(z)=1）、有益组LD（S_i(z)=z）、有害组DL（S_i(z)=1-z）。只有永远存活的LL组，其因果参数是有意义的，因为对于其他组来说，两个潜在结果至少有一个是无定义的。因此，我们关心永远存活组的平均因果作用 SACE=E[Y_i(1)-Y_i(0)|G=LL]。

类似地，为了识别存活组平均因果作用，需要做出额外的假设：关于S和Y的可忽略性假设、单调性假设、排他性约束假设、替代相关性假设。通过工具变量对人群分层，进而使用参数模型估计出存活组的因果作用[7, 8]。在单调性假设下，有害组DL组被排除了。如果要放宽单调性假设，可将其替换为随机单调性假设，也就是允许DL组的存在，但需要事先给定S_i(1)、S_i(0)和LL组之间的概率关系。

周晓华和他的同事在国际上率先提出了用于超过三个组别且存在死亡截断的多处理随机临床试验的统计方法[9]。此外，周晓华和他的同事还发展了新的推断方法来检验总体治疗效果，并且证明了该方法在大样本下的收敛性，完善了大样本下该方法的统计理论。周晓华和他的同事还提出了适用于结果是二分类和连续型变量的情形下，在非参数和半参数模型中识别感兴趣因果参数SACE的方法。证明了SACE在部分正则假设下可识别的数学性质，同时提出当违背部分假设时，减少估计偏差的统计方法和理论。

今年，周晓华团队在《中国疾病预防控制中心周报》上的一篇文章回顾了用瑞德西韦治疗COVID-19的两个临床试验，这两个试验都涉及到死亡截断问题。例如，当以病情恢复时间为结局，则在恢复前被死亡截断的个体没有结局定义。如果死亡率与基础疾病存在正相关性，则把死亡个体当作在终点处删失的做法低估了真实的治疗效应。作者通过模拟实验表明，把死亡当作删失的方法以及存活案例分析方法（无论是否按协变量分层）都会导致有偏的治疗效应估计[10]。

4. 中介分析

在一些问题中，干预和结局之间可能存在中介。用T表示治疗方案，M表示中介，Y表示结局。这时，从T到Y存在两条路径：一条是从T直接到Y，称为直接因果作用；另一条是从T经过M再到Y，称为间接因果作用。

例如，在关于大学生酗酒的研究中，受试者被随机分配到三种方案T之一：放松（对照）、酒精扫描+干预、酒精扫描+干预+情绪控制技术，中介M是情绪控制能力，结局Y是3个月后的抑郁得分。经典的中介分析方法需要很强的序列可忽略性假设，要求中介和结局之间不能存在混杂。周晓华和合作者放宽了序列可忽略性假设，提出的新方法允许中介和结局之间存在混杂，并且允许存在缺失数据[11]。

今年周晓华团队在Statistics in Medicine发表的一篇文章关注了用确定结果随机事件模型研究中介分析，对于二值变量的情形，在一些识别性条件下，作者给出了直接因果作用和间接因果作用的估计，允许序列可忽略性假设不成立。另外，当识别性条件被破坏时，作者给出了直接因果作用和间接因果作用的上下界[12]。

5. 个性化治疗

个性化治疗（或精准医疗）与因果推断有着密切的关系。精准医疗旨在根据患者的特征，选择恰当的治疗方案，达到治疗收益的最大化。但是，同一位患者一次只能接受一种治疗，我们无法同时观察到患者接受不同治疗的结果。用因果推断的语言说，设X表示患者的协变量，Z为治疗方案，可以取Z=0或1，用Y(z)表示患者接受治疗z后的潜在结果。精准医疗却需要根据协变量特异因果作用E[Y(1)-Y(0)|X]选择相应的治疗方案。具体而言，假设潜在结果数值越大表示患者收益越大，那么如果上面的协变量特异因果作用大于0，就选择Z=1，否则选择Z=0。这样，我们可以从理论上把精准医疗转化为因果推断问题。

给定潜在结果的单指标模型，在一些假设下就可以识别出协变量特异因果作用。比如，当我们比较中药和西药的疗效时，可能会发现当生物标志物在一定范围内(c1<X<c2)时，中药治疗比西药治疗更有效，否则西药治疗更有效。这样，我们能画出一条曲线，用于描述中药相对于西药的因果作用随生物标志物变化的情形，即协变量特异因果作用曲线，这条曲线在临床应用中为研究者提供了很大便利[13]。

在最优个性化治疗规则的因果作用估计方面，周晓华和他的同事首次提出使用生物标志物调整效应曲线（BATE）、协变量特征的治疗效果曲线（CSTE）来表示给定生物标志物水平下的条件平均处理效应，并为每个病人选择最优治疗方案，同时严格证明了新提出统计方法的数学性质。针对治疗结果为二分类变量时，周晓华和他的同事提出采用B样条方法估计 CSTE 曲线，采用Hotelling tube方法估计CSTE曲线置信带的新统计方法及理论，定义并推导出二值响应变量数据的CSTE曲线可表示为优势比的对数，同时利用B样条方法估计广义变系数模型的系数，并推导出该方法的大样本数学性质[14, 15]。

今年周晓华团队在 Journal of the American Statistical Association 新发表的论文把协变量特异的治疗效果曲线推广到含有高维协变量的情形，在保证模型灵活性的同时允许高维变量选择，完善了协变量特异的治疗效果曲线估计量的数学理论。作者分别针对治疗组和控制组构造新的半参数模型，并引入惩罚项以实现变量选择，基于向后样条拟合核的极值估计量研究收敛速度，构造出同时置信带，支持最大化因果作用的个性化治疗策略的选择和推断[16]。作者还用R软件包实现了文中所描述的方法[17]。

参考文献

[1] Keisuke, H. , Imbens, G. W. , Rubin, D. B. , & Xiao-Hua, Z. . (2000). Assessing the effect of an influenza vaccine in an encouragement design. Biostatistics.

[2] Hua, C. , Geng, Z. , Zhou, X. H. , Small, D. S. , & Jing, C. . (2009). Identifiability and estimation of causal effects in randomized trials with noncompliance and completely nonignorable missing data. Biometrics.

[3] Frangakis, C. E. , Rubin, D. B. , & Zhou, X. . (2002). Clustered encouragement designs with individual noncompliance: bayesian inference with randomization, and application to advance directive forms. Biostatistics.

[4] Taylor, L. , & Zhou, X. H. . (2010). Multiple imputation methods for treatment noncompliance and nonresponse in randomized clinical trials. Biometrics.

[5] Taylor, L. , & Zhou, X. H. . (2011). Methods for clustered encouragement design studies with noncompliance and missing data. Biostatistics.

[6] Qiu, Y. , Tao, J. , & Zhou, X. H. . (2021). Inference of heterogeneous treatment effects using observational data with high-dimensional covariates. Journal of the Royal Statistical Society: Series B.

[7] Ding, P. , Geng, Z. , Yan, W. , & Zhou, X. H. (2011). Identifiability and estimation of causal effects by principal stratification with outcomes truncated by death. Journal of the American Statistical Association.

[8] Wang, L. , Zhou, X. H. , & Richardson, T. S. . (2016). Identification and estimation of causal effects with outcomes truncated by death. Biometrika.

[9] Wang, L. , Richardson, T. S. , & Zhou, X. H. . (2017). Causal analysis of ordinal treatments and binary outcomes under truncation by death. Journal of the Royal Statistical Society: Series B (Statistical Methodology).

[10] Deng, Y. , & Zhou, X. H. . (2021). Caution About Truncation-By-Death in Clinical Trial Statistical Analysis: A Lesson from Remdesivir. China CDC Weekly.

[11] Zheng, C. , & Zhou, X. H. . (2015). Causal mediation analysis in the multilevel intervention and multicomponent mediator case. Journal of the Royal Statistical Society Series B, 77.

[12] Li, W. , Geng, Z. , & Zhou, X. H. . (2021). Causal mediation analysis with sure outcomes of random events model. Statistics in Medicine.

[13] 韩开山, 周晓华, & 刘宝芳. (2017). 结局变量为二值时利用cste曲线选择最优治疗方案. 中国科学:数学.

[14] Zhou, X. H. , & Ma, Y. B. . (2012). Bate curve in assessment of clinical utility of predictive biomarkers. Science China Mathematics.

[15] Ma, Y. , & Zhou, X. H. . (2017). Treatment selection in a randomized clinical trial via covariate-specific treatment effect curves. Statistical Methods in Medical Research.

[16] Guo, W. , Zhou, X. H., & Ma, S. (2020). Estimation of optimal individualized treatment rules using a covariate-specific treatment effect curve with high-dimensional covariates. Journal of the American Statistical Association.

[17] https://cran.r-project.org/web/packages/CSTE/index.html

（参考文献可上下滑动查看）

2021泛太平洋因果推断大会

2021泛太平洋因果推断大会是自2019年起由北京大学公共卫生学院生物统计系，北京大学北京国际数学研究中心发起的因果科学领域一年一度的盛会。本次盛会将于9月11日至12日线上举办，面向所有对因果科学感兴趣的相关研究人员，盛会聚集了海内外大批因果科学相关领域的学术大拿，汇聚一堂共同探讨因果科学的最新进展，探索新的科研方向。本次泛太平洋因果推断大会为了满足对因果科学感兴趣的相关研究者，我们采取了提名+邀请嘉宾的方式，吸引力来大批业界大拿的关注，也为因果科学的研究进程提供更加开放平台。