Parametric Bootstrap Confidence Intervals for the Multivariate Fay-Herriot Model

06/26/2020 ∙ by Takumi Saegusa, et al. ∙ 0

The multivariate Fay-Herriot model is quite effective in combining information through correlations among small area survey estimates of related variables or historical survey estimates of the same variable or both. Though the literature on small area estimation is already very rich, construction of second-order efficient confidence intervals from multivariate models have so far received very little attention. In this paper, we develop a parametric bootstrap method for constructing a second-order efficient confidence interval for a general linear combination of small area means using the multivariate Fay-Herriot normal model. The proposed parametric bootstrap method replaces difficult and tedious analytical derivations by the power of efficient algorithm and high speed computer. Moreover, the proposed method is more versatile than the analytical method because the parametric bootstrap method can be easily applied to any method of model parameter estimation and any specific structure of the variance-covariance matrix of the multivariate Fay-Herriot model avoiding all the cumbersome and time-consuming calculations required in the analytical method. We apply our proposed methodology in constructing confidence intervals for the median income of four-person families for the fifty states and the District of Columbia in the United States. Our data analysis demonstrates that the proposed parametric bootstrap method generally provides much shorter confidence intervals compared to the corresponding traditional direct method. Moreover, the confidence intervals obtained from the multivariate model is generally shorter than the corresponding univariate model indicating the potential advantage of exploiting correlations of median income of four-person families with median incomes of three and five person families.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

For the last few decades, there has been an increasing demand to produce reliable estimates for small geographic areas, commonly referred to as small areas, since such estimates are routinely used for fund allocation and regional planning. The primary data, usually a survey data, are usually too sparse to produce reliable direct small area estimates that use data from the small area under consideration. To improve upon direct estimates, different small area estimation techniques that use multi-level models to combine information from relevant auxiliary data have been proposed in the literature. The readers are referred to Jiang and Lahiri (2006) and Rao and Molina (2015) for a comprehensive review of small area estimation.

In estimating per-capita income of small places (population less than 1000), Fay and Herriot (1979)

proposed an empirical Bayes method to improve on direct survey-weighted estimates by borrowing strength from administrative data and survey estimates from a bigger area. Their method uses a two-level normal model in which the first level captures the variability of the survey estimates and the second level links the true small area means to aggregate statistics from administrative records and survey estimates for a bigger area. Researchers working on small area estimation have found the Fay-Herriot model useful in investigating various theoretical properties as well as implementing methodology in different applied problems when we do not have access to micro-data because of confidentiality and other reasons. For a review on the Fay-Herriot model and the related empirical best predictions, readers are referred to

Lahiri (2003b).

Following the pioneering paper by Fay and Herriot (1979), several multivariate extensions of the Fay-Herriot model have been considered to combine information from small area estimates of related variables or from past small area estimates of the same variable or both. They are essentially special cases of the general multivariate random effects or two-level multivariate model. In the context of estimating median income of four-person families for the fifty states and the District of Columbia (small areas), Fay (1987)

suggested a multivariate extension of the Fay-Herriot model, commonly referred to as the multivariate Fay-Herriot model, in order to borrow strength from the corresponding survey estimates of median income of three and five person families for the small areas. Alternatively, in Fay’s setting one could think of using survey estimates of median income for the three-person and five-person families as auxiliary variables in an univariate Fay-Herriot model. But, unlike the univariate Fay-Herriot model, the multivariate Fay-Herriot model incorporates sampling variance-covariance matrix of direct survey estimates of median income of the 3-, 4- and 5-person families for each small area. Moreover, the multivariate Fay-Herriot model borrows strengths through correlations of the components of area specific vector of random effects associated with the true median income of 3-, 4-, and 5- person families. Inferences on the four-person median income for the small areas drawn from the multivariate Fay-Herriot model are expected to be more efficient and reasonable when compared to the inferences drawn from an univariate Fay-Herriot model with survey estimates of the median income for three and five person families as auxiliary variables. This is because the univariate Fay-Herriot model would ignore the sampling variability of the survey estimates of median income for the three and five person families in the small areas.

