Fast Robust (M/MM-Estimator) Linear Regression (C++)

Fits a robust linear regression by iteratively reweighted least squares (IRLS), minimizing \(\sum_i \rho(r_i / \hat\sigma)\) for residuals \(r_i = y_i - x_i^\top\beta\) and a fixed robustness scale \(\hat\sigma\), rather than ordinary least squares' \(\sum_i r_i^2\). method = "M" uses Huber's weight function \(w(u) = 1\) for \(|u| \le c\) and \(w(u) = c/|u|\) otherwise (c, default 1.345, tuned for 95% efficiency under normality); any other value of method (including the default, "MM") uses Tukey's bisquare weight \(w(u) = (1 - (u/c_b)^2)^2\) for \(|u| \le c_b\) and \(0\) otherwise, with \(c_b = 4.685\) hardcoded (not settable through this exported wrapper, though the internal fitter accepts it). The scale \(\hat\sigma\) is fixed once at the start as the normalized median absolute deviation of the OLS residuals, \(\hat\sigma = \mathrm{median}(|r_i|) / 0.6745\) (the internal fitter also accepts a caller-supplied fixed scale, but this wrapper always estimates it). Each IRLS iteration re-weights and re-solves the weighted normal equations \(X^\top W X\, \beta = X^\top W y\) via Eigen::LDLT; convergence is declared when the relative change in \(\beta\) falls below tol. Column 1 of a real design typically holds the intercept, but no columns are treated specially except via fixed_idx.

Usage

fast_robust_regression_cpp(
  X,
  y,
  warm_start_beta = NULL,
  smart_cold_start = TRUE,
  method = "MM",
  j = 2L,
  c = 1.345,
  maxit = 50L,
  tol = 1e-07,
  fixed_idx = NULL,
  fixed_values = NULL,
  warm_start_weights = NULL,
  warm_start_fisher_info = NULL,
  estimate_only = FALSE
)

Arguments

X

A numeric matrix of predictors.

y

A numeric vector of responses.

warm_start_beta

Optional starting values for coefficients. If provided, smart_cold_start is ignored.

smart_cold_start

Logical. If TRUE (the default) and no warm_start_beta is supplied, use an OLS (QR) initial guess; see Details.

method

Robust estimation method: "M" for Huber weighting, anything else (default "MM") for Tukey bisquare weighting; see Details.

j

1-based index of the coefficient whose asymptotic variance to return in ssq_b_j.

c

Huber tuning constant (default 1.345; only used when method = "M").

maxit

Maximum number of IRLS iterations.

tol

Relative parameter-change convergence tolerance.

fixed_idx

Optional indices of fixed parameters.

fixed_values

Optional values for fixed parameters.

warm_start_weights

Optional initial working weights for the first IRLS iteration.

warm_start_fisher_info

Optional initial curvature (\(X^\top W X\)) matrix for the first IRLS iteration.

estimate_only

If TRUE, skip the post-fit asymptotic-variance computation and return only coefficients, scale, converged, and iterations.

Value

If estimate_only = TRUE: a list with coefficients, scale (the fixed MAD-based robustness scale \(\hat\sigma\)), converged, iterations. Otherwise, additionally: ssq_b_j and fisher_information (the final IRLS \(X^\top W X\) curvature matrix). ssq_b_j is computed only if the fit converged or ran the full maxit iterations, as the standard M-estimator asymptotic variance \(\widehat{\mathrm{Var}}(\hat\beta_j) = \left(\frac{n}{n-p}\right) \frac{\sum_i \psi(r_i)^2}{n\,\bar\psi'^2} \, [(X^\top X)^{-1}]_{jj}\), where \(\psi\) is the derivative of \(\rho\) (i.e. \(\psi(r) = w(r/\hat\sigma)\,r\)) and \(\bar\psi'\) is the mean of \(\psi'\) across observations, matching the classical Huber (1981) sandwich-free M-estimator variance formula; NA if j indexes a fixed coefficient or the fit neither converged nor exhausted maxit.

Fixed parameters, warm starts

fixed_idx and fixed_values optionally hold a subset of coefficients fixed at caller-supplied constant values (subtracted out of y as an offset) rather than estimated. warm_start_beta supplies a starting coefficient vector directly; otherwise, if smart_cold_start = TRUE (the default), an ordinary QR least-squares fit seeds the start (and, when variance will later be requested via j, also caches the QR-based \([(X^\top X)^{-1}]_{jj}\) entry for reuse in the variance formula below). warm_start_weights seeds the IRLS weights for the first iteration only (skipping that iteration's Huber/bisquare weight computation); warm_start_fisher_info similarly seeds the first iteration's \(X^\top W X\) curvature matrix.