第十二章：M-估计量与 Z-估计量（二）

在上一部分中，我们讨论了标准的 M-估计量和 Z-估计量的相合性与渐近正态性。在实际的统计建模中，我们往往还会遇到包含滋扰参数 (Nuisance Parameter) 或者需要通过迭代算法（如牛顿-拉夫逊法）来求解的情况。本部分将深入探讨带有估计参数的 Z-估计量以及单步估计量 (One-step Estimator)。

8. 包含估计滋扰参数的估计量 (Estimated Parameter)

假设样本 \(X_1, \dots, X_n \sim i.i.d.\ F_{(\theta, \eta)}\)，其中：

\(\theta\) 是我们感兴趣的目标参数 (Parameter of interest)。
\(\eta\) 是未知的滋扰参数 (Nuisance parameter)。

通常，我们需要在 Z-估计方程中代入 \(\eta\) 的一个估计量（记为 \(\hat{\eta}_n\)）：

\[ P_n \psi_{(\theta, \eta)} \Rightarrow P_n \psi_{(\theta, \hat{\eta}_n)} = \frac{1}{n} \sum_{i=1}^n \psi_{(\theta, \hat{\eta}_n)}(X_i) = 0 \]

这本质上是一个两步过程 (2-step procedure)：首先估计 \(\eta\)，然后将其作为“已知”参数代入方程求解 \(\theta\)。

8.1 动机示例：响应变量缺失的回归 (Missing Response Regression)

考虑回归模型：

\[ y_i = m_{\theta_0}(x_i) + \epsilon_i, \quad E(\epsilon_i | x_i) = 0 \]

假设响应变量 \(y_i\) 存在缺失情况。我们定义观测指示变量： \(R_i = 1\) 如果 \(y_i\) 被观测到；\(R_i = 0\) 如果 \(y_i\) 缺失。

随机缺失假设 (Missing at Random, MAR): 假设缺失机制满足：

\[ P(R_i = 1 | x_i, y_i) = P(R_i = 1 | x_i) = \omega_{\eta_0}(x_i) \]

其中 \(\omega_{\eta_0}(\cdot)\) 被称为缺失倾向得分函数 (Missing propensity function)。在这里我们假设它是一个参数化的二分类回归模型。

MAR 的含义：给定 \(X_i\)，缺失指示变量 \(R_i\) 和真实响应 \(y_i\) 是条件独立的。
这种缺失机制被称为可忽略的随机缺失 (Ignorable missing at random)，因为只要我们观测到 \(x_i\)，缺失的 \(y_i\) 对推断的干扰就是“可忽略的”。

如何估计目标参数 \(\theta\)？

方法 1：仅对观测到数据的样本 (即 \(R_i = 1\)) 进行最小二乘估计 (LSE)

目标函数为：

\[ LS_n(\theta) = \sum_{i=1}^n R_i (y_i - m_\theta(x_i))^2 \]

对其求导得到得分方程：

\[ \frac{\partial LS_n(\theta)}{\partial \theta} = -2 \sum_{i=1}^n R_i (y_i - m_\theta(x_i)) \frac{\partial m_\theta(x_i)}{\partial \theta} = 0 \quad \text{--- (8.13)} \]

在真实参数 \(\theta_0\) 处，利用 \(E(\epsilon_i | x_i) = 0\) 且在 MAR 假设下：

\[ E \left\{ R_i (y_i - m_{\theta_0}(x_i)) \frac{\partial m_{\theta_0}(x_i)}{\partial \theta} \right\} = E \left\{ \epsilon_i \frac{\partial m_{\theta_0}(x_i)}{\partial \theta} \omega_{\eta_0}(x_i) \right\} = 0 \]

因此，在满足一定正则条件的情况下，求解 (8.13) 得到的最小二乘估计量是相合的且渐近正态 (AN) 的。

方法 2：逆概率加权估计 (Inverse Probability Weighted Estimation, IPW)

通过使用真实的缺失倾向 \(\omega_{\eta_0}(x_i)\) 对 (8.13) 进行加权，我们得到 IPW 得分方程：

