Vignette: Basic Robust Estimators

Outline

In this vignette, we discuss how to use robust estimators of location (and scale). The estimators are organized by i) estimating method, 2) population characteristic and 3) type of implemented function.

Preparations

First, we load the namespace of the robsurvey package and attach it to the search path.

> library("robsurvey", quietly = TRUE)

The argument quietly = TRUE suppresses the start-up message in the call of library("robsurvey").

1 LOS (Length-of-stay) Hospital Data

The losdata are a simple random sample without replacement of n = 70 patients from a (fictive) hospital population of N = 2 479 patients in inpatient treatment.^{Note 2} Next, we attach the data to the search path.

> data("losdata")
> attach(losdata)

The first 3 rows of the data are:

> head(losdata, 3)
  los   weight  fpc
1  10 34.91549 2479
2   7 34.91549 2479
3  21 34.91549 2479

where

• los length-of-stay in hospital (days)
• weight sampling weight
• fpc population size (finite population correction)

We consider estimating average length-of-stay in hospital (LOS, days).

> library("survey")

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

2 Trimming

> weighted_mean_trimmed(los, weight, LB = 0, UB = 0.95)
[1] 9.323529

> m <- weighted_mean_trimmed(los, weight, LB = 0, UB = 0.95, info = TRUE)
> names(m)
[1] "characteristic" "estimator"      "estimate"       "variance"      
[5] "residuals"      "model"          "design"         "call"          

> m <- svymean_trimmed(~los, dn, LB = 0, UB = 0.95)
> m
     mean    SE
los 9.324 1.064

> coef(m)
     los 
9.323529 
> vcov(m)
    Variance
los 1.131988
> SE(m)
[1] 1.063949

> summary(m)
Weighted trimmed estimator (0, 0.95) of the population mean 

     mean    SE
los 9.324 1.064

Sampling design:
Independent Sampling design
svydesign(ids = ~1, fpc = ~fpc, weights = ~weight, data = losdata, 
    calibrate.formula = ~1)

3 Winsorization

> weighted_mean_winsorized(los, weight, LB = 0, UB = 0.95)
[1] 10.40845

> weighted_mean_k_winsorized(los, weight, k = 1)
[1] 11.40845

4 Weight Reduction Methods

Winsorization and trimming act directly on the values of an estimator. Other estimation methods reduce the sampling weight of potential outliers instead. A hybrid method of winsorization and weight reduction to treat influential observations has been proposed by Dalén (1987). An observation y_i is called influential if its expanded value, w_iy_i, is exceedingly large. Let c > 0 denote a winsorization or censoring cutoff value. Dalén’s estimator "Z2" and "Z3" of the population y-total are given by ∑_i ∈ s[w_iy_i]_∘^c, where ∘ is a placeholder for "Z2" or "Z3" and $$ \begin{align*} [w_i y_i]_{Z2}^c = \begin{cases} w_i y_i & \text{if} \quad w_i y_i \leq c, \\ c & \text{otherwise}, \end{cases} &\qquad \text{and} \qquad [w_i y_i]_{Z3}^c = \begin{cases} w_i y_i & \text{if} \quad w_i y_ \leq c, \\ c + (y_i - c/w_i) & \text{otherwise}. \end{cases} \end{align*} $$

Estimator "Z2" censors the terms w_iy_i at c. In estimator "Z3", observations y_i such that w_iy_i > c contribute to the estimated total only with c plus the excess over the cutoff, (w_iy_i − c). Note that the excess over the threshold has a weight of 1.0 (Lee, 1995). An estimator of the population y-mean obtains by dividing the estimator of the estimated y-total by the (estimated) population size.

From a practical point of view, the choice of constant c in Dalén’s estimators is rather tricky because we cannot only derive c from a large order statistic, say y_(k), k < n (like for trimming). Instead, the corresponding weight w_(k) needs to be taken into account.

Good to know.

It is helpful to plot w_iy_i (weight times los) against y_i (los). The censoring constant c = 1500 (see dotted horizontal line) is such that the two largest (w_iy_i)’s are censored to 1 500.

> weighted_mean_dalen(los, weight, censoring = 1500)
2 of 71 observations censored
[1] 10.73129

5 M-Estimation

> weighted_mean_huber(los, weight, type = "rwm", k = 8)
[1] 11.17228

> huber2(los, weight, k = 8)
[1] 13.02817

> m <- svymean_huber(~los, dn, type = "rwm", k = 8)
> m
     mean    SE
los 11.17 1.328

> summary(m)
Huber M-estimator (type = rwm) of the population mean 

     mean    SE
los 11.17 1.328

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

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

Sampling design:
Independent Sampling design
svydesign(ids = ~1, fpc = ~fpc, weights = ~weight, data = losdata, 
    calibrate.formula = ~1)

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

> m <- svymean_huber(~los, dn, type = "rwm", k = 8)

> mer(m)
Search interval: [1, 10]
Minimum found for k =  10 
Rel. efficiency gain: 35%
     mean    SE
los 11.46 1.478

6 Utility Functions

> weighted_quantile(los, weight, probs = c(0.1, 0.9))
10% 90% 
  3  22 
> weighted_median(los, weight)
50% 
  8 

> weighted_mad(los, weight)
[1] 5.930408

> weighted_IQR(los, weight)
[1] 6.6717

Notes

¹ All bare-bone methods can be called with the argument info = TRUE to return a list with the following entries: characteristic (e.g., mean), estimator (e.g., trimmed estimator), estimate (numerical value), variance (by default: NA), robust (list of arguments that specify robustness), residuals (numerical vector), model (list of data used for estimation), design (by default: NA), call.

