Vignette: Robust Horvitz-Thompson Estimator

1.1 Survey design object

In order to use the survey methods (not bare-bone methods), we must attach the survey package (Lumley, 2010, 2021) to the search path

> library("survey")

and specify a survey or sampling design object

> dn <- svydesign(ids = ~ID, strata = ~strat, fpc = ~fpc, weights = ~weight,
+                 data = workplace, calibrate.formula = ~-1 + strat)

Note. Since version 4.2, the survey package allows the definition of pre-calibrated weights (see argument calibrate.formula of the function svydesign()). This vignette uses this functionality (in some places). If you have installed an earlier version of the survey package, this vignette will automatically fall back to calling svydesign() without the calibration specification. See vignette Pre-calibrated weights of the survey package to learn more.

1.2 Exploring the data

To get a first impression of the distribution of payroll, we examine two (design-weighted) boxplots of payroll (on raw and logarithmic scale). The boxplots are obtained using function survey::svyboxplot.

From the boxplot with payroll on raw scale, we recognise that the sample distribution of payroll is skewed to the right; the boxplot on logarithmic scale demonstrates that log-transform is not sufficient to turn the skewed distribution into a symmetric distribution. The outliers need not be errors. Following Chambers (1986), we distinguish representative outliers from non-representative outliers (→ see vignette “Basic Robust Estimators” for an introduction to the notion of non-/ representative outliers).

The outliers visible in the boxplot refer to a few large workplaces. Moreover, we assume that these outliers represent other workplaces in the population that are similar in value (i.e., representative outliers).

2.1 Bare-bone methods

The following bare-bone estimating methods are available:

• weighted_mean_huber()
• weighted_total_huber()
• weighted_mean_tukey()
• weighted_total_tukey()

The functions with postfix _tukey are M-estimators with the Tukey biweight ψ-function. The Huber RHT M-estimator of the payroll total can be computed with

> weighted_total_huber(payroll, weight, k = 8, type = "rht")
[1] 15587090084

Note that we must specify type = "rht" for the RHT [the case type = "rhj" is discussed in the vignette “Basic Robust Estimators”]. Here, we have chosen the robustness tuning constant k = 8.

2.2 Survey methods

The following survey method are available;

• svymean_huber()
• svytotal_huber()
• svymean_tukey()
• svytotal_tukey()

The survey method of the RHT (and its standard error) is

> m <- svytotal_huber(~payroll, dn, k = 8, type = "rht")
> m
            total        SE
payroll 1.559e+10 1.182e+09

The summary() method summarizes the most important facts about the estimate.

> summary(m)
Huber M-estimator (type = rht) of the population total 

            total        SE
payroll 1.559e+10 1.182e+09

Robustness:
  Psi-function: with k = 8 
  mean of robustness weights: 0.9917 

Algorithm performance:
  converged in 4 iterations
  with residual scale (weighted MAD): 89474

Sampling design:
Stratified Independent Sampling design
svydesign(ids = ~ID, strata = ~strat, fpc = ~fpc, weights = ~weight, 
    data = workplace, calibrate.formula = ~-1 + strat)

The estimated location, variance, and standard error can be extracted from object m with the following commands.

> coef(m)
    payroll 
15587090084 
> vcov(m)
            Variance
payroll 1.398281e+18
> SE(m)
[1] 1182489162

For M-estimators, the estimated scale (weighted MAD) can be extracted with the scale() function.

> scale(m)
[1] 89474.01

Additional utility functions are:

• residuals() to extract the residuals
• fitted() to extract the fitted values under the model in use
• robweights() to extract the robustness weights

In the following figure, the robustness weights of object m are plotted against the residuals. The Huber RHT M-estimator downweights observations at both ends of the residuals’ distribution.

> plot(residuals(m), robweights(m))

2.3 Adaptive estimation

An adaptive M-estimator of the total (or mean) is defined by letting the data chose the tuning constant k. This approach is available for the RHT estimator → see vignette “Basic Robust Estimators”, Chap. 5.3 on M-estimators.