\[ \sum_{i=1}^n \frac{R_i (y_i - m_\theta(x_i)) \frac{\partial m_\theta(x_i)}{\partial \theta}}{\omega_{\eta_0}(x_i)} = 0 \quad \text{--- (8.14)} \]

取期望同样可以验证无偏性：

\[ E \left\{ \sum_{i=1}^n \frac{R_i (y_i - m_\theta(x_i)) \frac{\partial m_\theta(x_i)}{\partial \theta}}{\omega_{\eta_0}(x_i)} \right\} = E \left\{ \epsilon_i \frac{\partial m_{\theta_0}(x_i)}{\partial \theta} \right\} = 0 \]

问题： 忽略缺失值的估计量 (8.13) 和 IPW 估计量 (8.14)，哪一个更有效（即渐近方差更小）？

8.2 滋扰参数的二分类似然估计

在实际中，\(\eta_0\) 通常是未知的。我们可以通过最大化二分类似然函数来估计它：

\[ L_n(\eta) = \prod_{i=1}^n \omega_\eta^{R_i}(x_i) (1 - \omega_\eta(x_i))^{1 - R_i} \]

对数似然函数为：

\[ l_n(\eta) = \sum_{i=1}^n \{ R_i \log \omega_\eta(x_i) + (1 - R_i) \log(1 - \omega_\eta(x_i)) \} \]

求导并令其为 0，得到 \(\hat{\eta}\)：

\[ \frac{\partial l_n(\eta)}{\partial \eta} = \sum_{i=1}^n \left\{ \frac{R_i}{\omega_\eta(x_i)} - \frac{1 - R_i}{1 - \omega_\eta(x_i)} \right\} \frac{\partial \omega_\eta(x_i)}{\partial \eta} = 0 \Rightarrow \hat{\eta} \]

此时，(8.14) 变为带有估计滋扰参数的方程：

\[ \sum_{i=1}^n \frac{R_i (y_i - m_\theta(x_i)) \frac{\partial m_\theta(x_i)}{\partial \theta}}{\omega_{\hat{\eta}}(x_i)} = 0 \quad \text{--- (8.15)} \]

反直觉的结论： Chen, Leung and Qin (2008) 证明，基于估计参数 \(\hat{\eta}\) 的方程 (8.15) 求解出的 \(\theta\)，比基于真实已知参数 \(\eta_0\) 的方程 (8.14) 求解出的 \(\theta\) 渐近更有效！

8.3 缺失倾向的非参数估计

假设参数化的模型 \(P(R_i = 1 | x_i) = \omega_{\eta_0}(x_i)\) 过于严格。我们可以考虑非参数形式：

\[ P(R_i = 1 | x_i, y_i) = P(R_i = 1 | x_i) = \omega(x_i) \]

缺失倾向 \(\omega(\cdot)\) 可以通过核平滑 (Kernel smoothing) 方法（如 Nadaraya-Watson 估计量）进行估计：

\[ \hat{\omega}_h(x) = \frac{\sum_{i=1}^n K\left(\frac{x - x_i}{h}\right) R_i}{\sum_{i=1}^n K\left(\frac{x - x_i}{h}\right)} \]

其中 \(K\) 是核函数（对称的概率密度函数），\(h\) 是平滑带宽 (smoothing bandwidth)。要求当 \(n \rightarrow \infty\) 时，\(h \rightarrow 0\) 且 \(nh \rightarrow \infty\)。其期望本质上是样本点处权重的加权平均：

\[ E(\hat{\omega}_h(x)) = E \left( \sum_{i=1}^n E \left( \frac{K\left(\frac{x - x_i}{h}\right)}{\sum K\left(\frac{x - x_j}{h}\right)} R_i \Big| x_1, \dots, x_n \right) \right) = E \left( \sum_{i=1}^n \frac{K\left(\frac{x - x_i}{h}\right)}{\sum K\left(\frac{x - x_j}{h}\right)} \omega(x_i) \right) \]

