Processing math: 17%

数学物理学报, 2025, 45(3): 992-1012

响应变量缺失下条件平均处理效应的k近邻核估计

曾华俊1,2, 明瑞星1,2, 苏培娟1,2, 黄绍航1,2, 肖敏,1,3,*

1浙江工商大学统计数据工程技术与应用协同创新中心 杭州 310018

2浙江工商大学统计与数学学院 杭州 310018

3浙江工商大学经济运行态势预警与模拟推演实验室 杭州 310018

k-Nearest Neighbor Kernel Estimation of Conditional Average Treatment Effect with Missing Response Variables

Zeng Huajun1,2, Ming Ruixing1,2, Su Peijuan1,2, Huang Shaohang1,2, Xiao Min,1,3,*

1Collaborative Innovation Center of Statistical Data Engineering Technology & Application, Zhejiang Gongshang University, Hangzhou 310018

2School of Statistics and Mathematics, Zhejiang Gongshang University, Hangzhou 310018

3Economic Forecasting and Policy Simulation Laboratory, Zhejiang Gongshang University, Hangzhou 310018

通讯作者: 肖敏,Email: xiaomin@zjgsu.edu.cn

收稿日期: 2024-06-18   修回日期: 2024-09-13  

基金资助: 浙江省社会科学规划课题(22GXSZ001Z)
数字+学科建设项目(SZJ2022B004)
浙江省教育厅一般项目(Y202353084)
浙江省登峰学科
Summit Advancement Disciplines of Zhejiang Province(Zhejiang Gongshang University-Statistics)
浙江省省属高校基本科研业务费专项资金(XT202302)

Received: 2024-06-18   Revised: 2024-09-13  

Fund supported: Zhejiang Provincial Philosophy and Social Sciences Planning Project(22GXSZ001Z)
Digital Science and Engineering Construction Project(SZJ2022B004)
Zhejiang Provincial Department of Education General Project(Y202353084)(Zhejiang Gongshang University-Statistics)
Fundamental Research Funds for the Provincial Universities of Zhejiang(XT202302)

摘要

基于 Neyman-Rubin 潜在结果框架, 构建 k 近邻核估计量来测度响应变量随机缺失情形下的条件平均处理效应 (conditional average treatment effect, CATE), 旨在评估不同处理方式对个体的影响. 证明了 k 近邻核估计量的几乎完全收敛性和渐近正态性. 数值模拟表明 k 近邻核估计量的表现优良, 利用真实数据进行实证分析, 实证结果显示 k 近邻核估计量具有较小的平均绝对偏差和均方根误差.

关键词: 条件平均处理效应; 随机缺失; k 近邻核估计量; 渐近正态性

Abstract

Under the Neyman-Rubin potential outcome framework, we construct a k-nearest neighbor kernel estimator to measure the conditional average treatment effect in the case of random missing response variables, aiming to evaluate the impact of different treatments on individuals. The paper proves the almost complete convergence and the asymptotic normality of the estimator. The numerical simulation shows that the k-nearest neighbor kernel estimator performs well. The real-world data is used for empirical analysis, and the empirical results show that mean absolute error and root mean square error of the k-nearest neighbor kernel estimator are smaller.

Keywords: conditional average treatment effect; random missing; k-nearest neighbor kernel estimator; asymptotic normality

