Chapter 12: M-Estimators and Z-Estimators (Part II)

In the previous part, we discussed the consistency and asymptotic normality of standard M-estimators and Z-estimators. In practical statistical modeling, we often encounter situations involving nuisance parameters or requiring iterative algorithms (such as the Newton-Raphson method) to obtain solutions. This part delves into Z-estimators with estimated parameters and one-step estimators.

8. Estimators with Estimated Nuisance Parameters

Suppose the samples \(X_1, \dots, X_n \sim i.i.d.\ F_{(\theta, \eta)}\), where:

\(\theta\) is the parameter of interest.
\(\eta\) is an unknown nuisance parameter.

Typically, we need to substitute an estimator of \(\eta\) (denoted \(\hat{\eta}_n\)) into the Z-estimating equation:

\[ 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 \]

This is essentially a two-step procedure: first estimate \(\eta\), then treat it as a "known" parameter and solve for \(\theta\) in the equation.

8.1 Motivating Example: Missing Response Regression

Consider the regression model:

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

Assume the response \(y_i\) is subject to missingness. We define an indicator variable: \(R_i = 1\) if \(y_i\) is observed; \(R_i = 0\) if \(y_i\) is missing.

Missing at Random (MAR) Assumption: Assume the missingness mechanism satisfies:

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

where \(\omega_{\eta_0}(\cdot)\) is called the missing propensity function. Here we assume it follows a parametric binary regression model.

Meaning of MAR: Given \(X_i\), the missing indicator \(R_i\) and the true response \(y_i\) are conditionally independent.
This missingness mechanism is called ignorable missing at random, because as long as we observe \(x_i\), the missing \(y_i\) does not introduce additional bias for inference.

How to estimate the target parameter \(\theta\)?

Method 1: Least squares estimation (LSE) using only the observed data (i.e., \(R_i = 1\))

The objective function is:

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

Taking derivatives gives the score equation:

\[ \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)} \]

At the true parameter \(\theta_0\), using \(E(\epsilon_i | x_i) = 0\) and under the MAR assumption:

\[ 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 \]

Thus, under certain regularity conditions, the least squares estimator obtained by solving (8.13) is consistent and asymptotically normal (AN).

Method 2: Inverse Probability Weighted Estimation (IPW)

By weighting (8.13) with the true missing propensity \(\omega_{\eta_0}(x_i)\), we obtain the IPW score equation:

\[ \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)} \]

Taking expectations also verifies unbiasedness:

\[ 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 \]

Question: Which estimator is more efficient (i.e., has smaller asymptotic variance): the estimator ignoring missing values (8.13) or the IPW estimator (8.14)?

8.2 Binary Likelihood Estimation of the Nuisance Parameter

In practice, \(\eta_0\) is usually unknown. We can estimate it by maximizing the binary likelihood function:

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

The log-likelihood function is:

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

Differentiating and setting to zero yields \(\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} \]

Then, (8.14) becomes an equation with an estimated nuisance parameter:

\[ \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)} \]

Counterintuitive conclusion: Chen, Leung and Qin (2008) proved that the \(\theta\) obtained from equation (8.15) using the estimated parameter \(\hat{\eta}\) is asymptotically more efficient than the \(\theta\) obtained from equation (8.14) using the true known parameter \(\eta_0\)!

8.3 Nonparametric Estimation of the Missing Propensity

Assuming the parametric model \(P(R_i = 1 | x_i) = \omega_{\eta_0}(x_i)\) may be too restrictive. We can consider a nonparametric form:

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

The missing propensity \(\omega(\cdot)\) can be estimated via kernel smoothing methods (e.g., the Nadaraya-Watson estimator):

\[ \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)} \]

where \(K\) is a kernel function (symmetric probability density function) and \(h\) is the smoothing bandwidth. We require that as \(n \rightarrow \infty\), \(h \rightarrow 0\) and \(nh \rightarrow \infty\). Its expectation is essentially a weighted average of the weights at the sample points:

\[ 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) \]

It can be shown that if we choose an appropriate bandwidth \(h \approx O(n^{-\frac{1}{4+d}})\) (where \(d\) is the dimension of \(x\)), then:

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

9. Asymptotic Theorem for Z-Estimators with Estimated Parameters

Theorem 8.13 (Asymptotic Expansion for Two-Step Z-Estimators)