可以证明，如果选择合适的带宽 \(h \approx O(n^{-\frac{1}{4+d}})\)（其中 \(d\) 为 \(x\) 的维度），那么：

\[ \sqrt{nh^d} (\hat{\omega}_h(x) - \omega(x)) \rightarrow N(\mu, v^2) \quad (x \in \mathbb{R}^d) \]

9. 包含估计参数的 Z-估计量渐近定理

定理 8.13 (两步 Z-估计量的渐近展开)

对于参数空间 \(\mathbb{R}^k\) 的开子集中的每一个 \(\theta\)，以及度量空间中的每一个 \(\eta\)，设 \(x \rightarrow \phi_{\theta, \eta}(x)\) 是一个 \(\mathbb{R}^k\)-值的可测函数，满足：函数族 \(\{\phi_{\theta, \eta} : \|\theta - \theta_0\| < \delta, d(\eta, \eta_0) < \delta\}\) 对于某个 \(\delta > 0\) 是 Donsker 类；且当 \((\theta, \eta) \rightarrow (\theta_0, \eta_0)\) 时，\(P\|\phi_{\theta, \eta} - \phi_{\theta_0, \eta_0}\|^2 \rightarrow 0\)。

假设 \(P\phi_{\theta_0, \eta_0} = 0\)，且映射 \(\theta \rightarrow P\phi_{\theta, \eta}\) 在 \(\theta_0\) 处可微（在 \(\eta_0\) 的邻域内关于 \(\eta\) 一致成立），其导数矩阵 \(V_{\theta_0, \eta}\) 非奇异，并且满足当 \(\eta \rightarrow \eta_0\) 时 \(V_{\theta_0, \eta} \rightarrow V_{\theta_0, \eta_0}\)。

此外，如果存在相合的估计量 \((\hat{\theta}, \hat{\eta}) \xrightarrow{p} (\theta_0, \eta_0)\) 且满足 \(\sqrt{n}\mathbb{P}_n\phi_{\hat{\theta}, \hat{\eta}} = o_p(1)\)，那么：

\[ \sqrt{n}(\hat{\theta}_n - \theta_0) = -V_{\theta_0, \eta_0}^{-1} \sqrt{n} P\psi_{\theta_0, \hat{\eta}_n} - V_{\theta_0, \eta_0}^{-1} \mathbb{G}_n\psi_{\theta_0, \eta_0} + o_P(1 + \sqrt{n}\|P\psi_{\theta_0, \hat{\eta}_n}\|) \]

关于定理 8.13 的注记：

等式右侧的第一个漂移项 (drift term) \(-V_{\theta_0, \eta_0}^{-1} \sqrt{n} P\psi_{\theta_0, \hat{\eta}_n}\) 是由于我们额外估计了滋扰参数 \(\eta\) 而产生的。
\(\eta\) 这里可以是无穷维的函数参数（不仅仅是标量或向量参数）。
在一定的平滑性条件下，这个多出来的漂移项可以通过 Delta 方法进行泰勒展开近似：

\[ P\psi_{\theta_0, \hat{\eta}_n} = P\psi_{\theta_0, \eta_0} + \frac{\partial P\psi_{\theta_0, \eta_0}}{\partial \eta} (\hat{\eta}_n - \eta_0) + o_p(\hat{\eta}_n - \eta_0) \]
正交性 (Orthogonality)：如果目标参数和滋扰参数满足（近似）正交条件，即交叉偏导数为 0（\(\frac{\partial P\psi_{\theta_0, \eta_0}}{\partial \eta} = 0\)），那么估计滋扰参数 \(\hat{\eta}_n\) 的误差将不会影响目标参数 \(\hat{\theta}_n\) 在一阶极限下的渐近分布（漂移项消失）。

9.1 示例：对称位置模型 (Symmetric Location)

假设 \(X_1, \dots, X_n \sim i.i.d.\ F\)，且分布 \(F\) 关于真实位置参数 \(\theta_0\) 对称。我们要估计 \(\theta_0\)，但尺度参数 (scale parameter) \(\sigma\) 是未知的。设 \(x \rightarrow \psi(x)\) 是一个反对称的奇函数。

