An Application of a Small Area Procedure with Correlation Between Measurement Error and Sampling Error to the Conservation Effects Assessment Project

Abstract

County level estimates of mean sheet and rill erosion from the Conservation Effects Assessment Project (CEAP) are useful for program development and evaluation. Since county sample sizes in the CEAP survey are insufficient to support reliable direct estimators, small area estimation procedures are needed. The quantity of water runoff is a useful covariate but is unavailable for the full population. We use an estimate of mean runoff from the CEAP survey as a covariate in a small area model with sheet and rill erosion as the response. As the runoff and sheet and rill erosion are estimators from the same survey, the measurement error in the covariate is important as is the correlation between the measurement error and the sampling error. We conduct a detailed investigation of small area estimation in the presence of a correlation between the measurement error in the covariate and the sampling error in the response. In simulations, the proposed predictor is superior to small area predictors that assume the response and covariate are uncorrelated or that ignore the measurement error entirely.

1. Introduction

The Conservation Effects Assessment Project (CEAP) is a program comprised of several surveys that are intended to evaluate the environmental impacts of agricultural production. We consider data from a CEAP survey of cropland that was conducted over the period 2003 to 2006. An important variable collected in CEAP is sheet and rill erosion (soil loss due to the flow of water). County estimates of sheet and rill erosion can improve the efficiency of allocation of resources for conservation efforts. Sample sizes in CEAP are too small to support reliable direct county estimates. Past analyses have explored a variety of issues that arise in the context of small area estimation using CEAP data (Berg and Chandra 2014; Berg and Lee 2019; Erciulescu and Fuller 2016; Lyu et al. 2020).

Traditional small area estimation procedures utilize population-level auxiliary information from censuses or administrative databases (Jiang and Lahiri 2006; Pfeffermann 2013; Rao and Molina 2015). A critical assumption underlying the seminal Fay and Herriot (1979) predictor is that one can condition on the observed value of the covariate. As discussed in Lyu et al. (2020), the task of obtaining covariates that are related to sheet and rill erosion and are known for the full population of cropland of interest is difficult. Use of variables collected in the CEAP survey as covariates is therefore desirable. We use an estimate of mean water runoff from the CEAP survey as a covariate in an area level model with sheet and rill erosion as the response. As the covariate and response are both estimates from the CEAP survey, the analysis should recognize not only the sampling error in the covariate but also the correlation between the covariate and the response.

If the covariate is an estimator from a sample survey, naive application of standard Fay-Herriot procedures can lead to erroneous inferences (Arima et al. 2017; Bell et al. 2019). A widely used technique to model sampling error in the covariate is to employ either a structural (Arima et al. 2012; Ghosh et al. 2006; Torabi 2012; Torabi et al. 2009) or a functional (Datta et al. 2010; Ghosh and Sinha 2007; Torabi 2011) measurement error model. In the structural model, the latent covariate is stochastic, while the functional model treats the unobserved covariate as fixed (Carroll et al. 2006; Fuller 2009).

Ybarra and Lohr (2008) develops predictors for an area level model in which a covariate is subject to functional measurement error. Arima et al. (2015), Arima et al. (2017), and Burgard et al. (2021) extend Ybarra and Lohr (2008) to bivariate and Bayesian frameworks. Burgard et al. (2020) conducts an analysis of the model of Ybarra and Lohr (2008) under the stronger assumption that the random terms have normal distributions. Bell et al. (2019) compares the properties of functional and structural measurement error models to the naive Fay-Herriot model. Mosaferi et al. (2023) extends Ybarra and Lohr (2008) to lognormal data. All of these works assume that the measurement error in the covariate and the sampling error in the response are uncorrelated.

Models in which the measurement error and the sampling error are correlated have received little attention in the small area estimation literature. Franco and Bell (2022) defines a bivariate model that is equivalent to a structural model with a correlation between the measurement error and the sampling error. We adopt the functional modeling approach, which, unlike the structural model, requires no assumptions about the distribution of the latent covariate. Kim et al. (2015) permits a correlation between the covariate and response but conceptualizes the parameter of interest as the unobserved value of a covariate that is measured with error. Ybarra (2003) generalizes the functional measurement error model of Ybarra and Lohr (2008) to allow for a correlation between the measurement error and the sampling error. Burgard et al. (2022) uses likelihood-based arguments to derive a predictor for a model in which the measurement error and the sampling error are correlated.

We conduct a thorough analysis of a model in which the measurement error and the sampling error are correlated. Our work expands on Ybarra (2003) in several dimensions. We conduct extensive simulation studies in a framework where the measurement error in the covariate is correlated with the sampling error in the response. We rigorously discuss the theoretical properties of the predictors with estimated parameters. Further, we provide comprehensive software at https://github.com/emilyjb/SAE-Correlated-Errors/.

The rest of this manuscript is organized as follows. We derive a predictor using properties of the bivariate normal distribution in Section 2.1. We propose estimators of the fixed model parameters in Section 2.2, and we study the theoretical properties of the proposed estimators. Then, we derive the mean squared prediction error (MSPE) of the proposed predictor. We also conduct extensive simulations to assess the properties of the proposed procedures in Section 3. We apply the proposed method to data from the CEAP survey in Section 4. We conclude in Section 5.

2. Model and Predictor

We define an area-level model, where the measurement error in the covariate is correlated with the sampling error in the response. We denote the true, unknown value of the covariate by $x_i={\left(x_{i1},\dots ,x_{ip}\right)}^{\prime }$ for $i=1,\dots ,n$ , where $n$ is the total number of small areas. We do not observe $x_i$ directly. Instead, we observe a contaminated version of $x_i$ denoted by $W_i$ . In the common situation, $W_i$ represents an estimator of $x_i$ obtained from a survey. The measurement error is functional instead of structural because $x_i$ is regarded as a fixed quantity.

The parameter of interest is

$${\theta }_i={\beta }_0+{x^{\prime }}_i{\beta }_1+b_i,$$

where $b_i\stackrel{iid}{\sim }N(0,{\sigma }_b^2)$ for $i=1,\dots ,n$ . We let $Y_i$ denote an estimator of ${\theta }_i$ . The variables representing observable quantities are then $\left\{{\left({W_i}^{\prime },Y_i\right)}^{\prime }:i=1,\dots ,n\right\}$ . We define a model for ${\left({W_i}^{\prime },Y_i\right)}^{\prime }$ as,

$$\begin{array}{c}{Y}_i={\beta }_0+{x^{\prime }}_i{\beta }_1+b_i+e_i,\\ W_i=x_i+u_i,\end{array}$$ (1)

and ${\left({u^{\prime }}_i,e_i\right)}^{\prime }\stackrel{ind}{\sim }MVN\left(0,{\Psi }_i\right),$ where $\left({u^{\prime }}_i,e_i\right)$ is independent of $b_i$ for $i=1,\dots ,n$ . We assume that ${\Psi }_i$ is known.

Typically, ${\Psi }_i$ is the design-variance of ${\left({W^{\prime }}_i,Y_i\right)}^{\prime }$ , and ${\Psi }_i$ is constructed from the unit-level data using standard procedures for complex surveys. We parametrize ${\Psi }_i$ as

$${\Psi }_i=\left(\begin{array}{cc}{\Psi }_{uui}& {\Psi }_{uei}\\ {{\Psi }^{\prime }}_{uei}& {\psi }_{eei}\end{array}\right).$$

