Bayes, buttressed by design-based ideas, is the best overarching paradigm for sample survey inference
Section 4. Examples

• Previous
• Next

A variety of examples are offered, admittedly somewhat slanted towards my own work with colleagues. I avoid topics such as small area or time series estimation, because the need for modeling is well established there.

Example 1 continued. Normal random sample. A common critique of model-based methods is “if I base an inference on a model and the model is wrong, the inference must be wrong. Because (to paraphrase George Box), all models are wrong, therefore all model-based inference is wrong. So I prefer design-based inference, which does not require modeling assumptions”. The reasoning is plausible but it’s not that simple. The validity of model-based methods depends both on the design and on the degree and nature of model misspecification; the fact that design-based methods do not overtly depend on models does not mean they are necessarily superior to methods that do.

Suppose a simple random sample of size $n$ is taken from a continuous distribution with unknown mean $\mu$ and standard deviation $\sigma .$ Consider three interval estimates for the population mean:

• (i) Interval A is the standard design-based 95% confidence interval, equation (3.1).
• (ii) Interval B replaces the 97.5^th normal percentile in equation (3.1), namely 1.96, by the 97.5^th percentile of the $t$ distribution with $n-1$ df.
• (iii) Interval C is the 95% posterior credible interval based on a normal model with the Jeffreys’ prior distribution $p\left(\mu ,\sigma \right)\propto 1/\sigma .$

From a frequentist perspective, which of these intervals is the best? Interval A makes no distributional assumption on $Y,$ and B makes a normal assumption $–$ does that mean that A is superior? If $n$ is large then the intervals are essentially the same, and if $n$ is small then B is arguably better than A even if the data are not normal, because it is reflecting uncertainty about the variance.

In an informal survey, two thirds of a recent class of our well-trained Ph.D. students preferred B to C, on the grounds that it avoids the choice of prior distribution and hence makes fewer assumptions. But B and C are equally good or bad, because they are the same procedure! That this interval is an exact 95% confidence interval under normality, and is a 95% credible interval for the stated choice of prior, are just two properties of the procedure. Judging a method by its overt assumptions is an over-simplification.

Example 2. Simple random sample with a lower bound on the variance (Little, 2006). Suppose in the previous example that $n=7,$ $\overline{y}=1,s=1.$ The standard $t$ interval (ignoring finite population corrections) is

$$\overline{y}\pm t_{6,975}s/\sqrt{n}=1\pm 2.447/\sqrt{7}=1\pm 0.925(4.1)$$

if in fact we know that $\sigma =$ 1.5, a better interval is

$$\overline{y}\pm z_{0.975}\sigma /\sqrt{n}=1\pm 1.96\times 1.5/\sqrt{7}=1\pm 1.11,(4.2)$$

the wider interval reflecting the fact that $\sigma$ is greater than the particular value of $s$ for this sample. The t interval (4.1) has exact confidence coverage, but, given what we know about $\sigma ,$ it is the wrong inference for this specific data set: we should not pick it over (4.2) because it is narrower!

Now suppose we know that $\sigma >$ 1.5, because of some unaccounted source of additional variation. The Bayesian approach then incorporates this information into the prior distribution for $\sigma ,$ resulting in a credible interval that is wider than (4.2). What is the frequentist answer? The t confidence interval still has exactly nominal coverage in repeated sampling, but it is clearly too narrow as an inference for the observed dataset, because it fails to reflect what is known about $\sigma .$ The interval (4.2) is anticonservative, despite the fact that it is wider than (4.1) for the realized data set $–$ a property that can happen for confidence intervals but cannot happen for credible intervals under a specific model. Asymptotic frequentist methods are no help here, so what is the alternative to Bayes in this example?

A related example arises in one group random-effects analysis of variance, when the least squares estimate of the between-group variance is negative $–$ a Bayesian analysis addresses this with a prior distribution on the between-group variance that does not allow negative values. Random-effects and mixed-effects models are important to handle clustering in surveys, and Bayesian methods are better than ML in this setting.

Example 3. Post-stratification on a categorical covariate. Prediction (as in modeling) is a more reliable general approach to inference than weighting (as in design-based inference). I illustrate this general statement with the simple example of post-stratification on a single variable $Z.$ A more complex example $–$ Example 7 $–$ is given later.