PDF (1399KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

曾华俊, 明瑞星, 苏培娟, 黄绍航, 肖敏. 响应变量缺失下条件平均处理效应的k近邻核估计[J]. 数学物理学报, 2025, 45(3): 992-1012

Zeng Huajun, Ming Ruixing, Su Peijuan, Huang Shaohang, Xiao Min. k-Nearest Neighbor Kernel Estimation of Conditional Average Treatment Effect with Missing Response Variables[J]. Acta Mathematica Scientia, 2025, 45(3): 992-1012

1 引言

因果推断是统计学研究的核心领域之一, 其中处理效应分析近年来受到越来越多学者的关注. 因果推断的处理效应研究旨在评估某个处理 (例如介入, 治疗, 政策干预等) 对特定结果变量的因果效应. 通常, 学者们采用 Neyman-Rubin 潜在结果框架, 通过估计平均处理效应 (average treatment effect, ATE) 来探究处理变量对结果变量的影响[1,2]. 当处理效应存在异质性时, ATE 模型仅关注整个群体的平均效应, 忽略了不同个体对处理的反应差异. 为考虑个体处理效应的异质性, 学者们构建条件平均处理效应 (conditional average treatment effect, CATE) 来测度不同处理方式对个体的影响[3,4], 并被广泛地应用在精准医疗, 经济, 社会学等多个领域[5,6,7].

在日益重视个体差异和个性化处理的背景下, 为更好地适应个体特征的多样性和非线性关系, 提供更加灵活和稳健的估计, 非参数方法成为 CATE 估计的一种有力工具. 常见的 CATE 估计的非参数方法有 Nadaraya-Watson (NW) 核估计方法, k 近邻核估计法 (k-nearest neighbor kernel estimator, kNN), GBDT (gradient boosting decision trees), XGBoost (extreme gradient boosting) 和随机森林等. NW 核估计方法使用核函数对附近所有观测数据进行加权, 在处理效应研究中可以提供较为平滑的预测结果. Abrevaya 等[8]首次提出将亚组水平上的处理效应定义为条件平均处理效应 (CATE), 并采用 NW 核方法来估计 CATE, 证明了估计量的一致收敛性和渐近正态性. Li 等[9]利用 NW 核估计方法来估计处理效应的异质性, 并在理论上详细讨论了四种不同结构的模型下估计量的渐近正态性. 然而, NW 核估计方法的性能受带宽参数的选择影响, 因此, 选择合适的带宽是一个巨大的挑战. k 近邻核估计法同 NW 核方法一样都是基于核函数的非参数方法, 该方法是基于最近的 k 个邻居数据点进行加权平均, 通常使用欧氏距离或其他距离度量来衡量观测点之间的相似性. Kudraszow 和 Vieu[10] 利用 k 近邻核方法来估计核估计的带宽, 使得带宽可以根据协变量来自动选择, 证明了所得估计量的几乎一致收敛性. Ling 等[11]在 Kudraszow 和 Vieu[10] 的基础上, 将该方法推广到完整数据下协变量具有 α 混合下的函数型数据的 k 近邻核估计方法, 并证明了估计量的完全收敛性及其收敛速率. 与传统的 NW 核估计方法相比, k 近邻核估计法通过交叉验证寻找最优的离散参数 k, 可以避免手动调整带宽参数的主观性, 提供更客观和准确的参数选择.

需要说明的是, 虽然 NW 核方法和 k 近邻核估计法在处理异质性强, 非平稳的数据时具有优势, 但是, 当数据维数很大或数据集很大时, NW 核方法和 k 近邻核估计法的预测能力较 GBDT, XGBoost, 随机森林和深度学习等机器学习方法不一定具有优势. GBDT[12], XGBoost[13] 和随机森林[14]都是基于机器学习中集成学习的方法, 通过训练弱决策树模型来提高预测性能. 因果森林在随机森林的基础上引入了因果推断的思想, 是近年来条件平均处理效应研究的热门方法之一[15,16]. 上述非参数方法都能灵活地适应复杂的数据结构和特征关系, 捕捉到不同处理对个体的效应差异. 相较于 NW 核方法和 k 近邻核估计法, 尽管 GBDT, XGBoost 和随机森林对于高维数据和大规模数据具有较好的性能和效率, 但在处理数据时容易过拟合, 可解释性较差, 理论结果也较为薄弱. 其它机器学习方法, 例如深度学习, 也有类似弱点, 在此不作过多说明.

现实中, 在很多研究中存在响应变量数据缺失情况[17,18]. 数据缺失机制一般分为三种类型: 完全随机缺失 (missing completely at random, MCAR), 随机缺失 (missing at random, MAR) 和非随机缺失 (missing not at random, MNAR). 关于处理效应研究与缺失数据方法结合的最新进展, 见文献 [19,20,21]. 由于 MAR 假设情形下的缺失更为普遍, 本文采用 k 近邻核方法构建响应变量 MAR 缺失情形下的 k 近邻核估计量. 据笔者所知, 现有的文献尚未发现采用 k 近邻核方法研究缺失响应变量的 CATE, 从这种意义上来讲, 本文填补了对缺失情形下 CATE 的非参数估计的空缺. 本文证明 k 近邻核估计量是几乎完全收敛的, 并给出其收敛速度, 该结果可视为将 Kudraszow 和Vieu[10] 的结果从完整数据的均值估计情形推广到缺失数据的 CATE 估计情形. 此外, 本文还证明 k 近邻核估计量的渐近正态性. 本文的理论结果完善了响应变量 MAR 情形下的 CATE 估计量的大样本性质, 为理解和推动缺失情形下 CATE 的非参数研究发挥重要作用. 通过数值模拟, 本文将 k 近邻核估计量与 GBDT, XGBoost 进行对比, 说明响应变量 MAR 情形下 k 近邻核估计量的表现优良. 最后, 利用真实数据进行实证分析, 实证结果与模拟结果基本一致.

文章其余部分的结构如下. 第 2 节介绍响应变量随机缺失下 CATE 模型, 并构造出 CATE 估计量. 第 3 节展示估计量的几乎完全收敛性和渐近正态性. 第 4 节采用数值模拟来探索 k 近邻核估计量在样本量和缺失度不同情况下的表现, 并与 GBDT, XGBoost 两种方法进行比较. 第 5 节进行实证分析, 研究 k 近邻核估计量在真实数据中的表现. 第 6 节给出总结与展望. 附录给出几乎完全收敛性和渐近正态性的证明过程.

2 模型与估计

本文采用 Neyman-Rubin 潜在结果框架估计响应变量 MAR 情形下的条件平均处理效应. 对于个体 i{1,2,,n}, 假设其潜在结果不受其他个体是否分配处理的影响[22]. 考虑二值处理分配变量 Wi, Wi=1 表示个体 i 分配处理 (处理组), Wi=0 表示个体 i 不分配处理 (对照组). 设 (Yi(1),Yi(0)) 为个体 i 的潜在结果 (或响应变量), 其中 Yi(1) 对应个体 i 分配处理的结果变量, Yi(0) 对应个体 i 不分配处理的结果变量. 又设 δi(Wi) 表示 Yi(Wi) 的二值缺失指示变量, 当 Yi(Wi) 缺失时, δi(Wi)=0, 反之 δi(Wi)=1.Yi(Wi) 不存在缺失, 且可以观测到个体 i 的两种潜在结果, 此时个体处理效应 (individual treatment effect, ITE) 为个体 i 的两个潜在结果之差, 即 Yi(1)Yi(0). 注意到对于每个个体 i, 不可能同时观测到 Yi(1)Yi(0), 因此可将完整数据下的实际观测结果 Yobsi 表示为 Yobsi=WiYi(1)+(1Wi)Yi(0). 在假设 Yi(Wi) 存在缺失的条件下, 实际观测结果为 Yi=Wiδi(Wi)Yi(1)+(1Wi)δi(Wi)Yi(0).

给定可观测数据集 Di={Yi,Xi,Wi,δi(Wi),i=1,,n}, 其中协变量 XiF, 通常 FRd. 因果推断中, 协变量 Xi 是用来控制混杂因素的关键变量, 关于混杂因素的判定和控制可参见文献 [23,24]. 设处理组样本量为 n1, 对照组样本量为 n0, 这里 n0=nn1. 本文关注响应变量在 MAR 情形下的条件平均处理效应 (CATE)

τ(x)=E[Y(1)Y(0)X=x].
(2.1)

为估计 (2.1) 式, 需对处理分配机制和响应变量 MAR 机制作出必要的假设. 处理分配机制满足以下可忽略性假设[3].

(HT1 ) 给定 X, 有 Y(0), Y(1)W 相互独立, 即 Y(0),Y(1)WX;

(HT2 ) 给定 X, 记 π:=π(X)=P(W=1X) 是倾向得分, π 满足 0<π<1.

由 (HT1) 和 (HT2) 知, P(W=1Y(W),X=x)=P(W=1X=x).pπ:=P(W=1), 易知 n1 约为 npπ, n0 约为 n(1pπ).

在响应变量MAR机制中, 缺失指示变量满足以下假设[25].

(HM1 ) 给定 XW, 变量 Y(W)δ(W) 是独立的, 即 δ(W)Y(W)(X,W);

(HM2 ) 给定 XW, 记 θ:=θ(W,X)=P(δ(W)X,W), θ 满足 0<θ<1.

由 (HM1)-(HM2) 知, P(δ(W)=1Y(W),X=x,W=w)=P(δ(W)=1X=x,W=w).pθ(w):=P(δ(W)=1W=w), 易知 Y(1) 的可观测样本量约为 n1pθ(1), Y(0) 的可观测样本量约为 n0pθ(0).

根据假设 (HT1), (HT2), (HM1) 和 (HM2), 容易算得

τ(x)=μ1(x)μ0(x),
(2.2)

其中, μw(x)=E[YW=w,X=x,δ(W)=1],w=0,1. μ1(x)μ0(x) 分别是响应变量 MAR 情形下处理组和对照组的条件平均潜在结果, 其 k 近邻核估计量为

ˆμw,kNN(x)=ni=1I(Wi=w)K(Hnw,kw(x)1d(x,Xi))δi(Wi)Yi(Wi)ni=1I(Wi=w)K(Hnw,kw(x)1d(x,Xi))δi(Wi),xF,w=0,1,
(2.3)

其中, K()[] 是核函数; d(x,Xi) 是半度量; B(x,hw) 是与半度量 d 相关的拓扑中心为 x, 半径为 hw 的球体, 即 B(x,hw)={uFd(u,x)hw},w=0,1; Hnw,kw(x)

Hnw,kw(x)=min{hwR+:ni=1I(Wi=w)IB(x,hw)(Xi)=kw},w=0,1.
(2.4)

在 (2.4) 式中, k1k0 表示处理组和对照组中用于估计的最近邻样本量. 本文采用 5 折交叉验证来获得最优的 k1k0 值. 具体过程如下. 首先, 将数据集随机均等地分成 5 个子集, 依次选择 1 个子集作为验证集, 其余 4 个子集作为训练集. 接着, 将 k0k1 的候选范围分别设定为 [2,1.5×I(Wi=0)I(δ(Wi)=1)n][2,1.5×I(Wi=1)I(δ(Wi)=1)n] 的整数部分. 此设定能够根据数据集的大小自动调整 k 值的上限, 使其适应不同样本量的情况. 然后, 在每一折交叉验证中, 使用 4 个子集作为训练集训练模型, 然后在剩余的验证集上计算均方误差. 重复此过程, 遍历所有的 k0k1 组合, 计算每次验证集上的均方误差. 最后, 计算不同 k0k1 组合 5 次交叉验证的平均均方误差, 选择使平均均方误差最小的一组 k0k1 作为最优值, 从而确保模型的最优性能.

这样, τ(x)k 近邻核估计量 ˆτkNN(x)

ˆτkNN(x)=ˆμ1,kNN(x)ˆμ0,kNN(x).
(2.5)

3 渐近结果

本节研究 CATE 的 k 近邻核估计量的几乎完全收敛性和渐近正态性. 首先, 给出几乎完全收敛性和 Kolmogorov 熵的定义.

定义 3.1{Tn(x);n1} 是定义在集合 F 上的函数列, T(x) 也是定义在集合 F 上的函数. 若存在随机变量 un, 满足

supxFTn(x)T(x)∣=Oa.co.(un),

其中, 符号 “a.co” 表示几乎完全收敛性 (almost complete convergence), 则称 \{T_n(x); n \geq 1\}\mathcal{F} 上几乎完全收敛到 T(x), 且收敛速率为 u_n.

定义 3.2S_{w,\mathcal{F}},w \in \{0,1\}, 是 \mathcal{F} 的子集, 其 Kolmogorov 熵为

\psi_ {S_{w,\mathcal{F}}}(\varepsilon_w)=\mathrm{log}(N_{\varepsilon_w}(S_{w,\mathcal{F}})),

其中, N_ {\varepsilon_w}(S_{w,\mathcal{F}}) 表示 \mathcal{F} 中能覆盖集合 S_{w,\mathcal{F}}, 半径为 \varepsilon_w 的最小开球数.

在介绍本文的主要结果前, 先给出一些必要的假设条件.

假设 3.1 对任意 \varepsilon_w>0, 存在非负函数 \phi(\cdot), 满足 \phi(0)=0, \lim\limits_ {\varepsilon_w\rightarrow0}\phi(\varepsilon_w)=0, 且存在 \xi_0>0 和常数 C, 使得对任意的 \xi \in (0,\xi_0), \phi^{\prime}(\xi)<C. 此外, 存在 \xi\geq0, 使得对任意 z\in[0,1], 满足 \lim\limits_ {\varepsilon_w\rightarrow0}\bigl(\phi(u\varepsilon_w)/\phi(\varepsilon_w)\bigr)=z^\xi.

假设 3.2 对任意 \varepsilon_w>0, 记 \varphi_x(\varepsilon_w):=\mathbb{P}(X\in B(x,\varepsilon_w)). 假设 \varphi_x(\varepsilon_w)0 的邻域内连续, 且 \varphi_x(\varepsilon_w)>0, \varphi_x(0)=0. 此外, 存在函数 \phi(\cdot) (满足假设 3.1), 有界的正函数 f(x), 常数 a_w>0, 使得当 {\varepsilon_w}\rightarrow0 时, 有\sup\limits_{x\in S_{w,\mathcal{F}}}\bigl| \varphi_x(\varepsilon_w)/\phi(\varepsilon_w)-f(x)\bigr|=O\bigl({\varepsilon_w}^{a_w}\bigr).

假设 3.3\mu_w(x)S_{w,\mathcal{F}}b_w 阶有界的 Lipschitz 函数, 即存在常数 b_w>0 和常数 C, 对任意 x_1,x_2\in\ S_{w,\mathcal{F}}, 有 \bigl| \mu_w(x_1)-\mu_w(x_2)\bigr| \leq Cd^{b_w}(x_1,x_2).

假设 3.4 对任意 l\geq2, 存在 S_{w,\mathcal{F}} 上连续的有界函数 \sigma_{w,l}(x), 使得

\mathbb{E}\left(\mid Y(W)\mid^l\mid W=w, X=x\right)<\sigma_{w,l}(x).

假设 3.5 核函数 K 是非负有界的单调不增函数, 且在支撑集 [0,1] 上 Lipschitz 连续. 如果设 K(1)=0, 那么 K(t) 的导数 K^{\prime}(t) 必满足 -\infty< C \leq K^{\prime}(t) \leq \tilde{C} <0, 其中符号 C\tilde{C} 表示常数.

假设 3.6 存在 \lambda>1, 使得 S_{w,\mathcal{F}} 的 Kolmogorov 的 \varepsilon_w 熵满足

\sum\limits_ {n_w=1}^\infty\exp\Biggl((1-\lambda)\psi_{S_{w,\mathcal{F}}}\biggl(\frac{\log {n_w}}{n_w}\biggr)\Biggr)<\infty, w=0,1.

假设 3.7 对任意 i, j \in \{1,\cdots,n\}, i\neq j, 记 g(u,v)=\bigl[(Y_i(1)-\mu_1(u))(Y_j(0)-\mu_0(v))\mid X_i=u,X_j=v\bigr],u,v \in \mathcal{F}. 假设 g(u,v)(x,x) 的某一邻域内连续, 即当 h_1\rightarrow0, h_0\rightarrow0 时, 有

\sup\limits_{u\in B(x,h_1),v\in B(x,h_0)}\bigl| g(u,v)-g(x,x)\bigr|=o(1).

假设 3.8g_w(u)=\mathrm{Var}\left[Y(W)\mid W=w,X=u\right],u \in \mathcal{F}. 假设 g_w(u)x 的某一邻域内连续, 即当 h_w\rightarrow0 时, 有

\sup\limits_{u\in B(x,h_w)}\bigl| g_w(u)-g_w(x)\bigr|=o(1), w=0,1.

假设 3.9 对任意 i, j \in \{1,\cdots,n\}, i\neq j, 在 X_i=x 给定的情况下, Y_i(1)Y_j(0) 的相关系数为 \rho(x)\in[-1,1], 此时 g(x,x)=\rho(x) \sqrt{g_1(x)g_0(x)}.

注 3.1 假设 3.1 在零点邻域内可将函数 \phi(\cdot) 视为 Lipschitz 函数. 假设 3.2 是对函数空间拓扑结构的一般表示, 通过函数 \phi(\cdot) 来控制函数变量在小球上的概率浓度. 假设 3.3-3.5 是标准的非参数技术设置. 假设 3.6 是关于 S_{w,\mathcal{F}} 的 Kolmogorov 的 \varepsilon_w 熵的条件. 假设3.1-3.6 的详细说明可参见文献 [10,26]. 假设 3.7 是给定 X_iX_j 下对 Y_i(1)Y_j(0) 协方差的连续性假设, 假设 3.8 是给定 X_i 下对 Y_i(W_i) 条件方差的连续性假设, 假设 3.9 是对 Y_i(1)Y_j(0) 相关关系的假设.

在上述假设条件的基础上, 还需引入一些关键的记号. 记

K_{i, w}(x)=K\left(H_{n_w,k_w}(x)^{-1}d\left(x, X_{i}\right)\right), x\in S_{w,\mathcal{F}}, w=0,1.

B_{n_1}(x)B_{n_0}(x) 分别表示 \mu_1(x)\mu_0(x) 的估计偏差, 易得

B_{n_w}(x)=\frac{\mathbb{E}[\tilde{\mu}_{2n_w}(x)]}{\mathbb{E}[\tilde{\mu}_{1n_w}(x)]}-\mu_w(x), w=0,1,

其中

\tilde{\mu}_{1 n_w}(x)=\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w) K_{i, w}(x)\delta_{i}(W_i)}{n_w \mathbb{E} \left[K_{i, w}(x)\right]},\qquad\qquad\ w=0,1,
(3.1)
\tilde{\mu}_{2 n_w}(x)=\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w) K_{i, w}(x)\delta_{i}(W_i)Y_i(W_i)}{n_w \mathbb{E} \left[K_{i, w}(x)\right]}, w=0,1.
(3.2)

令 CATE 估计偏差 \tilde{\tau}(x)

\tilde{\tau}(x)=\sqrt{n}(\hat{\tau}_{kNN}(x)-\tau(x)-(B_{n_1}(x)-B_{n_0}(x))).

下面给出本文的主要结果.

定理 3.1 [几乎完全收敛性] 在假设 3.1-3.6 和响应变量 MAR 情形下, 对足够大的 n, 及任意满足条件当 n\rightarrow\infty 时, {k_w/n_w}\rightarrow0 的正实数列 k_w, 有

\frac{(\log n_w)^2}{k_w}<\psi_{S_{w,\mathcal{F}}}(\frac{\log n_w}{n_w})<\frac{k_w}{\log n_w},

此时,

\sup\limits_{x\in S_{w,\mathcal{F}}}\left|\hat{\tau}_{kNN}(x)-\tau(x)\right|=O_{\rm a.co}\left(u_{n_1}+u_{n_0}\right),
(3.3)

其中

u_{n_w}=\left(\bigg(\phi^{-1}\left(\frac{k_w}{n_w}\right)\bigg)^{b_w}+\sqrt{\frac{\psi_{S_{w,\mathcal{F}}}(\log n_w/n_w)}{k_w}}\right).

定理 3.2 [渐近正态性] 在假设 3.2, 3.4,3.7-3.9 和响应变量 MAR 情形下, 当 n\rightarrow\infty 时, 有

\frac{\tilde{\tau}(x)}{\sqrt{\sigma_0^2+\sigma_1^2}} \stackrel{L}{\rightarrow} N\left(0, 1\right),
(3.4)

其中

\sigma_1^2=\frac{g_1(x)\mathbb{E}[K_{1,1}^2(x)]}{p_{\pi}\mathbb{E}^2[K_{1,1}(x)]}, \sigma_0^2=\frac{g_0(x)\mathbb{E}[K_{n,0}^2(x)]}{(1-p_{\pi})\mathbb{E}^2[K_{n,0}(x)]}.

4 数值模拟

本节采用数值模拟实验来探索 k 近邻核估计量在样本量不同, 缺失度不同时的表现, 并与机器学习中的 GBDT, XGBoost 作比较. 4.1 节主要介绍 k 近邻核估计的模拟过程和参数设定, 4.2 节展示模拟结果, 并对结果进行分析.