The component ${\Psi }_{uei}$ captures the correlation between the measurement error and the sampling error. We denote the fixed vector of parameters which needs to be estimated as $\omega ={\left({\beta }_0,{{\beta }^{\prime }}_1,{\sigma }_b^2\right)}^{\prime }$ . The objective is to predict the small area parameter ${\theta }_i$ .

• Remark 1: The model in (1) has strong connections to other models in the small area estimation literature. The model is identical to the model of Ybarra (2003). If ${\Psi }_{uei}=0$ , then the model in (1) simplifies to the functional measurement error models of Ybarra and Lohr (2008) and Burgard et al. (2020). The structure of our model is also similar to the one given by Kim et al. (2015). The parameter of interest in our framework is the conditional mean of the response denoted by ${\theta }_i$ , which differs from the parameter of interest in Kim et al. (2015), the unobserved covariate, $x_i$ . In small area estimation, the parameter of interest is usually the mean of a response variable. Thus, we think that our formulation is more useful to practitioners than that of Kim et al. (2015).

• Remark 2: In many situations, as in the CEAP data analysis of Section 4, the observed covariate is an estimator from a sample survey. In this case, the error in the observed covariate is a sampling error. Because the model (1) has the form of a measurement error model, we refer to the random term $u_i$ as a measurement error instead of a sampling error. This terminology is common in the small area estimation literature (Ybarra and Lohr 2008).

2.1. Predictor as a Function of the True Model Parameters

We first define a predictor as a function of the unknown $\omega$ . This predictor is also defined in Ybarra (2003). We provide a derivation that differs slightly from that of Ybarra (2003). One can express the parameter of interest as ${\theta }_i=Y_i-e_i$ . We define a predictor of the parameter of interest as $Y_i-{\widehat{e}}_i$ , where ${\widehat{e}}_i$ is an appropriately defined predictor of $e_i$ . We now proceed to develop a form for ${\widehat{e}}_i$ . As in Fuller (2009) and Ybarra and Lohr (2008), define

$$v_i=Y_i-{\beta }_0-{{\beta }^{\prime }}_1{W}_i=b_i+e_i-{{\beta }^{\prime }}_1{u}_i.$$

Then, using properties of the bivariate normal distribution (as explained in Appendix A),

$$e_i|v_i\sim N\left(\frac{{\psi }_{eei}-{{\beta }^{\prime }}_1{\Psi }_{uei}}{{\sigma }_b^2+{\sigma }_{\delta i}^2}{v}_i,{\psi }_{eei}-\frac{{\left({\psi }_{eei}-{{\beta }^{\prime }}_1{\Psi }_{uei}\right)}^2}{{\sigma }_b^2+{\sigma }_{\delta i}^2}\right),$$ (2)

where ${\sigma }_{\delta i}^2={{\beta }^{\prime }}_1{\Psi }_{uui}{\beta }_1+{\psi }_{eei}-2{{\beta }^{\prime }}_1{\Psi }_{uei}$ . A predictor of ${\theta }_i$ is then

$${\tilde{\theta }}_i={\tilde{\theta }}_i(\omega )=Y_i-{\widehat{e}}_i(\omega ),$$ (3)

where

$${\widehat{e}}_i(\omega )=E(e_i|v_i)=\frac{{\psi }_{eei}-{{\beta }^{\prime }}_1{\Psi }_{uei}}{{\sigma }_b^2+{\sigma }_{\delta i}^2}{v}_i.$$

The MSPE of ${\tilde{\theta }}_i(\omega )$ is

$$\begin{array}{c}MSPE\left({\tilde{\theta }}_i(\omega )\right)=E\left[{\left(Y_i-e_i-\left(Y_i-{\widehat{e}}_i(\omega )\right)\right)}^2\right]\\ =V\left(e_i|v_i\right)\\ ={\psi }_{eei}-\frac{{\left({\psi }_{eei}-{{\beta }^{\prime }}_1{\Psi }_{uei}\right)}^2}{{\sigma }_b^2+{\sigma }_{\delta i}^2}:=M_{1i}(\omega ).\end{array}$$ (4)

In (4), $V\left(e_i|v_i\right)=E\left[V\left(e_i|v_i\right)\right]$ because $V\left(e_i|v_i\right)$ does not depend on any random variables.

• Remark 3: We use the properties of the bivariate normal distribution to derive the predictor (3). A different way to develop a predictor is to find the convex combination of $Y_i$ and ${\beta }_0+{{\beta }^{\prime }}_1{W}_i$ that minimizes the MSPE. Ybarra (2003) demonstrates that the predictor (3) is the minimum MSPE convex combination of $Y_i$ and ${\beta }_0+{{\beta }^{\prime }}_1{W}_i$ under moment conditions that do not require normality.

• Remark 4: Burgard et al. (2022) define a predictor for a generalization of the model (1). The predictor proposed in (3) differs from that of Burgard et al. (2022). In the supplement, we provide empirical evidence that the predictor (3) is more efficient than the predictor of Burgard et al. (2022).

2.2. Estimation of Parameters

We require estimators of ${\beta }_0$ , ${\beta }_1$ , and ${\sigma }_b^2$ . We estimate the regression coefficients by matching the sample moments with their theoretical expectations. Define ${\zeta }_1=n^{-1}{\Sigma }_{i=1}^n{Y}_i{W}_i-n^{-1}{\Sigma }_{i=1}^n{\Psi }_{uei},{\zeta }_2=n^{-1}{\Sigma }_{i=1}^n{Y}_i,{\zeta }_3=n^{-1}{\Sigma }_{i=1}^n{W}_i,and{\zeta }_4=n^{-1}{\Sigma }_{i=1}^n{W}_i{W^{\prime }}_i$ $-n^{-1}{\Sigma }_{i=1}^n{\Psi }_{uui}.$ Then, $E[{\zeta }_2]={\beta }_0+{{\beta }^{\prime }}_{1}E[{\zeta }_3]$ and $E[{\zeta }_1]={\beta }_{0}E[{\zeta }_3]+E[{\zeta }_4]{\beta }_1.$ We define an estimator of $({\beta }_0,{{\beta }_1}^{\prime })$ . as

$$\left(\begin{array}{c}{\widehat{\beta }}_0\\ {\widehat{\beta }}_1\end{array}\right)={\left(\begin{array}{cc}1& {{\zeta }^{\prime }}_3\\ {\zeta }_3& {\zeta }_4\end{array}\right)}^{-1}\left(\begin{array}{c}{\zeta }_2\\ {\zeta }_1\end{array}\right).$$

Theorem 1 of Appendix B states that the estimator of the regression coefficients is consistent, where we simplify and only consider a univariate covariate. We outline the proof of Theorem 1 in Appendix B and provide further details in the Supplemental Material.

The estimator of ${\sigma }_b^2$ defined in Ybarra and Lohr (2008) can be easily adjusted to the case of correlated errors, as in Ybarra (2003). Following Ybarra and Lohr (2008), one can define an estimator of ${\sigma }_b^2$ as

$${\widehat{\sigma }}_{b,YL}^2=\frac{1}{n-p-1}\sum _{i=1}^n\left\{{\left(Y_i-{\widehat{\beta }}_0-{\widehat{\beta }}_1{W}_i\right)}^2-{\psi }_{eei}-{{\widehat{\beta }}^{\prime }}_1{\Psi }_{uui}{\widehat{\beta }}_1+2{{\widehat{\beta }}^{\prime }}_1{\Psi }_{uei}\right\}.$$ (5)

A drawback of the estimator (5) is that it can be negative. Thus, we use a profile likelihood.

We define the profile likelihood for estimating ${\sigma }_b^2$ by