For every \(\theta\) in an open subset of \(\mathbb{R}^k\) and every \(\eta\) in a metric space, let \(x \rightarrow \phi_{\theta, \eta}(x)\) be a \(\mathbb{R}^k\)-valued measurable function, such that: the class of functions \(\{\phi_{\theta, \eta} : \|\theta - \theta_0\| < \delta, d(\eta, \eta_0) < \delta\}\) is Donsker for some \(\delta > 0\); and \(P\|\phi_{\theta, \eta} - \phi_{\theta_0, \eta_0}\|^2 \rightarrow 0\) as \((\theta, \eta) \rightarrow (\theta_0, \eta_0)\).

Assume \(P\phi_{\theta_0, \eta_0} = 0\), and the map \(\theta \rightarrow P\phi_{\theta, \eta}\) is differentiable at \(\theta_0\) (uniformly in \(\eta\) in a neighborhood of \(\eta_0\)), with a nonsingular derivative matrix \(V_{\theta_0, \eta}\), and \(V_{\theta_0, \eta} \rightarrow V_{\theta_0, \eta_0}\) as \(\eta \rightarrow \eta_0\).

Furthermore, if there exist consistent estimators \((\hat{\theta}, \hat{\eta}) \xrightarrow{p} (\theta_0, \eta_0)\) satisfying \(\sqrt{n}\mathbb{P}_n\phi_{\hat{\theta}, \hat{\eta}} = o_p(1)\), then:

\[ \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}\|) \]

Remarks on Theorem 8.13:

The drift term on the right-hand side, \(-V_{\theta_0, \eta_0}^{-1} \sqrt{n} P\psi_{\theta_0, \hat{\eta}_n}\), arises because we have additionally estimated the nuisance parameter \(\eta\).
Here \(\eta\) can be an infinite-dimensional functional parameter (not just a scalar or vector parameter).
Under certain smoothness conditions, this extra drift term can be approximated by a Taylor expansion using the Delta method:

\[ 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: If the target parameter and the nuisance parameter satisfy an (approximate) orthogonality condition, i.e., the cross partial derivative is 0 (\(\frac{\partial P\psi_{\theta_0, \eta_0}}{\partial \eta} = 0\)), then the error in estimating the nuisance parameter \(\hat{\eta}_n\) will not affect the first-order asymptotic distribution of \(\hat{\theta}_n\) (the drift term vanishes).

9.1 Example: Symmetric Location Model

Suppose \(X_1, \dots, X_n \sim i.i.d.\ F\), and the distribution \(F\) is symmetric about the true location parameter \(\theta_0\). We want to estimate \(\theta_0\), but the scale parameter \(\sigma\) is unknown. Let \(x \rightarrow \psi(x)\) be an antisymmetric odd function.

Consider the following Z-estimating equation (with \(\hat{\sigma}\) being a consistent estimator of \(\sigma\)):

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

Since \(F(\cdot)\) is symmetric about \(\theta_0\) and \(\psi(\cdot)\) is an odd function, under the true parameter, the expectation is zero for any scale estimator \(\hat{\sigma}\):

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

Therefore, by Theorem 8.13, since the drift term is zero, we obtain:

\[ \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) \]

This means that the estimation procedure is asymptotically equivalent to using the true \(\eta_0\). In this model, the estimation error of \(\hat{\sigma}\) does not appear in the first-order asymptotic expansion.

10. One-step Estimator

In many complex models, the exact root of the Z-estimating equation \(\Psi_n(\theta) = 0\) can be very difficult to obtain (e.g., multiple local extrema, high nonlinearity).

We can improve the estimation by performing one optimization step based on an initial estimator \(\tilde{\theta}_n\). This is exactly the idea of the Newton-Raphson algorithm.

Consider a first-order Taylor expansion around \(\tilde{\theta}_n\):

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

We obtain a one-step update of \(\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)} \]

Usually, we can choose a computationally easy but inefficient \(\sqrt{n}\)-consistent estimator as the initial value \(\tilde{\theta}_n\) (e.g., method of moments estimator). We will see that the one-step estimator \(\hat{\theta}_n\) achieves exactly the same asymptotic efficiency as the fully iterated Z-estimator that would otherwise be difficult to compute.

10.1 Regularity Conditions

Uniform convergence condition: For any constant \(M\) and a given nonsingular matrix \(\dot{\Psi}_0\), we have:

\[ \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)} \]

where \(\dot{\Psi}_0\) is the limit matrix of the derivative \(\dot{\Psi}_n(\theta)\) as \(\theta \rightarrow \theta_0\).