4.1 模拟过程与参数设定

首先介绍 k 近邻核估计的数值模拟过程, 再分别对 k 近邻核估计, GBDT, XGBoost 进行参数设定.

(一) k 近邻核估计的模拟过程

步骤一 生成协变量和响应变量. 给定处理组样本量 n_1, 生成 n_1 个处理组协变量 X_i, 响应变量为 Y_i(1); 给定对照组样本量 n_0, 生成 n_0 个对照组的协变量 X_i, 对应的响应变量为 Y_i(0). 表 1 是协变量和响应变量的具体生成方式. 表 1 中的 \varepsilon_i 满足 \varepsilon_i \sim N(0,1), i \in \{1,2,\cdots,n\}, 这里 n=n_1+n_0. 在不同生成方式下, \tau(X_i) 为真实的条件平均处理效应 (CATE).

表 1   模拟数据生成方式

新窗口打开| 下载CSV


步骤二 生成缺失变量. 给定处理组缺失度 p_{\theta}(1) 和对照组缺失度 p_{\theta}(0), 生成缺失指示变量 \delta_i(1)\delta_i(0), 其中 \delta_i(1) \sim b(1,p_{\theta}(1)), \delta_i(0) \sim b(1,p_{\theta}(0)).

步骤三 选择 k_1k_0 值. 以均方误差最小化为标准, 用 5 折交叉验证方法来选择最优的 k_0k_1. 选择 Epanechnikov 核函数 K(t)=\frac{3} {4}(1-t^2)\mathbb{I}(\mid t\mid \leq 1) (满足假设 3.5), 距离为 d(x,u)=\sqrt{(x-u)^2}.

步骤四 计算评价标准. 选定 k_0k_1 后, 利用 X-Learner 计算 \hat{\tau}_{kNN}(x), 并与真实值 {\tau}(x) 作对比, 通过 100 次模拟求出平均绝对偏差 (MAE) 和均方根误差 (RMSE),

\mathrm{MAE} = \frac{\sum\limits_{i=1}^{100} \mid \hat{\tau}_{kNN}(x)-\tau(x)\mid}{100},\quad \mathrm{RMSE} = \sqrt{\frac{\sum\limits_{i=1}^{100}\bigl(\hat{\tau}_{kNN}(x)-\tau(x)\bigr)^2}{100}}.

(二) 参数设定

在模拟过程中, 数据生成分布和缺失指示变量分布是随机生成的. 为评估 k 近邻核估计量在不同数据分布, 样本量及缺失度设置下的表现, 需要设定的参数为处理组样本量 n_1 和缺失度 p_{\theta}(1), 对照组样本量 n_0 和缺失度 p_{\theta}(0). 模拟设置 n_1 取 100, 200, 400, 600, p_{\theta}(1) 取 0.1, 0.2, 0.3, 0.5, 0.7, 0.9, n_0 取 100, 200, 400, 600, p_{\theta}(0) 取 0.1, 0.2, 0.3, 0.5, 0.7, 0.9. 通过这些参数设定, 分别讨论 k 近邻核估计量在一维正态, 二维正态, 一维重尾, 二维重尾时的模拟结果, 能够系统地评估 k 近邻核估计量在不同数据分布, 样本量和缺失度下的稳定性和精确性, 特别是在高缺失度和复杂分布情况下的表现.

为减少 GBDT, XGBoost 方法之间的差异, 调整 GBDT 和 XGBoost 的参数: 迭代次数为 100, 最大深度为 6, 学习率为 0.1, 其他参数设置为默认值. 这些参数能够在模型复杂性, 计算资源和预测性能之间取得良好的平衡, 确保模型在训练过程中有效收敛并避免过拟合. 在此基础上, 根据 X-Learner 预测 CATE, 并与真实值比较求出 GBDT, XGBoost 的 MAE 和 RMSE.

4.2 模拟结果与分析

通过模拟, 可以得到 k 近邻核估计量, GBDT, XGBoost 在各种情况下的 MAE 和 RMSE. 由于情况复杂和篇幅有限, 本节仅选取数据生成方式为一维正态分布的部分数据来具体说明 k 近邻核估计量, GBDT, XGBoost 的表现, 并进一步探索 k 近邻核估计量在样本量以及缺失度不同时的表现趋势.

(一) 样本量的影响

分析 k 近邻核估计量, GBDT, XGBoost 在样本量不同时的表现. 完整的数值模拟结果表明样本量不同的情况下, k 近邻核估计量的表现绝大部分优于 GBDT 和 XGBoost 的表现. 为具体说明, 选取 n_0=100, p_ {\theta}(0)=0.1, 数据生成方式为一维正态分布的数据, 比较 k 近邻核估计量, GBDT, XGBoost 在 n_1=100,200,400,600 时的 MAE 和 RMSE. 结果显示, 绝大部分情况下 k 近邻核估计量的 MAE 和 RMSE 低于 GBDT 和 XGBoost 的 MAE 和 RMSE, 该结果说明在样本量不同的情况下, 与 GBDT 和 XGBoost 相比, k 近邻核估计具有较好的表现.

为探索 k 近邻核估计量在不同样本量下的表现, 得出 k 近邻核估计量的 MAE 和 RMSE 随 p_{\theta}(1) 变化的折线图. 由图 1 可知, 随着 n_1 变大, 不同 p_{\theta}(1)k 近邻核估计量的 MAE 和 RMSE 都呈下降趋势. 此外, 不同样本量下 k 近邻核估计量的 MAE 基本随 p_{\theta}(1) 的增大而减小, 仅在 p_{\theta}(1)=0.9 时 MAE 有细微偏差, 而 RMSE 都随着 p_{\theta}(1) 的增大而增大. 该结果说明样本量越大, 采用 k 近邻核方法估计响应变量随机缺失下的条件平均处理效应越为准确.

图1

图1   (a) n_0=100, p_{\theta}(0)=0.1, 数据生成方式为一维正态分布时, 不同样本量下 k 近邻核估计量随 p_{\theta}(1) 变化的 MAE; (b) n_0=100, p_{\theta}(0)=0.1, 数据生成方式为一维正态分布时, 不同样本量下 k 近邻核估计量随 p_{\theta}(1) 变化的 RMSE


(二) 缺失度的影响

分析 k 近邻核估计量, GBDT, XGBoost 在缺失度不同时的表现. 完整的数值模拟结果表明与 GBDT 和 XGBoost 相比, 缺失度不同时 k 近邻核估计量的表现与样本量不同时的表现基本一致, 仅在缺失度较大时 k 近邻核估计量的表现较差. 为具体说明, 选取 n_0=100, n_1=100, 数据生成方式为一维正态分布的数据, 比较 k 近邻核估计量, GBDT 和 XGBoost 三种方法在 p_ {\theta}(0)=0.1,0.2,0.3,0.5,0.7,0.9 时的 MAE 和 RMSE. 结果显示不同缺失度下 k 近邻核估计量的 MAE 明显低于 GBDT, XGBoost 的 MAE, 而 k 近邻核估计量的 RMSE 在 p_{\theta}(1)=0.7 之前基本低于 GBDT 和 XGBoost 方法, 但在 p_{\theta}(1)=0.7p_ {\theta}(1)=0.9 时多高于 GBDT 和 XGBoost, 表明 k 近邻核估计量在缺失度不太大时表现优于 GBDT 和 XGBoost.

为探索 k 近邻核估计量在不同缺失度下的表现, 图 2 表示不同缺失度 (p_{\theta}(0))k 近邻核估计量的 MAE 和 RMSE 随 p_ {\theta}(1) 变化的折线图. 通过观察得知, 随着 p_ {\theta}(0) 变大, 不同 p_{\theta}(1)k 近邻核估计量的 MAE 呈下降趋势, 且 p_ {\theta}(0) 越大, MAE 的下降趋势越为显著, 而 RMSE 绝大部分会随着 p_{\theta}(0) 增大而增大, 且 p_ {\theta}(0) 较大时, RMSE 的增长趋势较为明显. 该结果进一步表明, k 近邻核估计量具有低偏差的特点, 且缺失度非常大时, k 近邻核估计量的表现越差.

图2

图2   (a) n_0=100, n_1=100, 数据生成方式为一维正态分布时, 不同缺失度下 k 近邻核估计量随 p_{\theta}(1) 变化的 MAE; (b) n_0=100, n_1=100, 数据生成方式为一维正态分布时, 不同缺失度下 k 近邻核估计量随 p_{\theta}(1) 变化的 RMSE


综上, 通过比较样本量和缺失度不同情形下的 MAE 和 RMSE 可知, 较于 GBDT 和 XGBoost 两种机器学习方法, 采用 k 近邻核方法来估计响应变量 MAR 情形下的 CATE 具有明显的优良性. 此外, 样本量越大, k 近邻核估计量的表现越好, 但在缺失度非常大时 k 近邻核估计量的 RMSE 较大. 其他模拟数据与该结论呈现近似的走势和形态, 在此不再赘述.

5 实证分析

本节利用真实数据研究响应变量 MAR 情形下, 协变量分别为一维和二维时 k 近邻核估计量, GBDT 和 XGBoost 的表现. 本文使用来自婴儿健康与发展项目 (Infant Health and Development Program, IHDP) 的半合成数据集[27]. 该数据集收集了 1985-1988 年美国早产儿随机试验的 747 条完整记录, 其中 treatment 是处理分配指标, treatment = 1 表示处理组, treatment = 0 表示对照组; \rm y\_factual 是完整的潜在结果, 不存在缺失情况; x1-x25 表示协变量, x1-x6 为连续变量, x7-x25 为离散变量, 涉及母亲和婴儿的各项特征; mu0 和 mu1 分别是合成的对照组和处理组的条件平均潜在结果.

在使用 IHDP 数据集研究之前, 需考虑协变量选取和 \rm y\_factual 缺失数据的生成. 根据图 3, 从连续变量 x1-x6 之中选择协变量, 这些变量都近似服从标准正态分布. 因此, 在本文中选择 x1 作为一维协变量, 选择 x1 和 x2 作为二维协变量. 由于 \rm y\_factual 是完整的数据, 因此本文通过随机删除一定比例的响应变量数据来模拟实证分析中可能遇到的随机缺失情形. 具体来说, 设定缺失度, 标记出缺失的响应变量并将其删除, 同时生成缺失指示变量以满足响应变量 MAR 条件, 这里设置缺失度为 0.1, 0.2, 0.3, 0.5, 0.7. 需要说明的是, 当缺失度为 0.9 时, 缺失的数据量非常大, 生成的缺失样本可能无法获得有效的信息, 因此不做考虑. 此外, k 近邻核估计量中核函数的选取, k_0k_1 的交叉验证方式, GBDT 和 XGBoost 的参数设置都与模拟一致, 这里不再重复. 不同的是求解 MAE 和 RMSE 时, 将 IHDP 数据集中预先给定的 mu0 和 mu1 作为真实值, 并与 k 近邻核估计量, GBDT 和 XGBoost 的预测值作对比.

图3

图3   IHDP 数据集中 x1-x6 的直方图和密度曲线


本文将生成的缺失数据按照 7 : 3 的比例划分为训练集和测试集, 分别得出不同缺失度下 k 近邻核估计量, GBDT, XGBoost 三种方法的 MAE 和 RMSE. 表 2 是不同缺失度下 IHDP 数据的 MAE 和 RMSE 汇总表.表 2 可知, 无论是一维协变量还是二维协变量, 不同缺失度情况下 k 近邻核估计量的 MAE 和 RMSE 都小于 GBDT 和 XGBoost 的 MAE 和 RMSE, 这表明响应变量 MAR 情形下 k 近邻核估计量的表现更优. 此外, 随着缺失度的增加, k 近邻核估计量的 MAE 逐渐减小, 这与模拟结果一致. 另外, k 近邻核估计量的 RMSE 也随着缺失度的增加而减小, 这可能是因为在缺失响应变量的 IHDP 数据集中, 较大的缺失度减少了邻近样本点的数量, 从而减小估计方差, 使得 k 近邻核估计量的 RMSE 减小.

表 2   不同缺失度下 IHDP 数据的 MAE 和 RMSE 汇总表

新窗口打开| 下载CSV


6 总结与展望

