Two-Stage M-Estimation

• Two-stage M-estimation is a sequential statistical framework that leverages a preliminary estimator to localize the parameter space and improve efficiency in the subsequent estimation stage.
• It decomposes complex estimation tasks into practical subproblems—such as adaptive design, nuisance plug-in, and reweighting—to enhance computational tractability and accuracy.
• The method effectively addresses challenges like non-smoothness and high dimensionality by integrating preliminary bias corrections and variance adjustments across stages.

Two-stage M-estimation denotes a broad class of procedures in which a preliminary stage constructs an estimator, design, nuisance proxy, or low-dimensional representation, and a subsequent stage solves an M-estimation problem conditional on that output. In the settings covered on arXiv, this architecture appears in adaptive multistage sampling, covariate-shift correction, control-function estimation for endogenous Tobit models, composite likelihood for multivariate probit models, robust high-dimensional variable selection, simulation-based decision theory, and self-supervised pre-training followed by downstream fitting. The common statistical feature is that the second-stage criterion is random not only because of sampling noise, but also because it inherits the uncertainty, geometry, and sometimes the symmetries of the first stage (Mallik et al., 2014, Zhang et al., 30 Jun 2025, Tinati et al., 29 Mar 2026).

1. General formulation

In classical M-estimation, one minimizes an empirical criterion such as $\sum_{i=1}^N \rho(y_i;\theta)$, or equivalently solves an estimating equation of the form $\sum_{i=1}^N \psi(V_i;\theta)=0$. Two-stage variants retain that basic structure, but replace a single criterion by a sequential composition. One common form is adaptive multistage design: a pilot estimator

$$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$

is computed from a first batch, and then a second-stage design $P_2(\cdot\mid \hat\theta_{n_1})$ is chosen on a localized region, after which the final estimator is either a weighted pooled argmin or the argmin of a second-stage criterion alone (Mallik et al., 2014).

A second form is nuisance plug-in. Here the target parameter $\theta$ enters a population objective through an unknown function $g_0$, so the oracle problem is

$$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$

A two-stage estimator first builds $\hat g$ nonparametrically and then solves

$$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$

This structure underlies the two-stage maximum score framework and many semiparametric plug-in estimators (Gao et al., 2020).

A third form uses weighted estimating equations. Under two-phase multiwave sampling, a valid M-estimator solves

$\sum_{i=1}^N W_i\,\psi(V_i;\theta)=0,$

where the multiwave inverse-probability weights $\sum_{i=1}^N \psi(V_i;\theta)=0$0 correct for adaptive observation of expensive measurements. In that setting the weights themselves are products of the first-stage design decisions, so the stage distinction is built directly into the estimating equation (Kluger et al., 18 Feb 2026).

A fourth form appears in representation learning. Pre-training solves

$\sum_{i=1}^N \psi(V_i;\theta)=0$1

and downstream fitting then solves

$\sum_{i=1}^N \psi(V_i;\theta)=0$2

This is explicitly analyzed as two-stage M-estimation in self-supervised pre-training, where the first-stage parameter may be identifiable only up to a symmetry group (Tinati et al., 29 Mar 2026).

The resulting variants differ in what is carried from stage 1 to stage 2: a localized design, a nuisance estimate, an influence-weight, a control-function residual, or a learned representation.

Setting	Stage 1	Stage 2
Multistage sampling	Pilot argmin from $\sum_{i=1}^N \psi(V_i;\theta)=0$3	Localized design and second-stage argmin
Covariate shift	Estimate $\sum_{i=1}^N \psi(V_i;\theta)=0$4 under $\sum_{i=1}^N \psi(V_i;\theta)=0$5	Likelihood-ratio reweighting or truncated debiasing
Endogenous Tobit	OLS for $\sum_{i=1}^N \psi(V_i;\theta)=0$6 and residuals $\sum_{i=1}^N \psi(V_i;\theta)=0$7	Tobit M-estimation for $\sum_{i=1}^N \psi(V_i;\theta)=0$8
Multivariate probit	$\sum_{i=1}^N \psi(V_i;\theta)=0$9 univariate probits	$$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$0 bivariate pairwise fits
Bi-level selection	Group-penalized robust M-estimation	Hard-thresholding within groups
Pre-training/fine-tuning	Minimize $$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$1	Minimize downstream loss with $$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$2