Consider an equal probability sample with a single categorical post-stratifying variable $Z,$ for which known population counts $N_h$ are available for each post-stratum $h,$ $h=1,2,\dots ,H.$ Let ${\overline{y}}_h$ be the sample mean in post-stratum $h,$ based on sample size $n_h,$ and $n={\displaystyle {\sum }_{h=1}^H{n}_h},$ $N={\displaystyle {\sum }_{h=1}^H{N}_h}.$ The standard estimate of the population mean is the post-stratified mean:

$${\overline{y}}_{\text{PS}}={\displaystyle \sum _{h=1}^H{P}_h{\overline{y}}_h},(4.3)$$

where $P_h=N_h/N$ and $N$ is the population size. This can be viewed as a weighted mean

$${\overline{y}}_{\text{PS}}={\displaystyle \sum _{h=1}^H{\displaystyle \sum _{i=1}^{n_h}{w}_i{y}_i}},$$

where $w_i=N_h/(N{n}_h)$ is the post-stratification weight for sampled units in poststratum $h.$ These weights can be very large in post-strata with small sample sizes $n_h,$ which, unlike stratified sampling based on $Z,$ are not under the control of the sampler. These large weights can lead to excessive variability in ${\overline{y}}_{\text{PS}}.$ In fact, strictly speaking, ${\overline{y}}_{\text{PS}}$ does not have a distribution in repeated sampling, because with positive probability the sample sizes in some post-strata may be zero. This remains true if the post-strata are modified to ensure that the post-strata sample counts are all positive for the observed sample, for example by pooling adjacent strata.

The standard design-based approach to excessive variability of ${\overline{y}}_{\text{PS}}$ is to modify the weights $\left\{w_i\right\},$ for example by trimming the large ones. However, from a prediction perspective, this is misguided. The problem is not the weights $–$ the population proportions $\left\{P_h\right\}$ in each post-stratum are known, after all $–$ the problem is that sparsity of sample in some post-strata renders the estimates $\left\{{\overline{y}}_h\right\}$ unreliable. It is the estimates in sparse post-strata that need to be modified, not the weights $\left\{w_i\right\}$ attached to sampled units. The principled way to modify $\left\{{\overline{y}}_h\right\}$ is to assume a model relating $Y$ and $Z.$ The design-based approach, by avoiding such a model, leads to the wrong principle $–$ modifying the weights rather than the predictions of non-sampled values.

A related point: a common design-based approximation (Kish, 1992) measures the proportionate increase in variance from weighting as $\text{1}\text{+}\text{cv}\left(w_i\right),$ where $\text{cv}\left(w_i\right)$ is the coefficient of variation of the weights; trimming the weights reduces $\text{cv}\left(w_i\right)$ and hence this proportionate increase. However, this rule of thumb is only valid when $Y$ and $Z$ are unrelated, in which case post-stratification is useless. If $Y$ and $Z$ are related and the sample size is not too small, the variance of the post-stratified mean is smaller, not larger, than the variance of the unweighted sample mean (Holt and Smith, 1979; Little and Vartivarian, 2005). The rule of thumb fails because the relationship between $Y$ and $Z$ is not modeled.

What is the Bayesian approach to excessive variability of ${\overline{y}}_{\text{PS}}?$ The latter is the posterior mean for the stratified normal model

$$\left(y_{hi}|{\mu }_h,{\sigma }^2\right)\stackrel{\text{iid}}{{\displaystyle ~}}G\left({\mu }_h,{\sigma }^2\right),(4.4)$$

$$p\left({\mu }_h,\sigma \right)\propto 1/\sigma ,(4.5)$$

where $G\left(a,b^2\right)$ denotes the normal (Gaussian) distribution with mean $a,$ variance $b^2,$ and equation (4.5) is the Jeffreys’ prior distribution on the mean and standard deviation in each post-stratum. Given a sparse sample in some post-strata, the prior distribution needs to be modified to allow borrowing of strength from other strata. One approach is to assume the normal random-effects model

$$p\left({\mu }_h|\mu ,\tau ,\sigma \right)\stackrel{\text{iid}}{{\displaystyle ~}}G\left(\mu ,{\tau }^2\right),p\left(\mu ,\sigma ,\tau \right)\propto {\sigma }^{-1},(4.6)$$

which treats the stratum means as random effects. The variances in each post-stratum might also be treated as distinct random effects and assigned a prior distribution, rather than pooled. The posterior mean of $\overline{Y}$ for the prior distribution (4.6) moves the weight $w_i$ of sampled units in post-stratum $h$ towards one, with a degree of shrinkage that depends on the relative size of estimates of $\sigma$ and $\tau$ (Lazzeroni and Little, 1998).