In estimating median income of four-person families for the fifty states and the District of Columbia, Datta et al. (1991) used a bivariate Fay-Herriot model with a general structure for the variance-covariance matrix of the vector of area specific random effects. However, in many small area applications, structured variance-covariance matrices for the vector of area specific random effects arise naturally. For example, in order to combine information from the related past data, Rao and Yu (1994) proposed a stationary time series cross-sectional model while Datta et al. (2002) proposed a random walk time series and cross-sectional model. Although the time series cross-sectional models can be viewed as special cases of the multivariate Fay-Herriot model, one can achieve greater efficiency in estimating the unknown variance covariance matrix by reducing the number of parameters in the variance-covariance matrix through time series cross-sectional models.

Empirical best linear unbiased predictions and associated uncertainty measures for multivariate Fay-Herriot models with or without structured variance-covariance matrices for the vector of random effects have been adequately studied; see, e.g., Datta et al. (1991), Rao and Yu (1994), Benavent and Morales (2016), Datta et al. (2002), and others. However, the problem of constructing second-order efficient confidence intervals for the multivariate Fay-Herriot model, i.e., confidence intervals with coverage error being the number of small areas, received very little attention. Datta et al. (2002) obtained a second-order efficient confidence interval for a small area mean using an analytical method. To this end, they first obtained the exact expression for the term of order

in a higher order expansion of coverage probability of a normality-based empirical Bayes confidence interval, originally proposed by

Cox (1975), and then, using the term in the expansion, suggested an adjustment to the normal percentile in order to lower the coverage error to . The approach of Ito and Kubokawa (2021) in obtaining a second-order efficient confidence region for the vector of means for each area is essentially a multivariate generalization of Datta et al. (2002)

. However, their results are specifically designed for the multivariate Fay-Herriot model with an unstructured variance-covariance matrix when a method-of-moment estimator of the variance-covariance matrix of the vector of random effect is used. The derivation of the second-order efficient confidence intervals by the analytical method of

Datta et al. (2002) or Ito and Kubokawa (2021) is cumbersome and one needs to go through the such derivation each time one changes the model (say, a multivariate Fay-Herriot model with model variance-covariance structure suggested by the time series cross-sectional model of Rao and Yu (1994) or Datta et al. (2002)) or estimation method for the model parameters.

Parametric bootstrap method for obtaining second-order unbiased mean squared error estimation was first proposed by Butar and Lahiri (2002)

. Construction of the second-order efficient confidence interval based on the empirical best linear predictor of a small area parameter for a general linear mixed model was proposed by

Chatterjee et al. (2008). For parametric bootstrap confidence intervals for the univariate Fay-Herriot model, see Lahiri (2003a) and Li and Lahiri (2010). In this paper, we develop a parametric bootstrap method for obtaining second-order efficient confidence intervals for small area parameters from a multivariate Fay-Herriot model. Compared to the analytical method, our parametric bootstrap approach for constructing second-order confidence intervals for small area parameters is versatile and theoretically complete because our method applies to any variance estimator with minimal assumptions and theoretical justification is directly provided to the proposed method.

In section 2, we describe the multivariate model, associated estimation of the model parameters, and the proposed parametric confidence interval for a linear combination of small area means. We present our data analysis in section 3. An outline of the technical proof of our main result is deferred to the Appendix.

2 Parametric Bootstrap Confidence Intervals for the Multivariate Fay-Herriot Model

2.1 Multivariate Fay-Herriot model

Let and be a vector of characteristics of interest and a vector of direct survey estimates of for area , respectively, where is the number of small areas. The multivariate Fay-Herriot model (Fay, 1987; Benavent and Morales, 2016) is given by


where is a matrix of known explanatory variables; and are vectors of area specific sampling errors and random effects, respectively; and are all independent with and , being the known sampling variance-covariance matrix of .

We assume that , the variance-covariance matrix of the random effects , depends on unknown parameters with . For the small area application considered by Datta et al. (1991), is an unstructured variance-covariance matrix with and . For the stationary time series cross-sectional model of Rao and Yu (1994), is a structured variance-covariance matrix with as the number of time points, and For the random walk time series cross-sectional model of Datta et al. (2002), is a structured variance-covariance matrix with as the number of time points, and Let be a vector of all the unknown parameters.