$$L({\sigma }_b^2|{\widehat{\beta }}_0,{\widehat{\beta }}_1)=\left[\prod _{i=1}^n\left(\frac{1}{\sqrt{2\pi ({\sigma }_b^2+{\widehat{\sigma }}_{\delta i}^2)}}\right)\right]\exp \left(-0.5\sum _{i=1}^n\frac{{(Y_i-{\widehat{\beta }}_0-{{\widehat{\beta }}^{\prime }}_1{W}_i)}^2}{{\sigma }_b^2+{\widehat{\sigma }}_{\delta i}^2}\right),$$ (6)

where ${\widehat{\sigma }}_{\delta i}^2={\psi }_{eei}+{\widehat{{\beta }^{\prime }}}_1{\Psi }_{uui}{\widehat{\beta }}_1-2{\widehat{{\beta }^{\prime }}}_1{\Psi }_{uei}$ . The estimator of unknown parameter ${\sigma }_b^2$ is

$${\widehat{\sigma }}_b^2={\mathrm{argmax}}_{{\sigma }_b^2\ge 0}L({\sigma }_b^2|{\widehat{\beta }}_0,{\widehat{\beta }}_1),$$ (7)

where the maximization is over the parameter space for ${\sigma }_b^2$ . The profile likelihood estimator is similar to a maximum likelihood estimator (MLE) in the sense that the estimator does not account for the loss of degrees of freedom from estimating regression coefficients. Another possibility is to construct a restricted ML (REML)-type estimator to improve upon the properties of ${\widehat{\sigma }}_b^2$ , and this is a possible future research direction. We give theoretical consideration to the properties of ${\widehat{\sigma }}_b^2$ for a univariate covariate. Theorem 2 in Appendix B states that ${\widehat{\sigma }}_b^2$ is a consistent estimator of ${\sigma }_b^2$ . We present a proof of Theorem 2 in the Supplement.

• Remark 5: For the estimation of ${\sigma }_b^2$ , Li and Lahiri (2010) proposed an adjusted likelihood. We tried this technique and found that the resulting estimator can have a positive bias in simulations. We therefore use the profile likelihood function (6) for our simulations and data analysis. The benefits of our proposed estimator are that the estimator is simple and has tractable theoretical properties, as we discuss in Appendix B.

• Remark 6: It may be observed that we use a likelihood-based estimator of ${\sigma }_b^2$ and a moment-based estimator of the regression coefficients. An alternative is to develop a likelihood-based estimator of the regression coefficients, along the lines of Burgard et al. (2020). We prefer the moment-based estimator for two main reasons. First, the estimator can be calculated in one step, enabling a computationally simple procedure. Second, the moment-based estimator is robust to the assumption of normality, as we demonstrate through the simulation study. We prefer the profile-likelihood estimator of ${\sigma }_b^2$ over the moment-based estimator for the purpose of avoiding negative estimates.

2.3. Predictors with Estimated Parameters

We evaluate the predictor (3) at the estimator of $\omega$ . We define the predictor as

$${\widehat{\theta }}_i={\tilde{\theta }}_i(\widehat{\omega })=Y_i-{\widehat{e}}_i(\widehat{\omega }),$$ (8)

where $\widehat{\omega }=({\widehat{\beta }}_0,{\widehat{{\beta }^{\prime }}}_1,{\widehat{\sigma }}_b^2{)}^{\prime }$ . The vector of estimated parameters $\widehat{\omega }$ is obtained using the procedure of Section 2.2. The MSPE of ${\widehat{\theta }}_i$ decomposes into a sum of three terms as

$$\begin{array}{c}MSPE({\widehat{\theta }}_i)=E\left[{\left(Y_i-{\widehat{e}}_i(\widehat{\omega })-(Y_i-e_i)\right)}^2\right]\\ =E\left[{\left({\widehat{e}}_i(\widehat{\omega })-e_i\right)}^2\right]\\ =E\left[{\left({\widehat{e}}_i(\omega )-e_i\right)}^2\right]+E\left[{\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)}^2\right]\\ +2E\left[\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)\left({\widehat{e}}_i(\omega )-e_i\right)\right]\\ =V(e_i|v_i)+E\left[{\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)}^2\right]+2E\left[\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)\left({\widehat{e}}_i(\omega )-e_i\right)\right]\\ =M_{1i}(\omega )+M_{2i}+2E\left[\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)\left({\widehat{e}}_i(\omega )-e_i\right)\right].\end{array}$$

The first term, $M_{1i}(\omega )=V(e_i|v_i)=E\left[V(e_i|v_i)\right]$ , is defined in (4). The second term, $M_{2i}=E\left[{\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)}^2\right]$ , accounts for the variance of $\widehat{\omega }$ . Consider the cross term defined as $2E\left[\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)\left({\widehat{e}}_i(\omega )-e_i\right)\right]$ . Suppose $\widehat{\omega }$ is independent of $e_i$ given $v_i$ . Then,

$$\begin{array}{c}MSPE({\widehat{\theta }}_i)=E\left[{\left(Y_i-{\widehat{e}}_i(\widehat{\omega })-(Y_i-e_i)\right)}^2\right]\\ =E\left[{\left({\widehat{e}}_i(\widehat{\omega })-e_i\right)}^2\right]\\ =E\left[{\left({\widehat{e}}_i(\omega )-e_i\right)}^2\right]+E\left[{\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)}^2\right]\\ +2E\left[\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)\left({\widehat{e}}_i(\omega )-e_i\right)\right]\\ =V(e_i|v_i)+E\left[{\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)}^2\right]+2E\left[\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)\left({\widehat{e}}_i(\omega )-e_i\right)\right]\\ =M_{1i}(\omega )+M_{2i}+2E\left[\left({\widehat{e}}_i(\widehat{\omega })-{\widehat{e}}_i(\omega )\right)\left({\widehat{e}}_i(\omega )-e_i\right)\right].\end{array}$$

The final equality holds because ${\widehat{e}}_i(\omega )=E[e_i|v_i]$ , and the preceding equality holds by the assumption that $\widehat{\omega }$ is independent of $e_i$ given $v_i$ . The MSPE of ${\widehat{\theta }}_i$ then decomposes into a sum of two terms as

$$MSPE({\widehat{\theta }}_i)=M_{1i}(\omega )+M_{2i}.$$

We use a plug-in estimator of $M_{1i}(\omega )$ defined as

$${\widehat{M}}_{1i}=M_{1i}(\widehat{\omega }),$$

where $M_{1i}(\omega )$ is defined in (4). We use the jackknife technique to estimate $M_{2i}$ as well as the bias of ${\widehat{M}}_{1i}$ for $M_{1i}(\omega )$ . Let ${\widehat{\omega }}^{(k)}$ denote the estimator of $\omega$ with area $k$ omitted. The jackknife estimator of $M_{2i}$ is defined as

$${\widehat{M}}_{2i,JK}=\sum _{k=1}^n{\left({\widehat{e}}_i\left({\widehat{\omega }}^{(k)}\right)-n^{-1}\sum _{k=1}^n{\widehat{e}}_i\left({\widehat{\omega }}^{(k)}\right)\right)}^2.$$

The jackknife estimate of the bias of the estimator of $M_{1i}$ is

$${\widehat{b}}_{i,JK}=n^{-1}\sum _{k=1}^n{M}_{1i}\left({\widehat{\omega }}^{(k)}\right)-{\widehat{M}}_{1i}(\widehat{\omega }).$$

The estimator of the MSPE is then defined as

$${\widehat{MSPE}}_i={\widehat{M}}_{1i}+{\widehat{M}}_{2i,JK}-{\widehat{b}}_{i,JK}.$$ (9)

