Relative performance of methods based on model-assisted survey regression estimation: A simulation study
Section 2. Model-assisted estimation under probability sampling

• Previous
• Next

2.1   GREG estimators

Consider the estimation of a finite population total $t_y={\displaystyle {\sum }_{i\in U}{y}_i},$ where $U=\left\{1,\dots ,N\right\}$ is the set of units of the finite population and $y_i$ is the value of the survey variable of interest for the unit $i\in U.$ Let $s\subset U$ be a sample selected according to a sampling design $p\left(\mathrm{.}\right),$ where $p\left(s\right)$ is the probability of selecting $s.$ For $i \in U,$ let ${\pi }_i=\text{Pr}\left[i\in s\right]$ denote the first-order inclusion probabilities of the design. We assume ${\pi }_i>0$ for all $i\in U.$ Additionally, assume $d$ auxiliary variables, $x_i={\left(x_{i1},x_{i2},\dots ,x_{id}\right)}^T$ are known for each $i \in U.$ A standard approach is to use the Horvitz-Thompson estimator

$${\widehat{t}}_{y,\text{HT}}={\displaystyle \sum _{i\in s}\frac{y_i}{{\pi }_i}}={\displaystyle \sum _{i\in s}{d}_i{y}_i}$$

where $d_i={\pi }_i^{-1}$ denotes design weights. Under this strictly design-based framework, the auxiliary data do not impact the form of the estimator but can impact the design weights, $d_i,$ through the specification of the sampling design.

One strategy to use auxiliary data in estimation is to employ a model-assisted estimator of $t_y$ by specifying a working model for the mean of $y$ given $x$ and use this model to predict $y$ values. Specifying a linear regression working model leads to the generalized regression (GREG) estimator (Cassel, Särndal and Wretman, 1976). The GREG estimator typically has smaller variance than the Horvitz-Thompson estimator if the working model has some predictor power for $y.$ Here, we consider the GREG estimator under a linear regression working model

$$y_i=x_i^T\beta +{\epsilon }_i(2.1)$$

with $\beta ={\left({\beta }_0,{\beta }_1,\dots ,{\beta }_p\right)}^T,$ ${\epsilon }_i$ independent and identically distributed with mean zero and variance ${\sigma }^2$ and $x_i={\left(1, x_{i1},\dots ,x_{ip}\right)}^T.$ The GREG estimator is given by

$${\widehat{t}}_{y,\text{GREG}}={\displaystyle \sum _{i\in s}\frac{y_i-x_i^T{\widehat{\beta }}_s}{{\pi }_i}}+{\displaystyle \sum _{i\in U}{x}_i^T{\widehat{\beta }}_s}(2.2)$$

with the regression coefficients $\beta$ estimated as

$${\widehat{\beta }}_s=\underset{\beta }{\text{argmin}}{\left(Y_s-X_s\beta \right)}^T{\Pi }_s^{-1}\left(Y_s-X_s\beta \right)={\left(X_s^T{\Pi }_s^{-1}{X}_s\right)}^{-1}{X}_s^T{\Pi }_s^{-1}{Y}_s,(2.3)$$                                                                             

where $X_s$ is a $n\times \left(p+1\right)$ matrix, $Y_s$ is a $n$ -vector and ${\text{Π}}_s$ is an $n\times n$ diagonal matrix of first-order inclusion probabilities for the sampled units.

The GREG estimator can also be written as a weighted sum of the variable of interest, $y,$ yielding regression weights that are independent of $y$ and, therefore, can be applied to any study variable, $y:$

$${\widehat{t}}_{y,\text{GREG}}={\displaystyle \sum _{i\in s}\left[1+{\left(t_x-{\widehat{t}}_{x\text{,}\text{HT}}\right)}^T{\left({\displaystyle \sum _{k\in s}{x}_k{x}_k^T{d}_k}\right)}^{-1}{x}_i\right]d_i{y}_i}={\displaystyle \sum _{i\in s}{w}_i{y}_i}(2.4)$$