For unified representations over areas, we define , , and define and in the same way as . Then, the model can be expressed as

where with and with . With this notation, we can write . In this paper, we are interested in constructing confidence intervals for , where is a -dimensional vector of known constants. For example, if we let , is the first characteristics in the first area, and also can be the difference of characteristics in different areas by setting appropriately.

Under the model (1), the best linear unbiased predictor of with known parameters is given by

which shrinks toward the regression part . Note that each element in depends not only on the corresponding observation but also other observations in the same area when or have non-zero off-diagonal elements. Exploiting the information on the correlation structure, the best linear unbiased predictor would be able to provide more accurate estimates of than simple applications of the univariate FH models to each element. In fact, it will be numerically shown that such advantage is inherited to interval lengths of confidence intervals. The multivariate Fay-Herriot model provides more efficient confidence intervals than the univariate Fay-Herriot model by borrowing information from related components.

2.2 Estimation of model parameters

Because the best linear unbiased predictor depends on unknown parameters, statistical inference on

is carried out via the empirical best linear unbiased estimator given by

where and are some estimators of and . We estimate by the generalized least squares estimator

once in is obtained. There are several different methodology to estimate (e.g. the restricted maximum likelihood estimator (Benavent and Morales, 2016) and moment-based estimators (Ito and Kubokawa, 2021)), but the proposed method to construct the empirical Bayes confidence interval does not depend on a specific variance estimator. For the data analysis below, we adopt the maximum likelihood estimator that maximizes

by the EM algorithm. Note that this method estimates and simultaneously and automatically yields the generalized least squares estimator as the maximum likelihood estimator.

2.3 Confidence intervals via parametric bootstrap

We describe our methodology to construct the empirical Bayes confidence interval for . To motivate our method, we first consider a traditional approach to interval estimation. The key observation for this approach is that the conditional distribution of under the model (1) is given by , where


follows the standard normal distribution, one can find

such that for a fixed . Because the resultant interval for contains unknown parameters and , the traditional approach replaces these parameters by their consistent estimators and to obtain the confidence interval for . Though this interval has a correct coverage asymptotically, it tends to be too short or too long in practice. This undesirable phenomenon is due to the reliance on the rather crude approximation of the standard normal distribution by , which yields the coverage error of . Because and must be estimated, the issue of the asymptotic approximation is not avoidable. Instead, we consider the distribution of from the beginning and consider a method to precisely approximate it. We achieve this goal through the parametric bootstrap.

We construct the bootstrap sample in a prametric way as follows. First we independently generate and . Because and in the model (1), we construct

The resultant bootstrap sample is . To approximate , we compute with . Bootstrap estimators and of and is obtained in the same way as and by replacing the original sample by the bootstrap sample . For example, one can compute the bootstrap maximum likelihood estimator by maximizing

as in our data analysis below. Once is computed, we plug this in to obtain and .

The conditional distribution of

given the data is the parametric bootstrap approximation of the distribution of

. Because the random variable

can be generated as described above, one can find the quantity that satisfies as precisely as possible. Because parametric bootstrap provides a precise approximation, is expected to yield a similar probability for . The proposed parametric bootstrap confidence interval is

The following theorem states that the proposed empirical Bayes confidence interval achieves correct coverage asymptotically with error .

Theorem 1.

We assume the following conditions:

  • The matrix is of full rank satisfying .

  • is a strictly positive definite matrix satisfying where be the Frobenius norm.

  • There exists positive constants and such that the sampling variance-covariance matrix satisfies for .

Let . Suppose satisfies


3 Application

In this section, we use old data used earlier by Datta et al. (1991)

to compare three different confidence interval methods: direct method, parametric bootstrap confidence interval methods – one based on an univariate Fay-Herriot model and the other based on multivariate Fay-Herriot model. The data contain direct survey estimates of median income of 3-, 4- and 5-person families and their associated standard errors for the fifty states and the District of Columbia during years 1979-88. In addition, data contain census median income of 3-, 4-, and 5-person families obtained from the 1970 and 1980 decennial censuses. The U.S. Department of Health and Human Services (HHS) administers a program of energy assistance to low-income families. Eligibility for the program is determined by a formula where the most important variable is an estimate of the current median income for four-person families by states.

