Safe (Generalized) Huber M-Estimator of Location

A (generalized) Huber M-estimator of location with MAD scale that handles the degenerate case of zero scale gracefully, where MASS::huber() would return an error.

Usage

huberM(
  x,
  conf.level = NA,
  sides = c("two.sided", "left", "right"),
  method = c("wald", "boot"),
  k = 1.345,
  mu = NULL,
  s = NULL,
  weights = NULL,
  na.rm = FALSE,
  ...
)

Arguments

x: numeric vector of data values
conf.level: confidence level of the interval. If set to NA (the default), only the point estimate is returned.
sides: character string specifying the sidedness of the confidence interval (one of "two.sided" (default), "left" or "right"). See ConfidenceIntervals(). The location is unbounded, so the open side is $\pm$Inf. One-sided intervals require conf.level > 0.5.
method: confidence interval method: "wald" (default) or "boot".
k: positive tuning constant; the algorithm winsorizes at k standard deviations. Default is 1.345.
mu: initial location estimate. NULL (default) uses median(x), computed after na.rm is applied.
s: scale estimate held constant through the iterations. NULL (default) uses mad(x, center = mu), computed after na.rm is applied.
weights: optional non-negative frequency weights of the same length as x. The starting values then become the weighted median and the weighted MAD (both via medianX(), the MAD with the usual constant 1.4826), and $n$ in the Wald interval is sum(weights), so integer weights give exactly the result of the replicated data. The bootstrap resamples observations together with their weights.
na.rm: logical; whether to remove missing values before computation; default is FALSE
...: further arguments passed to the bootstrap engine when method = "boot": R, type, parallel, and ncpus; see Details

Value

if conf.level = NA, a numeric scalar. Otherwise a named numeric vector with elements:

est: location estimate from Huber's M-estimator
lci: lower confidence interval bound
uci: upper confidence interval bound

Details

Wald interval (method = "wald")

Uses an empirical sandwich standard error for the Huber M-estimator: $$ \hat\mu \;\pm\; t_{1-\alpha/2,\,n-1} \hat s \sqrt{\frac{\hat\tau}{n}}, $$ where $$ \hat\tau = \frac{n \sum_{i=1}^n \psi_k(r_i)^2} {\left\{\sum_{i=1}^n I(|r_i| \leq k)\right\}^2}, \qquad r_i = \frac{x_i-\hat\mu}{\hat s}. $$ Here, $\psi_k$ denotes the Huber score function. No additional small-sample adjustment is applied beyond use of the t-quantile.

Bootstrap interval (method = "boot")

The statistic $\hat\mu$ is resampled $R$ times. Note that mu and s are fixed at their initial values (computed from the full data before resampling) and are not re-estimated on each resample. The bootstrap therefore targets the variability of the location estimator with fixed scale initialization, not a fully re-estimated robust estimator. Bootstrap arguments are passed through ... and extracted via .extractBootArgs():

R: number of bootstrap replicates (default 999)
type: confidence interval type: "perc", "basic", "norm" or "bca" (default)
parallel: parallelization mode: "no", "multicore", or "snow" (default "no")
ncpus: number of CPUs for parallel bootstrap (default getOption("boot.ncpus", 1L))

The original internal estimator is accessible as DescToolsX:::.huberM.

Note

Adapted from code by Martin Maechler to conform to package standards

Random number generation

method = "boot" resamples and therefore advances R's global random number generator. Call base::set.seed() beforehand for reproducible intervals.

References

Huber, P. J. and Ronchetti, E. M. (2009). Robust Statistics, 2nd ed. Wiley.

Examples

huberM(c(1:9, 1000))
#> [1] 5.553915
mad(c(1:9, 1000))
#> [1] 3.7065

set.seed(7)
x <- c(round(rnorm(1000), 1), round(rnorm(50, mean = 10, sd = 10)))

huberM(x, conf.level = 0.95)
#>         est         lci         uci 
#>  0.05600116 -0.01020176  0.12220408 
huberM(x, conf.level = 0.95, method = "boot", R = 499, type = "bca")
#>         est         lci         uci 
#>  0.05600116 -0.01018489  0.12185434 

# degenerate case: scale zero
huberM(rep(9, 100))
#> Warning: scale 's' is zero -- returning initial 'mu'
#> [1] 9