svyreg_tukeyM {robsurvey}

Tukey Biweight Robust Survey Regression M- and GM-Estimator

Description

svyreg_tukeyM and svyreg_tukeyGM compute, respectively, a survey weighted M- and GM-estimator of regression using the biweight Tukey psi-function.

Usage

svyreg_tukeyM(formula, design, k, var = NULL, na.rm = FALSE, verbose = TRUE,
              ...)
svyreg_tukeyGM(formula, design, k, type = c("Mallows", "Schweppe"),
               xwgt, var = NULL, na.rm = 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).

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_tukeyM and svyreg_tukeyGM 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_tukeyM 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.

See Also

Overview (of all implemented functions)

summary, coef, residuals, fitted, SE and vcov

plot for regression diagnostic plot methods.

Other robust estimating methods svyreg_huberM and svyreg_huberGM

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 Tukey bisquare psi-function
m <- svyreg_tukeyM(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))