本文基于 Neyman-Rubin 潜在结果框架, 构建了响应变量随机缺失情形下的条件平均处理效应的 k 近邻核估计量, 证明了该估计量的几乎完全收敛性和渐近正态性. 数值模拟上, 选取样本量分别为 100, 200, 400, 600, 缺失度分别为 0.1, 0.2, 0.3, 0.5, 0.7, 0.9, 数据生成分布分别为一维正态, 二维正态, 一维重尾, 二维重尾的数据, 以平均绝对偏差和均方根误差作为衡量标准, 模拟结果显示绝大部分情况下响应变量随机缺失情形下 k 近邻核估计量的表现优于 GBDT 和 XGBoost. 实证上, 利用 IHDP 数据集比较 k 近邻核估计量, GBDT 和 XGBoost 的平均绝对偏差和均方根误差, 说明 k 近邻核估计量都具有较小的平均绝对偏差和均方根误差.

本文提出了响应变量随机缺失情形下的 CATE 的估计方法. 然而, 该方法没有考虑其他缺失情况和缺失机制, 如协变量的随机缺失和响应变量的非随机缺失等. 此外, 在模拟中需要选择两个不同的 k_0k_1, 导致模拟时间增加, 这也是该估计量在实际应用中的缺陷. 从理论框架上看, 虽然 k 近邻核估计量具有较好的理论基础, 但与现有的机器学习方法相比可能仍有一定的局限性. 相比之下, 随机森林等机器学习方法相对更加灵活, 可能在估计缺失响应变量的 CATE 时表现较优. 关于 k 近邻核估计量与其他机器学习方法的比较, 本文不作详细比较. 后续工作将对上述三点进行深入探究. 由于篇幅有限, 本文不再过多讨论.

附录

先介绍一些符号. 设 (\mathcal{F},\mathscr{A}) 是一般可测空间, (\mathbb{R},\mathscr{B}(\mathbb{R})) 是 Borel 可测空间, (A_{i,w},B_{i,w})\in \mathscr{A}\times\mathscr{B}(\mathbb{R}) 是随机序偶, i=1,2,\cdots,n, w=0,1. 又设 G_w 是定义在 (\mathbb{R}, S_{w,\mathcal{F}}\times \mathcal{F}) 上的非负实值函数, 且 G_w 关于第一个变元是单调不减的, 这里 S_ {w,\mathcal{F}}\subseteq\mathcal{F}, w=0,1.

为证明定理 3.1, 需要引理 A.1 和引理 A.2, 其中引理 A.1 是引理 A.2 的基础, 利用引理 A.2 可以直接证明定理 3.1. 此外, 为了叙述的简洁, 定理 3.1 中关于几乎完全收敛性 (“a.co”) 的极限都是在 n 趋向于无穷大时求得.

事实上, 对于响应变量 MAR 情形下条件平均潜在结果 \mu_{w}(x) 的估计, k 近邻核方法可视为对传统核方法在邻域上的扩展, 可以通过寻找最优 k 值来避免手动调整带宽 h_{n_w}(x), 以达到与传统核估计量相同的最优收敛速度. 为得到 k 近邻核估计量的几乎完全收敛性, 先给出响应变量 MAR 情形下的传统核估计量, 定义为

\hat{\mu}_{w}(x)=\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)K(h_{n_w}(x)^{-1}d(x,X_i))\delta_i(W_i)Y_i(W_i)}{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)K(h_{n_w}(x)^{-1}d(x,X_i))\delta_i(W_i)}, \forall x \in \mathcal{F}, w= 0,1,

其中, h_{n_w}(x) 是正实数序列, 满足 \lim\limits_{n_w\rightarrow\infty}h_{n_w}(x)=0. 利用该估计量容易得到引理 A.2 中的收敛速度.

下面的引理 A.1 证明了响应变量 MAR 情形下 \mu_{w}(x) 的传统核估计量的几乎完全收敛性.

引理 A.1 在假设3.1-3.6 和响应变量 MAR 情形下, 当 n 足够大时, 存在与 x 无关的序列 h_{n_w}, w=0,1, 使得 h_ {n_w}(x) 满足

0<Ch_{n_w}\leq\inf\limits_{x\in S_{w,\mathcal{F}}}h_{n_w}(x)\leq\sup\limits_{x\in S_{w,\mathcal{F}}}h_{n_w}(x)\leq\tilde{C}h_{n_w}<\infty.

此外, 当 n\rightarrow\infty 时, h_{n_w}\rightarrow0, 且对足够大的 n, 有

\frac{\mathrm{log}^2 n_w}{n_w\phi(h_{n_w})}<\psi_{S_{w,\mathcal{F}}}\bigg(\frac{\mathrm{log} n_w}{n_w}\bigg)<\frac{n_w\phi(h_{n_w})}{\mathrm{log} n_w}.
(A.1)

此时

\begin{matrix} \sup\limits_{x\in S_{w,\mathcal{F}}}\mid \hat{\mu}_{w}(x)-\mu_w(x)\mid=O\left(p_{\theta}(w)h_{n_w}^{b_w}\right)+O_{\rm a.co}\left(p_{\theta}(w) \sqrt{\frac{\psi_{S_{w, \mathcal{F}}}(\mathrm{log} n_w/n_w)}{n_{w} \phi\left(h_{n_w}\right)}}\right). \end{matrix}
(A.2)

注意到

\begin{align*} \hat{\mu}_{w}(x)-\mu_w(x)&=\frac{1}{\hat{\mu}_{1 n_w}(x)}\bigl[\hat{\mu}_{2n_w}(x)-\mathbb{E}\left[\hat{\mu}_{2n_w}(x)\right]\bigr]\\ &~~~+\frac{1}{\hat{\mu}_{1n_w}(x)}\bigl[\mathbb{E} \left[\hat{\mu}_{2n_w}(x)\right]-p_{\theta}(w)\mu_w(x)\bigr] +\frac{\mu_w(x)}{\hat{\mu}_{1n_w}(x)}\bigl[p_{\theta}(w)-\hat{\mu}_{1 n_w}(x)\bigr], \end{align*}

其中

\begin{matrix} \hat{\mu}_{1 n_w}(x)=\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)\left[K\left(\frac{d\left(x, X_{i}\right)}{h_{n_w}(x)}\right)\right]\delta_{i}(W_i)}{n_w \mathbb{E} \left[K\left(\frac{d\left(x, X_{i}\right)}{h_{n_w}(x)}\right)\right]}, \end{matrix}
(A.3)
\begin{matrix} \hat{\mu}_{2 n_w}(x)=\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)\left[K\left(\frac{d\left(x, X_{i}\right)}{h_{n_w}(x)}\right)\right]\delta_{i}(W_i)Y_i(W_i)}{n_w \mathbb{E} \left[K\left(\frac{d\left(x, X_{i}\right)}{h_{n_w}(x)}\right)\right]}. \end{matrix}
(A.4)

这样证明 (A.2) 式等价于证明以下三式成立

\begin{matrix} &\sup _{x \in S_{w, \mathcal{F}}}\bigl|\mathbb{E}[ \hat{\mu}_{2 n_w}(x)]-p_{\theta}(w) \mu_{w}(x)\bigr|=O\left(p_{\theta}(w)h_{n_w}^{b_w}\right), \end{matrix}
(A.5)
\begin{matrix} &\sup _{x \in S_{w, \mathcal{F}}}\bigl|\hat{\mu}_{1 n_w}(x)-p_{\theta}(w)\bigr|=O_{\rm a.co}\left(p_{\theta}(w) \sqrt{\frac{\psi_{S_{w, \mathcal{F}}}(\mathrm{log} n_w/n_w)}{n_{w} \phi\left(h_{n_w}\right)}}\right),\end{matrix}
(A.6)
\begin{matrix} &\sup _{x \in S_{w, \mathcal{F}}}\bigl|\hat{\mu}_{2 n_w}(x)-\mathbb{E}[\hat{\mu}_{2 n_w}(x)]\bigr|=O_{\rm a.co}\left(p_{\theta}(w) \sqrt{\frac{\psi_{S_{w, \mathcal{F}}}(\mathrm{log} n_w/n_w)}{n_{w} \phi\left(h_{n_w}\right)}}\right). \end{matrix}
(A.7)

(A.5)-(A.7) 式的证明类似于文献 [定理 2] 的证明, 与之不同的地方是这里考虑响应变量 MAR 情形下的 \hat{\mu}_{w}(x), 其证明是平凡的, 此处省略证明细节.

注 A.1 不等式 (A.1) 通过控制 S_{w, \mathcal{F}} 的熵来处理拓扑问题. 对于一个半径不太大的小球, 需满足 \psi_{S_{w, \mathcal{F}}}(\varepsilon_w) 大小适中, 且满足在 n\rightarrow\infty 时, \psi_{S_{w, \mathcal{F}}}(\varepsilon_w)/(n_w\phi(h_w))\rightarrow0. 在一般情况下, \mathrm{log}^2n_w=O\left(n_w\phi(h_w)\right). 此外, 需要说明的是, (2.3) 式中的 H_{n_w,k_w}(x)h_ {n_w}(x) 的子序列.

在介绍引理 A.2 前, 引入一个新的记号. 对于任意 x\in S_{w,\mathcal{F}}, {n_w}\in N_{\varepsilon_w}(S_{w,\mathcal{F}}), 记

q_{{n_w},x}(t)=\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)B_{i,w}G_w(t,(x,A_{i,w}))}{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)G_w(t,(x,A_{i,w}))}, w=0,1.

引理 A.2 假设 \left(D_{n_w}(x)\right)_{{n_w}\in N_{\varepsilon_w}(S_{w,\mathcal{F}})} 是定义在 S_ {w,\mathcal{F}} 上的实随机变量, q:S_{w,\mathcal{F}}\rightarrow \mathbb{R} 是非随机函数, 且满足 \sup\limits_{x\in S_{w,\mathcal{F}}}\mid q(x)\mid <\infty. 如果存在非负递减序列 (u_{n_w})_{{n_w}\in N_{\varepsilon_w}(S_{w,\mathcal{F}})} 满足 \lim\limits_ {{n_w}\rightarrow \infty }u_{n_w}=0, 那么对任意递增序列 \beta_{n_w}\in(0,1), 都有 \beta_{n_w}-1=O\bigl(p_{\theta}(w)u_{n_w}\bigr). 进一步, 如果存在两列实随机变量D_{n_w}^{-}(\beta_{n_w},x)_{{n_w}\in N_{\varepsilon_w}(S_{w,\mathcal{F}})}D_{n_w}^{+}(\beta_{n_w},x)_{{n_w}\in N_{\varepsilon_w}(S_{w,\mathcal{F}})} 满足条件

\begin{align*} &({\bf H1})\quad \forall {n_w}\in N_{\varepsilon_w}(S_{w,\mathcal{F}}), \forall x\in S_{w, \mathcal{F}},D_{n_w}^{-}(\beta_{n_w},x)\leq D_{n_w}^{+}(\beta_{n_w},x);\\ &({\bf H2})\quad \lim_{{n_w \to \infty}}\mathbb{I}_{\{D_{n_w}^{-}(\beta_{n_w},x)\leq D_{n_w}(x)\leq D_{n_w}^{+}(\beta_{n_w},x), \forall x\in S_{w, \mathcal{F}}\}} = 1;\\ &({\bf H3})\quad \sup\limits_{x\in S_{w, \mathcal{F}}}\left|\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)G_w\left(D_{n_w}^{-}(\beta_{n_w},x),(x,A_{i,w})\right)}{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)G_w\left(D_{n_w}^{+}(\beta_{n_w},x),(x,A_{i,w})\right)}-\beta_{n_w}\right|=O_{\rm a.co}\left(p_{\theta}(w)u_{n_w}\right);\\ &({\bf H4})\quad \sup\limits_{x\in S_{w, \mathcal{F}}}\bigl| q_{{n_w},x}\left(D_{n_w}^{-}(\beta_{n_w},x)\right)-q_w(x)\bigr|=O_{\rm a.co}\left(p_{\theta}(w)u_{n_w}\right);\\ &({\bf H5})\quad \sup\limits_{x\in S_{w, \mathcal{F}}}\bigl| q_{{n_w},x}\left(D_{n_w}^{+}(\beta_{n_w},x)\right)-q_w(x)\bigr|=O_{\rm a.co}\left(p_{\theta}(w)u_{n_w}\right), \end{align*}

那么

