Vignette: Robust Survey Regression

1.2 Goal

The goal is to regress the county-specific size of the farm population in 1990 (variable farmpop) on a set of explanatory variables. The simplest population regression model is

$$\begin{equation} \xi: \quad \mathrm{farmpop}_i = b_0 + b_1 \cdot \mathrm{numfarm}_i + \sqrt{v_i} E_i, \qquad i \in U, \end{equation}$$

where U is the set of labels of all 3 141 counties in the U.S. (population), b₀, b₁ ∈ ℝ are unknown regression coefficients, v_i are known constants (i.e., known up to a constant multiplicative factor, σ), and E_i are regression errors (random variables) with mean zero and unit variance such that E_i and E_j are conditionally independent given numfarm_i, i ∈ U.

Another candidate model is

$$\begin{equation} \xi_{log}: \quad \log(\mathrm{farmpop}_i) = b_0 + b_1 \cdot \log(\mathrm{numfarm}_i) + \sqrt{v_i} E_i, \qquad i \in U. \end{equation}$$

1.3 Sampling design

The counties dataset is of size n = 100. In what follows, we restrict attention to the subset of the 98 counties with farmpop_i > 0. Hence, we define the sampling design dn.

> dn <- svydesign(ids = ~1, fpc = ~fpc, weights = ~weights,
+                 data = counties[counties$farmpop > 0, ],
+                 calibrate.formula = ~1)

1.4 Exploring the data

The figures below show scatterplots of farmpop plotted against numfarm (in raw and log scale). The scatterplot in the left panel indicates that the variance of farmpop increases with larger values of numfarm (heteroscedasticity). The same graph but in log-log scale (see right panel) shows an outlier, which is located far from the bulk of the data.

2.1 Heteroscedasticity

In the above scatterplot, we observe that the variance of farmpop increases with larger values of numfarm (heteroscedasticity). We conjecture that the v_i’s in model ξ are proportional to the square root of the variable numfarm. Thus, we add variable vi to the design object dn with the update() command

> dn <- update(dn, vi = sqrt(numfarm))

The argument var of svyreg() is now used to specify the heteroscedastic variance (it can be defined as a formula, i.e. ~vi, or the variable name in quotation marks, "vi").

> svyreg(farmpop ~ -1 + numfarm, dn, var = ~vi, na.rm = TRUE)

Weighted least squares

Call:
svyreg(formula = farmpop ~ -1 + numfarm, design = dn, var = ~vi, 
    na.rm = TRUE)

Coefficients:
numfarm  
  1.775  

Scale estimate: 99.65  

Note that we have dropped the regression intercept.

2.2 Regression M-estimator

We continue to assume the heteroscedastic model ξ, but now we compute a weighed regression M-estimator. Two types of estimators are available:

• svyreg_huberM() M-estimator with Huber or asymmetric Huber ψ-function (the latter obtains by specifying argument asym = TRUE);
• svyreg_tukeyM() M-estimator with Tukey biweight (bisquare) ψ-function.

The tuning constant of both ψ-functions is called k. The Huber M-estimator of regression with k = 3 is (note that we use na.rm = TRUE to remove observations with missing values)

> m <- svyreg_huberM(farmpop ~ -1 + numfarm, dn, var = ~vi, k = 1.3,
+                    na.rm = TRUE)
> m

Survey regression M-estimator (Huber psi, k = 1.3)

Call:
svyreg_huberM(formula = farmpop ~ -1 + numfarm, design = dn, 
    k = 1.3, var = ~vi, na.rm = TRUE)

IRWLS converged in 6 iterations

Coefficients:
numfarm  
  1.669  

Scale estimate: 83.35 (weighted MAD) 

and the summary obtains by

> summary(m)

Call:
svyreg_huberM(formula = farmpop ~ -1 + numfarm, design = dn, 
    k = 1.3, var = ~vi, na.rm = TRUE)

Residuals:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-165.951  -64.126  -10.095    6.848   47.253  458.062 

Coefficients:
        Estimate Std. Error t value Pr(>|t|)    
numfarm  1.66859    0.09758    17.1   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 83.35 on 3077.18 degrees of freedom

Robustness weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2369  1.0000  1.0000  0.9294  1.0000  1.0000 

The output of the summary method is almost identical with the one of an lm() call. The paragraph on robustness weights is a numerical summary of the M-estimator’s robustness weights, which can be extracted from the estimated model with robweights().

The subsequent figure shows the graph for plot(residuals(m), robweights(m)). From the graph, we can see by how much the residuals have been downweighted.

2.3 Regression GM-estimator

We want to fit the population regression model

$$ \begin{equation*} \log(\mathrm{farmpop}_i) = b_0 + b_1 \log(\mathrm{numfarm}_i) + b_2 \log(\mathrm{totpop}_i) + b_3 \log(\mathrm{farmacre}_i) + \sqrt{v_i} E_i, \qquad i \in U, \end{equation*} $$