• Remark 7: The simplifying assumption that $e_i$ is independent of $\widehat{\omega }$ given $v_i$ facilitates construction of a simple MSPE estimator. The simulation studies presented in Section 3 verify that the MSPE estimator, constructed under this assumption, has good properties.

• Remark 8: Alternatives to the jackknife variance estimator are Taylor linearization and the bootstrap. For this model, Taylor linearization is possible, but the operations are tedious. We prefer the jackknife relative to Taylor linearization for simplicity of implementation.

• Remark 9: We use the assumption of normality when formulating the predictor and when proving that the parameter estimators are consistent. The procedures, however, do not rely heavily on the normality assumption. The development of the predictor in Ybarra (2003) as the optimal convex combination between the direct estimator and ${\beta }_0+{{\beta }^{\prime }}_1{W}_i$ does not require normality. The estimators of regression coefficients remain consistent under suitable assumptions on the fourth moments. We study the robustness of the prediction procedure to departures from normality through simulations.

3. Simulations

We conduct simulations with two goals. The first is to understand the effect of the nature of ${\Psi }_i$ on the properties of the predictor. The second is to assess effects of departures from normality. We first simulate data from normal distributions, as specified in model (1). For the simulations of Section 3.1, we use unequal ${\Psi }_i$ . For the simulations of Section 3.2, we use a constant value of ${\Psi }_i=\Psi$ for $i=1,\dots ,n$ . We generate data from $t$ distributions in Section 3.3. For both normal and $t$ distributions, we use a univariate covariate so that ${\Psi }_i$ is a 2 × 2 matrix, and ${\beta }_1$ is a scalar. We generate a fixed set of $x_i$ as independent chi-square random variables with five degrees of freedom. We set $({\beta }_0,{\beta }_1,{\sigma }_b^2)=(1,2,0.36)$ .

As one of the simulation objectives is to understand the effects of the form of ${\Psi }_i$ on the properties of the predictor, we define a general form for ${\Psi }_i$ by

$${\Psi }_j=(0.75+0.25j{)}^2\mathrm{diag}(\sqrt{a},\sqrt{b})\left(\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right)\mathrm{diag}(\sqrt{a},\sqrt{b}),$$ (10)

for $j=1,2,3,4$ . We conduct simulations with equal and uneuqal ${\Psi }_i$ . For the simulations with unequal ${\Psi }_i$ , we set one-fourth of the $\{{\Psi }_i:i=1,\dots ,n\}$ equal to ${\Psi }_j$ for $j=1,2,3,4$ . Areas $(j-1)n/4+1,\dots ,(jn)/4$ are assigned ${\Psi }_j$ . For the simulations with equal ${\Psi }_i$ , we set ${\Psi }_i={\Psi }_1$ for $i=1,\dots ,n$ . The configurations are chosen to reflect a range of conditions.

We define eight simulation configurations by four combinations of $(a,b,\rho )$ and two sample sizes. First, we set $(a,b,\rho )=(0.25,0.75,0.2)$ . For this configuration, the measurement error variance is smaller than the sampling error variance, and the correlation between the measurement error and the sampling error is 0.2. We then increase the correlation and set $(a,b,\rho )=(0.25,0.75,0.8)$ . Next, we reverse $a$ and $b$ so that the measurement error variance exceeds the sampling error variance. For the third and fourth choices of ${\Psi }_i$ , we define $(a,b,\rho )=(0.75,0.25,0.2)$ and $(a,b,\rho )=(0.75,0.25,0.8)$ . For each choice of ${\Psi }_i$ defined by $(a,b,\rho )$ , we use two sample sizes of $n=100$ and $n=500$ . For each of the eight configurations, we conduct a Monte Carlo (MC) simulation with a MC sample size of 1,000.

We refer to the procedure proposed in Section 2 as ME-Cor. We compare the proposed procedure to two primary competitors. One competitor is the approach of Ybarra and Lohr (2008), which assumes that $Cor(u_i,e_i)=0$ . We abbreviate the Ybarra and Lohr (2008) procedure as “YL”. We implement estimation and prediction for the model of Ybarra and Lohr (2008) using the R package saeME. The other competitor is the standard estimator and predictor for the traditional Fay and Herriot (1979) model, abbreviated as “FH”. This competitor is of practical interest because naive application of the Fay-Herriot model when the covariate is from a sample survey is tempting for its simplicity. We implement estimation, prediction, and MSPE estimation for the standard Fay-Herriot model using the R package SAE.

We do not include the predictor outlined in Ybarra (2003) in the simulations for two main reasons. One is that the predictor of Ybarra (2003) is not fully developed for the case of a correlation between the measurement error and the sampling error. The other is that we do not view the procedure of Ybarra (2003) as a competitor to our approach. Instead, our objective is to build on the predictor of Ybarra (2003) and study its properties in more detail.

We also do not compare our predictor to a predictor for a bivariate model in which the covariate is included as a second response variable. As discussed in Section 1, Franco and Bell (2022) considers a bivariate modeling approach. Their approach is equivalent to a structural measurement error model with no covariates. We prefer to remain in the framework of functional measurement error. Therefore, we do not include a comparison to a predictor based on a bivariate model.

3.1. Normal Distributions, Unequal ${\Psi }_i$

We first simulate data with normally distributed random components, as specified in model (1). We use the unequal ${\Psi }_i$ defined in (10). We compare the alternative estimators and predictors for the case of normal distributions and unequal ${\Psi }_i$ in Tables 1 and 2. Table 1 contains the MC means and standard deviations of the alternative estimators of the fixed model parameters. The properties of the estimators are similar across the eight configurations. At the sample size of $n=500$ , the proposed estimator (ME-Cor) is approximately unbiased for the intercept and slope. For $n=500$ , the MC bias of the proposed estimator of ${\sigma }_b^2$ is below 0.01 for each combination of $(a,b,\rho )$ .

Table 1.

MC Means and Standard Deviations of Estimators of Fixed Model Parameters. True Values Are $({\beta }_0,{\beta }_1,{\sigma }_b^2)=(1,2,0.36)$ . Normal Distributions; Unequal ${\Psi }_i$ .

	$a$	$b$	$\rho$	$n$	ME-Cor		YL		FH
					MC mean	MC SD	MC mean	MC SD	MC mean	MC SD
β0				100	0.987	0.298	0.942	0.263	1.266	0.260
β1	0.250	0.750	0.200	100	2.002	0.054	2.011	0.046	1.944	0.045
${\sigma }_b^2$				100	0.351	0.236	0.083	0.172	1.050	0.291
β0				500	0.993	0.122	0.953	0.112	1.218	0.109
β1	0.250	0.750	0.200	500	2.001	0.021	2.009	0.019	1.956	0.018
${\sigma }_b^2$				500	0.357	0.108	0.019	0.053	1.045	0.125
β0				100	0.996	0.176	0.840	0.170	1.071	0.164
β1	0.250	0.750	0.800	100	2.001	0.027	2.028	0.026	1.988	0.025
${\sigma }_b^2$				100	0.344	0.106	0.000	0.000	0.058	0.073
β0				500	1.001	0.077	0.810	0.079	1.084	0.075
${\beta }_1$	0.250	0.750	0.800	500	2.000	0.013	2.038	0.013	1.983	0.012
${\sigma }_b^2$				500	0.359	0.047	0.000	0.000	0.044	0.037
β0				100	0.968	0.470	0.912	0.431	1.972	0.363
β1	0.750	0.250	0.200	100	2.008	0.085	2.018	0.075	1.803	0.062
${\sigma }_b^2$				100	0.377	0.375	0.175	0.335	3.410	0.568
β0				500	1.002	0.186	0.957	0.163	1.833	0.151
β1	0.750	0.250	0.200	500	2.000	0.031	2.009	0.027	1.838	0.024
${\sigma }_b^2$				500	0.361	0.191	0.067	0.135	3.502	0.257
β0				100	0.991	0.352	0.778	0.325	1.758	0.272
β1	0.750	0.250	0.800	100	2.004	0.064	2.048	0.060	1.844	0.047
${\sigma }_b^2$				100	0.350	0.282	0.001	0.018	2.174	0.392
β0				500	0.996	0.149	0.819	0.138	1.708	0.119
β1	0.750	0.250	0.800	500	2.000	0.026	2.036	0.024	1.856	0.020
${\sigma }_b^2$				500	0.360	0.135	0.000	0.000	2.205	0.169