² We have constructed the losdata as a showcase; though, the LOS measurements are real data that we have taken from the 201 observations in Ruffieux et al. (2000). Our losdata are a sample of size n = 70 from the 201 original observations.

³ The function huber2() is is similar to MASS::hubers() (Venables and Ripley, 2002). It differs from the implementation in MASS in that it allows for weights and is initialized by the (normalized) weighted interquartile range (IQR) not the median absolute deviations (MAD).

Bibliographical notes

The paper of Chambers (1986) is the landmark paper about outliers in finite population sampling. Lee (1995) and Beaumont and Rivest (2009) are a good starting point to learn about robustness in finite population sampling.

Trimming and winsorization are discussed in Lee (1995) and Beaumont and Rivest (2009). The variance estimators of the weighted trimmed and winsorized estimators are straightforward adaptions of the classical estimators; see Huber and Ronchetti (2009, Chap. 3.3) or Serfling (1980, Chap. 8). A rigorous treatment in the context of finite population sampling can be found in Shao (1994).

Rao (1971) was among the first to propose weight reduction. Consider a sample of size n, and suppose that the ith observation is an outlier. He suggested to reduce the outlier’s sampling weight w_i to one, and redistribute the weight difference w_i − 1 among the remaining observations. As a result, observation i does not represent other values like it. Dalén’s estimator offers a more general notion of weight reducution; see Dalén (1987) and also Chen et al. (2017).

In the context of finite population sampling, M-estimators were first studied by Chambers (1986). He investigated robust methods in the model- or prediction based framework of Royall and Cumberland (1981). Model-assisted estimators were introduced (for ratio estimation) by Gwet and Rivest (1992) and studied by Lee (1995), and Hulliger (1995, 1999, 2005). A recent comprehensive treatment can be found in Beaumont and Rivest (2009).

References

BEAUMONT, J.-F. AND RIVEST, L.-P. (2009). Dealing with outliers in survey data. In: Sample Surveys: Theory, Methods and Inference ed. by Pfeffermann, D. and Rao, C. R. Volume 29A of Handbook of Statistics, Amsterdam: Elsevier, Chap. 11, 247–280. DOI: 10.1016/S0169-7161(08)00011-4

CHAMBERS, R. (1986). Outlier Robust Finite Population Estimation. Journal of the American Statistical Association 81, 1063–1069, DOI: 10.1080/01621459.1986.10478374.

DALEN, J. (1987). Practical Estimators of a Population Total Which Reduce the Impact of Large Observations. Research report, Statistics Sweden, Stockholm.

CHEN, Q., ELLIOTT, M. R., HAZIZA, D., YANG, Y., GHOSH, M., LITTLE, R. J. A., SEDRANSK, J. AND THOMPSON, M. (2017). Approaches to Improving Survey-Weighted Estimates. Statistical Science 32, 227–248, DOI: 10.1214/17-STS609.

GWET, J.-P. AND RIVEST, L.-P. (1992). Outlier Resistant Alternatives to the Ratio Estimator. Journal of the American Statistical Association 87, 1174–1182, DOI: 10.1080/01621459.1992. 10476275.

HUBER, P. J. AND RONCHETTI, E. (2009). Robust Statistics, New York: John Wiley & Sons, 2nd edition, DOI: 10.1002/9780470434697.

HULLIGER, B. (2006). Horvitz–Thompson Estimators, Robustified. In: Encyclopedia of Statistical Sciences ed. by Kotz, S. Volume 5, Hoboken (NJ): John Wiley & Sons, 2nd edition, DOI: 10.1002/0471667196.ess1066.pub2.

HULLIGER, B. (1999). Simple and robust estimators for sampling. In: Proceedings of the Survey Research Methods Section, American Statistical Association, 54–63, American Statistical Association.

HULLIGER, B. (1995). Outlier Robust Horvitz–Thompson Estimators. Survey Methodology 21, 79–87.

LEE, H. (1995). Outliers in business surveys. In: Business survey methods ed. by Cox, B. G., Binder, D. A., Chinnappa, B. N., Christianson, A., Colledge, M. J. and Kott, P. S. New York: John Wiley & Sons, Chap. 26, 503–526, DOI: 10.1002/9781118150504.ch26.

LUMLEY, T. (2021). survey: analysis of complex survey samples. R package version 4.0, URL https://CRAN.R-project.org/package=survey.

LUMLEY, T. (2010). Complex Surveys: A Guide to Analysis Using R: A Guide to Analysis Using R, Hoboken (NJ): John Wiley & Sons.

RAO, J. N. K. (1971). Some Aspects of Statistical Inference in Problems of Sampling from Finite Populations. In: Foundations of Statistical Inference ed. by Godambe, V. P. and Sprott, D. A. Toronto: Holt, Rinehart, and Winston, 171–202.

ROYALL, R. M. AND CUMBERLAND, W. G. (1981). An Empirical Study of the Ratio Estimator and Estimators of its Variance. Journal of the American Statistical Association 76, 66–82, DOI: 10.2307/2287043.

RUFFIEUX, C., PACCAUD, F. AND MARAZZI, A. (2000). Comparing rules for truncating hospital length of stay. Casemix Quarterly 2, 3–11.

SHAO, J. (1994). L-Statistics in complex survey problems. The Annals of Statistics 22, 946–967, DOI: 10.1214/aos/1176325505.

SERFLING, R. J. (1980). Approximation theorems of mathematical statistics, New York: John Wiley & Sons. DOI: 10.1002/9780470316481.

VENABLES, W. N. AND RIPLEY, B. D. (2002). Modern Applied Statistics with S, New York: Springer-Verlag, 4th edition, DOI: 10.1007/978-0-387-21706-2.