The prior distribution (4.6) makes the non-trivial assumption that the post-stratum means are exchangeable. It can be relaxed by restricting the random effects model to a subset of post-strata with small sample counts; or the constant mean $\mu$ in (4.6) might be replaced by a regression on known post-stratum characteristics $C_h,$ as in:

$$p\left({\mu }_h|{\beta }_0,{\beta }_1,\tau ,\sigma ,c_h\right)\stackrel{\text{ind}}{{\displaystyle ~}}G\left({\beta }_0+{\beta }_1{c}_h,{\tau }^2\right),p\left({\beta }_0,{\beta }_1,\sigma ,\tau \right)\propto {\sigma }^{-1},(4.7)$$

which limits the exchangeability assumption to the errors in the regression of ${\mu }_h$ on $c_h.$ For extensive generalizations of this basic example, see Gelman and Little (1997), Elliott and Little (2000), Elliott (2007), Gelman (2007) and Si, Trangucci, Gabry and Gelman (2020).

Example 4. Regression estimator of the mean, given a population auxiliary variable. If the auxiliary variable $Z$ in the previous example is continuous, a common way to incorporate it in the inference is via the regression estimate of the mean:

$$q={\overline{y}}_{\text{REG}}=\overline{y}+{\widehat{\beta }}_1\left(\overline{Z}-\overline{z}\right),$$

where ${\widehat{\beta }}_1$ is the least squares estimate of the slope of $Y$ on $Z$ in the sample, and $\overline{z}$ and $\overline{Z}$ are respectively the sample and population mean of $Z.$ In a simulation study of five real populations, Royall and Cumberland (1981, 1985) assess inferences centered at $q,$ with (a) the standard design-based standard error based on simple random sampling, namely:

where $s_{Y.Z}^2$ is the sample residual variance and $f=n/N$ is the sampling fraction; and (b) the prediction standard error based on the normal linear regression model with constant variance, namely:

The design-based confidence intervals $I_{0.95D}$ exhibit very poor conditional confidence coverage when the observed $\overline{z}$ deviates substantially from $\overline{Z}.$ The model-based confidence interval $I_{0.95M}$ takes into account this lack of balance with respect to $\overline{z},$ but is vulnerable to model misspecification, specifically lack of linearity in the relationship between $Y$ and $Z$ or non-constant residual variance. Robust estimates of standard error, based on the sandwich estimator or the jackknife, yield intervals with better conditional coverage properties, although still sometimes deviating from nominal coverage levels. An alternative approach is Bayesian inference based on a more flexible model relating $Y$ to $Z,$ such as the penalized spline model:

$$\begin{array}{l}\left(y_i|z_i,\beta ,{\sigma }^2\right)\stackrel{\text{ind}}{{\displaystyle ~}}G\left(\text{spline}\left(z_i,\beta \right),{\sigma }^2{z}_i^{\alpha }\right),\hfill \\ \text{spline}\left(z_i,\beta \right)={\beta }_0+{\displaystyle \sum _{j=1}^p{\beta }_j{z}_i^j}+{\displaystyle \sum _{\text{ℓ}=1}^m{\beta }_{\text{ℓ}+p}{\left(z_i-{\kappa }_{\text{ℓ}}\right)}_{+}^p},\hfill \\ \left({\beta }_{l+p}|\tau \right)\stackrel{\text{iid}}{{\displaystyle ~}}N\left(0,{\tau }^2\right),l=1,\dots ,m;p\left({\beta }_0,\dots ,{\beta }_p,\alpha ,\sigma ,\tau \right)\propto 1/\sigma ,0<\alpha <2,(4.8)\hfill \end{array}$$

where the constants ${\kappa }_1<\dots <{\kappa }_m$ are selected fixed knots, and ${(u)}_{+}^p=u^p$ if $u>0$ and 0, otherwise (see, for example Ruppert, Wand and Carroll, 2003). The parameter $\alpha$ allows for a variety of common forms of heteroskedasticity. The Bayesian standard errors then reflect imbalance in distribution of $Z$ in the sample and population, and the flexibility of the model limits bias from model misspecification.