\sup\limits_{x\in S_{w, \mathcal{F}}}\big| q_{{n_w},x}(D_{n_w}(x))-q_w(x)\big|=O_{\rm a.co}\left(p_{\theta}(w)u_{n_w}\right).
(A.8)

类似于文献 [引理 3] 的证明, 与之不同的地方是这里的收敛速度为 p_{\theta}(w)u_{n_w}, 此处省略证明细节.

引理 A.1 说明将 k_w 的取值设置为 n_{w} \phi\left(h_{n_w}\right), 响应变量 MAR 情形下 CATE 的 k 近邻核估计量能够达到与传统核估计量相同的的最优收敛速度. 引理 A.2 说明可以直接从核估计量的几乎完全收敛性出发, 推导出 k 近邻核估计量在同样条件下也具有类似的几乎完全收敛性. 基于上述两个引理, 下面给出响应变量 MAR 情形下 CATE 的 k 近邻核估计量几乎完全收敛性的证明过程.

定理 3.1 的证明 本定理的证明需利用引理 A.2 的结论, 故需要对引理 A.2 的五个条件进行验证. 沿用引理 A.2 的记号, 取 A_{i,w}=X_i, B_{i,w}=Y_i(W_i)\delta_i(W_i), G_w(t,(x,A_{i,w}))=K(t^{-1}d(x,A_{i,w})), D_{n_w}(x)=H_{n_w,k_w}(x); 又取 \beta_{n_w}\in (0,1) 是单调递增数列, 满足 \beta_{n_w}-1=O\left(p_{\theta}(w)u_{n_w}\right), 其中

u_{n_w}=\left(\bigg(\phi^{-1}\left(\frac{k_w}{n_w}\right)\bigg)^{b_w}+\sqrt{\frac{\psi_{S_{w,\mathcal{F}}}(\log n_w/n_w)}{k_w}}\right);

再取 D_{n_w}^{-}(\beta_{n_w},x)=H_{n_w,k_w}^{-}(x), D_{n_w}^{+}(\beta_{n_w},x)=H_{n_w,k_w}^{+}(x), h_{n_w}=\phi^{-1}\left(\frac{k_w}{n_w}\right), \varphi_x\left(D_{n_w}^{-}(\beta_{n_w},x)\right)=\frac{\sqrt{\beta_{n_w}}k_w}{n_w}, \varphi_x\left(D_{n_w}^{+}(\beta_{n_w},x)\right)=\frac{k_w}{\sqrt{\beta_{n_w}}n_w}. 下面验证引理 A.2 中的条件 (H1)-(H5).

(i) 条件 (H1) 是显然成立的.

(ii) 验证 (H2). 设 x_1,x_2,\cdots,x_{N_{\varepsilon_w}(S_{w,\mathcal{F}})}S_{w, \mathcal{F}}\varepsilon_w 网, 其中 \varepsilon_w=\log n_w/n_w. 对于任意 \zeta>0, 有

\begin{align*} &\hspace{0.4cm}\mathbb{P}\left(\left| \mathbb{I}_{\{D_{n_w}^{-}(\beta_{n_w},x)\leq D_{n_w}(x)\leq D_{n_w}^{+}(\beta_{n_w},x),\forall x\in S_{w,\mathcal{F}}\}}-1\right|>\zeta\right) \nonumber \\ &\leq \mathbb{P}\left(\inf\limits_{x\in S_{w,\mathcal{F}}}\left(D_{n_w}(x)-D_{n_w}^{-}(\beta_{n_w},x)\right)<0\right)\nonumber +\mathbb{P}\left(\sup\limits_{x\in S_{w,\mathcal{F}}}\left(D_{n_w}(x)-D_{n_w}^{+}(\beta_{n_w},x)\right)>0\right)\nonumber\\ &\leq N_{\varepsilon_w}(S_{w,\mathcal{F}})\max\limits_{1\leq j\leq N_{\varepsilon_w}(S_{w,\mathcal{F}})}\mathbb{P}\left(\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)\mathbb{I}_{B(x_j,H_{n_w,k_w}^{-}(x_j)+2\varepsilon_w)(x_i)}\geq k_w\right)\nonumber\\ &\hspace{0.5cm}+N_{\varepsilon_w}(S_{w,\mathcal{F}})\max\limits_{1\leq j\leq N_{\varepsilon_w}(S_{w,\mathcal{F}})}\mathbb{P}\left(\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)\mathbb{I}_{B\left(x_j,H_{n_w,k_w}^{+}(x_j)-2\varepsilon_w\right)(x_i)}\leq k_w\right), \end{align*}

再利用文献 [28] 中的 Chernoff 不等式, 可得

\begin{align*} &\hspace{0.4cm}\mathbb{P}\bigg(\Bigl| \mathbb{I}_{\{D_{n_w}^{-}(\beta_{n_w},x)<D_{n_w}(x)<D_{n_w}^{+}(\beta_{n_w},x),\forall x\in S_{w,\mathcal{F}}\}}-1\Bigr|>\zeta\bigg) \nonumber \\ &\leq N_{\varepsilon_w}(S_{w,\mathcal{F}})\left({\rm e}^{-C{k_w}(1-\beta_{n_w})^2}+{\rm e}^{-C{k_w}\beta_{n_w}}\right)\\ &\leq N_{\varepsilon_w}(S_{w,\mathcal{F}})\left({\rm e}^{-\log[C\frac{n_w}{k_w}\phi(h_{n_w}^{-}+2\varepsilon_w)\exp(1-\beta_{n_w})]}\right)^{-k_w}\nonumber\\ &\hspace{0.4cm}+N_{\varepsilon_w}(S_{w,\mathcal{F}})\left({\rm e}^{\frac{\beta_{n_w}}{2}\left(1-C\frac{n_w}{k_w}\phi(h_{n_w}^{+}-2\varepsilon_w)\right)^2}\right)^{-k_w}. \end{align*}

又因为当 \varepsilon_w\rightarrow0 时, \psi_ {S_{w,\mathcal{F}}}(\varepsilon_w)/k_w\rightarrow0, 且注意到 \psi_ {S_{w,\mathcal{F}}}(\varepsilon_w)=\mathrm{log}(N_{\varepsilon_w}(S_{w,\mathcal{F}})), 易有

\mathbb{P}\bigg(\Bigl|\mathbb{I}_{\{D_{n_w}^{-}(\beta_{n_w},x)\leq D_{n_w}(x)\leq D_{n_w}^{+}(\beta_{n_w},x),\forall x\in S_{w,\mathcal{F}}\}}-1\Bigr|>\zeta\bigg)\leq {\rm e}^{(1-\lambda)\psi_{S_{w,\mathcal{F}}}(\varepsilon_w)}.

最后利用假设 3.6, 知条件 (H2) 成立.

(iii) 验证 (H3). 记

\begin{align*} &l_w^{*}\left(x,D_{n_w}(\beta_{n_w},x)\right)=\mathbb{E}\Bigl[K\left(H_{n_w,k_w}^{-1}(x)d(x,X_1)\right)\Bigr], x\in S_{w,\mathcal{F}}, w=0,1;\\ &m_w^{*}(x,D_{n_w}(\beta_{n_w},x))=\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)K\left(H_{n_w,k_w}^{-1}(x)d(x,X_i)\right)}{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)K\left(H_{n_w,k_w}^{-1}(x)d(x,X_i)\right)\delta_i(W_i)}, x\in S_{w,\mathcal{F}}, w=0,1. \end{align*}

经过简单的计算, 有

\begin{align*} \left|\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)G\left(D_{n_w}^{-}(\beta_{n_w},x),(x,A_i)\right)}{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)G\left(D_{n_w}^{+}(\beta_{n_w},x),(x,A_i)\right)}-\beta_{n_w}\right| \leq\mid m_1\mid \mid l_1\mid \left(\mid l_2\mid+\mid l_3-m_1^{-1}\mid\right), \end{align*}

其中

\begin{align*} &m_1=\frac{m_w^{*}\left(x,D_{n_w}^{-}(\beta_{n_w},x)\right)}{m_w^{*}\left(x,D_{n_w}^{+}(\beta_{n_w},x)\right)}, l_1=\frac{l_w^{*}\left(x,D_{n_w}^{-}(\beta_{n_w},x)\right)}{l_w^{*}\left(x,D_{n_w}^{+}(\beta_{n_w},x)\right)},\nonumber\\ &l_2=\frac{\hat{\mu}_{1n_w}{\left(x,D_{n_w}^{-}(\beta_{n_w},x)\right)}}{{\hat{\mu}_{1n_w}}\left(x,D_{n_w}^{+}(\beta_{n_w},x)\right)}-1, l_3=\frac{l_w^{*}\left(x,D_{n_w}^{+}(\beta_{n_w},x)\right)}{l_w^{*}\left(x,D_{n_w}^{-}(\beta_{n_w},x)\right)}\beta_{n_w}. \end{align*}

n\rightarrow\infty 时, 有 \mid m_1\mid\rightarrow1. 根据假设 3.5 知, 存在常数 C, 使得 \sup\limits_{x\in S_{w,\mathcal{F}}}\mid l_1\mid<C. 此外, 根据 (A.6) 式易得

\begin{align*} \sup\limits_{x\in S_{w,\mathcal{F}}}\mid l_2\mid &\leq\frac{\sup\limits_{x\in S_{w,\mathcal{F}}}\Bigl| \hat{\mu}_{1n_w}\left(x,H_{n_w,k_w}^{-}(x)\right)-p_{\theta}(w)\Bigr|+ \sup\limits_{x\in S_{w,\mathcal{F}}}\Bigl|\hat{\mu}_{1n_w}\left(x,H_{n_w,k_w}^{+}(x)\right)-p_{\theta}(w)\Bigr|}{\inf\limits_{x\in S_{w,\mathcal{F}}}\Bigl|\hat{\mu}_{1n_w}\left(x,H_{n_w,k_w}^{+}(x)\right)\Bigr|} \nonumber \\ &=O_{\rm a.co.}\bigg(p_{\theta}(w)\sqrt{\frac{\psi_{S_{w,\mathcal{F}}}(\varepsilon_w)}{k_w}}\bigg). \end{align*}

根据文献 [引理 1] 和假设3.1-3.2 得, 存在 \gamma>0, 使得

\begin{align*} l_w^{*}\left(x,H_{n_w,k_w}(x)\right)=\gamma\varphi_x\bigl(H_{n_w,k_w}(x)\bigr)+O\bigl(\phi(h_{n_w})h_{n_w}^{b_w}\bigr), x\in S_{w,\mathcal{F}}. \end{align*}

显然 \varphi_x\left(D_{n_w}^{-}(\beta_{n_w},x)\right)/\varphi_x\left(D_{n_w}^{+}(\beta_{n_w},x)\right)=\beta_{n_w}.

\begin{align*} \sup\limits_{x\in S_{w,\mathcal{F}}}\left| l_3-m_1^{-1}\right| =O\Biggl(p_{\theta}(w)\biggl(\phi^{-1}\left(\frac{k_w}{n_w}\right)\biggr)^{b_w}\Biggr). \end{align*}

于是条件 (H3) 成立.

(iv) 验证 (H4) 和 (H5). 取 q_{n_w,x}(D_{n_w}(x))=\tilde{\mu}_{w,kNN}(x), q_w(x)=\mu_w(x)p_{\theta}(w), 其中

\begin{align*} \tilde{\mu}_{w,kNN}(x)=\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)K\bigl(H_{n_w,k_w}(x)^{-1}d(x,X_i)\bigr)\delta_i(W_i)Y_i(W_i)}{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)K\bigl(H_{n_w,k_w}(x)^{-1}d(x,X_i)\bigr)}. \end{align*}

由于 H_{n_w,k_w}^{-}(x), H_{n_w,k_w}^{+}(x)H_{n_w,k_w}(x) 同时满足引理 A.1 的条件, 因此条件 (H4) 和条件 (H5) 成立.

由于条件 (H1)-(H5) 成立, 应用引理 A.2 可得

\begin{matrix} \sup\limits_{x\in S_{w,\mathcal{F}}}\bigl|\tilde{\mu}_{w,kNN}(x)-\mu_w(x)p_{\theta}(w)\bigr|=O_{\rm a.co}\bigl(p_{\theta}(w)u_{n_w}\bigr), w=0,1. \end{matrix}
(A.9)

根据 (2.3) 和 (3.1) 式知