我们考虑如下的 Z-估计方程（其中 \(\hat{\sigma}\) 是 \(\sigma\) 的某种相合估计）：

\[ \frac{1}{n} \sum_{i=1}^n \psi\left( \frac{X_i - \theta}{\hat{\sigma}} \right) = 0 \]

由于 \(F(\cdot)\) 关于 \(\theta_0\) 对称，且 \(\psi(\cdot)\) 是奇函数，在真实的参数下，对于任意的尺度估计 \(\hat{\sigma}\)，期望恒为 0：

\[ P\psi_{\theta_0, \hat{\sigma}} = \int \psi\left( \frac{x - \theta_0}{\hat{\sigma}} \right) dF(x) = 0, \quad \forall \hat{\sigma} \]

因此，根据定理 8.13，由于漂移项为 0，我们可以得到：

\[ \sqrt{n}(\hat{\theta}_n - \theta_0) = -V_{\theta_0, \eta_0}^{-1} \frac{1}{\sqrt{n}} (\mathbb{P}_n\psi_{\theta_0, \eta_0} - P\psi_{\theta_0, \eta_0}) + o_p(1) \]

这意味着该估计过程渐近等效于使用了真实的 \(\eta_0\)。在这个模型中，\(\hat{\sigma}\) 的估计误差没有在一阶渐近展开中表现出来。

10. 单步估计量 (One-step Estimator)

在很多复杂的模型中，Z-估计方程 \(\Psi_n(\theta) = 0\) 的精确根可能非常难以求解（例如存在多个局部极值、高度非线性等）。

我们可以通过在某个初始估计量 (initial estimator) \(\tilde{\theta}_n\) 的基础上进行一步优化来改进估计结果。这正是基于牛顿-拉夫逊算法 (Newton-Raphson Algorithm) 的思想。

考虑在 \(\tilde{\theta}_n\) 处进行一阶泰勒展开：

\[ \Psi_n(\tilde{\theta}_n) + \dot{\Psi}_n(\tilde{\theta}_n)(\theta - \tilde{\theta}_n) = 0 \quad \text{--- (8.16)} \]

我们可以得到对 \(\tilde{\theta}_n\) 的单步更新：

\[ \hat{\theta}_n = \tilde{\theta}_n - \dot{\Psi}_n^{-1}(\tilde{\theta}_n) \Psi_n(\tilde{\theta}_n) \quad \text{--- (8.17)} \]

通常我们可以选择一个容易计算但不是最有效 (inefficient) 的 \(\sqrt{n}\) 相合估计量作为初始值 \(\tilde{\theta}_n\)（例如矩估计 ME）。我们将看到，单步估计量 \(\hat{\theta}_n\) 能够完美达到原本难以求解的全 Z-估计量同样的渐近有效性。

10.1 正则条件

一致收敛条件： 对于任意常数 \(M\) 和给定的非奇异矩阵 \(\dot{\Psi}_0\)，满足：

\[ \sup_{\sqrt{n}\|\theta - \theta_0\| < M} \|\sqrt{n}(\Psi_n(\theta) - \Psi_n(\theta_0)) - \dot{\Psi}_0 \sqrt{n}(\theta - \theta_0)\| \xrightarrow{P} 0 \quad \text{--- (8.18)} \]

其中 \(\dot{\Psi}_0\) 是导数 \(\dot{\Psi}_n(\theta)\) 当 \(\theta \rightarrow \theta_0\) 时的极限矩阵。

在实践中，我们可以将 (8.17) 中的经验导数替换为任何相合的估计矩阵 \(\dot{\Psi}_{n,0}\)（满足 \(\dot{\Psi}_{n,0} \xrightarrow{p} \dot{\Psi}_0\)）：

\[ \hat{\theta}_n = \tilde{\theta}_n - \dot{\Psi}_{n,0}^{-1} \Psi_n(\tilde{\theta}_n) \quad \text{--- (8.19)} \]