In practice, we can replace the empirical derivative in (8.17) with any consistent estimate matrix \(\dot{\Psi}_{n,0}\) (satisfying \(\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)} \]

Theorem 8.14 (Asymptotic Equivalence of One-step Estimators)

Assume the central limit theorem holds: \(\sqrt{n}\Psi_n(\theta_0) \xrightarrow{d} Z\), and condition (8.18) above is satisfied. For any \(\sqrt{n}\)-consistent initial estimator \(\tilde{\theta}_n\) and a derivative estimator satisfying \(\dot{\Psi}_{n,0} \xrightarrow{p} \dot{\Psi}_0\), the one-step estimator \(\hat{\theta}_n\) constructed by (8.19) satisfies:

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

This means the one-step estimator has exactly the same asymptotic distribution as the true Z-estimator.

Proof of Theorem 8.14 (click to expand)

Starting from the one-step update formula (8.19), multiply both sides by \(\sqrt{n}\dot{\Psi}_{n,0}\) to obtain:

\[ \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) \]

Rewrite by subtracting \(\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) \]

Since the initial estimator is \(\sqrt{n}\)-consistent, i.e., \(\sqrt{n}(\tilde{\theta}_n - \theta_0) = O_p(1)\), we can choose a large enough \(M\) such that \(P(\sqrt{n}\|\tilde{\theta}_n - \theta_0\| < M)\) is arbitrarily close to 1. Therefore, by condition (8.18), on the neighborhood where this probability tends to 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) \]

Substituting this expression into the previous equation:

\[ \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) \]

Factor and combine terms:

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

Since \(\dot{\Psi}_{n,0} - \dot{\Psi}_0 = o_p(1)\) and \(\sqrt{n}(\tilde{\theta}_n - \theta_0) = O_p(1)\), their product is \(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) \]

Finally, applying Slutsky's lemma (multiply on the left by \(\dot{\Psi}_{n,0}^{-1}\), which converges in probability to \(\dot{\Psi}_0^{-1}\)), we directly obtain:

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

The proof is complete. \(\square\)

Additional remarks: If \(\Psi_n(\theta) = \mathbb{P}_n \psi_\theta\), then under the conditions of Theorem 8.11 (standard Z-estimator asymptotic normality), condition (8.18) is naturally satisfied, with \(\dot{\Psi}_0 = V_{\theta_0}\); under the conditions of van der Vaart (1998) Theorem 5.41, this condition also holds, with \(\dot{\Psi}_0 = P\dot{\psi}_{\theta_0}\).

10.2 Example: One-step Estimator for the Cauchy Location Model

Let \(X_1, X_2, \dots, X_n \sim i.i.d.\ F\), with probability density function given by the Cauchy distribution:

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

The score function is:

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

For large \(x\), the score function \(\dot{l}_\theta(x) \sim \frac{1}{x}\), so it has finite Fisher information: \(E(\dot{l}_\theta^2(X)) < \infty\). Since this equation has multiple roots, directly finding the global maximum likelihood estimator (MLE) via numerical optimization is very difficult.

A safer strategy is to use the one-step estimator: Take the sample median as the initial estimator \(\tilde{\theta}_n\) (the median is a stable \(\sqrt{n}\)-consistent estimator), then perform a single Newton step on the score equation near it. This ensures the asymptotic efficiency of the estimator (equivalent to the MLE) while avoiding the computational pitfalls of multiple roots.

10.3 Example: One-step Estimator for a Mixture Model

Suppose the observed data come from a mixture of two location-scale models:

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

where:

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

Here, \((\mu, \nu)\) are location parameters, \((\sigma, \tau)\) are scale parameters. The parameter vector is \(\theta = (\mu, \sigma, \nu, \tau, p)\). The mixture density is:

\[ 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} \]

The likelihood function is:

\[ 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\} \]

If we set the parameter \(\mu = X_1\) and let \(\sigma \rightarrow 0\) (i.e., let the first distribution collapse onto the first data point), one term in the likelihood tends to infinity. This causes the entire likelihood function to become unbounded:

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

This is known as the unbounded likelihood problem. Standard MLE cannot be used directly.

One-step strategy: Since the likelihood function does not have a finite global maximum, we start from another stable estimator. First, obtain method of moments estimators \(\hat{\theta}_{me}\) by matching the first five moments. This estimator is consistent and does not depend on the shape of the likelihood. Then, using \(\hat{\theta}_{me}\) as the starting point, perform a single-step maximization or a Newton step nearby to obtain an asymptotically efficient and statistically well-behaved optimal estimator.