Let and denote the true median income of 3-, 4- and 5-person families, respectively, for , where is the number of states and the District of Columbia in the United States. Let and be the corresponding direct survey estimates. Our primary interest is the four-persons family median income, , and we consider estimating the parameter by borrowing strength from not only area specific auxiliary variables but also from the direct survey estimates of median income for the 3- and 5-person families. As for the area specific auxiliary variables, we consider the median income data obtained from the most recent decennial census and an ’adjusted’ census median income obtained by multiplying the most recent census median income by the ratio of per-capita income of the current year to the most recent decennial census year. The per-capita income information is available from administrative records maintained by the Bureau of Economic Analysis (BEA). Then, the covariate matrix is a matrix given by

where and denote the census data and adjusted census in the th area for three-person , four-family and five-person family median incomes.

We first applied the multivariate Fay-Herriot model (MFH) given in (1) separately to the survey data in each year from 1981 to 1988, where we used median income of 1979 obtained from the 1980 decennial census data as auxiliary variables. For comparison, we also applied the univariate Fay-Herriot (UFH) model only to the four-person family income data with the corresponding census data as auxiliary variables. We found that the maximum likelihood estimates of the random effects variance in the UFH model were 0 in 1982, 1983 and 1986, in which confidence intervals of cannot be obtained. On the other hand, we observed that the MFH model produces positive definite estimates for in all the years, and correlations are quite high in some years. This indicates that the random effects variance in can be stably estimated by borrowing strength from other information such as and through the MFH model (1). For illustration, we focus on the results in 1984 and 1987 in which the estimated correlation matrices are given by

respectively. Note that the correlations are quite high in 1987 while relatively small in 1984.

Based on 1000 bootstrap replications, we computed 95 confidence intervals of under both MFH and UFH models. We also computed 95 confidence intervals based on the direct estimator (denoted by DIR), given by , where is the upper quantile of the standard normal distribution, and is the -element of . In Figure 1, we present the differences in lengths of 95 confidence interval based on the MFH model, the UFH model and the DIR method, where the states are arranged in the ascending order of sampling variances. Negative values of the difference indicate that the lengths of confidence intervals from the MFH model are shorter than those from the UFH model or the DIR method. We also reported summary values of area-wise confidence intervals in Table 1. Comparing the inverval from the MFH model and the DIR method, the difference tends to be larger as the sampling variance increases. This is reasonable because we can improve the accuracy of inference on parameters in areas with large sampling variance by borrowing strength through the model. Comparing intervals from the MFH and UFH models, two lengths are comparable in 1984 possibly because the correlations among three median incomes are not so strong. The advantage of borrowing strength from the other incomes can be limited. On the other hand, in 1987, the MFH model produces shorter confidence intervals than the UFH model in almost all the areas due to the high correlations.

We next investigated the performance of the confidence intervals using the census data in 1979 as if they were true values. We applied the MFH and UFH methods to the survey data in 1979 using 1969 census data as covariates. In this case, we applied the UFH method to not only but also and . We first found that in the UFH model both maximum likelihood and restricted maximum likelihood estimates of the random effects variances are zero (for and ) or very small (for ). Based on 1000 bootstrap replications and the maximum likelihood method, we obtained confidence intervals of with , in the MFH and UFH models. We calculated mean and median lengths of area-wise confidence intervals, denoted by Len1 and Len2, respectively. We also computed the empirical coverage rate (CR) by considering 1979 census data as true values. The results are reported in Table 2. Although the UFH model provides shorter confidence intervals than the MFH model, the empirical coverage rate is quite low compared with the nominal level . This suggests that the confidence intervals in the UFH model are too liberal in this case, possibly because of the small estimates of random effects variance. On the other hand, the MFH model provides reasonable confidence intervals. Their coverage rates are quite high and their lengths are much shorter than the those in the direct method.