2. Statistical rationale and efficiency mechanisms

The principal motivation for two-stage M-estimation is that a first-stage estimate can make the second-stage problem more informative, more localized, or more computationally tractable than a one-stage analogue. In multistage sampling problems, the stage-1 estimator identifies a “zoomed-in” region, and the second-stage design concentrates sampling budget there. For change-point estimation with signal size $$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$3, $$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$4, the two-stage design yields

$$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$5

an accelerated rate relative to the one-stage $$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$6. The same paper reports analogous acceleration from $$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$7 to $$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$8 in monotone inverse problems (Mallik et al., 2014).

An allied mechanism appears in one-step refinement for Markov sequences. There, a relatively short learning sample of size $$\hat\theta_{n_1}=\arg\min_{\theta\in\Theta}\mathcal M_{1,n}(\theta)$$9 with $P_2(\cdot\mid \hat\theta_{n_1})$0 produces a preliminary $P_2(\cdot\mid \hat\theta_{n_1})$1-consistent estimator, and the remaining observations are used in a one-step MLE update: $P_2(\cdot\mid \hat\theta_{n_1})$2 Under the stated regularity conditions, $P_2(\cdot\mid \hat\theta_{n_1})$3 is asymptotically equivalent to the full MLE,

$P_2(\cdot\mid \hat\theta_{n_1})$4

so the second stage recovers full asymptotic efficiency while only requiring a short pilot interval (Kutoyants et al., 2016).

In computational statistics, the gain can be algorithmic rather than rate-theoretic. For multivariate probit models, a two-stage composite likelihood replaces a single $P_2(\cdot\mid \hat\theta_{n_1})$5-variate integral by $P_2(\cdot\mid \hat\theta_{n_1})$6 univariate and $P_2(\cdot\mid \hat\theta_{n_1})$7 bivariate subproblems. Stage 1 fits the marginal regression parameters $P_2(\cdot\mid \hat\theta_{n_1})$8; stage 2 fits the pairwise correlations $P_2(\cdot\mid \hat\theta_{n_1})$9 conditional on the first-stage margins. The procedure is described as “trivially parallelizable and much faster than full-dimensional GHK or MC,” while empirical studies on the Six Cities and MEPS data are reported to show minor loss of efficiency, “<5 %,” relative to full ML (Ting et al., 2020).

These examples show that two-stage M-estimation is not tied to a single source of gain. Depending on the problem, stage 1 may concentrate design mass, generate a near-efficient starting value, or decompose a numerically difficult criterion into manageable pieces.

3. Dependence across stages and asymptotic theory

The central technical difficulty is that the second-stage criterion is usually conditionally i.i.d. only after conditioning on the first-stage output. In the generic multistage framework, the proofs proceed by conditioning on the event that the pilot estimate lies in a high-probability set, treating $\theta$0 as fixed, applying i.i.d. empirical-process theory to the second-stage sample, and then integrating out the randomness of the pilot. The key conditions are identifiability and curvature,

$\theta$1

together with a modulus bound for the second-stage empirical process. If $\theta$2 is non-increasing for some $\theta$3 and $\theta$4 satisfies

$\theta$5

then

$\theta$6

After localization, one studies

$\theta$7

and obtains a weak limit $\theta$8, typically a Gaussian process with deterministic drift, followed by an argmin theorem (Mallik et al., 2014).

In two-stage control-function Tobit estimation, the first-stage OLS fit for the endogenous regressor contributes directly to second-stage uncertainty. The joint asymptotic normality of $\theta$9 has block matrices