Suppose that the target quantity is not the population mean of $Y,$ but the least squares slope of $Y$ on $Z$ in the population. A robust approach is to impute the non-sampled values of $Y$ using the model (4.8), and then estimate the slope of $Y$ on $Z$ as the least squares slope estimated on the filled-in population data. Uncertainty can be propagated by multiple imputation (Rubin, 1987), a method founded on Bayesian ideas. In this context, Little (2004) distinguishes between the “target model” that determines the target population quantity of interest, here the linear regression of $Y$ on $X,$ and the “working model” (4.8) that is the basis for inference, and is used to predict survey variables for the non-sampled and nonresponding units in the population. Distinguishing between these two models provides for a robust form of Bayesian survey inference.

Szpiro, Rice and Lumley (2010) apply a similar idea in a superpopulation regression setting, and Little (2019) argues that this is more straightforward than changing the interpretation of the estimand, the approach adopted by Buja, Berk, Brown, George, Pitkin, Zhan and Zhang (2019). Szpiro et al. (2010) show that the approach provides a Bayesian interpretation of the sandwich estimator of variance in regression, which is asymptotically equivalent to sample reuse estimates of variance like the bootstrap or jackknife, which are commonly applied in sample survey settings.

Example 5. Inference for samples with unequal probabilities of selection. For designs with unequal selection probabilities, classic papers critiquing the modeling approach to surveys (Kish and Frankel, 1974; Hansen, Madow and Tepping, 1983) do not include the selection probabilities in the model, yielding inferences that are vulnerable to model misspecification. The selection probabilities play an important role in robust model-based inference, but as model covariates rather than as sampling weights.

Consider, for example, inference about a population mean $\overline{Y}.$ If $y_i$ is the value of a survey variable $Y$ and ${\pi }_i$ is the selection probability for unit $i,$ the usual design-based estimator of the mean of $Y$ weights sampled units (say units $i=1,\dots ,n)$ by the inverse of ${\pi }_i$ resulting in the Horvitz-Thompson estimate (Horvitz and Thompson, 1952)

$${\overline{y}}_{\text{HT}}=N^{-1}{\displaystyle \sum _{i=1}^n{y}_i/{\pi }_i},(4.9)$$

if the population size $N$ is known, or the Hájek (1971) estimate

$${\overline{y}}_{\text{HK}}=\left({\displaystyle \sum _{i=1}^n{y}_i/{\pi }_i}\right)/\left({\displaystyle \sum _{i=1}^{n}1/{\pi }_i}\right),(4.10)$$

if $N$ is estimated. Weighting is a “one size fits all” approach, with sampled units receiving the same weight irrespective of the relationship between ${\pi }_i$ and $y_i.$ This is potentially inefficient $–$ if the relationship is weak, weighting simply reduces the precision of the estimate without a compensating reduction in bias.

The modeling approach incorporates the selection probabilities by regressing $y_i$ on ${\pi }_i.$ The strength of relationship between $y_i$ and ${\pi }_i$ then moderates how the selection probability affects the estimator $–$ if the relationship is weak, the regression coefficient of ${\pi }_i$ is small and the sampling weight has little influence. This results in more efficient estimates.

A linear regression of $y_i$ on ${\pi }_i$ is vulnerable to bias if the linearity is misspecified, but the impact of misspecification can be reduced by choosing a model that results in a design-consistent estimate. Many models satisfy this requirement; see for example Firth and Bennett (1998).

In stratified sampling with stratifying variable $Z,$ the natural regression model includes covariates that are dummy variables for the strata. The resulting estimate of the population mean of a continuous survey variable $Y$ is the stratified mean, which has the same form as equation (4.3) but with $Z$ forming strata rather than post-strata. Weighting by the inverse of the selection probability in each stratum and dummy variable regression both yield the estimator (4.3) in this situation. The Bayesian approaches to reducing variance described in example 3 tend to have less pay-off in the case of stratified sampling than for post-stratification, because the sampler has control over the sample sizes in each stratum.

In probability proportional to size (PPS) sampling, the covariate $Z$ is the size of the unit and ${\pi }_i=\min \left(c{z}_i,1\right).$ Now $Z$ is a continuous variable, and weighting and regression may yield different answers. The approaches can be unified by considering models that yield the design-weighted estimates when used to predict the non-sampled units. In particular, ignoring finite population corrections, the HT estimate (4.9) is the posterior mean for the “Horvitz Thompson model”:

$$\left(y_i|z_i,\beta ,{\sigma }^2\right)\stackrel{\text{ind}}{{\displaystyle ~}}G\left(\beta z_i,{\sigma }^2{z}_i^2\right),p\left(\beta ,\sigma \right)\propto 1/\sigma ,(4.11)$$