Table 2.

MC MSPE’s of Alternative Predictors (MC MSPE of Predictor) and MC Means of Estimated MSPE’s (MC Mean Est. MSPE) of ME-Cor and FH Predictors. Normal Distributions; Unequal ${\Psi }_i$ .

		MC MSPE of predictor				MC mean est. MSPE
$(a,b,\rho )$	$n$	Direct	ME-Cor	YL	FH	ME-Cor-MSPE	FH-MSPE
(0.250, 0.750, 0.200)	100	1.010	0.748	0.759	0.809	0.744	0.499
	500	1.007	0.741	0.758	0.804	0.739	0.487
(0.250, 0.750, 0.800)	100	1.010	1.002	1.105	1.584	1.001	0.087
	500	1.007	1.000	1.112	1.614	1.001	0.048
(0.750, 0.250, 0.200)	100	0.337	0.335	0.345	0.357	0.334	0.302
	500	0.336	0.334	0.345	0.356	0.333	0.301
(0.750, 0.250, 0.800)	100	0.334	0.214	0.440	0.559	0.215	0.285
	500	0.336	0.213	0.443	0.563	0.212	0.285

The proposed estimator of ${\sigma }_b^2$ usually has a negative bias when $n=100$ . A negative small-sample bias for the estimator of ${\sigma }_b^2$ is not surprising because the estimator does not incorporate a correction for the loss of degrees of freedom associated with estimating ${\beta }_0$ and ${\beta }_1$ . An exception to the negative bias occurs when the measurement error exceeds the sampling error and when the correlation between $u_i$ and $e_i$ is only 0.2. For this configuration, the estimator of ${\sigma }_b^2$ has a positive bias at $n=100$ . Increasing the sample size to $n=500$ rectifies the bias of the estimator of ${\sigma }_b^2$ for the configuration with $(a,b,\rho )=(0.75,0.25,0.2)$ . When $(a,b,\rho )=(0.75,0.25,0.2)$ , the distribution of the estimator of ${\sigma }_b^2$ for $n=100$ is highly skewed right and has extreme values. When the sample size increases to $n=500$ , the distribution of the estimator of ${\sigma }_b^2$ for this configuration is more unimodal and symmetric.

The presence of a nontrivial correlation between $u_i$ and $e_i$ causes the YL estimator to have a negative bias for the intercept and a positive bias for the slope. The YL estimator of ${\sigma }_b^2$ has a severe negative bias in the presence of a nonzero correlation between $u_i$ and $e_i$ . The R function FHme, used to implement the YL procedure, applies a lower bound of zero to the estimator of ${\sigma }_b^2$ , and it is apparent from the results that many of the estimates reach the lower bound of zero.

The measurement error attenuates the FH estimator of the slope toward zero and leads to a positive bias in the estimator of the intercept. The FH estimator of ${\sigma }_b^2$ usually has a positive bias because the FH estimator of ${\sigma }_b^2$ includes part of the measurement error. For configurations with $a<b$ and $\rho =0.8$ , the covariance structure causes the FH estimator of ${\sigma }_b^2$ to have a negative bias.

Table 2 summarizes the empirical properties of the alternative predictors and MSPE estimators. The columns under the heading “MC MSPE of Predictor” contain the average MC MSPE’s of the alternative predictors, where the average is across areas. The columns under the heading “MC Mean Est. MSPE” contain the average MC means of the MSPE estimators for the ME-Cor and FH procedures. The column labeled “Direct” indicates the average MC MSPE of the direct estimator, $Y_i$ .

The YL predictor is superior to the FH predictor but inferior to the ME-Cor predictor. For all but the configuration with $(a,b,\rho )=(0.25,0.75,0.2),$ the YL predictor has MC MSPE exceeding that of the direct estimator. The results for the YL predictor demonstrate the importance of accounting for a correlation between $u_i$ and $e_i$ .

The properties of the FH predictor depend on the structure of ${\Psi }_i$ and are similar for $n=100$ and $n=500$ . When $\rho =0.2$ and $a<b$ (measurement error variance is smaller than the sampling variance), the FH predictor is superior to the direct estimator but inferior to the ME-Cor predictor. For all other configurations, the MC MSPE of the FH predictor exceeds the MC MSPE of the direct estimator. This empirical finding echoes a theoretical result in Ybarra and Lohr (2008) that the Fay-Herriot predictor can have MSPE greater than the variance of the direct estimator if the covariate is measured with error.

The efficiency of the FH predictor relative to the direct estimator is best when $a<b$ and $\rho =0.2$ . This is reasonable because this configuration most closely approximates a situation where the covariate is measured without error. The MC means of the estimators of ${\sigma }_b^2$ in Table 1 provide insight into why the FH predictor is less efficient than the direct estimator for $(a,b,\rho )\ne (0.25,0.75,0.2)$ . The MC means of the estimators of ${\sigma }_b^2$ in Table 1 reveal the impacts of measurement error on the shrinkage parameters for the FH predictor, where the shrinkage parameter is defined as ${\widehat{\sigma }}_b^2/({\widehat{\sigma }}_b^2+{\psi }_{eei})$ . For the scenarios with $(a,b,\rho )=(0.25,0.75,0.8)$ , the estimator of ${\sigma }_b^2$ has a severe negative bias, leading to considerable over-shrinkage toward a covariate that is itself measured with error. For scenarios with $a>b$ , the positive bias of the estimator of ${\sigma }_b^2$ is overwhelming, and the FH estimator does not exhibit enough shrinkage toward the estimated regression line.

An important implication of measurement error is that the Fay-Herriot MSPE estimator (FH-MSPE) has a severe negative bias for the MSPE of the Fay-Herriot predictor (FH). When $a<b$ and $\rho =0.8$ , the MC mean of the FH MSPE estimator is more than an order of magnitude lower than the MC MSPE of the FH predictor. A risk of naive application of Fay-Herriot procedures in the presence of measurement error is under-estimation of the MSPE.

The ME-Cor predictor has smaller MC MSPE than the alternatives considered for all configurations. When $a<b$ , the gain in efficiency from the ME-Cor predictor relative to the direct estimator is greater for $\rho =0.2$ than for $\rho =0.8$ . When $a>b$ , the opposite pattern holds, as the ratio of the MC MSPE of the direct estimator to the MC MSPE of the ME-Cor predictor is greater for $\rho =0.8$ than for $\rho =0.2$ . The ME-Cor procedure renders only trivial improvements in efficiency over the direct estimator for configurations with $(a,b,\rho )=(0.25,0.75,0.8)$ or $(a,b,\rho )=(0.75,0.25,0.2)$ . Increasing $n$ from 100 to 500 has little effect on the properties of the ME-Cor predictor.

