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 multilevel 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 percapita income of small places (population less than 1000), Fay and Herriot (1979)
proposed an empirical Bayes method to improve on direct surveyweighted estimates by borrowing strength from administrative data and survey estimates from a bigger area. Their method uses a twolevel 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 FayHerriot model useful in investigating various theoretical properties as well as implementing methodology in different applied problems when we do not have access to microdata because of confidentiality and other reasons. For a review on the FayHerriot 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 FayHerriot 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 twolevel multivariate model. In the context of estimating median income of fourperson families for the fifty states and the District of Columbia (small areas), Fay (1987)
suggested a multivariate extension of the FayHerriot model, commonly referred to as the multivariate FayHerriot 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 threeperson and fiveperson families as auxiliary variables in an univariate FayHerriot model. But, unlike the univariate FayHerriot model, the multivariate FayHerriot model incorporates sampling variancecovariance matrix of direct survey estimates of median income of the 3, 4 and 5person families for each small area. Moreover, the multivariate FayHerriot 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 fourperson median income for the small areas drawn from the multivariate FayHerriot model are expected to be more efficient and reasonable when compared to the inferences drawn from an univariate FayHerriot model with survey estimates of the median income for three and five person families as auxiliary variables. This is because the univariate FayHerriot 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 fourperson families for the fifty states and the District of Columbia, Datta et al. (1991) used a bivariate FayHerriot model with a general structure for the variancecovariance matrix of the vector of area specific random effects. However, in many small area applications, structured variancecovariance 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 crosssectional model while Datta et al. (2002) proposed a random walk time series and crosssectional model. Although the time series crosssectional models can be viewed as special cases of the multivariate FayHerriot model, one can achieve greater efficiency in estimating the unknown variance covariance matrix by reducing the number of parameters in the variancecovariance matrix through time series crosssectional models.
Empirical best linear unbiased predictions and associated uncertainty measures for multivariate FayHerriot models with or without structured variancecovariance 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 secondorder efficient confidence intervals for the multivariate FayHerriot model, i.e., confidence intervals with coverage error being the number of small areas, received very little attention. Datta et al. (2002) obtained a secondorder 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 normalitybased 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 secondorder 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 FayHerriot model with an unstructured variancecovariance matrix when a methodofmoment estimator of the variancecovariance matrix of the vector of random effect is used. The derivation of the secondorder 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 FayHerriot model with model variancecovariance structure suggested by the time series crosssectional model of Rao and Yu (1994) or Datta et al. (2002)) or estimation method for the model parameters.Parametric bootstrap method for obtaining secondorder unbiased mean squared error estimation was first proposed by Butar and Lahiri (2002)
. Construction of the secondorder 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 FayHerriot model, see Lahiri (2003a) and Li and Lahiri (2010). In this paper, we develop a parametric bootstrap method for obtaining secondorder efficient confidence intervals for small area parameters from a multivariate FayHerriot model. Compared to the analytical method, our parametric bootstrap approach for constructing secondorder 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 FayHerriot Model
2.1 Multivariate FayHerriot 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 FayHerriot model (Fay, 1987; Benavent and Morales, 2016) is given by
(1) 
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 variancecovariance matrix of .
We assume that , the variancecovariance matrix of the random effects , depends on unknown parameters with . For the small area application considered by Datta et al. (1991), is an unstructured variancecovariance matrix with and . For the stationary time series crosssectional model of Rao and Yu (1994), is a structured variancecovariance matrix with as the number of time points, and For the random walk time series crosssectional model of Datta et al. (2002), is a structured variancecovariance 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 nonzero offdiagonal 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 FayHerriot model provides more efficient confidence intervals than the univariate FayHerriot 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 momentbased 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
Since
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 isThe 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 variancecovariance matrix satisfies for .
Let . Suppose satisfies
Then
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 FayHerriot model and the other based on multivariate FayHerriot model. The data contain direct survey estimates of median income of 3, 4 and 5person families and their associated standard errors for the fifty states and the District of Columbia during years 197988. In addition, data contain census median income of 3, 4, and 5person 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 lowincome families. Eligibility for the program is determined by a formula where the most important variable is an estimate of the current median income for fourperson families by states.
Let and denote the true median income of 3, 4 and 5person 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 fourpersons 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 5person 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 percapita income of the current year to the most recent decennial census year. The percapita 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 threeperson , fourfamily and fiveperson family median incomes.
We first applied the multivariate FayHerriot 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 FayHerriot (UFH) model only to the fourperson 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 areawise 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 areawise 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.
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 
threepersons family  fourpersons family  fivepersons 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 
4 Concluding Remarks
In this paper, we proposed the parametric bootstrap method for computing a secondorder accurate confidence interval of small area parameters from the multivariate FayHerriot 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 nonpositive definiteness of the estimated variancecovariance matrices. Singularities of estimated variancecovariance 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).
Acknowledgement
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 SES1758808.
Appendix
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 variancecovariance 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 thatBecause 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
(2) 
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 CauchySchwartz inequality thatBecause 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).
References
 Benavent and Morales (2016) Benavent, R. and D. Morales (2016). Multivariate fayherriot 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 smallarea 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 crosssectional 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 fayherriot 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 smallarea estimation. Statistical Science 18, 199–210.
 Lahiri (2003b) Lahiri, P. (2003b). A review of empirical best linear unbiased prediction for the fayherriot smallarea 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 crosssectional data. Canadian Journal of Statistics 22, 511–528.
Comments
There are no comments yet.