svyreg_huberM: Huber Robust Survey Regression M- and GM-Estimator

Description

svyreg_huberM and svyreg_huberGM compute, respectively, a survey weighted M- and GM-estimator of regression using the Huber psi-function.

Usage

svyreg_huberM(formula, design, k, var = NULL, na.rm = FALSE, asym = FALSE,
              verbose = TRUE, ...)
svyreg_huberGM(formula, design, k, type = c("Mallows", "Schweppe"),
               xwgt, var = NULL, na.rm = FALSE, asym = FALSE, verbose = TRUE,
               ...)

Value

Object of class svyreg.rob

Arguments

formula: a [formula] object (i.e., symbolic description of the model)
design: an object of class survey.design; see svydesign.
k: [double] robustness tuning constant (\(0 < k \leq \infty\)).
var: a one-sided [formula] object or variable name ([character]) that defines the heteroscedastic variance or [NULL] indicating homoscedastic variance (default: NULL).
na.rm: [logical] indicating whether NA values should be removed before the computation proceeds (default: FALSE).
asym: [logical] toggle for asymmetric Huber psi-function (default: FALSE).
verbose: [logical] indicating whether additional information is printed to the console (default: TRUE).
type: [character] "Mallows" or "Schweppe".
xwgt: [numerical vector] or [NULL] of weights in the design space (default: NULL); xwgt is only relevant for type = "Mallows" or type = "Schweppe".
...: additional arguments passed to the method (e.g., maxit: maxit number of iterations, etc.).

Failure of convergence

By default, the method assumes a maximum number of maxit = 100 iterations and a numerical tolerance criterion to stop the iterations of tol = 1e-05. If the algorithm fails to converge, you may consider changing the default values; see svyreg_control.

Details

Package survey must be attached to the search path in order to use the functions (see library or require).

svyreg_huberM and svyreg_huberGM compute, respectively, M- and GM-estimates of regression by iteratively re-weighted least squares (IRWLS). The estimate of regression scale is (by default) computed as the (normalized) weighted median of absolute deviations from the weighted median (MAD; see weighted_mad) for each IRWLS iteration. If the weighted MAD is zero (or nearly so), the scale is computed as the (normalized) weighted interquartile range (IQR).

M-estimator

The regression M-estimator svyreg_huberM is robust against residual outliers (granted that the tuning constant k is chosen appropriately).

GM-estimator

Function svyreg_huberGM implements the Mallows and Schweppe regression GM-estimator (see argument type). The regression GM-estimators are robust against residual outliers and outliers in the model's design space (leverage observations; see argument xwgt).

Numerical optimization

See svyreg_control.

Models

Models for svyreg_rob are specified symbolically. A typical model has the form response ~ terms, where response is the (numeric) response vector and terms is a series of terms which specifies a linear predictor for response; see formula and lm.

A formula has an implied intercept term. To remove this use either y ~ x - 1 or y ~ 0 + x; see formula for more details of allowed formulae.

Examples

Run this code

head(workplace)

library(survey)
# Survey design for stratified simple random sampling without replacement
dn <- if (packageVersion("survey") >= "4.2") {
        # survey design with pre-calibrated weights
        svydesign(ids = ~ID, strata = ~strat, fpc = ~fpc, weights = ~weight,
                  data = workplace, calibrate.formula = ~-1 + strat)
    } else {
        # legacy mode
        svydesign(ids = ~ID, strata = ~strat, fpc = ~fpc, weights = ~weight,
                  data = workplace)
    }

# Compute regression M-estimate with Huber psi-function
m <- svyreg_huberM(payroll ~ employment, dn, k = 8)

# Regression inference
summary(m)

# Extract the coefficients
coef(m)

# Extract variance/ covariance matrix
vcov(m)

# Diagnostic plots (e.g., standardized residuals against fitted values)
plot(m, which = 1L)

# Plot of the robustness weights of the M-estimate against its residuals
plot(residuals(m), robweights(m))

Run the code above in your browser using DataLab