\begin{align*} \tilde{\mu}_{w,kNN}(x)=\hat{\mu}_{w,kNN}(x)\tilde{\mu}_{1n_w}(x)\frac{n_w\mathbb{E}\Bigl[K\bigl(H_{n_w,k_w}(x)^{-1}d(x,X_i)\bigr)\Bigr]}{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)K\bigl(H_{n_w,k_w}(x)^{-1}d(x,X_i)\bigr)}, w=0,1. \end{align*}

由 (A.6) 式得

\begin{matrix} \sup\limits_{x\in S_{w,\mathcal{F}}}\bigl|\hat{\mu}_{w,kNN}(x)-\mu_w(x)\bigr|=O_{\rm a.co}(u_{n_w}), w=0,1. \end{matrix}
(A.10)

注意到

\begin{align*} \sup\limits_{x\in S_{w,\mathcal{F}}}\bigl|\hat{\tau}_{kNN}(x)-\tau(x)\bigr|\leq\sup\limits_{x\in S_{w,\mathcal{F}}}\bigl|\hat{\mu}_{1,kNN}(x)-\mu_1(x)\bigr|+\sup\limits_{x\in S_{w,\mathcal{F}}}\bigl|\hat{\mu}_{0,kNN}(x)-\mu_0(x)\bigr|, \end{align*}

所以

\begin{align*}\sup\limits_{x\in S_{w,\mathcal{F}}}\bigl|\hat{\tau}_{kNN}(x)-\tau(x)\bigr|=O_{\rm a.co}\left(u_{n_1}+u_{n_0}\right).\end{align*}

定理 3.1 得证.

在证明定理 3.2 之前, 先对 \tilde{\tau}(x) 进行分析. 由 (2.3), (3.1) 和 (3.2) 式可以得到

\begin{align*} \hat{\mu}_{w,kNN}(x)-\mu_w(x)-B_{n_w}(x)&=\frac{\bigl(\tilde{\mu}_{2n_w}(x)-\mathbb{E}[\tilde{\mu}_{2n_w}(x)]\bigr)-\mu_w(x)\bigl(\tilde{\mu}_{1n_w}(x)-\mathbb{E}[\tilde{\mu}_{1n_w}(x)]\bigr)}{\tilde{\mu}_{1n_w}(x)}\nonumber\\ &\hspace{0.5cm}+\frac{B_{n_w}(x)\bigl(\mathbb{E}[\tilde{\mu}_{1n_w}(x)]-\tilde{\mu}_{1n_w}(x)\bigr)}{\tilde{\mu}_{1n_w}(x)}, w=0,1. \end{align*}

计算可得

\begin{matrix} &~~~~\hat{\mu}_{w,kNN}(x)-\mu_w(x)-B_{n_w}(x)\nonumber\\ &=\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)Z_{n_w,i}(x)}{n_w\mathbb{E}[K_{i,w}(x)]\tilde{\mu}_{1n_w}(x)}+\frac{B_{n_w}(x)\bigl(\mathbb{E}[\tilde{\mu}_{1n_w}(x)]-\tilde{\mu}_{1n_w}(x)\bigr)}{\tilde{\mu}_{1n_w}(x)},w=0,1, \end{matrix}
(A.11)

其中

\begin{matrix} Z_{n_w,i}(x)=\delta_i(w)K_{i,w}(x)\bigl(Y_i(w)-\mu_w(x)\bigr)-\mathbb{E}\Bigl[\delta_i(w)K_{i,w}(x)\bigl(Y_i(w)-\mu_w(x)\bigr)\Bigr]. \end{matrix}
(A.12)

由于 B_{n_w}(x)=o(1), w=0,1, 故

\begin{align*} \tilde{\tau}(x)=\sqrt{n}\left(\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)Z_{n_1,i}(x)\bigl(1+o_p(1)\bigr)}{n_1\mathbb{E}[K_{1,1}(x)]\tilde{\mu}_{1n_1}(x)}-\frac{\sum\limits_{i=1}^{n}\mathbb{I}(W_i=w)Z_{n_0,i}(x)Q_{n_0}(x)\bigl(1+o_p(1)\bigr)}{n_0\mathbb{E}[K_{n,0}(x)]\tilde{\mu}_{1n_0}(x)}\right). \end{align*}

假设处理组个体 m (m=1,2,\cdots,n_1) 的协变量为 X_m, 潜在结果为 Y_m(1); 对照组个体 j (j=n_1+1,n_1+2,\cdots,n) 的协变量为 X_j, 潜在结果为 Y_j(0).B(x,H_{n_1,k_1})B(x,H_{n_0,k_0}) 两个小球来估计 Y(1)Y(0) 的带宽, 所以 Y(1)Y(0) 相关的点需要满足 X_m 在小球 B(x,H_{n_1,k_1}) 中, X_j 在小球 B(x,H_{n_0,k_0}) 中. 不妨假设 Y(1)Y(0) 相关的 k_1X_m 下标为 m=1,2,\cdots,k_1, Y(1)Y(0) 相关的 k_0X_j 下标为 j=n-k_0+1,\cdots,n, Y(1)Y(0) 不相关的 n_1-k_1X_m 的下标为 m=k_1+1,k_1+2,\cdots,n_1, n_0-k_0X_j 下标为 j=n_1+1,\cdots,n-k_0. 这样, 可将 \tilde{\tau}(x) 分成独立的三个部分, 即

\begin{matrix} \tilde{\tau}(x)=\sqrt{n}\left(J_0+J_1+J_{01}\right), \end{matrix}
(A.13)

其中

\begin{align*} &J_0=\frac{\sum\limits_{j=n_1+1}^{n-k_0}Z_{n_0,j}(x)\bigl(1+o_p(1)\bigr)}{{n_0}\mathbb{E}[K_{n,0}(x)]\tilde{\mu}_{1n_0}(x)}, J_1=\frac{\sum\limits_{m=k_1+1}^{n_1}Z_{n_1,m}(x)\bigl(1+o_p(1)\bigr)}{{n_1}\mathbb{E}[K_{1,1}(x)]\tilde{\mu}_{1n_1}(x)},\nonumber\\ &J_{01}=\frac{\sum\limits_{m=1}^{k_1}Z_{n_1,m}(x)\bigl(1+o_p(1)\bigr)}{{n_1}\mathbb{E}[K_{1,1}(x)]\tilde{\mu}_{1n_1}(x)}+\frac{\sum\limits_{j=n-k_0+1}^{n}Z_{n_0,j}(x)\bigl(1+o_p(1)\bigr)}{{n_0}\mathbb{E}[K_{n,0}(x)]\tilde{\mu}_{1n_0}(x)}. \end{align*}

至此, 我们已将 \tilde{\tau}(x) 分成独立不同分布的和的形式, 后文将利用林德伯格中心极限定理来证明定理 3.2. 显然, 当 n \to \infty 时, \tilde{\tau}(x) 的期望收敛到 0. 因此, 在证明中只需计算 \tilde{\tau}(x) 的渐近方差以及验证林德伯格条件即可, 详见文献 [30].

定理 3.2 的证明 先计算 \tilde{\tau}(x) 的方差. 由 (A.13) 式得 {\rm Var}[\tilde{\tau}(x)]=n(T_0+T_1+T_{01}), 其中

\begin{align*} &T_0={\rm Var}\left[\frac{\sum\limits_{j=n_1+1}^{n-k_0}Z_{n_0,j}(x)}{{n_0}\mathbb{E}[K_{n,0}(x)]\tilde{\mu}_{1n_0}(x)}\right], T_1={\rm Var}\left[\frac{\sum\limits_{m=k_1+1}^{n_1}Z_{n_1,m}(x)}{{n_1}\mathbb{E}[K_{1,1}(x)]\tilde{\mu}_{1n_1}(x)}\right],\\ &T_{01}={\rm Var}\left[\frac{\sum\limits_{m=1}^{k_1}Z_{n_1,m}(x)}{{n_1}\mathbb{E}[K_{1,1}(x)]\tilde{\mu}_{1n_1}(x)}+\frac{\sum\limits_{j=n-k_0+1}^{n}Z_{n_0,j}(x)}{{n_0}\mathbb{E}[K_{n,0}(x)]\tilde{\mu}_{1n_0}(x)}\right]. \end{align*}

下面分别计算 T_0, T_1T_{01}. 先计算 T_0. 易得

T_0=\frac{n_0-k_0}{{{n_0}^2}\tilde{\mu}^2_{1n_0}(x)\mathbb{E}^2[K_{n,0}(x)]}{\rm Var}[Z_{n_0,1}(x)].

v_{n_w}(x)=\mathbb{E}\Bigl[\delta_i(w)K_{i,w}(x)\bigl(Y_i(w)-\mu_w(x)\bigr)\Bigr], w=0,1, 根据 (A.12) 式知

\begin{align*} {\rm Var}[Z_{n_0,1}(x)]&=\mathbb{E}\Bigl[\delta_n^2(0)K^2_{n,0}(x)\bigl(Y_1(0)-\mu_0(x)\bigr)^2\Bigr]-v_{n_0}^2(x)\\ &=p_{\theta}^2(0)\mathbb{E}\Bigl[K_{n,0}^2(x)\bigl(Y_1(0)-\mu_0(X_1)\bigr)^2\Bigr]\\ &~~~+p_{\theta}^2(0)\mathbb{E}\Bigl[K_{n,0}^2(x)\bigl(\mu_0(X_1)-\mu_0(x)\bigr)^2\Bigr]-v_{n_0}^2(x)\\ &=I_{1}+I_{2}-v_{n_0}^2(x). \end{align*}

根据假设 3.3 得

\begin{align*} v_{n_0}(x)&\leq \sup\limits_{x_i\in B(x,h_{n_0})}\bigl |\mu_0(x_i)-\mu_0(x)\bigr|\mathbb{E}\left[\delta_i(0)K_{i,0}(x)\right]\leq Ch_{n_0}^{b_0}p_{\theta}(0)\mathbb{E}[K_{n,0}(x)], \end{align*}

故而有

\begin{matrix} v_{n_0}^2(x)\leq Ch_{n_0}^{2b_0}p_{\theta}^2(0)\mathbb{E}^2[K_{n,0}(x)]. \end{matrix}
(A.14)

对于 I_{1}, 根据假设 3.8,

\begin{align*} I_{1}&=p_{\theta}^2(0)\mathbb{E}\Bigl[K_{n,0}^2(x)\bigl(Y_1(0)-\mu_0(X_1)\bigr)^2\Bigr]\\ &=p_{\theta}^2(0)\mathbb{E}\bigl[K_{n,0}^2(x)g_2(0,X_1)\bigr]\\ &=p_{\theta}^2(0)\Bigl(g_2(0,x)\mathbb{E}[K_{n,0}^2(x)]+\mathbb{E}\bigl[K_{n,0}^2(x)\bigl(g_2(0,X_1)-g_2(0,x)\bigr)\bigr]\Bigr), \end{align*}

易知

\begin{align*} \Big|\mathbb{E}\bigl[K_{n,0}^2(x)\bigl(g_2(0,X_1)-g_2(0,x)\bigr)\bigr]\Big|\leq o(1)\mathbb{E}[K_{n,0}^2(x)], \end{align*}

故有

\begin{matrix} I_{1}=p_{\theta}^2(0)g_2(0,x)\bigl(1+o(1)\bigr)\mathbb{E}[K_{n,0}^2(x)]. \end{matrix}
(A.15)

对于 I_{2}, 根据假设 3.3

\begin{matrix} I_{2}&=p_{\theta}^2(0)\mathbb{E}\Bigl[K_{n,0}^2(x)\bigl(\mu_0(X_1)-\mu_0(x)\bigr)^2\Bigr]\nonumber\\ &\leq\sup\limits_{x_1\in B(x,h_{n_0})}\bigl(\mu_0(x_1)-\mu_0(x)\bigl)^2p_{\theta}^2(0)\mathbb{E}[K_{n,0}^2(x)]\nonumber\\ &\leq Ch_{n_0}^{2b_0}p_{\theta}^2(0)\mathbb{E}[K_{n,0}^2(x)]. \end{matrix}
(A.16)

由 (A.14), (A.15) 和 (A.16) 式可得

\begin{align*} {\rm Var}[Z_{n_0,1}(x)]&=I_{1}+I_{2}-v_{n_0}^2(x)\\ &=p_{\theta}^2(0)g_2(0,x)\bigl(1+o(1)\bigr)\mathbb{E}[K_{n,0}^2(x)]+Ch_{n_0}^{2b_0}p_{\theta}^2(0)\mathbb{E}[K_{n,0}^2(x)], \end{align*}