Figure 1: Differences of 95 confidence interval lengths based on the multivariate Fay-Herriot (MFH), univariate Fay-Herriot (UFH) and the naive method with direct estimator (DIR) in 1984 and 1987 surveys. The states are arranged in the ascending order with respect to the sampling variances.
Year Method min 25 Median Mean 75 max
MFH 4.12 5.55 6.13 6.02 6.38 8.28
1984 UFH 4.21 5.78 6.42 6.32 6.72 8.44
DIR 3.96 6.66 7.75 8.00 8.84 21.30
MFH 4.54 6.08 6.78 6.74 7.23 12.87
1988 UFH 5.86 7.60 8.57 8.39 8.95 11.85
DIR 6.18 8.23 10.71 10.84 12.15 32.39
Table 1: Summary measures of interval lengths of three methods.
three-persons family four-persons family five-persons income
Method CR Len1 Len2 CR Len1 Len2 CR Len1 Len2
MFH 100 3.21 3.13 96.1 2.40 2.29 100 4.52 4.36
UFH 74.5 1.77 1.69 - - - - - -
DIR 86.3 5.57 5.41 86.3 5.88 5.91 86.3 9.25 8.98
Table 2: Performance of confidence intervals of MFH, UFH and DIR. The results of UFH in four- and five-persons family incomes do not exist because of zero estimates of the random effects variance.

4 Concluding Remarks

In this paper, we proposed the parametric bootstrap method for computing a second-order accurate confidence interval of small area parameters from the multivariate Fay-Herriot model. The proposed parametric bootstrap method is easy to implement and is widely applicable to many variance estimators with minimal assumptions as seen in Theorem 1. This advantage forms a sharp contrast to the analytical calibration proposed by Datta et al. (2002) and Ito and Kubokawa (2021) where a different estimation method of model parameters requires cumbersome derivations of the correction terms and tedious checking of assumptions. We demonstrated the superior performance of the proposed methodology over the univariate and direct methods in the family income data. Better coverage and generally shorter length of the proposed interval is due to the effective use of the correlation structure in the same area and direct approximation of the distribution through parametric bootstrap.

There are several future directions to extend the proposed methodology. An immediate extension is to construct the confidence region of a vector of small area parameters studied by Ito and Kubokawa (2021). In the current paper, we focused on the linear combination of small area parameters because the primary interest lies in the single parameter or difference in two parameters in practice. When more than three parameters are of interest, our simple and versatile parametric bootstrap method is expected to be a powerful alternative to the analytical calibration. Though methodology itself is exactly the same as in the current paper, theoretical justification of parametric bootstrap is a challenging problem because it involves multivariate integrals. Another direction is to extend the parametric bootstrap to a more general multivariate linear mixed models. Because theoretical arguments by Chatterjee et al. (2008) for the general univariate case is similar to that in this paper, its multivariate extension can be done in a similar way. Another interesting question is to address the issue of non-positive definiteness of the estimated variance-covariance matrices. Singularities of estimated variance-covariance matrices may occur both in the original estimate and the bootstrap estimate. Because this issue compromises the validity of the parametric bootstrap procedure, it is important to develop reasonable adjustment methods possibly motivated from existing approaches in the univariate situation (e.g. Li and Lahiri, 2010).


The second author’s research was partially supported by Japan Society for Promotion of Science (KAKENHI) 18K12757. The third author’s research was partially supported by the U.S. National Science Foundation Grants SES-1758808.


For notational simplicity, we suppress the dependence on . For example, we write and for and .

Proof of Theorem

Recall that the conditional distribution of given is the multivariate normal distribution with mean and the variance-covariance matrix where . Let . The conditional distribution of given is then the normal distribution with mean and variance . Let and

be the cumulative distribution function and density function for the standard normal random variable. Define

. It follows that

Because for and is uniformly bounded,

for some constant . Thus, the evaluation of reduces to the evaluation of , and . In particular, if we obtain , and , then it follows that by Jensen’s inequality so that


where is a smooth function of . Because we consider the parametric bootstrap, the mathematical argument leading to the last display similarly yields

with some appropriate modifications. In the following, we provide a sketch of the proof of (2) by verifying . Once we obtain this result, proving the statement on the confidence interval is straightforward as in Chatterjee et al. (2008).

To analyze the moment of , first consider the element of . We have

where is a diagonal matrix with 1 in the th element with and 0 otherwise. Let . Thus, we can write

We evaluate moments of . Clearly, . For the second moment, a general term of the matrix is