定理 8.14 (单步估计量的渐近等效性)

假设中心极限定理成立 \(\sqrt{n}\Psi_n(\theta_0) \xrightarrow{d} Z\)，并且满足上述条件 (8.18)。对于任意 \(\sqrt{n}\)-相合的初始估计量 \(\tilde{\theta}_n\) 以及满足 \(\dot{\Psi}_{n,0} \xrightarrow{p} \dot{\Psi}_0\) 的导数估计，通过 (8.19) 构造的单步估计量 \(\hat{\theta}_n\) 满足：

\[ \sqrt{n}(\hat{\theta}_n - \theta_0) = -\dot{\Psi}_0^{-1} \sqrt{n} \Psi_n(\theta_0) + o_p(1) \quad \text{--- (8.20)} \]

这意味着单步估计量具有与真正的 Z-估计量完全相同的渐近分布。

定理 8.14 的证明（点击展开）

由单步更新公式 (8.19) 出发，两边同乘 \(\sqrt{n}\dot{\Psi}_{n,0}\) 得到：

\[ \dot{\Psi}_{n,0} \sqrt{n}(\hat{\theta}_n - \theta_0) = \dot{\Psi}_{n,0} \sqrt{n}(\tilde{\theta}_n - \theta_0) - \sqrt{n}\Psi_n(\tilde{\theta}_n) \]

将其改写并减去 \(\Psi_n(\theta_0)\) 的项：

\[ = \dot{\Psi}_{n,0} \sqrt{n}(\tilde{\theta}_n - \theta_0) - \sqrt{n}\{\Psi_n(\tilde{\theta}_n) - \Psi_n(\theta_0)\} - \sqrt{n}\Psi_n(\theta_0) \]

因为初始估计量是 \(\sqrt{n}\)-相合的，即 \(\sqrt{n}(\tilde{\theta}_n - \theta_0) = O_p(1)\)。这意味着通过选择足够大的 \(M\)，概率 \(P(\sqrt{n}\|\tilde{\theta}_n - \theta_0\| < M)\) 可以无限逼近于 1。因此，根据条件 (8.18)，在以概率趋于 1 的邻域内有：

\[ \sqrt{n}\{\Psi_n(\tilde{\theta}_n) - \Psi_n(\theta_0)\} = \dot{\Psi}_0 \sqrt{n}(\tilde{\theta}_n - \theta_0) + o_p(1) \]

将此表达式代入之前的等式：

\[ \dot{\Psi}_{n,0} \sqrt{n}(\hat{\theta}_n - \theta_0) = \dot{\Psi}_{n,0} \sqrt{n}(\tilde{\theta}_n - \theta_0) - \dot{\Psi}_0 \sqrt{n}(\tilde{\theta}_n - \theta_0) + o_p(1) - \sqrt{n}\Psi_n(\theta_0) \]

提取公因式并合并：

\[ = (\dot{\Psi}_{n,0} - \dot{\Psi}_0) \sqrt{n}(\tilde{\theta}_n - \theta_0) - \sqrt{n}\Psi_n(\theta_0) + o_p(1) \]

已知 \(\dot{\Psi}_{n,0} - \dot{\Psi}_0 = o_p(1)\)，且 \(\sqrt{n}(\tilde{\theta}_n - \theta_0) = O_p(1)\)，它们的乘积为 \(o_p(1)\)：

\[ = o_p(1) \cdot O_p(1) - \sqrt{n}\Psi_n(\theta_0) + o_p(1) = -\sqrt{n}\Psi_n(\theta_0) + o_p(1) \]

最后，应用 Slutsky 引理（左侧再乘上 \(\dot{\Psi}_{n,0}^{-1}\)，由于 \(\dot{\Psi}_{n,0}^{-1} \xrightarrow{p} \dot{\Psi}_0^{-1}\)），直接推导出：

\[ \sqrt{n}(\hat{\theta}_n - \theta_0) = -\dot{\Psi}_0^{-1} \sqrt{n} \Psi_n(\theta_0) + o_p(1) \]