where $t_x$ is the known population total vector of the covariates $x$ and ${\widehat{t}}_{x,\text{HT}}$ is the Horvitz-Thompson estimator vector of the covariate population totals $t_x={\displaystyle {\sum }_{i\in U}{x}_i}.$ The regression weights, $w_i,$ are termed calibration weights because they satisfy the calibration constraint ${\displaystyle {\sum }_{i\in s}{w}_i{x}_i}={\displaystyle {\sum }_{i\in U}{x}_i}.$ The calibration weight $w_i$ does not depend on the study variable $y_i.$ Note that the GREG estimator (2.4) can alternatively be expressed as

$${\widehat{t}}_{y,\text{GREG}}={\widehat{t}}_{y,\text{HT}}+{\left(t_x-{\widehat{t}}_{x,\text{HT}}\right)}^T{\widehat{\beta }}_s$$

which only requires known population totals $t_x.$ For the GREG estimator, the individual population values $x_i, i\in U$ are not needed.

If a variable selection procedure, such as a forward stepwise procedure, is implemented prior to fitting the linear regression model, then the calibration weights will depend on $y$ as the selected models may vary across study variables. This type of stepwise survey regression estimator is calibrated to the auxiliary variables selected by the variable selection procedure for a specific variable of interest, $y.$

Using a working linear regression model with many auxiliary variables, including interactions of categorical auxiliary variables, can produce substantially variable weights, and greatly increase the variance of the GREG estimator. Furthermore, some of the regression weights, $w_i,$ $i\in s,$ may be negative, thus losing the interpretation of a weight as the number of population units represented by the sampled unit.

2.2   Survey regression estimator with lasso

If the linear regression model in (2.1) is sparse, i.e., $p$ is large, and, say, only $p_0$ of the $p$ regression coefficients are nonzero, then the estimation of the zero coefficients in (2.3) leads to extra variation in the GREG estimator (2.2). In this case, model selection to remove extraneous variables could reduce the overall design variance of the GREG estimator, leading to more efficient estimates of finite population totals. The least absolute shrinkage and selection operator (lasso) method, developed by Tibshirani (1996), simultaneously performs model selection and coefficient estimation by shrinking some regression coefficients to zero. The lasso approach estimates coefficients by minimizing the sum of squared residuals subject to a penalty constraint on the sum of the absolute value of the regression coefficients.

McConville et al. (2017) proposed using survey-weight lasso estimated regression coefficients given by

$${\widehat{\beta }}_{s,L}=\underset{\beta }{\text{argmin}}{\left(Y_s-X_s\beta \right)}^T{\Pi }_s^{-1}\left(Y_s-X_s\beta \right) +\lambda {\displaystyle \sum _{j=1}^p\left|{\beta }_j\right|},$$

where $\lambda \ge 0.$ The lasso survey regression estimator for the total $t_y$ is then given by

$${\widehat{t}}_{y,\text{LASSO}}={\displaystyle \sum _{i\in s}\frac{y_i-x_i^T{\widehat{\beta }}_{s,L}}{{\pi }_i}}+{\displaystyle \sum _{i\in U}{x}_i^T{\widehat{\beta }}_{s,L}}.$$

The value of the penalty parameter $\lambda$ must be selected prior to obtaining the estimated coefficients. In general, this process of specifying hyperparameters prior to fitting the final model is called hyperparameter tuning. There are several potential selection criteria that can used to select the value of hyperparameters including Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC) or cross-validation. We used a version of cross-validation which incorporates the design weights in our simulation study; see McConville (2011) for discussion of the selection of the penalty parameter for survey-weighted lasso coefficient estimates.

2.3   Survey regression estimator with adaptive lasso

An issue with the use of the lasso criterion is that by shrinking the regression coefficients towards zero it yields biased estimates for regression coefficients that are far from zero. Under the adaptive lasso criterion (Zou, 2006), the coefficients in the $l_1$ penalty are weighted by the inverse of a root- $n$ consistent estimator of $\beta .$ Therefore, the bias for large coefficients tends to be smaller.

McConville et al. (2017) considered an adaptive lasso survey regression estimator

