A measurement error Rao–Yu model for regional prevalence estimation over time using uncertain data obtained from dependent survey estimates

The assessment of prevalence on regional levels is an important element of public health reporting. Since regional prevalence is rarely collected in registers, corresponding ﬁgures are often estimated via small area estimation using suitable health data. However, such data are frequently subject to uncertainty as values have been estimated from surveys. In that case, the method for prevalence estimation must explicitly account for data uncertainty to allow for reliable results. This can be achieved via measurement error models that introduce distribution assumptions on the noisy data. However, these methods usually require target and explanatory variable errors to be independent. This does not hold when data for both have been estimated from the same survey, which is sometimes the case in ofﬁcial statistics. If not accounted for, prevalence estimates can be severely biased. We propose a new measurement error model for regional prevalence estimation that is suitable for settings where target and explanatory variable errors are dependent. We derive empirical best predictors and demonstrate mean-squared error estimation. A maximum likelihood approach for model parameter estimation is presented. Simulation experiments are conducted to prove the effectiveness of the method. An application to regional hypertension prevalence estimation in Germany is provided.


Introduction
In recent years, the demand for detailed population health indicators has increased considerably. Policy-makers and researchers have realized the importance of monitoring disease frequencies on geographic regions and demographic groups over time. These insights mark important elements of public health reporting and are key to the provision of comprehensive health care programs. However, regional prevalence is rarely collected in registers due to confidentiality restrictions. Instead, such figures are often estimated via methods of small area estimation (SAE) using suitable health data. See Rao and Molina (2015) and Jiang and Lahiri (2006), or Pfeffermann (2002) and Pfeffermann (2013) for comprehensive overviews on SAE. Health data records are used as auxiliary variables in regression models that generate predictions for domain characteristics of interest. This approach allows for efficient estimates when the model has sufficient explanatory power (Münnich et al. 2019;Burgard et al. 2020b). Some recent papers in small area estimation for regional estimation are Trandafir et al. (2020), Bernal et al. (2020), Mills et al. (2020) and Saeedi et al. (2019), as well as Adin et al. (2019).
Due to the confidentiality of the data and the sensitive nature of health variables, area-level models are good options for predicting regional prevalence. Fay and Herriot (1979) proposed the original area-level model, which regresses regional direct estimates of target variables on covariates using cross-sectional data. This model was later extended and applied to multivariate data by Benavent and Morales (2016), Arima et al. (2017), Ubaidillah et al. (2019), Burgard et al. (2021) and Esteban et al. (2020). Since analyzing the development of domain characteristics over time is often relevant in empirical studies, generalizations of the area-level model to the temporal context soon appeared. Among them, the papers of Pfeffermann and Burck (1990), Rao and Yu (1994), Ghosh et al. (1996), Datta et al. (1999), You and Rao (2000) and Datta et al. (2002) provided generalizations of this model by considering time-series data across domains, multi-level geographic segmentation, Bayesian modeling and robust approaches. More recently, Esteban et al. (2012), Marhuenda et al. (2013), Morales et al. (2015), Jedrejczak and Kubacki (2017), Boubeta et al. (2017) and Benavent and Morales (2021) applied temporal area-level models to poverty mapping.
In SAE models, health variables are often subject to uncertainty. Their values are typically estimated from periodic large-scale health surveys, such as the National Health and Nutrition Examination Survey (NHANES) in the US (National Center for Health Statistics 2020) or Gesundheit in Deutschland aktuell (GEDA) in Germany (Robert Koch Institute 2012; Lange et al. 2015). Accordingly, the data records are not accurate, but associated with errors. In that case, the regression models must explicitly account for data uncertainty to allow for reliable results. These issues motivated researchers over the years to develop SAE methods that allow for uncertain data in both dependent and explanatory variables. Such SAE methods are often called measurement error models. The basic idea of these methods is to introduce distribution assumptions on the errors and make statistical inference under this premise. On that note, Ybarra and Lohr (2008), Burgard et al. (2020a) and Burgard et al. (2021) provided extensions of the area-level model that explicitly account for explanatory variable errors. A related approach was developed by Torabi et al. (2009) for unit-level observations. However, despite accounting for errors on both target and explanatory variables, these methods are still limited in their applicability. They require the target variable errors to be independent from the explanatory variable errors.
This condition is not fulfilled when some explanatory variables are estimated from the same survey as the target variable, which may occur in practice. That is to say, if we do not have enough data from administrative records, the alternative of using data from other surveys or from the same survey allows proposing new model-based estimates that may improve the direct estimates.
For prevalence estimation, the above-described scenario is relevant. Explanatory variables for health mapping are often health-related indicators as well. Since surveys like NHANES or GEDA are primary data sources for public health reporting in their respective countries, it may happen that both target variable and some explanatory variables are estimated from them. If the resulting dependencies between target and explanatory variable errors are ignored, prevalence estimates based on these records may be biased or inefficient.
We propose a new measurement error area-level model for regional prevalence estimation that allows for dependent errors on target and explanatory variables. With this feature, it is the first SAE method that is suitable for settings where data records for both model components have been estimated on a single survey. The model is an extension of the Rao-Yu model proposed by Rao and Yu (1994), which implies that it allows for either regional time-series data (domains and periods), or nested geographic segmentation (domains and subdomains). Thus, we refer to it as the measurement error Rao-Yu (MERY) model. We use it to estimate the development of regional prevalence over time based on uncertain data with dependent errors. Yet, we also show that the MERY model contains many existing measurement error SAE models as particular or limits cases, which makes it applicable to a wide range of empirical settings.
We first present a detailed description of the MERY model. A special emphasis is placed on the conditional component distributions, as the main innovation of the method is to allow for dependent errors on the variables. Next, the best predictors (BPs) and empirical best predictors (EBPs) of the domain characteristics under the model are derived. We subsequently address mean-squared error (MSE) estimation. The MSE of the BP is derived analytically. We further present a set of parametric bootstrap estimators for the MSE of the EBP. A maximum likelihood (ML) approach based on the Fisher scoring algorithm (Jennrich and Sampson 1976;Longford 1987) for model parameter estimation is stated. The methodology is tested in simulation experiments to demonstrate its effectiveness. An empirical application to temporal health mapping in Germany is provided. We estimate the hypertension prevalence for Nielsen regions crossed by age groups (domains), over three periods, using survey estimates from GEDA for the years 2009(Robert Koch Institute 2013, 2014a. We remark that the domains do not correspond to any geographic segmentation in administrative health records of German statistical offices. Therefore, it is not possible to apply the standard SAE methodology based on area-level models, where the auxiliary variables are taken from out-of-sample data sources. In addition, a quarter of the domains have a sample size of less than 271 and the minimum sample size is equal to 51. Therefore, the estimation of the hypertension prevalence by domains and time periods can be treated by using SAE statistical methodology.
We show that applying the new MERY model in this situation allows for considerable efficiency gains relative to direct estimation.
The remainder of the paper is organized as follows. Section 2 introduces the MERY model and derives the BPs. Section 3 states the MSE of the BP and gives parametric bootstrap estimators for the MSE of the EBP. Section 4 addresses model parameter estimation. Section 5 contains the simulation experiments. Section 6 presents the empirical application to temporal health mapping. Section 7 closes with some conclusive remarks. A supplemental material file with five appendices is provided. Appendix 1 contains the proof for the main theorem. Appendix 2 presents several MSE calculations. Appendix 3 derives the Fisher information matrix required for the ML approach for model parameter estimation. Finally, Appendices 4 and 5 show additional results for simulation and application.

Formulation
Let U be a finite population that is segmented into D disjoint sets U d , d = 1, . . . , D. Suppose that every set U d is further segmented into disjoint subsets U dt , d = 1, . . . , D, t = 1, . . . , m d . For ease of exposition, we assume that the number of subsets per set is constant with m d = T , d = 1, . . . , D. This is to say, U = ∪ D d=1 ∪ T t=1 U dt . In practice, the structure may correspond to either a population with a two-stage hierarchical geographic segmentation (domain-subdomain), or a population with geographic segmentation that is observed over several time periods (domain-period). This section introduces the MERY model by using the domain-subdomain terminology. The domain-period interpretation is straightforward. Let μ dt be a characteristic of interest in set d and subset t, for instance the mean or the total of some target variable. The objective is to estimate μ dt for all d = 1, . . . , D, t = 1, . . . , T . Suppose that for each μ dt a corresponding direct estimator y dt is available from a target survey S ⊂ U . By analogy with population structure, we assume that The MERY model is composed of three stages. In the first stage a model, called sampling model, is used to represent the sampling error of direct estimators. The sampling model indicates that direct estimators, y dt 's, are unbiased and can be expressed as where the sampling errors, e dt 's, are independent and normally distributed with known variances. This is to say, e dt ∼ N (0, σ 2 edt ), where σ 2 edt is known. In the second stage, called linking model, the true domain characteristics, μ dt 's, are assumed to vary linearly with a number p of domain-level auxiliary variables. The linking model is wherex 0,dt andx 1,dt are a 1 × p 0 and 1 × p 1 row vectors, respectively, containing the true aggregated (subdomain-level) values of p = p 0 + p 1 auxiliary variables for domain d. The terms β 0 and β 1 are the corresponding column vectors of regression coefficients. The u 1,d ' s and u 2,dt 's are model errors (domain and subdomain random effects), assumed to be independent and identically distributed (i.i.d.) from N (0, σ 2 1 ) and N (0, σ 2 2 ), respectively, with variance σ 2 1 and σ 2 2 unknown, and independent of the e dt 's. We assume thatx 0,dt is known and equal to x 0,dt , d = 1, . . . , D, t = 1, . . . , T . Further, we assume that thex 1,dt 's are unknown random vectors that are predicted from the same survey as μ dt or from different surveys. This is to say, the direct estimators y dt and x 1,dt , of μ dt andx 1,dt , respectively, can be calculated from the same or from a different survey sample. In the third stage, we consider the measurement error model where x 1,dt is a row vector containing the direct estimators ofx 1,dt . The term v dt is a column vector of random measurement errors, vdt is the known covariance matrix of v dt (typically taken as a design-based covariance matrix estimate) and v 11 , . . . , v DT , x 1,11 , . . . , x 1,DT are independent. Model (1) assumes that the vector of true aggregated valuesx 1,dt is random and it is the sum of two independent random vectors. The first one contains the predictions of the realizations ofx 1,dt and the second one is a random error term. Further, the model assumes that the joint distribution of the ( p 1 + 1) × 1 column vector ε dt = (e dt , v dt ) is multivariate normal with zero mean and covariance matrix where σ 12dt is a p 1 × 1 vector containing the covariances between e dt and the components of v dt . The covariances of σ 12dt , associated to auxiliary variables calculated from a different sample, are zero. The MERY model can be expressed as single model in the form (2) or in the simpler form y dt = x dt β + β 1 v dt + u 1,d + u 2,dt + e dt , d = 1, . . . , D, t = 1, . . . , T , where x dt = (x 0,dt , x 1,dt ) and β = (β 0 , β 1 ) . We finally assume that x 1,dt , u 1,d , u 2,dt , ε dt = (e dt , v dt ) , d = 1, . . . , D, t = 1, . . . , T , are mutually independent, but we only introduce inference procedures conditionally on (x 11 , . . . , x DT ). Under Model (2), it holds that var(β 1 v dt ) = β 1 vdt β 1 , cov(β 1 v dt , e dt ) = β 1 σ 12dt .
The variance of y dt , given x dt , is Thus, the distribution of y dt , given x dt , is y dt | x dt ∼ N x dt β, σ 2 ydt . Let us define the vectors where 1 T is a T × 1 vector with all elements equal to 1 and I T is the T × T identity matrix. At the domain level, Model (2) is Define M = DT , the vectors y = col 1≤d≤D ( y d ), and matrices V 1 = σ 2 1 diag In matrix form, Model (2) can be alternatively written as The variance matrix of y, given X, is Hereafter, we state some further conditional distributions (CDs) resulting under the model. They allow for a clearer picture of the dependencies between specific model components, as the main innovation of our contribution is to allow for SAE from dependent sample estimates. Recall that if y is normally distributed with partitioned mean vector μ = (μ s , μ r ) , where μ s and μ r are partitions of μ, as well as covariance matrix then the distribution of y r , given y s , is also multivariate normal, i.e., y r | y s ∼ N (μ r |s , V r |s ), with mean vector and covariance matrix At the subdomain level, we have the following CDs. The CD of e dt , given v dt , is N (μ e|vdt , σ 2 e|vdt ), where The CD of v dt , given e dt , is N p (μ v|edt , v|edt ), where The conditional expectation and variance of y dt , given x dt and v dt , are The CD of y dt , given The CD of y dt , given x dt and u 1,d , is The CD of y dt , given x dt and u 2,dt , is At the domain level, we have the following CDs. The CD of y d , given The CD of y d , given X d and u 1,d , is The CD of y d , given X d and u 2,d , is

Empirical best prediction
In this section, we obtain the EBP of the μ dt 's under the MERY model. For this, we first derive the respective BP under the assumption of known model parameters. Thereafter, the EBP is obtained from substituting the model parameters with consistent estimators. The BP is calculated according to the subsequent theorem.
The proof can be found in Appendix 1 of the supplemental material. Define the full . Please note that we describe how to obtainθ in Sect. 4. The MERY model is quite general. It allows for explanatory variables with no errors, errors that are independent of the sampling errors (out-of-sample errors) and errors that depend on the sampling errors (in-sample errors). If x dtk is measured without error, we have σ 2 dtk = var(x dtk ) = 0 and zero covariances. For the case that all explanatory variables are measured without error, that is, vdt = 0 p , which implies that ε dt = (e dt , 0) , d = 1, . . . , D, t = 1, . . . , T , thenμ ebp d corresponds to the Rao-Yu predictor proposed by Rao and Yu (1994). When x dtk is an out-of-sample direct estimator, then all its covariances with in-sample direct estimators are zero. For σ 2 1 = 0 and σ 12dt = 0, d = 1, . . . , D, t = 1, . . . , T , then μ ebp d is equivalent to the EBP under the measurement error Fay-Herriot model studied by Ybarra and Lohr (2008) as well as Burgard et al. (2020a).

Mean-squared error estimation
In this section, we elaborate on MSE estimation for the EBP under the MERY model. Let us first decompose the MSE into components. By summing and subtractingμ bp d , we obtain By taking expectations, we get In what follows, we present two parametric bootstrap estimators of M SE(μ ebp d ). The first is a standard estimator that approximates all four MSE components via the Monte Carlo method. The second estimator is a bias-corrected estimator that uses an analytical approximation for G 1d (θ). We will show via simulation in Sect. 5.4 that including the analytical formula improves MSE estimation over the standard estimator considerably.

Standard parametric bootstrap estimator
This section presents the standard estimator, denoted by mse * 1d , which approximates all four MSE components numerically. It is calculated with the subsequent algorithm.

Bias-corrected parametric bootstrap estimator
This section presents the bias-corrected estimator, denoted by mse * 2d , that uses an analytical approximation for G 1d (θ ). Observe that G 1d (θ) denotes the MSE of the BP in the course of the resampling process. Hence, we first derive the analytical expression for M SE(μ bp d |X d ) under the assumption of known model parameters θ . Afterward, we use this expression in another parametric bootstrap algorithm that yields us the bias-corrected estimator. Please note that we only present the final results of the corresponding mathematical developments. For further details regarding the derivation, see Appendix 2 of the supplemental material. Under model (2), the mean-squared error, where the first summand is and the third summand is , where the known model parameters θ are substituted by the estimateθ . The bias-corrected estimator is obtained as follows.
1. By using the data ( y d , which can be approximated by the Monte Carlo method, i.e., Observe that the difference between mse * 1d and mse * 2d lies within the approximation of G 1d (θ). While mse * 1d uses a purely numerical approximation over the Monte Carlo iterations, mse * 2d uses the bias-corrected term 2G 1d (θ )− H 1d * . Here, both components G 1d (θ ) and H 1d * are calculated based on the analytical formula for M SE(μ bp d |X d ). This allows for an efficiency advantage over the standard estimator.
Under the Fay-Herriot model, this fact was proved in Theorem 3 of González-Manteiga et al. (2008). These authors showed that the biases of the parametric bootstrap estimators mse * 1d and mse * 2d are of orders O(D −1/2 ) and o(D −1 ), respectively, under some model assumptions. The extension of their results to the MERY model is not straightforward, but similar asymptotic results may be expected under the MERY model after adding regularity conditions on the matrices vdt , d = 1, . . . , D, t = 1, . . . , T . Simulations of Sect. 5.4 cover the lack of theoretical results and gives some information about the behavior of mse * 1d and mse * 2d under the MERY model. Finally, note that the properties of the bootstrap estimators mse * 1d and mse * 2d rely on parametric distributional assumptions. Deviations from the normal distribution may produce biases in the MSE estimates. Fortunately, in the application to real data the sample sizes are not very small. Therefore, the asymptotic normality of the direct estimators, y dt and x 1,dt , does not allow important deviations from normality.

Model parameter estimation
This section presents a ML approach for model parameter estimation. It is based on the Fisher scoring algorithm (Jennrich and Sampson 1976;Longford 1987). Recall that For calculating V −1 d , we apply the inversion formula For calculating the derivatives of V d , we define the p × 1 vector For k = 1, . . . , p, the derivative of β, with respect to β k , is ∂β/∂β k = δ k . For k = 1, . . . , p 1 , the derivative of V d , with respect to β k , is ∂ V d /∂β k = 0. For k = p 1 + 1, . . . , p, we have the following derivatives In what follows, we apply the formulas The first partial derivatives of are The second partial derivatives and the elements of the Fisher information matrix can be found in Appendix 3 of the supplemental material. Let θ = (β , σ 2 1 , σ 2 2 ) denote the full parameter. In matrix form, the vector of scores and the Fisher information matrix are The updating formula of the Fisher scoring algorithm, for calculating the ML estimators of the MERY model parameters, is Please note that, depending on the data setting, Fisher scoring algorithms can be numerically instable and sensitive regarding starting values. Thus, we recommend using θ (0) =θ r y as seed, whereθ r y is the model parameter estimate obtained under the Rao-Yu model.

Simulation experiments
This section presents simulation experiments that are conducted in order to demonstrate the effectiveness of the method in a controlled environment. We investigate three aspects: (i) model parameter estimation, (ii) characteristic prediction, as well as (iii) mean-squared error estimation. Note that additional simulation results can be found in Appendix 4.

Results of model parameter estimation
This section studies the behavior of the ML fitting algorithm presented in Sect. 4. We use the root mean-squared error (RMSE) and the bias as performance measures. For τ ∈ {β 0 , β 1 , β 2 , σ 2 1 , σ 2 2 }, they are given by Since β 1 = (2, 2) , we can average the performance measures for the regression parameters. We start with the RMSE, which is displayed in In the B-and C-scenarios, the MERY model is clearly more efficient as it uses the additional information on the measurement error distributions. The biggest performance differences are in variance parameter estimation. Here, the RY model cannot distinguish between random effect variation and variation caused by the measurement errors. This leads to highly inefficient subdomain variance parameter estimates. The MERY model, on the other hand, shows a very stable estimation performance, as it can distinguish between the individual sources of variation.
We continue with the bias. The results are summarized in Table 3. For regression parameter estimation, the results are mixed. The MERY model is less biased with respect to intercept estimation, but often more biased in slope parameter estimation. However, the relative bias of slope parameter estimates is at most 2.6%, which is a reasonable magnitude under measurement error settings. Regarding variance parameter estimation, the MERY model is considerably less biased than the RY model. Especially on the subdomain level, which is also the level of the measurement errors, the RY model vastly overestimates the random effect variance. This is due to its incapability to distinguish between measurement errors and random effect variation. The MERY model, on the other hand, shows a much weaker tendency of overestimation because it uses the additional distribution information. The remaining bias is likely due to the likelihood function being flat under measurement errors.

Results of subdomain characteristic prediction
This section studies the performance of the EBPμ ebp dt based on the MERY model, derived in Sect. 2.2, as well as the EBLUPμ ebp0 dt based on the RY model. We consider RMSE, RRMSE, absolute bias and relative absolute bias as measures. Forμ For a compact presentation, we average the performance measures over the subdomains. The averaged results are presented in Table 4. We observe that the EBP based on the MERY model is superior to the EBLUP based on the RY model in all measures and scenarios. In the A-scenarios, the performance difference is very small, which is consistent with the theoretical developments from Sect. 2. In the B-and C-scenarios, the MERY model shows more distinct efficiency advantages. The inclusion of the additional information from the subdomain covariance matrices dt in combination with the more efficient model parameter estimates allows for more accurate results. Further, we see that the predictions under the MERY are also less biased than under the RY model in terms of the absolute bias.

Results of mean-squared error estimation
This section studies the parametric bootstrap procedures for MSE estimation in Sect. 3. In particular, we look at two estimators mse * 1 and mse * 2 , where the first is the standard parametric bootstrap estimator (Sect. 3.1) and the second is the bias-corrected estimator (Sect. 3.2). We consider bias, relative bias and RMSE as performance measures: Note that mse dt , d = 1, . . . , D, t = 1, . . . , T , is the MSE of the EBP over Monte Carlo iterations, which we take as approximation for the real value M SE(μ dt ). For a The results are summarized in Table 5. Note that the columns M SE and mse * correspond to the means of mse dt and mse * dt over all subdomains of the respective scenario. From Table 5, we observe that in the A-and B-scenarios, the performance of both MSE estimators is decent. The relative bias ranges between 0.1 and 7.2%, which is a reasonable magnitude for a mixed model with a multi-level random effect structure. However, we see that the bias of the MSE estimator mse * 2 is always smaller than for the standard estimator mse * 1 except for Scenario A.5. Further, by looking at the RMSE, it can be concluded that the second MSE estimator is more efficient than the first. In the C-scenarios, which resemble large measurement errors, the standard estimator mse * 1 become visibly biased. Its relative bias ranges from 10.0 to 24.6%. The bias-corrected estimator mse * 2 , on the other hand, is not as much affected by the large errors. Its relative bias only ranges from 0.6 1.9%. In addition to that, by looking at the RMSE, a significant efficiency advantage by the second estimator becomes evident. Accordingly, including the analytical approximation to the first MSE component improves the estimation accuracy considerably in these settings. Figure 1 shows the convergence behavior of mse * 2 . In particular, it shows the distribution of mse * 2dt (B)−mse * 2dt (500), d = 1, . . . , D, t = 1, . . . , T , when the parametric bootstrap estimate mse * 2dt (B) is calculated based on B = 1, . . . , 500 replicates. The light gray band displays the inner 90% of differences over all subdomains, while the gray band marks the inner 70%. We see that the mean stabilizes at around B = 150, implying that 150 bootstrap replicates are sufficient to obtain decent MSE estimates in terms of overall performance. However, if we require that the inner 90% differences vary less than 1% of M SE(μ), which would roughly be ±0.48 in Scenario 1, B = 300 is an appropriate choice.

Application
This section presents an empirical application to health mapping on the example of regional age-referenced hypertension prevalence in Germany. Note that further results are located in Appendix 5 of the supplemental material.

Data and specification
We use the health survey GEDA (Robert Koch Institute 2013, 2014a, b) from 2009 (n 1 = 21,262), 2010 (n 2 = 22,050) and 2012 (n 3 = 19,294) as data basis. It contains citizens of age 18+ from the German population that are interviewed via CATI in a representative telephone sample. The data sets contain health-related information of participants on the unit level. For further information on the survey as well as its sampling design and response rates, see Robert Koch Institute (2013) and Robert Koch Institute (2014a, b).
We define the domains as cross combinations of the seven Nielsen regions of Germany and the age groups I: 18-29, II: 30-39, III: 40-49, IV: 50-59, V: 60-69, VI: 70-79 and VII: 80+. Accordingly, we have D = 49 domains that are observed in T = 3, periods, which implies M = 147 domain-time sets in total. The Nielsen population structure does not correspond to the one used by the German statistical offices, and therefore, it is difficult to collect aggregated auxiliary variables to apply the standard SAE approach based on area-level models. However, we can obtain auxiliary data from the health survey GEDA and apply the statistical methodology based on the MERY model. For this sake, we calculate a set of Horvitz-Thompson estimators for the regional hypertension prevalence as target variable and other health indicators that we use as explanatory variables. Further, we exploit the fact that the hierarchical structure of the MERY model in Sect. 2.1 can be interpreted as a collection of domains that is observed over multiple time periods.
The objective is to estimate the hypertension prevalence per each domain d and time period t. That is to say, we define the number of individuals with high blood pressure as the model component μ dt , d = 1, . . . , D, t = 1, . . . , T , in the sense of Sect. 2.1, and divide by the known domain-time size N dt afterward. In order to predict the total of hypertensive people, we use the estimated totals of people with diabetes and increased blood lipids per domain and time period as covariates,x 1,dt = (x 1,1dt ,x 1,2dt ) in the notation of Sect. 2.1, since these diagnoses are well known to be associated with hypertension (see, e.g., Lastra et al. 2014). To determine whether a person suffers from any of these medical conditions in the sample, we rely on the 12-month prevalence profiles used by Robert Koch Institute (2012).
The main advantage of the new MERY model is that it allows for SAE based on uncertain data obtained from dependent survey estimates. In order to provide a suitable data basis for our method, we calculate Horvitz-Thompson estimators for the number of citizens per domain-time in the population that are associated with any of the three medical conditions. These results yield us the model components y dt and x 1,dt that are direct estimators of μ dt andx 1,dt , respectively. Note that the estimators are obtained from the same set of sample observations per domain and time period. In standard measurement error models for SAE, this would not be possible as y dt and x 1,dt are not allowed to have a covariances different from zero. The MERY model, on the other hand, allows for covariances between the direct estimators. In fact, one could argue that it even encourages dependency between the errors, as the covariances between y dt and x 1,dt mark additional information that can be used to improve the prediction. Recall that in Sect. 2.1, normality was assumed for these components. For the application, the applied Horvitz-Thompson estimators satisfy this assumption asymptotically, as demonstrated by Hájek (1960), Berger (1998) and Chen and Rao (2007). We further calculate direct estimators of the variances and covariances of y dt and x 1,dt . In this way, we construct the matrices dt and vdt for every domain and time period, in accordance with Sect. 2.1. This can be done using the unit-level observations in GEDA. We calculate direct estimators of the variances var(y dt ) and var(x 1, jdt ), j = 1, 2, based on the unit-level survey data. For j = 1, 2, the estimators of the covariances between the direct target and covariate estimators are obtained from Wood (2008) cov(y dt , x 1, jdt ) = 1 2 (var(y dt ) + var(x 1, jdt ) − var(y dt − x 1, jdt ) , The estimator of the covariance between the direct covariate estimators x 1,1dt and x 1,2dt is obtained analogously. Table 6 presents the quartiles, across domains and time periods, of the sample sizes n dt , the coefficients of variation, cv(y dt ), cv(x 1,1dt ) and cv(x 1,2dt ), of the direct estimators of the totals of the objective variable (hypertension) and of the two auxiliary variables (diabetes and lipids), and the corresponding correlation coefficients ρ(y dt , x 1,1dt ), ρ(y dt , x 1,2dt ), ρ(x 1,1dt , x 1,2dt ). This table shows that some domaintime sample sizes are small and that the direct estimators of the totals of the considered medical variables are correlated. The introduced MERY model is the first one that can deal with this situation.
The full specification of the MERY model is where y dt , x 1,1dt and x 1,2dt are the Horvitz-Thompson estimators of the number of individuals with high blood pressure, diabetes and increased blood lipids per domain d and time period t, respectively. We apply the ML algorithm from Sect. 4 for model parameter estimation. Further, we use the subsequent parametric bootstrap procedure in order to approximate the standard deviation of the model parameter estimates in order to evaluate their significations.
1. By using the data ( y d , 3. For τ ∈ {β 0 , β 1 , β 1 , σ 2 1 , σ 2 2 }, calculate the bootstrap estimator of sd(τ ); i.e., In accordance with Sect. 5.4, we rely on B = 300 bootstrap replicates. For significance testing, we use a t-test with test statistic t(τ ) =τ /sd * (τ ) that is asymptotically distributed as standard normal. The resulting values are located in the probability density function of the standard normal to find the respective p values. Confidence intervals for the model parameter estimates are calculated byτ ± t (m,1−α/2) sd * (τ ), where t (·) is the corresponding quantile of the t-distribution with m = DT − 3 degrees of freedom and significance level α. Please note that the mathematically correct test distribution would be the normal distribution, as maximum likelihood estimates are asymptotically normal distributed. However, the t-distribution provides slightly more conservative test decisions and confidence intervals. This seems reasonable to us in a small area setting, where asymptotic arguments, due to small sample sizes, are critical. In this particular application, the confidence intervals are almost the same, as the number of degrees of freedom is considerably large, such that the corresponding t-distribution already approximates the normal distribution quite well. Table 7 displays the results of model parameter estimation obtained from the ML algorithm as well as the parametric bootstrap procedure. For the sake of presentability, we show the random effect standard deviations,σ 1 andσ 2 , rather than the respective variances. It can be seen that all model parameters except for the intercept are significant on a 1% level. The intercept is not significantly different from zero. The signs of the regression parameters β 1 and β 2 are plausible, as diabetes and increased blood lipids are positively correlated with hypertension. Furthermore, the domain and domain-time-level random effects are clearly visible in the distribution of y dt .

Results
The estimated prevalence for 2009 over all domains is 26.1%. The estimated prevalence for 2010 is 26.4%, and for 2012, we have 27.1%. The estimated evolution of the hypertension prevalence over age groups for all three considered periods is displayed in Fig. 2. We see that the hypertension prevalence is low in the youngest cohorts (18-39). From age 40-79, a steep monotonous increase is evident for all three periods. For individuals of age 80+, the curve flattens on a high level in 2010 and 2012. An interesting observation is that in 2009, the prevalence for the eldest cohort is slightly lower than for the second eldest. A possible explanation is that for age 80+, individuals without hypertension have higher life expectancy than individuals with hypertension Empirical best predictionsμ dt /N dt and direct estimates y dt /N dt due to a lesser risk of suffering a stroke or other cardiovascular events (Kassai et al. 2005). Figure 3 shows the empirical best predictions dependent on the respective direct estimates. The left plot displays the point estimates. The right plot marks the comparison of the EBP RMSEs and the direct estimator standard deviations. We have used the bias-corrected parametric bootstrap estimator presented in Sect. 3.2, to estimate the EBP RMSEs of the domain-time characteristic predictions.
From the plot, we see that the direct estimates vary randomly around the prevalence predictions. This is in line with theory, as the EBP based on the MERY model exploits the design unbiasedness of the direct estimators with the error components being zero on expectation. Furthermore, it indicates that we do not have systematic prediction bias, which is further supported by the simulation results in Sect. 5.3. From the right plot, we see that the MERY predictions are always more efficient than the direct estimates. The EBP RMSEs are always smaller than the corresponding standard deviations of the direct estimators. Accordingly, we can conclude that the application of the MERY model allows for a considerable gain in accuracy. On the other hand, we also remark that most coefficients of variation of the direct estimators are below 20%, since the sample sizes in several domains are quite high. For some Statistical Offices, 20% marks a quality threshold that admits a publication of the results. However, the MERY still manages to improve these results uniformly, which is why its application is recommendable in this setting. Figure 4 displays the estimated regional hypertension prevalence over the seven Nielsen regions in percent for 2010. The figures are with respect to the total regional populations (without age referencing). We see that the highest prevalences are located in the east of Germany, which is the former territory of the German Democratic Republic. The lowest prevalences are located in Badem-Württemberg and Bavaria. The estimated distribution is plausible, as similar patterns have been found for hypertension for instance by Burgard et al. (2020b). Further, similar distributions are known for strongly correlated conditions, such as diabetes (Schipf et al. 2014).

Conclusions
This paper introduces a new small area model for the estimation of regional population characteristics that allows for explanatory variables with measurement errors that are possibly correlated with the errors of the target variable. The approach enriches the modeling process by admitting auxiliary information that comes from administrative records, from survey samples other than the target survey, or even from that same target survey. The MERY model is a generalization of the Rao-Yu model that can be applied to both temporal data (domains and periods) or geographically nested data (domains and subdomains). With these features, it is broadly applicable in practice and marks an important contribution to both the SAE and the measurement error model literature.
Concerning model parameter estimation, the paper proposes the maximum likelihood method by implementing a Fisher scoring algorithm. For predicting the measurement errors, the random effects, as well as the target population quantities, the paper derives best predictors and calculates the corresponding EBPs. After calculating the variances and covariances of the random effects and errors, the MSEs of the BPs are approximated and applied to construct parametric bootstrap procedures with bias correction term. Several simulation experiments are performed to study the behavior of the fitting algorithm, the EBPs and the parametric bootstrap estimators of the MSEs. In particular, the efficiency gain derived from the modeling of measurement errors and their correlations is shown and a recommendation about the number of bootstrap replicates is given. An interesting finding is that if the measurements errors are ignored, then the MSEs of predictors are underestimated, which leads to the predictions seeming more accurate than they actually are (Appendix 4 of supplemental material). The new statistical methodology is applied to the health survey GEDA from 2009, 2010 and 2012. The estimation target is the regional hypertension prevalence in Germany. In order to predict hypertension, the employed covariates, diabetes and increased blood lipids prevalence, are also taken from the health survey GEDA. The new MERY model is the only model in the statistical literature that admits this possibility. Finally, the obtained EBPs are shown to be more precise than the corresponding direct estimates.