and the Hájek estimate (4.10) is the posterior mean for the “Hájek model”:

$$\left(y_i|z_i,\beta ,{\sigma }^2\right)\stackrel{\text{ind}}{{\displaystyle ~}}G\left(\beta ,{\sigma }^2{z}_i\right),p\left(\beta ,\sigma \right)\propto 1/\sigma .(4.12)$$

These underlying models describe situations where the corresponding design-based estimates are optimal. However, they involve strong parametric assumptions. A robust Bayesian modeling approach embeds these models within a larger model, such as the penalized spline model (4.8), as proposed by Zheng and Little (2003). Zheng and Little (2005) and Chen, Elliott, Haziza, Yang, Ghosh, Little, Sedransk and Thompson (2017) provide simulation studies suggesting that this modeling approach can yield substantial gains over HT or Hájek estimation, both in terms of efficiency and closer to nominal (frequentist) confidence coverage in moderate samples $–$ remember, design-based results are asymptotic. The model (4.8) is readily expanded to include other auxiliary variables measured for all the population units, and the flexibility of small-sample inferences increased by including proper prior distributions for the model parameters.

Leon-Novelo and Savitsky (2019) consider Bayesian models for the joint distribution of $y_i$ and ${\pi }_i,$ calling models that fix ${\pi }_i$ a “plug-in” approach. However, assigning a distribution to ${\pi }_i$ seems to me both unnatural and unnecessary in this setting; if ${\pi }_i$ is recorded for all population units, it can be conditioned in the model, as in standard regression approaches that treat covariates as fixed. Even if values of ${\pi }_i$ for non-sampled units are not provided to analysts, their values can be predicted via Bayes Theorem, as in Zangeneh and Little (2015).

Example 6. Penalized spline of response propensity models. Consider unit nonresponse on a survey variable $Y,$ when a set of variables say $X_1,\dots ,X_p$ are observed for respondents and nonrespondents in the sample. The response propensity for unit $i,$ ${\theta }_i=\mathrm{Pr}\left(R_i=1|x_{i1},\dots ,x_{ip},\psi \right),$ where $\psi$ represents unknown parameters, plays an important role in both weighting and prediction approaches to survey nonresponse. In nonresponse weighting the sampling weight is multiplied by the nonresponse weight, as in

$$\begin{array}{ll}{w}_i\hfill & =1/\text{Pr}(\text{unit}i\text{is}\text{selected}\text{and}\text{responds})\hfill \\ \hfill & =(1/\mathrm{Pr}(i\text{selected}))\times (1/\text{Pr}( i\text{responds}|\text{selected})) \hfill \\ \hfill & \text{sampling}\text{weight}\text{×}\text{response}\text{weight}. \hfill \end{array}$$

Unlike the sampling weight, the response weight is unknown, and the definition needs to be clear on what the probability is conditioned. Little (2022) argues that it should condition on auxiliary and survey variables, but not on other variables that might affect it. As with sampling weights, weighting by the inverse of the estimated response propensity is a “one size fits all” approach that does not take into account the strength of relationship between $R$ and $Y.$ Penalized Spline of Propensity Prediction (PSPP, Zhang and Little, 2009) regresses $Y$ on a penalized spline of the estimated response propensity, with other variables that are observed for respondents and nonrespondents entering parametrically. The spline models a flexible relationship between the response propensity and $Y,$ providing robustness to model misspecification. The balancing property of the propensity score also provides a double robustness property, in that the parametric form of the regression on other predictors can be misspecified without bias, provided the penalized spline captures the relationship between the propensity and $Y.$ This method performs favorably with weighting methods in simulations (Zhang and Little, 2009; Yang and Little, 2015), and the fully Bayes version with prior distributions on the parameters propagates error in estimating the propensities.

Example 7. Unit nonresponse with poststratification. Another example where Bayesian and design-based approaches to inference differ concerns unit nonresponse with post-stratification. Consider the setting of Figure 4.1, where $Y$ is a survey variable subject to nonresponse, $X$ is a variable observed for all units in the sample, and $Z^{\text{aux}}=\left(Z_1,\dots ,Z_K\right)$ consists of a set of categorical auxiliary variables. The joint distribution of $Z^{\text{aux}}$ is observed for survey respondents, and the marginal distributions of $Z_k,$ $k=1,\dots ,K,$ are also observed for the population or sample from external sources. A distinctive feature is that the units in the auxiliary data are not linked with the units in the survey. This scenario occurs frequently in settings where post-stratification is used for nonresponse adjustment. Define $R_i$ as the response indicator for $\left(y_i,z_i\right)$ for sample unit $i,$ and suppose that

