Vignette: Robust Generalized Regression (GREG) and Ratio Prediction/ Estimation

2.1 Ratio predictor

Under the model, the least squares (LS) census estimator of θ is θ_N = ∑_i ∈ Uy_i/∑_i ∈ Ux_i. The sample weighted LS estimator is θ̂_n = ∑_i ∈ sw_iy_i/∑_i ∈ sw_ix_i, where w_i denotes the sampling weight. It is computed by

> rat <- svyratio_huber(~RMT85, ~P85, dn, k = Inf)
> rat

Survey ratio M-estimator (Huber psi, k = Inf)

Call:
svyratio_huber(numerator = ~RMT85, denominator = ~P85, design = dn, 
    k = Inf)

IRWLS converged in 1 iterations

Coefficients:
RMT85/P85  
    8.379  

Scale estimate: 4.891 (weighted MAD) 

Next, we predict the payroll total based on the estimated regression parameter θ̂_n (i.e., object rat) and the known population size in 1985 (total = 8339).

> tot <- svytotal_ratio(rat, total = 8339)
> tot
      total   SE
RMT85 69872 3214

The estimated mean square error of the ratio predictor is

> mse(tot) / 1e6
[1] 10.33016

which is equal to the estimated variance because the ratio predictor of the total is unbiased for the population total. The estimated total, standard error and variance covariance can be extracted by the functions coef, SE(), and vcov().

2.2 Robust ratio predictor

A robust ratio predictor (with Huber ψ-function and robustness tuning constant k = 20) can be computed by

> rat_rob <- svyratio_huber(~RMT85, ~P85, dn, k = 20)
> tot_rob <- svytotal_ratio(rat_rob, total = 8339)
> tot_rob
      total   SE
RMT85 68529 2550

In terms of the estimated standard error, this predictor is considerably more efficient than the (non-robust) ratio predictor. We come to the same conclusion if we consider the approximate mean square error

> mse(tot_rob) / 1e6
[1] 8.307178

In place of the Huber estimator, we can use the ratio M-estimator with Tukey biweight (bisquare) ψ-function; see svyratio_tukey(). The ratio estimator of the population of the mean is computed by svymean_ratio().

3.1 GREG

Consider the population regression model RMT85_i = θ₀ + P85_i ⋅ θ₁ + SS82_i ⋅ θ₂ + σE_i,   i ∈ U,

where CS82 is the number of Social-Democrat seats in municipal council. The E_i are independent and identically distributed random variables with zero mean and unit variance. The parameters θ₀, …, θ₂ and σ > 0 are unknown.

The weighted LS estimate is computed by

> wls <- svyreg(RMT85 ~ P85 + SS82, dn)
> wls

Weighted least squares

Call:
svyreg(formula = RMT85 ~ P85 + SS82, design = dn)

Coefficients:
(Intercept)          P85         SS82  
     61.756       11.590       -7.344  

Scale estimate: 161.2  

The summary method shows that variable CS82 contributes significantly to the explanation of the response variable.

> summary(wls)

Call:
svyreg(formula = RMT85 ~ P85 + SS82, design = dn)

Residuals:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-1066.08   -84.13   -25.84    -4.62    22.43  2001.02 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   61.756     36.108   1.710   0.0883 .  
P85           11.590      1.211   9.572   <2e-16 ***
SS82          -7.344      2.946  -2.493   0.0132 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 161.2 on 281 degrees of freedom

The diagnostic plots (below) indicate several issues. In particular, the fitted values are smaller than the responses (see “Response vs. Fitted values”). As a result, the fit tends to underestimate.

> plot(wls)

Letting the issues aside, we want to predict the 1985 revenues total. The population totals of P85 and SS82 are known quantities (from registers) and are passed to the function svytotal_reg() via the totals argument. Because the model includes a regression intercept, we must also specify the population size N. The total is predicted using

> tot <- svytotal_reg(wls, totals = c(P85 = 8339, SS82 = 6301), N = 284,
+                     type = "ADU")
> tot
      total   SE
RMT85 67914 2588

where type = "ADU" defines the “standard” GREG predictor, which is an asymptotically unbiased (ADU) estimator/ predictor, hence the name. By default, the argument check.names is set to TRUE in the call of svytotal_reg(). Thus, the names of arguments of totals are checked against the names of the estimated regression coefficients. If the names of totals are not specified, we call the function with check.names = FALSE. The estimated total, standard error and variance covariance can be extracted by the functions coef, SE(), and vcov().

The (estimated) approximate mean square error (MSE; which coincides with the estimated variance of the predictor) is

> mse(tot) / 1e6
[1] 6.699817

The estimated total, standard error and variance covariance can be extracted by the functions coef, SE(), and vcov().

3.2 Robust GREG

Consider the regression model from the last paragraph. Now, we compute a robust GREG predictor of the 1985 revenues total. The regression M-estimator with Tukey biweight (bisquare) ψ-function and the robustness tuning constant k = 15 is

> rob <- svyreg_tukeyM(RMT85 ~ P85 + SS82, dn, k = 15)
> rob

Survey regression M-estimator (Tukey psi, k = 15)

Call:
svyreg_tukeyM(formula = RMT85 ~ P85 + SS82, design = dn, k = 15)

IRWLS converged in 6 iterations

Coefficients:
(Intercept)          P85         SS82  
   -13.6992       8.3275      -0.2038  

Scale estimate: 16.46 (weighted MAD) 

The diagnostic plots look better than for the weighted LS estimator.

> plot(rob)

The robust GREG predictor of the 1985 revenues total is

> tot <- svytotal_reg(rob, totals = c(P85 = 8339, SS82 = 6301), N = 284,
+                     type = "huber", k = 50)
> tot
      total   SE
RMT85 66641 1300

where the prediction is based on the Huber ψ-function with tuning constant k = 50. The tuning constant for robust prediction should in general be larger than the one used for regression estimation. Observe that we have “mixed” the ψ-functions: Regression estimation based on the Tukey biweight ψ-function and prediction with the Huber ψ-function.

The (estimated) approximate MSE of the robust GREG predictor is

> mse(tot) / 1e6
[1] 3.308387

which is considerably smaller than the MSE of the “standard” GREG. The estimated total, standard error and variance covariance can be extracted by the functions coef, SE(), and vcov().

See the help file of svymean_reg() or svytotal_reg() to learn more.