The proposed MSPE estimator (ME-Cor-MSPE) is a good approximation for the MSPE of the ME-Cor predictor (ME-Cor). For each configuration, the average MC mean of the estimated MSPE for the ME-Cor predictor (ME-Cor-MSPE) is close to the average MC MSPE of the ME-Cor predictor (ME-Cor). The simulation results support the predictor and MSPE estimator proposed in Section 2.

3.2. Normal Distributions, Equal ${\Psi }_i=\Psi$

A special case in which the FH predictor retains reasonable properties occurs in the context of the structural model when the measurement error is uncorrelated with the sampling error and when the measurement error variance is constant (Bell et al. 2019). When the measurement error and sampling error are correlated, the naive Fay-Herriot predictor remains inappropriate, even if the measurement error variance is constant. To illustrate this point, we present simulation results with equal ${\Psi }_i=\Psi$ .

Tables 3 and 4 contain simulation results for ${\Psi }_i={\Psi }_1$ for $i=1,\dots ,n$ , where ${\Psi }_1$ is defined in (10). We continue to simulate the errors from normal distributions. The conclusions for equal ${\Psi }_i$ are the same as those for unequal ${\Psi }_i$ . As seen in Table 3, the FH and YL estimators of the fixed parameters remain biased when ${\Psi }_i=\Psi$ . In contrast, the MC means of the proposed estimators of the fixed parameters are close to the true parameter values. In Table 4, the proposed predictor has smaller MC MSPE than the alternatives for all configurations and sample sizes. The Fay-Herriot predictor remains inefficient when ${\Psi }_i=\Psi$ . Also, the FH estimator of the MSPE continues to have a negative bias in the case of equal ${\Psi }_i$ . The proposed MSPE estimator is nearly unbiased for the MSPE of the ME-Cor predictor.

Table 3.

MC Means and Standard Deviations of Estimators of Fixed Model Parameters. True Values Are $({\beta }_0,{\beta }_1,{\sigma }_b^2)=(1,2,0.36)$ . Normal Distributions, Equal ${\Psi }_i=\Psi$ .

	$a$	$b$	$\rho$	$n$	ME-Cor		YL		FH
					MC mean	MC SD	MC mean	MC SD	MC mean	MC SD
β0				100	0.994	0.252	0.949	0.253	1.202	0.243
β1	0.250	0.750	0.200	100	2.001	0.044	2.010	0.044	1.959	0.042
${\sigma }_b^2$				100	0.336	0.228	0.090	0.143	0.999	0.248
β0				500	0.997	0.111	0.953	0.111	1.204	0.107
β1	0.250	0.750	0.200	500	2.001	0.019	2.010	0.019	1.960	0.018
${\sigma }_b^2$				500	0.351	0.112	0.047	0.067	0.996	0.112
β0				100	0.996	0.154	0.826	0.156	1.070	0.151
β1	0.250	0.750	0.800	100	2.001	0.027	2.035	0.027	1.986	0.026
${\sigma }_b^2$				100	0.341	0.101	0.000	0.000	0.027	0.052
β0				500	0.999	0.077	0.798	0.078	1.086	0.075
β1	0.250	0.750	0.800	500	2.000	0.014	2.041	0.014	1.982	0.013
${\sigma }_b^2$				500	0.358	0.047	0.000	0.000	0.008	0.019
β0				100	0.962	0.377	0.914	0.380	1.697	0.317
β1	0.750	0.250	0.200	100	2.005	0.065	2.015	0.066	1.863	0.053
${\sigma }_b^2$				100	0.352	0.357	0.153	0.256	2.812	0.433
β0				500	0.988	0.148	0.948	0.149	1.604	0.129
β1	0.750	0.250	0.200	500	2.003	0.025	2.011	0.026	1.879	0.021
${\sigma }_b^2$				500	0.342	0.196	0.078	0.120	2.839	0.192
β0				100	0.983	0.286	0.825	0.294	1.475	0.251
β1	0.750	0.250	0.800	100	2.004	0.047	2.035	0.049	1.904	0.040
${\sigma }_b^2$				100	0.328	0.268	0.000	0.001	1.860	0.290
β0				500	0.987	0.127	0.822	0.131	1.502	0.110
β1	0.750	0.250	0.800	500	2.002	0.023	2.036	0.023	1.896	0.019
${\sigma }_b^2$				500	0.352	0.142	0.000	0.000	1.857	0.136

Table 4.

MC MSPE’s of Alternative Predictors (MC MSPE of Predictor), and MC Means of Estimated MSPE’s of ME-Cor and FH Predictors (MC Mean Est. MSPE). Normal Distributions, Equal ${\Psi }_i=\Psi$ .

$(a,b,\rho )$	$n$	MC MSPE of predictor				MC mean est. MSPE
		Direct	ME-Cor	YL	FH	ME-Cor-MSPE	FH-MSPE
(0.250, 0.750, 0.200)	100	0.745	0.564	0.575	0.580	0.564	0.441
	500	0.750	0.563	0.578	0.576	0.561	0.430
(0.250, 0.750, 0.800)	100	0.747	0.744	0.842	1.291	0.746	0.067
	500	0.752	0.747	0.847	1.329	0.745	0.017
(0.750, 0.250, 0.200)	100	0.249	0.247	0.255	0.255	0.248	0.230
	500	0.250	0.248	0.257	0.256	0.248	0.230
(0.750, 0.250, 0.800)	100	0.248	0.164	0.327	0.372	0.163	0.222
	500	0.250	0.162	0.330	0.376	0.162	0.221

3.3. Simulations with $t$ Distributions

We next simulate data from $t$ distributions. This allows us to assess robustness of the procedure to departures from normality. To simulate data from $t$ distributions, we first generate ${\tilde{u}}_i\stackrel{iid}{\sim }{t}_{(5)}/\sqrt{5/3}$ , ${\tilde{e}}_i\stackrel{iid}{\sim }{t}_{(5)}/\sqrt{5/3}$ , and ${\tilde{b}}_i\stackrel{iid}{\sim }{t}_{(5)}/\sqrt{5/3}$ . The notation $t_{(5)}/\sqrt{5/3}$ denotes a $t$ distribution with five degrees of freedom divided by $\sqrt{5/3}$ . The division by $\sqrt{5/3}$ standardizes the variables so that ${\tilde{u}}_i$ , ${\tilde{e}}_i$ , and ${\tilde{b}}_i$ have zero mean and unit variance. The ${\tilde{u}}_i$ , ${\tilde{e}}_i$ , and ${\tilde{b}}_i$ are mutually independent. We then define $(u_i,e_i{)}^{\prime }={\Psi }_i^{0.5}({\tilde{u}}_i,{\tilde{e}}_i{)}^{\prime }$ , where ${\Psi }_i^{0.5}$ is the square root matrix of ${\Psi }_i$ . We set $b_i={\sigma }_b{\tilde{b}}_i$ . We use the ${\Psi }_i$ defined in (10) with the four combinations of $(a,b,\rho )$ .

The results for the $t$ distribution are presented in Tables 5 and 6. Table 5 contains the MC means and standard deviations of the estimators of the fixed parameters when the random terms are generated from $t$ distributions. Table 6 contains the MC MSPE’s of the predictors as well as the MC means of the estimated MSPE’s for the proposed and Fay-Herriot predictors. The results for the $t$ distribution are largely similar to the results for the normal distribution.

Table 5.