$$\mathrm{Pr}\left(R_i\text{=}1|x_i,z_i,y_i,\psi \right)=\mathrm{Pr}\left(R_i\text{=}1|x_i,z_i,\psi \right),(4.13)$$

so that this probability does not depend on $y_i.$ Note that if this probability depends on $z_i,$ the response mechanism is missing not at random (MNAR) according to Rubin’s (1976) definition, because values of $z_i$ are missing for survey nonrespondents. Zangeneh and Little (2022) consider Bayes and ML estimation for data with this pattern. Considering for simplicity models that are i.i.d. over the units $i,$ the joint distribution of $X,$ $Z,$ $Y$ and $R$ is factored as

$$\begin{array}{ll}{f}_{X,Z,Y,R}\left(x_i,z_i,y_i,r_i|\theta ,\varphi \right)\hfill & =f_{Y|X,Z,R}\left(y_i|x_i,z_i,\theta ,R_i=r_i\right)f_{X,Z,R}\left(x_i,z_i,r_i|\varphi \right)\hfill \\ \hfill & \stackrel{\text{Under Eq}\text{. (4}\text{.13)}}{=}{f}_{Y|X,Z,R}\left(y_i|x_i,z_i,\theta \right)f_{X,Z,R}\left(x_i,z_i,r_i|\varphi \right),\hfill \end{array}$$

where $\theta$ and $\varphi$ are distinct parameters (Little and Rubin, 2019, Chapter 6). The parameters $\theta$ in this factorization can thus be estimated from the respondent survey data with $R_i=1,$ and parameters of the joint distribution of $X$ and $Z$ are estimated by combining the respondent and auxiliary data on those variables.

Two special cases considered by Zangeneh and Little (2022) are as follows:

1. No covariates $X$ and a single post-stratifier $Z.$ The missingness assumption in equation (4.13) then reduces to $\mathrm{Pr}\left(R_i=1|Z_i,Y_i,\psi \right)=\mathrm{Pr}\left(R_i=1|Z_i,\psi \right).$ The parameters of the conditional distribution of $Y$ given $Z$ are estimated from the survey respondents, and the parameters of the marginal distribution of $Z$ are estimated from the auxiliary data on $Z.$ In particular, if $Z$ is multinomial with $\mathrm{Pr}\left(z_i=j|\varphi \right)={\varphi }_j$ the ML estimate of ${\varphi }_j$ is simply the proportion of the auxiliary counts in post-stratum $j.$ If, in addition, $Y$ given $Z=j$ is assumed normal with mean ${\mu }_j$ and variance ${\sigma }_j^2,$ the resulting ML estimate of the population mean of $Y$ is simply the post-stratified mean given in equation (4.3). This simple example is interesting theoretically, because the response mechanism is MNAR but ignorable for likelihood inference; MAR is a sufficient but not always necessary condition for ignorability.
2. $X$ is a single categorical variable and $Z$ is a single categorical post-stratifier. Now the joint distribution of $X,$ $Z,$ and $R$ is not identified under a saturated model, and requires additional assumptions to identify the model. In particular, a constrained MNAR “RAKE” model assumes that the marginal distributions of $Z_1$ and $Z_2$ are different for respondents and nonrespondents, but the odds ratios of $Z_1$ and $Z_2$ are the same for respondents and nonrespondents. This yields a just-identified model. Raking the table of respondent counts of $Z_1$ and $Z_2$ to the auxiliary margins of $Z_1$ and $Z_2$ yields ML estimates of $\varphi$ under this RAKE model. (Little and Wu, 1991; Little, 1993). The post-stratified estimator of the mean of $Y$ is then

$${\overline{Y}}_{\text{rake}}={\displaystyle \sum _{j=1}^J{\displaystyle \sum _{k=1}^K{P}_{jk}\left(\widehat{\varphi }\right){\overline{Y}}_{jk}\left(\widehat{\theta }\right)}},$$

where $P_{jk}\left(\widehat{\varphi }\right)$ is the proportion of the population with $X=j,Z=k$ from raking the respondent counts to the margins, and ${\overline{Y}}_{jk}\left(\widehat{\theta }\right)$ is the estimated mean of $Y$ given $X=j,Z=k$ from the model for $Y$ given $X$ and $Z.$