Because and is fixed, we obtain . For the 8th moment, note that the 8th moment of the sum is the sum of the 8th moments up to constant. Thus we consider the 8th moment of where with . Let and

be a block matrix with blocks of zero matrices and one identity matrix such that

. Because , it follows from the Cauchy-Schwartz inequality that

Because is fixed and ,

for some constant . Since and is the sum of terms, the above expectation is . Hence we obtain .

To evaluate the 8th moment of the rest of terms in , we need to evaluate the moment of and To see this, we have, for example, that

Note that the evaluation of and involves asymptotic expansions of matrix entries. As pointed out by Chatterjee et al. (2008) on page 1240, this computation involves several hundreds pages of elementary calculations. In the end, both 8th moments reduce to the 8th moment of the Frobenius norm of which is . We omit these details and refer to Chatterjee et al. (2008).


  • Benavent and Morales (2016) Benavent, R. and D. Morales (2016). Multivariate fay-herriot models for small area estimation. Computational Statistics & Data Analysis 94, 372–390.
  • Butar and Lahiri (2002) Butar, F. and P. Lahiri (2002). On the measures of uncertainty of empirical bayes small-area estimators. Journal of Statistical Planning and Inference 112, 63–76.
  • Chatterjee et al. (2008) Chatterjee, S., P. Lahiri, and H. Li (2008). Parametric bootstrap approximation to the distribution of eblup and related prediction intervals in linear mixed models. The Annals of Statistics 36, 1221–1245.
  • Cox (1975) Cox, D. (1975). Prediction intervals and empirical bayes confidence intervals. Perspectives in Probability and Statistics. Papers in Honor of M. S. Bartlett (J. Gani, ed.), Applied Probability Trust, Univ. Sheffield, Sheffield. MR0403046, 47–55.
  • Datta et al. (1991) Datta, G., R. Fay, and M. Ghosh (1991). Hierarchical and empirical multivariate bayes analysis in small area estimation. Proceedings of 1991 ARC of Census Bureau, 63–79.
  • Datta et al. (2002) Datta, G., P. Lahiri, and T. Maiti (2002). Empirical bayes estimation of median income of four- person families by state using tie series and cross-sectional data. Journal of Statistical Planning and Inference 102, 83–97.
  • Datta et al. (2002) Datta, G. S., M. Ghosh, D. D. Smith, and P. Lahiri (2002). On an asymptotic theory of conditional and unconditional coverage probabilities of empirical Bayes confidence intervals. Scandinavian Journal of Statistics 29, 139–152.
  • Fay (1987) Fay, R. (1987). Application of multivariate regression to small domain estimation. Small Area Statistics (R. Platek, J.N.K. Rao, C.E. Sarndal, M.P. Singh, eds), 91–102.
  • Fay and Herriot (1979) Fay, R. and R. Herriot (1979). Estimators of income for small area places: an application of james–stein procedures to census. Journal of the American Statistical Association 74, 341–353.
  • Ito and Kubokawa (2021) Ito, T. and T. Kubokawa (2021). Corrected empirical bayes confidence region in a multivariate fay-herriot model. Journal of Statistical Planning and Inference 211, 12–32.
  • Jiang and Lahiri (2006) Jiang, J. and P. Lahiri (2006). Mixed model prediction and small area estimation, editor’s invited discussion paper. Test 15, 1–96.
  • Lahiri (2003a) Lahiri, P. (2003a). On the impact of bootstrap in survey sampling and small-area estimation. Statistical Science 18, 199–210.
  • Lahiri (2003b) Lahiri, P. (2003b). A review of empirical best linear unbiased prediction for the fay-herriot small-area model. The Philippine Statistician 52, 1–15.
  • Li and Lahiri (2010) Li, H. and P. Lahiri (2010). Adjusted maximum method for solving small area estimation problems.

    Journal of Multivariate Analysis

     101, 882–892.
  • Rao and Molina (2015) Rao, J. N. K. and I. Molina (2015). Small Area Estimation, 2nd Edition. Wiley.
  • Rao and Yu (1994) Rao, Y. and M. Yu (1994). Small area estimation by combining time series and cross-sectional data. Canadian Journal of Statistics 22, 511–528.