MC Means and Standard Deviations of Estimators of Fixed Model Parameters. True Values Are $({\beta }_0,{\beta }_1,{\sigma }_b^2)=(1,2,0.36)$ . $t$ Distributions; Unequal ${\Psi }_i$ .

	$a$	$b$	$\rho$	$n$	ME-Cor		YL		FH
					MC mean	MC SD	MC mean	MC SD	MC mean	MC SD
β0				100	0.990	0.298	0.954	0.272	1.204	0.268
β1	0.250	0.750	0.200	100	2.002	0.046	2.009	0.042	1.963	0.041
${\sigma }_b^2$				100	0.339	0.300	0.118	0.273	1.039	0.380
β0				500	1.004	0.134	0.959	0.123	1.245	0.120
β1	0.250	0.750	0.200	500	2.000	0.024	2.008	0.021	1.951	0.021
${\sigma }_b^2$				500	0.354	0.149	0.037	0.093	1.037	0.169
β0				100	1.000	0.185	0.806	0.182	1.084	0.174
β1	0.250	0.750	0.800	100	2.000	0.030	2.039	0.030	1.984	0.028
${\sigma }_b^2$				100	0.349	0.154	0.000	0.000	0.070	0.114
β0				500	1.001	0.080	0.817	0.080	1.082	0.077
β1	0.250	0.750	0.800	500	2.000	0.014	2.037	0.014	1.984	0.013
${\sigma }_b^2$				500	0.358	0.072	0.000	0.000	0.047	0.056
β0				100	0.945	0.504	0.903	0.446	1.928	0.390
β1	0.750	0.250	0.200	100	2.011	0.097	2.020	0.086	1.807	0.073
${\sigma }_b^2$				100	0.419	0.576	0.300	0.652	3.392	0.855
β0				500	0.992	0.237	0.943	0.205	1.950	0.187
β1	0.750	0.250	0.200	500	2.003	0.043	2.012	0.037	1.810	0.033
${\sigma }_b^2$				500	0.369	0.306	0.131	0.267	3.431	0.387
β0				100	0.975	0.404	0.797	0.343	1.750	0.311
β1	0.750	0.250	0.800	100	2.005	0.081	2.044	0.067	1.837	0.060
${\sigma }_b^2$				100	0.399	0.482	0.026	0.240	2.179	0.593
β0				500	0.992	0.175	0.799	0.156	1.701	0.140
β1	0.750	0.250	0.800	500	2.001	0.032	2.040	0.028	1.859	0.025
${\sigma }_b^2$				500	0.363	0.224	0.002	0.049	2.210	0.275

Table 6.

MC MSPE’s of Alternative Predictors (MC MSPE of Predictor), and MC Means of Estimated MSPE’s of ME-Cor and FH Predictors (MC Mean Est. MSPE). $t$ Distributions; Unequal ${\Psi }_i$ .

$(a,b,\rho )$	$n$	MC MSPE of predictor				MC mean est. MSPE
		Direct	ME-Cor	YL	FH	ME-Cor-MSPE	FH-MSPE
(0.250, 0.750, 0.200)	100	0.995	0.741	0.750	0.814	0.744	0.493
	500	1.004	0.738	0.754	0.802	0.739	0.484
(0.250, 0.750, 0.800)	100	0.989	0.983	1.096	1.559	1.001	0.094
	500	1.011	1.004	1.115	1.614	1.001	0.049
(0.750, 0.250, 0.200)	100	0.331	0.329	0.339	0.352	0.334	0.301
	500	0.334	0.332	0.342	0.353	0.333	0.301
(0.750, 0.250, 0.800)	100	0.338	0.219	0.444	0.562	0.218	0.284
	500	0.337	0.215	0.445	0.565	0.213	0.284

For configurations with $a<b$ , the estimator of ${\sigma }_b^2$ has a negative bias. The negative bias is expected because the objective function (6) does not account for the loss of degrees of freedom from estimating regression coefficients. When $a>b$ , the heavy tails of the $t$ distribution cause the estimator of ${\sigma }_b^2$ to have a positive bias. The positive bias is notable when $n=100$ . Because the random effects have $t$ distributions, the likelihood used to define the estimator of ${\sigma }_b^2$ is misspecified for this configuration. Nonetheless, increasing the sample size to $n=500$ markedly reduces the bias.

The positive bias of the estimator of ${\sigma }_b^2$ has minimal impacts on the properties of the predictor. As seen in Table 6, the proposed predictor has smaller MC MSPE than the alternatives. The proposed MSPE estimator is a good approximation for the MSPE of the predictor, even when the data are generated from $t$ distributions. The results for the $t$ distribution support the statement in remark 9 that the predictor and MSPE estimator are robust to the assumption of normality.

The FH and YL procedures remain inefficient when the data are generated from $t$ distributions. The bias in the estimators of fixed parameters from the FH and YL procedures is much more severe than the bias from the proposed procedure, even when the random terms are generated from $t$ distributions. The bias of the FH and YL estimators propagates into the predictor, resulting in high prediction MSPE’s for the FH and YL procedures.

3.4. Extended Simulations

We present extended simulation results in the Supplemental Material. First, we use a $t_{(3)}$ distribution instead of a $t_{(5)}$ distribution for the random terms. We consider the $t_{(3)}$ distribution because, unlike the $t_{(5)}$ distribution, the $t_{(3)}$ distribution does not have a finite fourth moment. We then consider a distribution where the random terms are distributed as centered and scaled ${\chi }_{(3)}^2$ random variables. We use the ${\chi }_{(3)}^2$ distribution because it is skewed.

These extended configurations allow us to assess the impacts of skewness and absence of fourth moments on the properties of the proposed procedure. The positive bias of the estimator of ${\sigma }_b^2$ is more severe for these configurations because the profile likelihood (6) is misspecified. The bias for the estimator of ${\sigma }_b^2$ has little impact on the other model parameters. The estimators of regression coefficients remain approximately unbiased under the $t_{(3)}$ and ${\chi }_{(3)}^2$ distributions. The predictors remain more efficient than the alternatives considered. Despite the bias for ${\sigma }_b^2$ , the proposed MSPE estimator continues to provide a reasonable approximation to the MSPE of the predictor. We refer the reader to the Supplemental Material for further detail.

We also validate the proposed procedure for a multivariate covariate in the Supplemental Material. We use two covariates, both of which are measured with error. The estimators of the fixed parameters remain nearly unbiased in the presence of a bivariate covariate. The proposed MSPE estimator is also nearly unbiased for the MSPE of the predictor.

4. CEAP Data Analysis

We apply the method proposed in Section 2 to predict mean log sheet and rill erosion in Iowa counties using CEAP data. Iowa has $n=99$ counties as small areas. In CEAP, sheet and rill erosion is measured using a computer model called RUSLE2. A variable that impacts the amount of sheet and rill erosion in a county is the quantity of water runoff. The mean runoff is unknown for the full population of cropland in Iowa. We use the sample mean of runoff obtained from the CEAP survey as the covariate for the small area model. The response is the log of the sample mean of RUSLE2.

We connect the context and notation of Section 2 to the CEAP data analysis. Let ${\tilde{Y}}_i$ denote the direct estimator of mean RUSLE2 erosion in county $i$ , where $i=1,\dots ,n$ . Let ${\tilde{W}}_i$ denote the direct estimator of the mean runoff for county $i$ . The direct estimators are defined in Appendix C. The unknown population mean runoff for county $i$ is denoted by ${\tilde{x}}_i$ . We define the small area model in the log scale. The response variable for the small area model is defined by $Y_i=\log ({\tilde{Y}}_i)$ . The covariate $W_i=\log ({\tilde{W}}_i)$ is a contaminated measurement of the log of the population mean runoff, defined as $x_i=\log ({\tilde{x}}_i)$ . Figure 1 contains a plot of $Y_i$ on the vertical axis against $W_i$ on the horizontal axis. The figure exhibits high variation in the association between runoff and RUSLE2. The variation may arise from inherent variability between the counties as well as measurement error in the covariate. The model (1) accounts for both of these sources of variation.

