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

svyreg_huberM {robsurvey}

R Documentation

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,
               ...)

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`).
`...`	additional arguments passed to the method (e.g., `maxit`: maxit number of iterations, etc.).
`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"`.

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.

Value

Object of class svyreg.rob

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.

Examples

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))

[Package robsurvey version 0.6 Index]