Some comments on this approach are as follows:

1. Note again this is a MNAR model, and it is weaker than the MAR model that assumes that response depends on $X$ but not on $Z.$ As a result, the method tends to be less biased than alternative methods that assume MAR for consistent estimates.
2. As was the case in example 4, the ML approach to this model does not involve modifying the weights applied to the data in cell $X=j,Z=k$ using arbitrary distance functions $–$ the ML estimates are fully efficient under the assumed model.
3. Sparse data in the cells can be addressed by adding proper prior distributions for the parameters and applying Bayesian methods. For example, for inference about the parameters $\theta$ of the distribution of $Y$ given $X=j,Z=k,$ one might assume a flat prior on the main effects of $X$ and $Z$ but a normal prior on the interactions, shrinking them towards zero. This achieves shrinkage of the cell mean ${\overline{y}}_{jk}$ towards the fitted means from an additive model. In frequentist terminology, this is a mixed ANOVA model with fixed main effects and random interactions.
4. The idea of raking to the $X$ and $Z$ margins is familiar, but note that if $X$ is a stratifier and $Z$ is a post-stratifier, the standard approach is to perform one iteration of raking, first matching the $X$ margin and then matching the $Z$ margin: under the assumed model, raking iteratively is the correct procedure. Raking is ML here, but Bayesian forms of raking also can be used to propagate parameter uncertainty.

Example 8. Proxy pattern-mixture analysis. Proxy pattern-mixture analysis is a method for assessing nonresponse bias for the mean of a survey variable $Y$ subject to nonresponse, when there is a set of covariates observed for nonrespondents and respondents. Historically the amount of missing data, as measured by the response rate, has been the most often-used metric for evaluating survey quality. However, response rates ignore the information contained in auxiliary covariates observed for nonrespondents. Methods based on the estimated probability of nonresponse such as the R indicator (Schouten, Cobben and Bethlehem, 2009) and q2 measure (Särndal and Lundström, 2010) do not take into account the strength of association of the survey variable of interest and the probability of response, which arguably should be factored into the assessment of bias. These measures assume MAR, but if the auxiliary variables are available for nonresponse adjustment, it is deviations from MAR that lead to bias.

Andridge and Little (2011) propose proxy pattern-mixture analysis (PPMA), a method based on a pattern-mixture model for nonresponse that combines in a simple and intuitive way the key features of nonresponse adjustment. Let $Y_i$ denote the value of a continuous survey outcome and $Z_i=\left(Z_{i1},Z_{i2},\dots ,Z_{ip}\right)$ denote the values of $p$ covariates for unit $i$ in the sample. Only $r$ of the $n$ sampled units respond, so observed data consist of $\left(Y_i,Z_i\right)$ for $i=1,\dots ,r$ and $Z_i$ for $i=r+1,\dots ,n.$ Let $R$ denote the response indicator, such that for unit $i$ $R_i=1$ if $Y_i$ is observed and $R_i=0$ if $Y_i$ is missing. To reduce dimensionality, we replace $Z$ by a single proxy variable $X$ that has the highest correlation with $Y$ in the respondent sample. This proxy variable can be estimated by regressing $Y$ on $Z$ using the respondent data, including important predictors of $Y,$ as well as interactions and nonlinear terms where appropriate. Specifically, we assume the regression model $E\left(Y|Z,R=1\right)={\alpha }_0+\alpha Z,$ and let $X=\alpha Z.$ The joint distribution of $X,$ $Y$ and $R$ using the following proxy pattern-mixture model, similar in form to that discussed in Little (1994):

$$\begin{array}{l}\left(X,Y|R=r\right)~G_2\left(({\mu }_x^{(r)},{\mu }_y^{(r)}),{\Sigma }^{(r)}\right)\\ R~\text{Bernoulli}\left(\pi \right)\\ {\Sigma }^{(r)}=\left(\begin{array}{cc}{\sigma }_{xx}^{(r)}& {\rho }^{(r)}\sqrt{{\sigma }_{xx}^{(r)}{\sigma }_{yy}^{(r)}}\\ {\rho }^{(r)}\sqrt{{\sigma }_{xx}^{(r)}{\sigma }_{yy}^{(r)}}& {\sigma }_{yy}^{(r)}\end{array}\right)\end{array}$$