Figure 1.

Scatterplot of $Y_i=\log ({\tilde{Y}}_i)$ on the y-axis against $W_i=\log ({\tilde{W}}_i)$ on the x-axis.

The model requires an estimate of ${\Psi }_i$ , the design variance of $(W_i,Y_i)$ . We explain how we estimate the design variance of $(W_i,Y_i)$ in Appendix C. The estimates of $Cor(u_i,e_i)$ range from −0.84 to 1, and the average correlation is 0.18 for the CEAP data. The estimate of the measurement error variance is uniformly smaller than the estimate of the variance of the sampling error for the CEAP data. The ratios of the estimates of $V(u_i)$ to the estimates of $V(e_i)$ range from 0.0002 to 0.3436, and the average ratio is 0.090. Although the measurement error is smaller than the sampling error, naive application of Fay-Herriot procedures has the potential to under-estimate the MSPE.

For this analysis, we treat the direct estimator of ${\Psi }_i$ as fixed. An alternative is to use a generalized variance function to smooth the estimator of ${\Psi }_i$ . For the purpose of this analysis, we think that use of the direct estimator is more illuminating because it enables us to demonstrate the impacts of the structure of ${\Psi }_i$ on the efficiency of the predictor.

We assume that the model (1) holds for $(W_i,Y_i)$ for $i=1,\dots ,n$ . The parameter of interest ${\theta }_i$ represents the mean log RUSLE2 erosion for county $i$ , where $i=1,\dots ,99$ . We use the procedure of Section 2.2 to estimate the parameters of model (1). We construct predictors of ${\theta }_i$ and MSPE estimators using the method of Section 2.3.

Table 7 contains the estimate of $\omega =({\beta }_0,{\beta }_1,{\sigma }_b^2)$ for the CEAP data. The standard errors are the square roots of the diagonal elements of the jackknife covariance matrix defined as

Table 7.

Estimates of Fixed Model Parameters and Associated Standard Errors for CEAP Data.

	${\beta }_0$	${\beta }_1$	${\sigma }_b^2$
$\widehat{\omega }$	−2.541	1.053	0.283
SE	0.250	0.185	0.041

$${\widehat{V}}_{JK}(\widehat{\omega })=\sum _{k=1}^n\left({\widehat{\omega }}^{(k)}-n^{-1}\sum _{k=1}^n{\widehat{\omega }}^{(k)}\right){\left({\widehat{\omega }}^{(k)}-n^{-1}\sum _{k=1}^n{\widehat{\omega }}^{(k)}\right)}^{\prime }.$$

The magnitude of the estimate of each parameter is more than double the corresponding standard error.

For the CEAP data analysis, several of the estimated MSPE’s are negative. We therefore apply a lower bound (LB) to the estimated MSPE for the CEAP analysis. We define the MSPE estimator for the CEAP study by

$${\widehat{MSPE}}_{i,LB}={\widehat{MSPE}}_{i}I[{\widehat{MSPE}}_i>0]+(M_{1i}(\widehat{\omega })+{\widehat{M}}_{2i,JK})I[{\widehat{MSPE}}_i\le 0],$$

where ${\widehat{MSPE}}_i$ is defined in (9).

Figure 2 contains a scatterplot illuminating the relationship between the efficiency gain from prediction and the components of the design covariance matrix. The $x$ –axis of the plot contains the ratios of the estimates of $Var(u_i)$ to the estimates of $Var(e_i)$ . The $y$ –axis depicts the estimates of $Cor(u_i,e_i)$ . The $z$ –axis has the relative mean square prediction error, defined as the ratio of the estimated MSPE of the predictor to the estimated sampling variance of the direct estimator. A ratio below one means that the predictor is more efficient than the direct estimator.

Figure 2.

Three-dimensional scatterplot with $Var(u_i)/Var(e_i)$ on the x-axis, $Cor(u_i,e_i)$ on the y-axis, and the ratio of the estimated MSPE to $Var(e_i)$ on the z-axis.

From Figure 2, it is apparent that an efficiency gain is attained for most counties. The relative mean square prediction errors range from 0.056 to 1.222, and the average is 0.781. The efficiency gains are often pronounced when $Cor(u_i,e_i)$ and $Var(u_i)/Var(e_i)$ are both small. For counties where the estimated MSPE of the predictor exceeds the estimated variance of the direct estimator, $Cor(u_i,e_i)$ tends to be relatively high. For instance, for the county with a relative mean square error of 1.222, $Cor(u_i,e_i)=0.946$ . This mirrors the simulation results for the configuration with small measurement error variance and high correlation. Overall, the estimated efficiency gains from small area modeling are substantial.

The left panel of Figure 3 contains a scatterplot of the predictors on the vertical axis against the direct estimators on the horizontal axis. The line in the plot is the 45-degree line through the origin. Small area prediction has the expected shrinkage effect. Prediction increases direct estimators that are unusually low and decreases direct estimators that are unusually high.

Figure 3.

Left: Scatterplot of predictors against direct estimators. Right: Scatterplot of root MSPE estimators against square roots of estimated variances of direct estimators.

The right panel of Figure 3 contains a plot of the square roots of the mean square prediction errors against the standard errors of the direct estimators. The line is again the 45-degree line through the origin. For most counties, prediction renders an efficiency gain relative to the direct estimator. The reduction in MSPE from prediction is often substantial. When the MSPE exceeds the estimated variance of the direct estimator, the estimated loss of efficiency is minimal.

5. Discussion

We conduct an extensive study of the properties of a small area predictor that recognizes a correlation between the measurement error in the covariate and the sampling error in the response. The simulation studies illustrate the dangers of naively applying the Fay-Herriot predictor when the covariate and response are estimators from the same survey. The Fay-Herriot predictor can be less efficient than the direct estimator, and the corresponding MSPE estimator can have a severe negative bias for the MSPE of the predictor. The problems with the Fay-Herriot predictor persist even when ${\Psi }_i=\Psi$ for all $i=1,\dots ,n$ .

The proposed predictor rectifies the problems with the Fay-Herriot predictor and is more efficient than other alternatives considered in the simulations. In both the simulations and the data analysis, the efficiency of the proposed method relative to the direct estimator depends on the nature of ${\Psi }_i$ . In the CEAP study, runoff varies less within a county than does sheet and rill erosion, and substantial gains in efficiency from small area prediction are observed for most counties. Without the methodology presented in this paper, the use of runoff as a covariate would be impossible. Our methodology is of general interest beyond the CEAP study. We provide a theoretically sound estimation procedure for use in conjunction with the simple practice of using estimators from related surveys as covariates and response variables. The methodology has potential use in a wide range of applications in the area of official statistics.

Supplemental Material

sj-pdf-1-jof-10.1177_0282423X241240835 – Supplemental material for An Application of a Small Area Procedure with Correlation Between Measurement Error and Sampling Error to the Conservation Effects Assessment Project

Supplemental material, sj-pdf-1-jof-10.1177_0282423X241240835 for An Application of a Small Area Procedure with Correlation Between Measurement Error and Sampling Error to the Conservation Effects Assessment Project by Emily Berg and Sepideh Mosaferi in Journal of Official Statistics