$g_0$0

and the marginal asymptotic variance of $g_0$1 contains the additional term

$g_0$2

where $g_0$3. The paper explicitly identifies this extra term as the effect of first-stage sampling variability (Shukla et al., 2023).

A closely related correction appears in two-stage composite likelihood. Writing the stacked estimating function as

$g_0$4

standard M-estimation yields

$g_0$5

but the stage-2 covariance must account for the fact that $g_0$6 is evaluated at $g_0$7 rather than $g_0$8. The paper therefore gives a Murphy–Topel-type robust covariance formula for the second-stage estimator (Ting et al., 2020).

In two-stage self-supervised learning, the dependence structure is more intricate because the first-stage representation may be identifiable only up to symmetry. With $g_0$9, the pre-training estimator satisfies

$$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$0

while the downstream least-squares coefficients satisfy a conditional $$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$1-CLT. The final downstream test risk then splits into the usual well-specified OLS variance term and an additional $$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$2-scaled term driven by pre-training fluctuations (Tinati et al., 29 Mar 2026).

4. Reweighting, debiasing, and robustness

A major strand of two-stage M-estimation uses stage 2 to debias or recalibrate a stage-1 estimator. Under covariate shift, the target is the moment

$$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$3

with source and target densities $$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$4 and $$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$5. In the idealized setting with known densities, the two-stage estimator first trains $$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$6 under $$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$7 and then forms the importance-weighted plug-in estimator

$$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$8

The paper establishes a minimax lower bound for the risk and shows that this two-stage estimator attains the minimax upper bound up to logarithmic factors (Zhang et al., 30 Jun 2025).

When $$\theta_0\in\arg\max_{\theta\in\Theta} Q(\theta,g_0), \qquad Q(\theta,g_0)=\mathbb E[m(W_i,\theta,g_0(Z_i))].$$9 and $\hat g$0 are unknown, density-ratio estimation can be unstable. The same work introduces a truncated, double-robust estimator

$\hat g$1

If the ratio tail obeys $\hat g$2 and $\hat g$3, then

$\hat g$4

Moreover, the estimator is double robust: it remains consistent if either $\hat g$5 or $\hat g$6 is consistent, with $\hat g$7 slowly (Zhang et al., 30 Jun 2025).

A parallel debiasing logic appears in two-phase multiwave sampling. The Multiwave Predict-Then-Debias estimator is

$\hat g$8

where $\hat g$9 uses cheap proxy data, $$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$0 applies inverse-probability weighting to the proxies, and $$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$1 uses the gold-standard measurements. The same estimator can be characterized as the minimizer of a debiased loss built from proxy predictions and a first-order correction term. Under the paper’s smoothness and overlap conditions, the estimator is asymptotically linear and asymptotically normal, and componentwise Wald intervals are asymptotically valid (Kluger et al., 18 Feb 2026).

These constructions clarify a recurring theme: stage 1 often supplies a biased or proxy-based surrogate, and stage 2 is designed not merely to refine it numerically but to restore the correct target through reweighting, residual correction, or inverse-probability debiasing.

5. Non-smooth, high-dimensional, and structured regimes

Two-stage M-estimation is especially useful when the second-stage objective is non-smooth or structurally constrained. In the two-stage maximum score estimator, the nuisance function

$$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$2

is estimated nonparametrically in stage 1, and stage 2 maximizes

$$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$3

The first stage acts as an “imperfect smoother” for the otherwise non-smooth criterion. The resulting rate depends on the covariate dimension $$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$4 and features phase transitions at $$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$5 and $$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$6. In particular,

$$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$7

with low-, intermediate-, and high-dimensional regimes having distinct optimal bandwidths and distinct limiting behavior, including non-Gaussian Gaussian-process argmax limits and bias-dominated limits (Gao et al., 2020).