where $N_2$ denotes the bivariate normal distribution. Our interest is bias in the marginal mean of $Y,$ which can be written as ${\mu }_y=\pi {\mu }_y^{(1)}+\left(1-\pi \right){\mu }_y^{(0)}.$ Andridge and Little (2011) assume that

$$\text{Pr}\left(R=1|Y,X\right)=f\left(X^{*}+\lambda Y\right),$$

for some unspecified function $f$ and known constant $\lambda .$ Here $X^{*}=X\sqrt{{\sigma }_{yy}^{(1)}/{\sigma }_{xx}^{(1)}}$ is the proxy $X$ scaled to have the same variance as $Y,$ which aids the interpretation of $\lambda$ by putting $X$ and $Y$ on the same scale. This mechanism is MAR when $\lambda =0,$ and deviates increasingly from MAR as $\lambda$ increases. With this assumption, the parameters are just identified and the ML estimate of the mean of $Y,$ averaging over patterns, is

$${\widehat{\mu }}_y={\overline{y}}_1+\left(\frac{n-r}{n}\right)\left(\frac{\lambda +\widehat{\rho }}{\lambda \widehat{\rho }+1}\right)\left({\overline{x}}_0-{\overline{x}}_1\right),$$

where ${\overline{x}}_1,{\overline{y}}_1$ are the respondent means of $X$ and $Y$ and ${\overline{x}}_0$ is the mean of $X$ for nonrespondents. The index of bias is then the adjustment of the sample mean ${\overline{y}}_1$ implied by this estimate, namely

$$\left(\frac{n-r}{n}\right)\left(\frac{\lambda +\widehat{\rho }}{\lambda \widehat{\rho }+1}\right)\left({\overline{x}}_0-{\overline{x}}_1\right),(4.14)$$

where $\widehat{\rho }$ is the respondent sample correlation. This adjustment incorporates three key factors that affect the potential bias in a simple and intuitive manner $–$ the nonresponse rate $\left(n-r\right)/n,$ the correlation $\widehat{\rho }$ between $X$ and $Y,$ and the deviation of ${\overline{x}}_0$ from ${\overline{x}}_1.$ In particular, the index increases with the nonresponse rate and the deviation of the mean of $X$ for respondents and nonrespondents.

There is no information about the parameter $\lambda$ in the data $–$ this is generally the case for methods that model deviations from MAR. Following Little (1994), Andridge and Little (2011) propose a sensitivity analysis, where estimates are generated for a range of values of $\lambda$ between 0 and infinity, specifically 0, 1 and infinity; the choice of 0 corresponds to MAR, the intermediate choice of 1 implies the bias of $Y$ is the same as the bias of the proxy variable $X^{*},$ and the choice of infinity is the most extreme deviation from MAR; estimates for this case have the highest variance. As $\lambda$ varies between 0 and infinity, the middle factor $\left(\lambda +\widehat{\rho }\right)/\left(\lambda \widehat{\rho }+1\right)$ varies between $\widehat{\rho }$ (when ${\widehat{\mu }}_y$ is the standard regression estimator of the mean) and $1/\widehat{\rho }$ (when ${\widehat{\mu }}_y$ is the inverse regression estimator proposed by Brown (1990). The sensitivity of the estimate to the choice of $\lambda$ is small when $\widehat{\rho }$ is close to 1, that is we have a strong proxy variable, and large when $\widehat{\rho }$ is close to 0, that is we have a weak proxy variable. So having auxiliary variables that are good predictors of the survey outcomes is crucial.

Bayesian versions of the proxy pattern-mixture model are readily developed, allowing for propagation of error in the model parameters. Also, the model can be applied to create multiple imputations of nonrespondent values, allowing the incorporation of complex sample design elements into the index.

More recently, the model has been used to develop indices of departure from random sampling, by applying it to model selection rather than nonresponse (Little, West, Boonstra and Hu, 2020; Boonstra, Little, West, Andridge and Alvaredo-Leiton, 2021). Refinements in this work are (a) to replace the parameter $\lambda$ by $\varphi =\lambda /\left(1+\lambda \right),$ a better parametrization because $\varphi$ ranges from 0 (MAR) to 1, and (b) (with some additional assumptions) to allow missingness also to depend on auxiliary variables orthogonal to $X.$ The method has also been extended to handle binary outcomes (Andridge, West, Little, Boonstra and Alvarado-Leiton, 2019) and indices of potential bias in regression coefficients (West, Little, Andridge, Boonstra, Ware, Pandit and Alvarado-Leiton, 2021).

• Previous
• Article main page
• Next