所以

\begin{matrix} T_0=\frac{(n_0-k_0)p_{\theta}^2(0)}{{n^2_0}\tilde{\mu}^2_{1n_0}(x)\mathbb{E}^2[K_{n,0}(x)]}g_2(0,x)\mathbb{E}[K_{n,0}^2(x)]+\frac{p_{\theta}^2(0)}{\tilde{\mu}^2_{1n_0}(x)}O_p(h_{n_0}^{2b_0}). \end{matrix}
(A.17)

再计算 T_{1}. 采取处理 T_{0} 时用到的相同技巧, 有

\begin{matrix} T_1=\frac{(n_1-k_1)p_{\theta}^2(1)}{{n^2_1}\tilde{\mu}^2_{1n_1}(x)\mathbb{E}^2[K_{1,1}(x)]}g_2(1,x)\mathbb{E}[K_{1,1}^2(x)]+\frac{p_{\theta}^2(1)}{\tilde{\mu}^2_{1n_1}(x)}O_p(h_{n_1}^{2b_1}). \end{matrix}
(A.18)

最后计算 T_{01}. 容易得到

\begin{align*} T_{01}&={\rm Var}\left[\frac{\sum\limits_{j=n-k_0+1}^{n}Z_{n_0,j}(x)}{{n_0}\tilde{\mu}_{1n_0}(x)\mathbb{E}[K_{n,0}(x)]}\right]+ {\rm Var}\left[\frac{\sum\limits_{m=1}^{k_1}Z_{n_1,m}(x)}{{n_1}\tilde{\mu}_{1n_1}(x)\mathbb{E}[K_{1,1}(x)]}\right]\\ &\hspace{0.4cm}+{\rm Cov}\left(\frac{\sum\limits_{m=1}^{k_1}Z_{n_1,m}(x)}{{n_1}\tilde{\mu}_{1n_1}(x)\mathbb{E}[K_{1,1}(x)]},\frac{\sum\limits_{j=n-k_0+1}^{n}Z_{n_0,j}(x)}{{n_0}\tilde{\mu}_{1n_0}(x)\mathbb{E}[K_{n,0}(x)]}\right) =\Gamma_{0}+\Gamma_{1}+\Gamma_{01}. \end{align*}

采用与处理 T_0, T_1 相同的技巧, 经过一系列常规的计算, 可以得到

\begin{align*} \Gamma_{0}=\frac{k_0p_{\theta}^2(0)}{{{n_0}^2}\tilde{\mu}^2_{1n_0}(x)\mathbb{E}^2[K_{n,0}(x)]}g_2(0,x)\mathbb{E}[K_{n,0}^2(x)]+\frac{p_{\theta}^2(0)}{\tilde{\mu}^2_{1n_0}(x)}O_p(h_{n_0}^{2b_0}),\\ \Gamma_{1}=\frac{k_1p_{\theta}^2(1)}{{{n_1}^2}\tilde{\mu}^2_{1n_1}(x)\mathbb{E}^2[K_{1,1}(x)]}g_2(1,x)\mathbb{E}[K_{1,1}^2(x)]+\frac{p_{\theta}^2(1)}{\tilde{\mu}^2_{1n_1}(x)}O_p(h_{n_1}^{2b_1}). \end{align*}

对于 \Gamma_{01}

\begin{align*} \Gamma_{01}=\frac{1}{{n_0}{n_1}\tilde{\mu}_{1n_0}(x)\tilde{\mu}_{1n_1}(x)\mathbb{E}[K_{1,1}(x)]\mathbb{E}[K_{n,0}(x)]}{\rm Cov}\left(\sum\limits_{m=1}^{k_1}Z_{n_1,m}(x),\sum\limits_{j=n-k_0+1}^{n}Z_{n_0,j}(x)\right). \end{align*}

根据 (A.12) 式, 有

\begin{align*} &\hspace{0.4cm}{\rm Cov}\left(\sum\limits_{m=1}^{k_1}Z_{n_1,m}(x),\sum\limits_{j=n-k_0+1}^{n}Z_{n_0,j}(x)\right)\\ &=\sum\limits_{m=1}^{k_1}\sum\limits_{j=n-k_0+1}^{n}p_{\theta}(1)p_{\theta}(0)Cov\Bigl(K_{m,1}(x)\bigl(Y_m(1)-\mu_1(x)\bigr),K_{j,0}(x)\bigl(Y_j(0)-\mu_0(x)\bigr)\Bigr)\\ &=\sum\limits_{m=1}^{k_1}\sum\limits_{j=n-k_0+1}^{n}p_{\theta}(1)p_{\theta}(0)\biggl(\mathbb{E}\Bigl[K_{m,1}(x)K_{j,0}(x)\bigl(Y_m(1)-\mu_1(x)\bigr)\bigl(Y_j(0)-\mu_0(x)\bigr)\Bigr]\\ &\hspace{0.5cm}-\mathbb{E}\Bigl[K_{m,1}(x)\bigl(Y_m(1)-\mu_1(x)\bigr)\Bigr]\mathbb{E}\Bigl[K_{j,0}(x)\bigl(Y_j(0)-\mu_0(x)\bigr)\Bigr]\biggr). \end{align*}

根据假设 3.3 可知

\begin{align*} \mathbb{E}\Bigl[K_{m,1}(x)\bigl(Y_m(1)-\mu_1(x)\bigr)\Bigr]\mathbb{E}\Bigl[K_{j,0}(x)\bigl(Y_j(0)-\mu_0(x)\bigr)\Bigr]\!=\!O_p(h_{n_1}^{b_1}h_{n_0}^{b_0})\mathbb{E}[K_{1,1}(x)]\mathbb{E}[K_{n,0}(x)]. \end{align*}

此外

\begin{align*} \mathbb{E}\Bigl[K_{m,1}(x)K_{j,0}(x)\bigl(Y_m(1)-\mu_1(x)\bigr)\bigl(Y_j(0)-\mu_0(x)\bigr)\Bigr]=V_{1}+V_{2}+V_{3}+V_{4}, \end{align*}

其中

\begin{align*} &V_{1}=\mathbb{E}\Bigl[K_{m,1}(x)K_{j,0}(x)\bigl(Y_m(1)-\mu_1(X_i)\bigr)\bigl(Y_j(0)-\mu_0(X_j)\bigr)\Bigr],\nonumber\\ &V_{2}=\mathbb{E}\Bigl[K_{m,1}(x)K_{j,0}(x)\bigl(Y_m(1)-\mu_1(X_i)\bigr)\bigl(\mu_0(X_j)-\mu_0(x)\bigr)\Bigr],\nonumber\\ &V_{3}=\mathbb{E}\Bigl[K_{m,1}(x)K_{j,0}(x)\bigl(Y_j(0)-\mu_0(X_j)\bigr)\bigl(\mu_1(X_i)-\mu_1(x)\bigr)\Bigr],\nonumber\\ &V_{4}=\mathbb{E}\Bigl[K_{m,1}(x)K_{j,0}(x)\bigl(\mu_1(X_m)-\mu_1(x)\bigr)\bigl(\mu_0(X_j)-\mu_0(x)\bigr)\Bigr]. \end{align*}

由假设 3.7 可得

\begin{align*} V_{1}&=\mathbb{E}\bigl[K_{m,1}(x)K_{j,0}(x)g(X_m,X_j)\bigr]\\ &=g(x,x)\mathbb{E}\bigl[K_{m,1}(x)K_{j,0}(x)\bigr]+\mathbb{E}\Bigl[K_{m,1}(x)K_{j,0}(x)\bigl(g(X_m,X_j)-g(x,x)\bigr)\Bigr]\\ &=g(x,x)\mathbb{E}\bigl[K_{m,1}(x)K_{j,0}(x)\bigr]+o(1)\mathbb{E}\bigl[K_{1,1}(x)K_{n,0}(x)\bigr], \end{align*}

V_{1}=\bigl(g(x,x)+o(1)\bigr)\mathbb{E}[K_{1,1}(x)K_{n,0}(x)]. 下面分析 V_{2}, V_{3}V_{4}, 根据假设 3.3 可得

\begin{align*} V_{2}=V_{3}=V_{4}=O_p(h_{n_1}^{b_1}h_{n_0}^{b_0})\mathbb{E}[K_{1,1}(x)]\mathbb{E}[K_{n,0}(x)]. \end{align*}

所以

\begin{align*} &~~~~\mathbb{E}\Bigl[K_{m,1}(x)K_{j,0}(x)\bigl(Y_m(1)-\mu_1(x)\bigr)\bigl(Y_j(0)-\mu_0(x)\bigr)\Bigr]\\ &=\bigl(g(x,x)+o(1)\bigr)\mathbb{E}[K_{1,1}(x)K_{n,0}(x)]+O_p(h_{n_1}^{b_1}h_{n_0}^{b_0})\mathbb{E}[K_{1,1}(x)]\mathbb{E}[K_{n,0}(x)]. \end{align*}

再根据假设 3.9 得

\begin{align*} \Gamma_{01}&=\frac{k_0k_1p_{\theta}(0)p_{\theta}(1)\mathbb{E}[K_{1,1}(x)K_{n,0}(x)]}{{n_0}{n_1}\tilde{\mu}_{1n_0}(x)\tilde{\mu}_{1n_1}(x)\mathbb{E}[K_{1,1}(x)]\mathbb{E}[K_{n,0}(x)]}\rho(x)\sqrt{g_1(x)g_0(x)}\\ &\hspace{0.5cm}+\frac{p_{\theta}(0)p_{\theta}(1)}{\tilde{\mu}_{1n_0}(x)\tilde{\mu}_{1n_1}(x)}O_p(h_{n_1}^{b_1}h_{n_0}^{b_0}). \end{align*}

易得

\begin{align*} T_0+\Gamma_{0}=\frac{p_{\theta}^2(0)g_0(x)\mathbb{E}[K_{n,0}^2(x)]}{n_0\tilde{\mu}^2_{1n_0}(x)\mathbb{E}^2[K_{n,0}(x)]}+\frac{p_{\theta}^2(0)}{\tilde{\mu}^2_{1n_0}(x)}O_p(h_{n_0}^{2b_0}),\\ T_1+\Gamma_{1}=\frac{p_{\theta}^2(1)g_1(x)\mathbb{E}[K_{1,1}^2(x)]}{n_1\tilde{\mu}^2_{1n_1}(x)\mathbb{E}^2[K_{1,1}(x)]}+\frac{p_{\theta}^2(1)}{\tilde{\mu}^2_{1n_1}(x)}O_p(h_{n_1}^{2b_1}). \end{align*}

利用上式可知

\begin{matrix} \frac{\Gamma_{01}}{T_0+\Gamma_{0}}=\frac{k_0k_1p_{\theta}(1)\tilde{\mu}_{1n_0}(x)\mathbb{E}[K_{1,1}(x)K_{n,0}(x)]\mathbb{E}[K_{1,0}(x)]}{n_1p_{\theta}(0)\tilde{\mu}_{1n_1}(x)\mathbb{E}[K_{n,0}^2(x)]\mathbb{E}[K_{1,1}(x)]}\rho(x)\sqrt{\frac{g_1(x)}{g_0(x)}}, \end{matrix}
(A.19)
\begin{matrix} \frac{\Gamma_{01}}{T_1+\Gamma_{1}}=\frac{k_0k_1p_{\theta}(0)\tilde{\mu}_{1n_1}(x)\mathbb{E}[K_{1,1}(x)K_{n,0}(x)]\mathbb{E}[K_{1,1}(x)]}{n_0p_{\theta}(1)\tilde{\mu}_{1n_0}(x)\mathbb{E}[K_{1,1}^2(x)]\mathbb{E}[K_{n,0}(x)]}\rho(x)\sqrt{\frac{g_0(x)}{g_1(x)}}. \end{matrix}
(A.20)

n\rightarrow\infty 时, 易知 k_1/n_1\rightarrow0, k_0/n_0\rightarrow0, h_{n_1}\rightarrow0h_{n_0}\rightarrow0. 此外, 根据 (A.6) 式可知, 当 n\rightarrow \infty 时, \tilde{\mu}_{1n_w}(x)\rightarrow p_{\theta}(w), w=0,1. 综合 (A.17)-(A.20) 式, 可得