$${\widehat{t}}_{y,\text{ALASSO}}={\displaystyle \sum _{i\in s}\frac{y_i-x_i^T{\widehat{\beta }}_{s,\text{AL}}}{{\pi }_i}}+{\displaystyle \sum _{i\in U}{x}_i^T{\widehat{\beta }}_{s,\text{AL}}},$$

where

$${\widehat{\beta }}_{s,\text{AL}}=\underset{\beta }{\text{argmin}}{\left(Y_s-X_s\beta \right)}^T{\Pi }_s^{-1}\left(Y_s-X_s\beta \right) +\lambda {\displaystyle \sum _{j=1}^p\frac{\left|{\beta }_j\right|}{\left|{\widehat{\beta }}_{sj}\right|}}$$

and ${\widehat{\beta }}_s$ is given by (2.3). The reliance of the adaptive lasso method on the standard weighted linear regression coefficient estimates, ${\widehat{\beta }}_s,$ leads to a loss of efficiency in settings when $p$ is large because the estimates ${\widehat{\beta }}_s$ tend to be very unstable.

2.4   Lasso calibration estimators

The lasso and adaptive lasso methods do not produce regression weights directly, as the estimators cannot be expressed as weighted combinations of the $y$ -values. McConville et al. (2017) developed lasso survey regression weights using a model calibration approach and a ridge regression approximation. These lasso regression weights depend on the variable of interest, $y.$

The lasso calibration estimator is calculated by regressing the variable of interest, $y_i,$ on an intercept and the lasso-fitted mean function $x_i^T{\widehat{\beta }}_{s,L}.$ The lasso calibration estimator can be written in the same form as (2.4), where $x_i$ is replaced by $x_i^{*}={\left(1,x_i^T{\widehat{\beta }}_{s,L}\right)}^T\text{}:$

$${\widehat{t}}_{y,\text{CLASSO}}={\displaystyle \sum _{i\in s}\left[1+{\left(t_{x^{*}}-{\widehat{t}}_{x^{*},\text{HT}}\right)}^T{\left({\displaystyle \sum _{k\in s}{x}_k^{*}{x}_k^{*T}{d}_k}\right)}^{-1}{x}_i^{*}\right]d_i{y}_i}.(2.5)$$

Similarly, the adaptive lasso calibration estimator is given by

$${\widehat{t}}_{y,\text{CALASSO}}={\displaystyle \sum _{i\in s}\left[1+{\left(t_{x^{**}}-{\widehat{t}}_{x^{**},\text{HT}}\right)}^T{\left({\displaystyle \sum _{k\in s}{x}_k^{**}{x}_k^{**T}{d}_k}\right)}^{-1}{x}_i^{**}\right]d_i{y}_i},$$

where the lasso-fitted mean for $x_i^{*}$ in (2.5) is replaced by the adaptive lasso fit, $x_i^{**}={\left(1,x_i^T{\widehat{\beta }}_{s,\text{AL}}\right)}^T.$ The weights for the lasso calibration estimators are calibrated to the population size $N$ and to the population total of the lasso-fitted mean functions.

2.5   Regression tree estimator

The GREG estimator can also be expressed as

$${\widehat{t}}_{y,r}={\displaystyle \sum _{i\in s}\frac{y_i-{\widehat{h}}_n\left(x_i\right)}{{\pi }_i}}+{\displaystyle \sum _{i\in U}{\widehat{h}}_n\left(x_i\right)},(2.6)$$

where ${\widehat{h}}_n\left(x_i\right)$ is an estimator of the mean function of $Y_i$ given $X_i=x_i,$ $h\left(x_i\right)=E\left(Y_i|X_i=x_i\right),$ based on the sample data $\left(y_i,x_i\right),i\in s.$ As an alternative to a linear regression model, McConville and Toth (2019) proposed estimating $h\left(x\right)$ with a regression tree model using the following algorithm:

• Let $k\left(n\right)$ be the minimum box size and $\alpha$ be a specified significance level.
• If the dataset contains at least $2k\left(n\right)$ observations then continue to step 3; otherwise, stop.
• Among the auxiliary variables $x_l, l=1,\dots ,d,$ choose a variable to split the data. The chosen $x_l$ is the variable that shows the largest significance difference after testing the null-hypothesis of homogeneous $E\left[y|x_l\right].$ If no variable leads to a significant difference, then stop.
• Split the data into two sets $S_L$ and $S_R$ by splitting based on the value of the selected variable $x_l$ that results in the largest decrease in the estimated mean square error, while satisfying the requirement that each subset contains at least $k\left(n\right)$ units.
• For each of the resulting subsets of the data, return to step 1.

The resulting regression tree model groups the categories of an auxiliary variable based on their relationship to the variable of interest and only includes auxiliary variables and interactions associated with this variable. Importantly, including a categorical variable does not require a split for each category, potentially reducing the model size substantially while still capturing important interactions.

After fitting a regression tree model, we obtain a set of boxes $Q_n=\left\{B_{n1},B_{n2},\dots ,B_{nq}\right\}$ which partition the data. Let $I\left(x_i\in B_{nk}\right)=1$ if $x_i\in B_{nk}$ and 0 otherwise, for $k=1,\dots ,q.$ This means that $I\left(x_i\in B_{nk}\right)=1$ for exactly one box $B_{nk}\in Q_n$ for every $i\in s.$ For every $x_i\in B_{nk},$ the estimator of $h\left(x_i\right)$ is given by

$${\tilde{h}}_n\left(x_i\right)=\widetilde{\#}{\left(B_{nk}\right)}^{-1}{\displaystyle \sum _{i\in s}{\pi }_i^{-1}{y}_{i}I\left(x_i\in B_{nk}\right)}={\tilde{\mu }}_{nk},(2.7)$$

where

$$\widetilde{\#}\left(B_{nk}\right)={\displaystyle \sum }_{i\in s}{\pi }_i^{-1}I\left(x_i\in B_{nk}\right)$$

is the HT estimator of the population size in box $B_{nk}.$ The regression tree estimator ${\widehat{t}}_{y,\text{TREE}}$ is obtained by inserting equation (2.7) into the generalized regression estimator, given in equation (2.6), leading to the post stratified estimator

$${\widehat{t}}_{y,\text{TREE}}={\displaystyle \sum _k{N}_k{\tilde{\mu }}_{nk}},$$

where $N_k$ is the number of units in $U$ that belong to box $k.$

Since ${\tilde{h}}_n\left(x_i\right)$ can be written as a linear regression estimator with $q$ indicator function covariates, the regression tree estimator is also a post-stratified estimator, where each box $B_{nk}$ represents a post-stratum. This implies that this estimator is calibrated to the population total of each box, providing a data-driven mechanism, dependent on $y,$ for selecting post-strata that ensures that none of them are empty. As a result, the regression weights are guaranteed to be non-negative. The weights produced by this estimation procedure depend on the variable of interest, $y.$ Therefore, unlike the GREG approach, a single set of generic weights to apply to all study variables is not available. Instead, a set of weights for each survey variable of interest is produced.

2.6   Variance estimation under stratified simple random sampling

Under stratified simple random sampling, a variance estimator of the model-assisted survey regression estimators described above is obtained by the Taylor linearization method and given by

$$\widehat{V}\left({\widehat{t}}_y\right) ={\displaystyle \sum _h\frac{N_h\left(N_h-n_h\right)}{n_h}\frac{1}{n_h-1}{\displaystyle \sum _{i\in s_h}{\left(e_{hi}-{\overline{e}}_h\right)}^2}},(2.8)$$

where $h$ indexes the strata, $N_h$ is the number of population units in stratum $h,$ $n_h$ is the number of sampled units $s_h$ in stratum $h,$ $e_{hi}=y_{hi}-{\widehat{h}}_n\left(x_{hi}\right)$ is the residual of sample unit $i$ in stratum $h$ under the regression model and ${\overline{e}}_h$ is the average residual in stratum $h.$

The variance estimators readily extend to more complex sampling designs, but for simplicity we have given the expression only for stratified simple random sampling which is used in the simulation study of Section 3.

• Previous
• Article main page
• Next