In high-dimensional bi-level variable selection, stage 1 solves a robust group-penalized M-estimation problem,

$$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$8

using losses such as Huber, Tukey’s biweight, or Cauchy and nonconvex group penalties such as group-MCP or SCAD. Stage 2 then applies elementwise hard-thresholding,

$$\hat\theta = \arg\max_{\theta\in\Theta} \frac1n\sum_{i=1}^n m\bigl(W_i,\theta,\hat g(Z_i)\bigr).$$9

to recover within-group sparsity. Under the stated restricted strong convexity and signal conditions, the paper proves local estimation consistency for the stage-1 estimator, a group-support oracle property, and bi-level consistency after hard-thresholding (Luo et al., 2019).

These examples show that two-stage M-estimation is not confined to smooth likelihood theory. It also organizes non-smooth, nonconvex, and high-dimensional procedures in which the first stage regularizes geometry or sparsity structure before a second-stage decision is taken.

6. Conceptual scope, symmetries, and recurring misconceptions

A broad decision-theoretic formulation makes clear that two-stage M-estimation is not reducible to sample splitting. When direct likelihood optimization is unavailable, one may restrict the decision rule to the composite form

$\sum_{i=1}^N W_i\,\psi(V_i;\theta)=0,$0

where $\sum_{i=1}^N W_i\,\psi(V_i;\theta)=0,$1 is a fixed compression map and $\sum_{i=1}^N W_i\,\psi(V_i;\theta)=0,$2 is optimized over a function class. The Bayesian version minimizes a Monte Carlo approximation to Bayes risk, while the minimax version solves a convex program that bounds worst-case empirical loss over simulated parameter draws. For i.i.d. data, the paper proposes empirical quantiles as the first-stage summaries and an affine second-stage map; under weak regularity, it states consistency, asymptotic normality, Monte Carlo approximation error of order $\sum_{i=1}^N W_i\,\psi(V_i;\theta)=0,$3, and statistical error of order $\sum_{i=1}^N W_i\,\psi(V_i;\theta)=0,$4 (Lakshminarayanan et al., 2022).

Another conceptual extension is symmetry. In self-supervised pre-training, the first-stage parameter can be invariant under a compact Lie group $\sum_{i=1}^N W_i\,\psi(V_i;\theta)=0,$5, so the true object is an orbit $\sum_{i=1}^N W_i\,\psi(V_i;\theta)=0,$6 rather than a single Euclidean vector. The theory therefore works on the quotient manifold, using local descriptor charts that are constant on orbits. If there is orthogonal equivariance,

$\sum_{i=1}^N W_i\,\psi(V_i;\theta)=0,$7

then the downstream hypothesis class is orbit-invariant and the minimum-norm downstream solution is well defined on the quotient (Tinati et al., 29 Mar 2026).

Several misconceptions are corrected by this literature. First, two-stage M-estimation is not only nuisance plug-in; it also includes adaptive design, one-step likelihood refinement, composite likelihood factorization, debiasing under two-phase sampling, and pre-training/fine-tuning pipelines. Second, the first stage is not automatically asymptotically negligible: it can generate additional variance terms, non-Gaussian limits, or bias-dominated behavior (Shukla et al., 2023, Gao et al., 2020, Tinati et al., 29 Mar 2026). Third, the second stage does not always stabilize the procedure by itself; under unknown covariate shift, importance weights may spike, and truncation is introduced precisely because ratio estimation can be unstable (Zhang et al., 30 Jun 2025).

Taken together, these results suggest a unifying interpretation. Two-stage M-estimation is a sequential statistical design principle in which the first stage reshapes the criterion seen by the second stage. That reshaping may localize the parameter space, alter the sampling law, replace latent nuisance functions by estimates, introduce proxy corrections, or quotient out symmetries. The success of the method depends on analyzing the induced stage dependence rather than treating the first stage as a preprocessing step external to inference.