with sample data by the weighted generalized M-estimator (GM) of regression. GM-estimators are robust against high leverage observations (i.e., outliers in the design space of the model) while M-estimators of regression are not. Two types of GM-estimators are available:

• svyreg_huberGM() with Huber or asymmetric Huber ψ-function (the latter obtains by specifying argument asym = TRUE);
• svyreg_tukeyGM() with Tukey biweight (bisquare) ψ-function.

> dn <- svydesign(ids = ~1, fpc = ~fpc, weights = ~weights,
+                 data = counties[counties$farmpop > 0, ])
> dn_exclude <- na.exclude(dn)

> f <- log(farmpop) ~ log(numfarm) + log(totpop) + log(farmacre)

> xmat <- model.matrix(f, dn_exclude$variables)[, -1]

> if (requireNamespace("wbacon", quietly = TRUE)) {
+     # package wbacon is available
+     mv <- wbacon::wBACON(xmat, weights = weights(dn_exclude))
+     # distances
+     dis <- wbacon::distance(mv)
+ } else {
+     # package wbacon is not available
+     center <- c(6.285968, 10.195002, 12.047715)
+     scatter <- matrix(0, 3, 3)
+     scatter[lower.tri(scatter, TRUE)] <- c(0.678646, 0.441020, 0.415634,
+                                            2.191174, -0.302097, 1.040932)
+     scatter <- scatter + t(scatter) - diag(3) * scatter
+     # distances
+     dis <- sqrt(mahalanobis(xmat, center, scatter))
+ }

> m <- svyreg_tukeyGM(f, dn_exclude, k = 4.6, xwgt = tukeyWgt(dis))
> summary(m)

Call:
svyreg_tukeyGM(formula = f, design = dn_exclude, k = 4.6, xwgt = tukeyWgt(dis))

Residuals:
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-1.004303 -0.352241 -0.052866  0.004878  0.360099  3.389738 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    1.62725    0.56957   2.857  0.00431 ** 
log(numfarm)   1.30079    0.07044  18.466  < 2e-16 ***
log(totpop)   -0.06550    0.03176  -2.062  0.03929 *  
log(farmacre) -0.20161    0.04617  -4.366 1.31e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5549 on 2979.95 degrees of freedom

Robustness weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.6270  0.7850  0.7114  0.8818  0.9959 

4.1 The MU284 population

The MU284 population of Särndal et al. (1992, Appendix B) is a dataset with observations on the 284 municipalities in Sweden in the late 1970s and early 1980s. It is available in the sampling package; see Tillé and Matei (2021). The population is divided into two parts based on 1975 population size (P75):

• the MU281 population, which consists of the 281 smallest municipalities;
• the MU3 population of the three biggest municipalities/ cities in Sweden (Stockholm, Göteborg, and Malmö).

The three biggest cities take exceedingly large values (representative outliers) on almost all of the variables. To account for this (at least to some extent), a stratified sample has been drawn from the MU284 population using a take-all stratum. The sample data, MU284strat, is of size n = 60 and consists of

• a stratified simple random sample (without replacement) from the MU281 population, where stratification is by geographic region (REG) with proportional sample size allocation;
• a take-all stratum that includes the three biggest cities/ municipalities (population M3).

The overall sampling fraction is 60/284 ≈ 21% (and thus not negligible). The data frame MU284strat includes the following variables.

We load the data and generate the sampling design.

> data(MU284strat)
> dn <- svydesign(ids = ~1, strata = ~ Stratum, fpc = ~fpc, weights = ~weights,
+                 data = MU284strat, calibrate.formula = ~-1 + Stratum)

4.2 Model fit and statistical inference

The Huber M-estimator of regression is

> m <- svyreg_huberM(RMT85 ~ REV84 + P85 + S82 + CS82, dn, k = 2)

The compound design- and model-based distribution of the estimated coefficients obtains by specifying the argument mode = "compound" in the call of the summary() method.

> summary(m, mode = "compound")

Call:
svyreg_huberM(formula = RMT85 ~ REV84 + P85 + S82 + CS82, design = dn, 
    k = 2)

Residuals:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
 -77.877  -18.469    5.044   66.250   17.635 2690.755 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 130.973688 106.868749   1.226    0.221    
REV84         0.005729   0.023422   0.245    0.807    
P85           9.499801   2.124650   4.471 1.13e-05 ***
S82          -3.845513   2.784454  -1.381    0.168    
CS82         -1.965082   5.292654  -0.371    0.711    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 23.96 on 279 degrees of freedom

Robustness weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.01781 1.00000 1.00000 0.95076 1.00000 1.00000 

We may call vcov(m, mode = "compound") to extract the estimated covariance matrix of the regression estimator under the compound design- and model-based distribution.