\begin{align*} \lim_{n \to \infty}{\rm Var}[\tilde{\tau}(x)]&=n\bigg(\frac{g_0(x)\mathbb{E}[K_{n,0}^2(x)]}{n_0\mathbb{E}^2[K_{n,0}(x)]} +O_p\left(h_{n_0}^{2b_0}\right)\bigl(1+o(1)\bigr)\bigg)\\ &~~~+\bigg(\frac{g_1(x)\mathbb{E}[K_{1,1}^2(x)]}{n_1\mathbb{E}^2[K_{1,1}(x)]}+O_p\left(h_{n_1}^{2b_1}\right)\bigg)\\ &=\frac{g_0(x)\mathbb{E}[K_{n,0}^2(x)]}{(1-{p_\pi})\mathbb{E}^2[K_{n,0}(x)]}+\frac{g_1(x)\mathbb{E}[K_{1,1}^2(x)]}{{p_\pi}\mathbb{E}^2[K_{1,1}(x)]} =\sigma_0^2+\sigma_1^2. \end{align*}

下面验证林德伯格条件. 记 A_n^2:=\sigma_0^2+\sigma_1^2. 不妨设

\begin{align*} D_l=\! \begin{cases} \frac{\sqrt{n}Z_{n_0,l+n_1}(x)\bigl(1+o_p(1)\bigr)}{n_0\mathbb{E}[K_{n,0}(x)]\tilde{\mu}_{1n_0}(x)}, \hskip 5.7cm l=1,2,\cdots,n_0-k_0,\\[3mm] \frac{\sqrt{n}Z_{n_1,l-n_0+k+0+k_1}(x)\bigl(1+o_p(1)\bigr)}{n_1\mathbb{E}[K_{1,1}(x)]\tilde{\mu}_{1n_1}(x)}, \hskip 2.6cm l=n_0-k_0+1,\cdots,n-k_0-k_1,\\[3mm] \frac{\sqrt{n}\sum\limits_{i=1}^{k_1}Z_{n_1,m}(x)\bigl(1+o_p(1)\bigl)}{n_1\mathbb{E}[K_{1,1}(x)]\tilde{\mu}_{1n_1}(x)}+\frac{\sqrt{n}\sum\limits_{j=n-k_0+1}^{n}Z_{n_0,j}(x)\bigl(1+o_p(1)\bigr)}{n_0\mathbb{E}[K_{n,0}(x)]\tilde{\mu}_{1n_0}(x)},~~ l=n-k_0-k_1+1. \end{cases} \end{align*}

D_l 的密度函数为 f_l(x), 又因为 \mathbb{E}[D_l]=0, 则林德伯格条件为

\begin{matrix} \forall \eta>0,\lim\limits_{(n-k_0-k_1+1)\rightarrow\infty}\frac{1}{\eta^2A_n^2}\sum_{l=1}^{n-k_0-k_1+1}\int_{\mid x\mid>\eta A_n}x^2f_l(x){\rm d}x=0. \end{matrix}
(A.21)

为证明 (A.21) 式成立, 将其分为三个部分, 即

\begin{align*} \sum_{l=1}^{n-k_0-k_1+1}\int_{\mid x\mid>\eta A_n}x^2f_l(x){\rm d}x&=\sum_{l=1}^{n_0-k_0}\int_{\mid x\mid>\eta A_n}x^2f_l(x){\rm d}x +\sum_{l=n_0-k_0+1}^{n-k_0-k_1}\int_{\mid x\mid>\eta A_n}x^2f_l(x){\rm d}x\\ &\hspace{0.5cm}+\int_{\mid x\mid>\eta A_n}x^2f_{n-k_0-k_1+1}(x){\rm d}x\\ &=R_1+R_2+R_3. \end{align*}

根据 (A.17), (A.18) 式和 T_{01} 计算过程, 便可证明 (A.21) 式成立, 因此满足林德伯格条件. 由林德伯格中心极限定理可得定理 3.2 成立.

声明 : 明瑞星与曾华俊在该项工作中贡献相同.

参考文献

Rosenbaum P R, Rubin D B.

The central role of the propensity score in observational studies for causal effects

Biometrika, 1983, 70(1): 41-55

[本文引用: 1]

Rosenbaum P R, Rubin D B.

Constructing a control group using multivariate matched sampling models that incorporate the propensity score

Amer Statist, 1985, 39(1): 33-38

[本文引用: 1]

Caron A, Baio G, Manolopoulou I.

Estimating individual treatment effects using non-parametric regression models: A review

J Roy Statist Soc Ser A, 2022, 185(3): 1115-1149

[本文引用: 2]

Yao L Y, Chu Z X, Li S, et al.

A survey on causal inference

ACM Transactions on Knowledge Discovery from Data, 2021, 15(5): 1-46

[本文引用: 1]

Sarraju A, Ward A, Li J, et al.

Personalizing cholesterol treatment recommendations for primary cardiovascular disease prevention

Sci Rep, 2022, 12(1): Article 23

[本文引用: 1]

Sorace M.

The ties that unbind: Intergovernmental decision rules and the policy-opinion link

J Eur Public Policy, 2023, 30(8): 1609-1632

[本文引用: 1]

Zhang Z W, Liu W, Zhang B, et al.

Causal inference with missing exposure information: Methods and applications to an obstetric study

Stat Methods Med Res, 2016, 25(5): 2053-2066

PMID:24318273      [本文引用: 1]

Causal inference in observational studies is frequently challenged by the occurrence of missing data, in addition to confounding. Motivated by the Consortium on Safe Labor, a large observational study of obstetric labor practice and birth outcomes, this article focuses on the problem of missing exposure information in a causal analysis of observational data. This problem can be approached from different angles (i.e. missing covariates and causal inference), and useful methods can be obtained by drawing upon the available techniques and insights in both areas. In this article, we describe and compare a collection of methods based on different modeling assumptions, under standard assumptions for missing data (i.e. missing-at-random and positivity) and for causal inference with complete data (i.e. no unmeasured confounding and another positivity assumption). These methods involve three models: one for treatment assignment, one for the dependence of outcome on treatment and covariates, and one for the missing data mechanism. In general, consistent estimation of causal quantities requires correct specification of at least two of the three models, although there may be some flexibility as to which two models need to be correct. Such flexibility is afforded by doubly robust estimators adapted from the missing covariates literature and the literature on causal inference with complete data, and by a newly developed triply robust estimator that is consistent if any two of the three models are correct. The methods are applied to the Consortium on Safe Labor data and compared in a simulation study mimicking the Consortium on Safe Labor.© The Author(s) 2013.

Abrevaya J, Hsu Y C, Lieli R P.

Estimating conditional average treatment effects

J Bus Econom Statist, 2015, 33(4): 485-505

[本文引用: 1]

Li L, Zhou N W, Zhu L X.

Outcome regression-based estimation of conditional average treatment effect

Ann Inst Statist Math, 2022, 74(5): 987-1041

[本文引用: 1]

Kudraszow N L, Vieu P.

Uniform consistency of kNN regressors for functional variables

Statist Probab Lett, 2013, 83(8): 1863-1870

[本文引用: 4]

Ling N X, Meng S Y, Vieu P.

Uniform consistency rate of kNN regression estimation for functional time series data

J Nonparametr Stat, 2019, 31(2): 451-468

[本文引用: 1]

Friedman J H.

Greedy function approximation: A gradient boosting machine

Ann Statist, 2001, 29: 1189-1232

[本文引用: 1]

Chen T Q, Guestrin C. Xgboost: A scalable tree boosting system//Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2016, 785-794

[本文引用: 1]

Breiman L.

Random forests

Mach Learn, 2001, 45(1): 5-32

[本文引用: 1]

Wager S, Athey S.

Estimation and inference of heterogeneous treatment effects using random forests

J Amer Statist Assoc, 2018, 113(523): 1228-1242

[本文引用: 1]

Athey S, Wager S.

Estimating treatment effects with causal forests: An application

Observ Stud, 2019, 5(2): 37-51

[本文引用: 1]

Little R J, D'agostino R, Cohen M L, et al.

The prevention and treatment of missing data in clinical trials

New Engl J Med, 2012, 367(14): 1355-1360

[本文引用: 1]

Ocampo A, Schmidli H, Quarg P, et al.

Identifying treatment effects using trimmed means when data are missing not at random

Pharm Stat, 2021, 20(6): 1265-1277

DOI:10.1002/pst.2147      PMID:34169641      [本文引用: 1]

Patients often discontinue from a clinical trial because their health condition is not improving or they cannot tolerate the assigned treatment. Consequently, the observed clinical outcomes in the trial are likely better on average than if every patient had completed the trial. If these differences between trial completers and non-completers cannot be explained by the observed data, then the study outcomes are missing not at random (MNAR). One way to overcome this problem-the trimmed means approach for missing data due to study discontinuation-sets missing values as the worst observed outcome and then trims away a fraction of the distribution from each treatment arm before calculating differences in treatment efficacy (Permutt T, Li F. Trimmed means for symptom trials with dropouts. Pharm Stat. 2017;16(1):20-28). In this paper, we derive sufficient and necessary conditions for when this approach can identify the average population treatment effect. Simulation studies show the trimmed means approach's ability to effectively estimate treatment efficacy when data are MNAR and missingness due to study discontinuation is strongly associated with an unfavorable outcome, but trimmed means fail when data are missing at random. If the reasons for study discontinuation in a clinical trial are known, analysts can improve estimates with a combination of multiple imputation and the trimmed means approach when the assumptions of each hold. We compare the methodology to existing approaches using data from a clinical trial for chronic pain. An R package trim implements the method. When the assumptions are justifiable, using trimmed means can help identify treatment effects notwithstanding MNAR data.© 2021 John Wiley & Sons Ltd.

Kuzmanovic M, Hatt T, Feuerriegel S. Estimating conditional average treatment effects with missing treatment information//Proceedings of the 26th International Conference on Artificial Intelligence and Statistics. Valencia: PMLR, 2023, 746-766

[本文引用: 1]

Stensrud M J, Robins J M, Sarvet A, et al.

Conditional separable effects

J Amer Statist Assoc, 2023, 118(544): 2671-2683

[本文引用: 1]

Zhao A Q, Peng D.

To adjust or not to adjust? Estimating the average treatment effect in randomized experiments with missing covariates

J Amer Statist Assoc, 2024, 119(545): 450-460

[本文引用: 1]

苗旺, 刘春辰, 耿直.

因果推断的统计方法

中国科学: 数学, 2018, 48(12): 1753-1778

[本文引用: 1]

Miao W, Liu C C, Geng Z.

Statistical approaches for causal inference (in Chinese)

Sci Sin Math, 2018, 48(12): 1753-1778 (in Chinese).

[本文引用: 1]

Greenland S, Pearl J, Robins J M.

Confounding and collapsibility in causal inference

Statist Sci, 1999, 14(1): 29-46

[本文引用: 1]

Geng Z, Guo J H, Fung W K.

Criteria for confounders in epidemiological studies

J R Stat Soc Ser B Stat Methodol, 2002, 64(1): 3-15

[本文引用: 1]

Han P S.

Multiply robust estimation in regression analysis with missing data

J Amer Statist Assoc, 2014, 109(507): 1159-1173

[本文引用: 1]

Ferraty F, Laksaci A, Tadj A, et al.

Rate of uniform consistency for nonparametric estimates with functional variables

J Statist Plann Inference, 2010, 140(2): 335-352

[本文引用: 1]

Hill J L.

Bayesian nonparametric modeling for causal inference

J Comput Graph Statist, 2011, 20(1): 217-240

[本文引用: 1]

Burba F, Ferraty F, Vieu P.

k-Nearest neighbour method in functional nonparametric regression

J Nonparametr Stat, 2009, 21(4): 453-469

[本文引用: 1]

Ezzahrioui M, Ould-Saïd E.

Asymptotic normality of a nonparametric estimator of the conditional mode function for functional data

J Nonparametr Stat, 2008, 20(1): 3-18

黄绍航. 随机缺失下条件平均处理效果的 k 近邻核估计. 杭州: 浙江工商大学, 2022

[本文引用: 1]

Huang S H. K-nearest Neighbor Kernel Estimation of Conditional Average Treatment Effect Under Random Missing Data. Hangzhou: Zhejiang Gongshang University, 2022

[本文引用: 1]

/