证明完毕。 \(\square\)

补充说明： 如果 \(\Psi_n(\theta) = \mathbb{P}_n \psi_\theta\)，那么在定理 8.11（标准的 Z-估计渐近正态性定理）的条件下，条件 (8.18) 是自然满足的，此时 \(\dot{\Psi}_0 = V_{\theta_0}\)；在 van der Vaart (1998) 的定理 5.41 条件下，该条件也成立，且 \(\dot{\Psi}_0 = P\dot{\psi}_{\theta_0}\)。

10.2 示例：柯西位置模型的单步估计量

设 \(X_1, X_2, \dots, X_n \sim i.i.d.\ F\)，其概率密度函数为 Cauchy 分布：

\[ p_\theta(x) = \frac{1}{\pi(1 + (x - \theta)^2)} \]

似然得分函数 (score function) 为：

\[ \dot{l}_\theta(x) = \frac{\dot{p}_\theta(x)}{p_\theta(x)} = \frac{2(x - \theta)}{1 + (x - \theta)^2} \]

对于较大的 \(x\) 而言，得分函数 \(\dot{l}_\theta(x) \sim \frac{1}{x}\)，因此它拥有有限的 Fisher 信息量：\(E(\dot{l}_\theta^2(X)) < \infty\)。由于该方程存在多个根 (multiple roots)，直接通过数值优化寻找全局最大似然估计 (MLE) 是非常困难的。

更稳妥的策略是使用单步估计量：以样本中位数作为初始估计量 \(\tilde{\theta}_n\)（中位数是稳定的 \(\sqrt{n}\)-相合估计），然后在其附近对似然得分方程执行一次单步更新。这样即保证了估计量的渐近效率（等同于MLE），又避开了多根寻优的计算陷阱。

10.3 示例：混合模型的单步估计量 (Mixture Model)

假设观测数据来自两个位置-尺度模型 (location-scale models) 的混合分布：

\[ X_1, X_2, \dots, X_n \sim i.i.d.\ pF_1 + (1-p)F_2 \]

其中：

\[ f_1(x) = f\left(\frac{x - \mu}{\sigma}\right)\frac{1}{\sigma} \quad \text{且} \quad f_2(x) = g\left(\frac{x - \nu}{\tau}\right)\frac{1}{\tau} \]

这里，\((\mu, \nu)\) 是位置参数，\((\sigma, \tau)\) 是尺度参数。参数向量为 \(\theta = (\mu, \sigma, \nu, \tau, p)\)。混合概率密度为：

\[ f(x) = p f\left(\frac{x - \mu}{\sigma}\right)\frac{1}{\sigma} + (1-p) g\left(\frac{x - \nu}{\tau}\right)\frac{1}{\tau} \]

似然函数为：

\[ L_n(\theta) = \prod_{i=1}^n \left\{ p f\left(\frac{X_i - \mu}{\sigma}\right)\frac{1}{\sigma} + (1-p) g\left(\frac{X_i - \nu}{\tau}\right)\frac{1}{\tau} \right\} \]

如果我们将参数设置为 \(\mu = X_1\) 并且让 \(\sigma \rightarrow 0\)（令第一个分布坍缩到第一个数据点上），似然函数中的一项趋于无穷大。这会导致整个似然函数无界：

\[ \sup L_n(\theta) = \infty \]

这被称为似然函数的无界性崩溃问题。标准的 MLE 无法直接使用。

单步策略 (One-step strategy): 由于似然函数整体不存在有限的全局极大值，我们可以从另一个稳定的估计量开始。首先，通过匹配前五阶矩的方法求解出矩估计量 (Method of Moment Estimators) \(\hat{\theta}_{me}\)，这个估计量不依赖于似然函数的形状，是相合的。然后，以 \(\hat{\theta}_{me}\) 为初始点，在附近进行单步最大化或利用牛顿步更新，从而获得渐近有效且统计性质良好的最优估计量。