Simultaneous inference for linear mixed model parameters with an application to small area estimation

TL;DR

This paper develops methods for valid simultaneous inference on parameters in linear mixed models, enabling joint analysis and comparisons of cluster-level estimates, with demonstrated effectiveness through simulations and a real data example.

Contribution

It introduces novel simultaneous prediction intervals and multiple testing procedures for mixed parameters in linear mixed models, using bootstrap methods.

Findings

Methods perform well in simulations

Effective in small area income estimation

Allows joint inference on cluster-level parameters

Abstract

Over the past decades, linear mixed models have attracted considerable attention in various fields of applied statistics. They are popular whenever clustered, hierarchical or longitudinal data are investigated. Nonetheless, statistical tools for valid simultaneous inference for mixed parameters are rare. This is surprising because one often faces inferential problems beyond the pointwise examination of fixed or mixed parameters. For example, there is an interest in a comparative analysis of cluster-level parameters or subject-specific estimates in studies with repeated measurements. We discuss methods for simultaneous inference assuming a linear mixed model. Specifically, we develop simultaneous prediction intervals as well as multiple testing procedures for mixed parameters. They are useful for joint considerations or comparisons of cluster-level parameters. We employ a consistent bootstrap approximation of the distribution of max-type statistic to construct our tools. The numerical performance of the developed methodology is studied in simulation experiments and illustrated in a data example on household incomes in small areas.

Figures

Tables

Table 1. Table 7: Empirical FWER under NERM. The nominal FWER is 0.05 0.05 0.05 .

	BS				BO
	$D=15$	$D=30$	$D=60$	$D=90$	$D=15$	$D=30$	$D=60$	$D=90$
${\sigma }_e^2=0.5$, ${\sigma }_u^2=1$	0.043	0.042	0.049	0.048	0.048	0.044	0.508	0.046
${\sigma }_e^2=1$, ${\sigma }_u^2=1$	0.039	0.037	0.045	0.042	0.042	0.042	0.046	0.044
${\sigma }_e^2=1$, ${\sigma }_u^2=0.5$	0.047	0.040	0.038	0.037	0.025	0.039	0.044	0.040

Table 2. Table 8: ECP (in %), WS and VS under NERM with t-distributed errors and random effects. The nominal coverage probability is 95 % percent 95 95\% .

		ECP				WS (VS)
	$D:n_d$	BS	MC	BE	BO	BS	MC	BE	BO
	15:5	93.7	87.9	77.1	92.4	2.325 (0.079)	2.325 (0.079)	2.145 (0.279)	2.484 (0.111)
$e_{dj}\sim t_6\left(1\right)$	30:10	92.3	90.9	87.4	92.0	1.883 (0.015)	1.883 (0.015)	1.800 (0.025)	1.905 (0.015)
$u_d\sim t_6\left(1\right)$	60:20	91.1	90.7	80.4	90.9	1.454 (0.002)	1.454 (0.002)	1.333 (0.008)	1.460 (0.002)
	90:30	91.4	91.3	72.2	91.5	1.238 (0.001)	1.238 (0.001)	1.096 (0.005)	1.241 (0.001)
	15:5	92.0	85.7	58.3	95.1	2.150 (0.082)	2.150 (0.082)	1.709 (0.330)	2.476 (0.145)
$e_{dj}\sim t_6\left(1\right)$	30:10	94.0	91.0	86.6	93.1	1.797 (0.011)	1.797 (0.011)	1.720 (0.023)	1.842 (0.012)
$u_d\sim N\left(0.5\right)$	60:20	93.5	93.7	81.6	93.8	1.420 (0.002)	1.420 (0.002)	1.293 (0.008)	1.429 (0.002)
	90:30	92.8	92.9	73.3	93.0	1.217 (0.001)	1.217 (0.001)	1.069 (0.005)	1.222 (0.001)
	15:5	95.0	92.0	89.6	93.4	1.752 (0.023)	1.752 (0.023)	1.732 (0.048)	1.795 (0.023)
$e_{dj}\sim N\left(0.5\right)$	30:10	95.2	94.9	89.3	95.1	1.367 (0.003)	1.367 (0.003)	1.287 (0.008)	1.376 (0.003)
$u_d\sim t_6\left(1\right)$	60:20	94.1	94.0	82.5	94.4	1.042 (0.000)	1.042 (0.000)	0.945 (0.004)	1.044 (0.000)
	90:30	95.4	95.3	76.6	95.5	0.883 (0.000)	0.883 (0.000)	0.774 (0.002)	0.885 (0.001)
	15:5	92.5	88.7	86.2	90.6	1.740 (0.043)	1.740 (0.043)	1.726 (0.071)	1.786 (0.048)
$e_{dj}\sim t_6\left(0.5\right)$	30:10	93.5	92.4	86.8	92.8	1.362 (0.007)	1.362 (0.007)	1.295 (0.011)	1.371 (0.007)
$u_d\sim t_6\left(1\right)$	60:20	92.4	92.8	82.4	93.0	1.040 (0.001)	1.040 (0.001)	0.952 (0.004)	1.043 (0.001)
	90:30	92.1	91.9	70.9	92.2	0.882 (0.000)	0.882 (0.000)	0.780 (0.003)	0.884 (0.000)

Table 3. Table 9: ECP (in %), WS and VS under NERM with chi-square distributed errors and random effects. The nominal coverage probability is 95 % percent 95 95\% .

		ECP				WS (VS)
	$D:n_d$	BS	MC	BE	BO	BS	MC	BE	BO
	15:5	94.0	90.6	87.7	91.9	2.452 (0.088)	2.452 (0.088)	2.438 (0.148)	2.516 (0.099)
$e_{dj}\sim {\chi }_5\left(1\right)$	30:10	90.9	90.1	85.2	90.3	1.934 (0.014)	1.934 (0.014)	1.839 (0.024)	1.946 (0.014)
$u_d\sim {\chi }_5\left(1\right)$	60:20	91.2	91.0	80.0	91.1	1.474 (0.002)	1.474 (0.002)	1.351 (0.009)	1.478 (0.002)
	90:30	92.2	91.9	75.1	92.0	1.248 (0.001)	1.248 (0.001)	1.105 (0.005)	1.251 (0.001)
	15:5	91.4	85.8	59.7	94.2	2.148 (0.092)	2.147 (0.092)	1.717 (0.360)	2.456 (0.121)
$e_{dj}\sim {\chi }_5\left(1\right)$	30:10	92.0	89.3	83.7	91.2	1.803 (0.011)	1.803 (0.011)	1.725 (0.023)	1.846 (0.011)
$u_d\sim N\left(0.5\right)$	60:20	92.3	91.7	81.3	92.3	1.419 (0.002)	1.419 (0.002)	1.293 (0.008)	1.428 (0.002)
	90:30	91.9	92.3	73.6	92.5	1.218 (0.001)	1.218 (0.001)	1.069 (0.005)	1.222 (0.001)
	15:5	95.2	93.1	92.0	93.5	1.800 (0.024)	1.800 (0.024)	1.785 (0.035)	1.823 (0.025)
$e_{dj}\sim N\left(0.5\right)$	30:10	95.0	94.1	89.1	94.5	1.386 (0.003)	1.386 (0.003)	1.303 (0.008)	1.392 (0.003)
$u_d\sim {\chi }_5\left(1\right)$	60:20	94.1	94.0	82.5	94.4	1.048 (0.000)	1.048 (0.000)	0.951 (0.004)	1.051 (0.000)
	90:30	95.4	95.3	76.6	95.5	0.886 (0.000)	0.886 (0.000)	0.777 (0.002)	0.888 (0.000)
	15:5	93.5	91.1	90.4	91.5	1.787 (0.049)	1.787 (0.049)	1.787 (0.066)	1.809 (0.052)
$e_{dj}\sim {\chi }_5\left(0.5\right)$	30:10	90.7	90.0	85.4	90.2	1.383 (0.007)	1.383 (0.007)	1.312 (0.012)	1.388 (0.007)
$u_d\sim {\chi }_5\left(1\right)$	60:20	90.7	90.7	79.4	90.9	1.049 (0.001)	1.049 (0.001)	0.960 (0.004)	1.051 (0.001)
	90:30	91.9	92.4	74.3	92.6	0.886 (0.000)	0.886 (0.000)	0.785 (0.003)	0.888 (0.000)

Equations

$${\mathbf{y}}_d={\mathbf{X}}_d\mathbf{\beta }+{\mathbf{Z}}_d{\mathbf{u}}_d+{\mathbf{e}}_d,d=1,\dots ,D,$$

$$y_{dj}={\mathbf{x}}_{dj}^t\mathbf{\beta }+u_d+e_{dj},d=1,\dots ,D,j=1,\dots ,n_d,$$

$$y_d={\mathbf{x}}_d^t\mathbf{\beta }+u_d+e_d,d=1,\dots ,D.$$

$${\mu }_d={\mathbf{k}}_d^t\mathbf{\beta }+{\mathbf{m}}_d^t{\mathbf{u}}_d,d=1,\dots ,D,$$

$${\hat{\mu }}_d:={\hat{\mu }}_d(\hat{\mathbf{\theta }})={\mathbf{k}}_d^t\hat{\mathbf{\beta }}+{\mathbf{m}}_d^t{\hat{\mathbf{u}}}_d,d=1,\dots ,D,$$

$$\mathrm{MSE}({\hat{\mu }}_d)=\mathrm{MSE}({\tilde{\mu }}_d)+\mathbb{E}{\left({\hat{\mu }}_d-{\tilde{\mu }}_d\right)}^2+2\mathbb{E}\left\{({\tilde{\mu }}_d-{\mu }_d)({\hat{\mu }}_d-{\tilde{\mu }}_d)\right\},$$

$${\mathbf{m}}_d^t({\mathbf{G}}_d-{\mathbf{G}}_d{\mathbf{Z}}_d^t{\mathbf{V}}_d^{-1}{\mathbf{Z}}_d{\mathbf{G}}_d){\mathbf{m}}_d+{\mathbf{b}}_d^t{\left(\sum _{d=1}^D{\mathbf{X}}_d^t{\mathbf{V}}_d^{-1}{\mathbf{X}}_d\right)}^{-1}{\mathbf{b}}_d=:g_{1d}(\mathbf{\theta })+g_{2d}(\mathbf{\theta }),$$

$$\alpha =P\left(\left|\frac{{\hat{\mu }}_d-{\mu }_d}{\hat{\sigma }({\hat{\mu }}_d)}\right|⩾c_{S_0}(1-\alpha )\text{for some }d\in [D]\right)=P\left(\underset{d=1,\dots ,D}{\max }\left|\frac{{\hat{\mu }}_d-{\mu }_d}{\hat{\sigma }({\hat{\mu }}_d)}\right|⩾c_{S_0}(1-\alpha )\right),$$

$$S_0:=\underset{d=1,\dots ,D}{\max }\left|S_{0d}\right|,\text{mbox}where{S}_{0d}=\frac{{\hat{\mu }}_d-{\mu }_d}{\hat{\sigma }({\hat{\mu }}_d)},c_{S_0}(1-\alpha ):=\inf \{t\in \mathbb{R}:P(S_0⩽t)⩾1-\alpha \}.$$

$${\mathcal{I}}_{1-\alpha }^S={\text{bigtimes}}_{d=1}^D{\mathcal{I}}_{d,1-\alpha }^S,\text{where}{\mathcal{I}}_{d,1-\alpha }^S=\left\{{\hat{\mu }}_d\pm c_{S_0}(1-\alpha )\hat{\sigma }({\hat{\mu }}_d)\right\},$$

$${\mathcal{I}}_{d,1-\alpha }^{CPI}=\left\{{\hat{\mu }}_d\pm c_d(1-\alpha )\times \hat{\sigma }({\hat{\mu }}_d)\right\}\forall d\in [D],$$

$$H_0:\mathbf{A\mu }=\mathbf{h}vs.H_1:\mathbf{A\mu }\ne \mathbf{h},$$

$$t_H:=\underset{d=1,\dots ,D^{\prime }}{\max }\left|t_{H_d}\right|,S_{H_0}:=\underset{d=1,\dots ,D}{\max }\left|S_{H_{0}d}\right|,t_{H_d}=\frac{{\hat{\mu }}_d^H-h_d}{\hat{\sigma }({\hat{\mu }}_d^H)}\text{and}{S}_{H_{0}d}=\frac{{\hat{\mu }}_d^H-{\mu }_d^H}{\hat{\sigma }({\hat{\mu }}_d^H)},$$

$$S_B^{\ast (b)}:=\underset{d=1,\dots ,D}{\max }\left|S_{Bd}^{\ast (b)}\right|,S_{Bd}^{\ast (b)}=\frac{{\hat{\mu }}_d^{\ast (b)}-{\mu }_d^{\ast (b)}}{{\hat{\sigma }}^{\ast }({\hat{\mu }}_d^{\ast (b)})},b=1,\dots ,B.$$

$${\mathcal{I}}_{1-\alpha }^{BS}={\text{bigtimes}}_{d=1}^D{\mathcal{I}}_{d,1-\alpha }^{BS},\text{where}{\mathcal{I}}_{d,1-\alpha }^{BS}=\left\{{\hat{\mu }}_d\pm c_{BS}(1-\alpha )\hat{\sigma }({\hat{\mu }}_d)\right\}.$$

$$\underset{q\in \mathbb{R}}{\sup }\left|P^{\ast }(S_B^{\ast }⩽q)-P(S_0⩽q)\right|=o_P(1).$$

$$P\left({\mu }_d\in {\mathcal{I}}_{1-\alpha }^{BS}\forall d\in [D]\right)\stackrel{D\to \infty }{\to }1-\alpha .$$

$$S_{B{H}_0}^{\ast (b)}:=\underset{d=1,\dots ,D}{\max }\left|S_{B{H}_{0}d}^{\ast (b)}\right|,S_{B{H}_{0}d}^{\ast (b)}=\frac{{\hat{\mu }}_d^{\ast H(b)}-{\mu }_d^{\ast H(b)}}{{\hat{\sigma }}^{\ast }({\hat{\mu }}_d^{\ast H(b)})},$$

$${\hat{\mathbf{\mu }}}^{\ast H(b)}=({\mathbf{a}}_1^t({\mathbf{k}}_1^t{\hat{\mathbf{\beta }}}^{\ast (b)}+{\mathbf{m}}_1^t{\hat{\mathbf{u}}}_1^{\ast (b)}),\dots ,{\mathbf{a}}_D^t({\mathbf{k}}_D^t{\hat{\mathbf{\beta }}}^{\ast (b)}+{\mathbf{m}}_D^t{\hat{\mathbf{u}}}_D^{\ast (b)}){)}^t:=\mathbf{A}{\hat{\mathbf{\mu }}}^{\ast (b)},$$

$$P\left({\mu }_d\in {\mathcal{I}}_{1-\alpha }^{B^{\prime }S}\forall d\in [D^{\prime }]\right)\stackrel{D^{\prime }\to \infty }{\to }1-\alpha .$$

$$\left[\begin{array}{cc}{\mathbf{X}}^t{\mathbf{R}}^{-1}\mathbf{X}& {\mathbf{X}}^t{\mathbf{R}}^{-1}\mathbf{Z}\\ {\mathbf{Z}}^t{\mathbf{R}}^{-1}\mathbf{X}& {\mathbf{Z}}^t{\mathbf{R}}^{-1}\mathbf{Z}+{\mathbf{G}}^{-1}\end{array}\right]\left[\begin{array}{c}\tilde{\mathbf{\beta }}\\ \tilde{\mathbf{u}}\end{array}\right]=\left[\begin{array}{c}{\mathbf{X}}^t{\mathbf{R}}^{-1}\mathbf{y}\\ {\mathbf{Z}}^t{\mathbf{R}}^{-1}\mathbf{y}\end{array}\right]$$

$$\mathbf{K}\tilde{\mathbf{\varphi }}={\mathbf{C}}^t{\mathbf{R}}^{-1}\mathbf{y},\text{where}\mathbf{K}={\mathbf{C}}^t{\mathbf{R}}^{-1}\mathbf{C}+{\mathbf{G}}^{+},{\mathbf{G}}^{+}=\left[\begin{array}{cc}{0}_{(p+1)\times (p+1)}& 0_{(p+1)\times D}\\ 0_{D\times (p+1)}& {\mathbf{G}}_{D\times D}^{-1}\end{array}\right],$$

$$Z=\frac{{\mathbf{c}}^t\left(\tilde{\mathbf{\varphi }}-\mathbf{\varphi }\right)}{\sqrt{\mathbb{V}\mathrm{ar}\left\{{\mathbf{c}}^t\left(\tilde{\mathbf{\varphi }}-\mathbf{\varphi }\right)\right\}}}=\frac{{\mathbf{c}}^t\left(\tilde{\mathbf{\varphi }}-\mathbf{\varphi }\right)}{\sqrt{{\sigma }_e^2{\mathbf{c}}^t{\left({\mathbf{C}}^t\mathbf{C}+{\sigma }_e^2{\mathbf{G}}^{+}\right)}^{-1}\mathbf{c}}}\sim \mathrm{N}(0,1).$$

$${\mathbf{l}}_M(\mathbf{x},\mathbf{\theta }):={\left({\mathbf{C}}^t\mathbf{C}+{\sigma }_e^2{\mathbf{G}}^{+}\right)}^{-1/2}\mathbf{c},$$

$$\alpha =P\left(\left|{\hat{\mu }}_d-{\mu }_d\right|⩾c_{VT}(1-\alpha ){\hat{\sigma }}_e\Vert {\mathbf{l}}_M(\mathbf{x},\hat{\mathbf{\theta }})\Vert ,\text{ for some }\mathbf{c}\in \mathcal{C}\right).$$

$$\alpha ⩽Q\{c_{VT}(1-\alpha ),{\eta }_0,{\kappa }_0,{\xi }_0\}+2\mathcal{E}P\left\{|t_v|>c_{VT}(1-\alpha ){\xi }_0\right\},$$

$$\left[\begin{array}{c}\hat{\mathbf{\beta }}-\mathbf{\beta }\\ \hat{\mathbf{u}}-\mathbf{u}\end{array}\right]\approx N\left\{0,{\left({\mathbf{C}}^t{\hat{\mathbf{R}}}^{-1}\mathbf{C}+{\hat{\mathbf{G}}}^{+}\right)}^{-1}\right\}.$$

$$S_0=\underset{d=1,\dots ,D}{\max }\left|S_{0d}\right|\approx \underset{d=1,\dots ,D}{\max }\frac{\left|{\bar{\mathbf{c}}}_d^t\left[\begin{array}{c}\hat{\mathbf{\beta }}-\mathbf{\beta }\\ {\hat{\mathbf{u}}}_d-{\mathbf{u}}_d\end{array}\right]\right|}{\hat{\sigma }({\hat{\mu }}_d)}=:\underset{d=1,\dots ,D}{\max }\left|S_{MCd}\right|=S_{MC},$$

$${\mathcal{I}}_{1-\alpha }^{MC}={\text{bigtimes}}_{d=1}^D{\mathcal{I}}_{d,1-\alpha }^{MC}\text{where}{\mathcal{I}}_{d,1-\alpha }^{MC}=\left\{{\hat{\mu }}_d\pm c_{MC}(1-\alpha )\hat{\sigma }({\hat{\mu }}_d)\right\}.$$

$${\mathcal{I}}_{1-\alpha }^{BO}={\text{bigtimes}}_{d=1}^D{\mathcal{I}}_{d,1-\alpha }^{BO},\text{where}{\mathcal{I}}_{d,1-\alpha }^{BO}=\left\{{\hat{\mu }}_d\pm c_{BO}(1-\alpha )\hat{\sigma }({\hat{\mu }}_d)\right\}.$$

Full text

Keywords: max-type statistic, mixed parameter, multiple testing, small area estimation, simultaneous confidence interval

1 Introduction

The family of linear mixed effects models (LMMs) developed by Henderson, (1950) has been extensively applied in the statistical analysis of clustered and longitudinal data (Jiang,, 2007; Verbeke and Molenberghs,, 2000) as well as for the treatment-level analysis in medicine (Francq et al.,, 2019). This modelling framework arises naturally in many fields such as environmental sciences, economics, medicine, etc. Under LMM one supposes that the extra between-cluster variation (or between-subject variation in longitudinal studies) is captured by cluster-specific random effects. Cluster-level parameters might be the most relevant part of the statistical analysis. In particular, they can be modelled by random effects themselves, or more frequently, by mixed effects which are often linear combinations of fixed and random effects. Mixed parameters are particularly appealing in, among others, animal husbandry, ecology and small area estimation (SAE) (see, for example, the monograph of Rao and Molina,, 2015, and the review of Tzavidis et al.,, 2018). The latter critically observed that although the resulting mixed parameter estimates "are a set of numbers of identical definition and simultaneous interest" and that one should thus consider a simultaneous rather than a point estimation problem, the topic of "ensemble properties of small area estimates (…) has been largely overlooked". They mention benchmarking and rank estimation as examples in this direction. This is in line with our observation regarding related literature on Bayesian hierarchical models in which the authors examine constrained Bayes and triple-goal estimation (Gosh,, 1992; Shen and Louis,, 1998). However, most of this literature hardly considers LMM. To the best of our knowledge, beyond such constrained and rank estimation, simultaneous inference for mixed parameters is still missing. This is surprising given the utility of such inference in applied domains, for example, within public health centres carrying out studies on demographic groups, or when statistical offices report to policy makers for resource distribution. The existing pointwise inference or joint estimation of mixed parameters is not less relevant or useful; nevertheless simultaneous inference would provide a framework for formulating statistically valid statements about a set of mixed parameters.

Numerous national and regional governments as well as international organisations conduct studies on socio-economic conditions in order to implement targeted policy interventions. The European Union, World Bank and statistical institutes regularly draft reports on economic development and poverty across countries, regions and provinces. When looking at such regional estimates, practitioners often aim to simultaneously assess and compare them. Nevertheless, existing methods are often not suitable to carry out such assessments or comparisons on a sound statistical basis. As soon as one begins to formulate a comparative statement about the situation in several areas simultaneously, the area-wise (or cluster-wise) analysis is rendered statistically invalid by an additional variability arising from the joint consideration. Consider cluster-wise prediction intervals (CPI) for mixed parameters; the coverage probabilities of $100(1-\alpha)$ intervals refer to the mean across all clusters. This implies that, by construction, about $100\alpha$ percent of the provided intervals (sometimes more) do not contain the true parameter. In other words, each time a statistical institute publishes its estimates for all areas with prediction intervals, the latter fail in at least $100\alpha\%$ to contain the true value. The same holds true for multiple comparisons via testing. The aim of this paper is to develop statistical tools that fill this gap. The investigation of such methods is not only theoretically appealing, but also relevant for practitioners.

We develop simultaneous prediction intervals (SPIs) and multiple testing (MT) procedures to disprove or support simultaneous hypotheses about certain characteristics. More specifically, our main proposal is to use a max-type statistic for a set of mixed parameters. We then employ a bootstrap procedure to consistently approximate the distribution of this statistic. The latter permits us to recover a critical value to construct an operational SPI or conduct MT procedures. Despite the unquestionable utility of such tools in the context of LMM, to the best of our knowledge, we are the first who investigate their theoretical and empirical properties. Furthermore, we compare the performance of our method with alternative simultaneous inference techniques that we adapt from regression and nonparametric curve estimation. Last but not least, we analytically derive simultaneous intervals based on the volume-of-tube formula of Weyl, (1939) to approximate the tail probabilities. Our mathematical investigation demonstrates that his technique is not operational in our context.

Our methods are different from those considered by Sun et al., (1999) or Maringwa et al., (2008) within the framework of longitudinal studies. They propose to apply, respectively, the volume-of-tube formula and Monte Carlo (MC) sampling to construct simultaneous bands for linear combinations of fixed effects only. In contrast, we investigate a more complex problem of examining mixed effects. Our proposal also differs from the derivation of Krivobokova et al., (2010) who employ a mixed model representation for penalized splines to construct uniform bands for one-dimensional regression curves. Contrary to us, the authors can use a simplified version of the volume-of-tube formula. Our results are distinct from those of Ganesh, (2009) who constructs simultaneous Bayesian credible intervals for a linear combination of area-level parameters under the model of Fay and Herriot, (1979). We consider a more general inferential problem in a broder class of LMM within the frequentist framework. Furthermore, our contribution to the area of MT is a practical methodology used under LMM for the first time. Employment of the max-type statistic might be considered as a complement to the study of Kramlinger et al., (2018) who examined chi-square statistics for constructing MT and confidence sets for mixed parameters. In the classical linear regression literature, max-type and chi-square statistics have been considered as complements, and are both well established in the practitioners’ toolbox. We believe that this is equally valid for mixed parameters. The former are more popular for SPI, whereas chi-square statistics are widely recognized for MT. Finally, Reluga et al., (2021) consider simultaneous inference for empirical best predictors under generalised linear mixed models, whereas our paper seeks to address simultaneous inference under LMM. Our study examines the statistical properties of SPI in contrast to the literature that investigates CPI. Starting from the work of Cox, (1975) and Morris, (1983), researchers proposed numerous methods based on analytical derivations (e.g. Basu et al.,, 2003; Kubokawa,, 2010; Yoshimori and Lahiri,, 2014) and resampling (e.g. Hall and Maiti,, 2006; Chatterjee et al.,, 2008). However, CPI and SPI serve different purposes and are not alternatives to each other.

In Section 2 we introduce the modelling framework and the parameter of interest. The construction of SPI and the MT procedure making use of the max-type statistic is outlined in Section 3. In Section 4 we introduce bootstrap-based SPI and MT and prove their consistency. Section 5 contains potential alternatives which we adapted to our setting. We investigate the finite sample performance of our method in Section 6, and apply it to study the household income in Galicia in Section 7. Section 8 contains final remarks and conclusions. Technical details are deferred to Appendix A and the supplementary material. The latter also includes the discussion of extensions of our method.

2 Linear mixed model inference

Consider a classical LMM formulation $\bm{y}=\bm{X\beta}+\bm{Zu}+\bm{e}$, where $\bm{X}$, $\bm{Z}$ are known, full column rank matrices for a fixed and a random part, $\bm{\beta}$ is a vector of fixed effects, $\bm{u}$ is a vector of random effects, and $\bm{e}$ denotes stochastic errors. It is common to assume $\bm{u}$ and $\bm{e}$ to be mutually independent with $\bm{u}\overset{ind}{\sim}N_{q}{(\bm{0},\bm{G})}$ and $\bm{e}\overset{ind}{\sim}N_{n}{(\bm{0},\bm{R})}$. More specifically, consider a LMM with a block diagonal covariance matrix (LMMb):

[TABLE]

where $n_{d}$ is the number of units in the $d^{th}$ cluster (or area), $\bm{y}_{d}\in\mathbb{R}^{n_{d}}$, $\bm{X}_{d}\in\mathbb{R}^{n_{d}\times(p+1)}$ and $\bm{Z}_{d}\in\mathbb{R}^{n_{d}\times q_{d}}$. Here, $D$ is the number of clusters, $\bm{\beta}\in\mathbb{R}^{p+1}$ an unknown vector of regression coefficients, $\bm{u}_{d}\overset{ind}{\sim}N_{q_{d}}{(\bm{0},\bm{G}_{d})}$ and $\bm{e}_{d}\overset{ind}{\sim}\prescript{}{n_{d}}{(\bm{0},\bm{R}_{d})}$, $n=\sum_{d=1}^{D}n_{d}$. We assume that $\bm{G}_{d}=\bm{G}_{d}(\bm{\theta})\in\mathbb{R}^{q_{d}\times q_{d}}$ and $\bm{R}_{d}=\bm{R}_{d}(\bm{\theta})\in\mathbb{R}^{n_{d}\times n_{d}}$ depend on variance parameters $\bm{\theta}=(\theta_{1},...,\theta_{h})^{t}$. LMM can be easily retrieved applying the notation of Prasad and Rao, (1990). Under this setup, suppose that the variance-covariance $\bm{V}$ is nonsingular $\forall\theta_{i}$, $i=1,\dots,h$ with $\mathbb{E}(\bm{y})=\bm{X\beta}$ and $\mathbb{V}\mathrm{ar}(\bm{y})=\bm{R}+\bm{ZG}\bm{Z}^{t}=\bm{V}(\bm{\theta})\eqqcolon\bm{V}$. Two important examples of LMM that are extensively used, especially in SAE, are the nested error regression model (NERM) of Battese et al., (1988), and the Fay-Herriot model (FHM) of Fay and Herriot, (1979). The former is defined as

[TABLE]

where $y_{dj}$ is the quantity of interest for the $j^{th}$ unit in the $d^{th}$ cluster, $\bm{x}_{dj}=(1,x_{dj1},\dots,x_{djp})^{t}$, $u_{d}\overset{iid}{\sim}N(0,\sigma^{2}_{u})$ and $e_{dj}\overset{iid}{\sim}N(0,\sigma^{2}_{e})$ for $d=1,\dots,D$, $j=1,\dots,n_{d}$. Here $\bm{y}_{d}=(y_{d1},\dots,y_{dn_{d}})$, $\bm{X}_{d}=\mathrm{col}_{1\leqslant j\leqslant n_{d}}\bm{x}^{t}_{dj}$, $q_{d}=1$, $Z_{d}=\bm{1}_{n_{d}}$ with $\bm{1}_{n_{d}}$ a $n_{d}$ vector of ones, $\bm{e}_{d}=(e_{d1},\dots,e_{dn_{d}})^{t}$, $\bm{\theta}=(\sigma_{e}^{2},\sigma_{u}^{2})^{t}$, $\bm{R}_{d}(\bm{\theta})=\sigma_{e}^{2}\bm{I}_{n_{d}}$ with $\bm{I}_{n_{d}}$ the $n_{d}\times n_{d}$ identity matrix and $\bm{G}_{d}(\bm{\theta})=\sigma_{u}^{2}$. In contrast, the FHM is often referred to as an area-level model and consists of two levels. The model at level 1, called the sampling model, assumes that the direct estimators $y_{d}$ of a cluster mean $\mu_{d}^{F}$ are design unbiased, and satisfy $y_{d}=\mu_{d}^{F}+e_{d}$, $e_{d}\overset{iid}{\sim}N(0,\sigma^{2}_{e_{d}})$, $d=1,\dots,D$. Under FHM, the sampling variance $\sigma^{2}_{e_{d}}=\mathbb{V}ar(y_{d}|\mu_{d}^{F})$ is supposed to be known for each cluster $d$. On the other hand, the linking model at level 2 is $\mu_{d}^{F}=\bm{x}^{t}_{d}\bm{\beta}+u_{d}$, $u_{d}\overset{iid}{\sim}N(0,\sigma^{2}_{u})$, $d=1,\dots,D$, where $\bm{x}_{d}=(1,x_{d1},\dots,x_{dp})^{t}$ is a $(p+1)$-vector of cluster-level auxiliary variables. Observe that the FHM can be rewritten as a LMM with $n_{d}=q_{d}=1$, $Z_{d}=1$, $\bm{\theta}=\sigma^{2}_{u}$, $\bm{R}_{d}(\sigma^{2}_{u})=\sigma^{2}_{e_{d}}$, that is

[TABLE]

Due to the data availability, the FHM is more frequently used in practice. While cluster-level information can be easily obtained (for example, using open access internet repositories), this is clearly not the case for units.

Assuming LMM, one is often interested in a simultaneous or comparative inference for general mixed parameters

[TABLE]

with $\bm{k}_{d}\in\mathbb{R}^{p+1}$ and $\bm{m}_{d}\in\mathbb{R}^{q_{d}}$ known. In our article $\mu_{d}$ is a cluster conditional mean, but other parameters can be explored as well. Henderson, (1975) developed the best linear unbiased predictor (BLUP) of a linear combination of random and fixed effects when $\bm{V}$ is completely known. Applying his idea one obtains BLUP estimator for (4) that is $\tilde{\mu}_{d}\coloneqq\tilde{\mu}_{d}(\bm{\theta})=\bm{k}^{t}_{d}\tilde{\bm{\beta}}+\bm{m}^{t}_{d}\tilde{\bm{u}}_{d}$, where $\bm{\theta}=(\theta_{1},\dots,\theta_{h})^{t}$, $\tilde{\bm{\beta}}=\tilde{\bm{\beta}}(\bm{\theta})=(\bm{X}^{t}\bm{V}^{-1}\bm{X})^{-1}\bm{X}^{t}\bm{V}^{-1}\bm{y}$, $\tilde{\bm{u}}_{d}=\tilde{\bm{u}}_{d}(\bm{\theta})=\bm{G}_{d}\bm{Z}^{t}_{d}\bm{V}_{d}^{-1}(\bm{y}_{d}-\bm{X}_{d}\tilde{\bm{\beta}})$. In practice $\bm{\theta}$ is usually unknown, hence one uses $\hat{\bm{\theta}}\coloneqq\hat{\bm{\theta}}(\bm{y})$ which yields the EBLUP

[TABLE]

with $\hat{\bm{\beta}}=\hat{\bm{\beta}}(\hat{\bm{\theta}})$, $\hat{\bm{u}}=\hat{\bm{u}}(\hat{\bm{\theta}})$ and $\hat{\bm{\theta}}=(\hat{\theta}_{1},\dots,\hat{\theta}_{h})^{t}$. Having assumed certain conditions on the distributions of random effects and errors, as well as the variance components $\bm{\theta}$ (see Appendix A.1), Kackar and Harville, (1981) proved that the two-stage procedure provides an unbiased estimator for $\mu_{d}$.

To construct a studentized max-type statistic, it is important to assess the variability of prediction. The most common measure of uncertainty is the mean squared error $\mathrm{MSE}(\hat{\mu}_{d})=\mathbb{E}(\hat{\mu}_{d}-\mu_{d})$. Here, $\mathbb{E}$ denotes the expectation with respect to model (1). We can decompose the MSE into

[TABLE]

where $\mathrm{MSE}(\tilde{\mu}_{d})$ accounts for the variability when the variance components $\bm{\theta}$ are known. Assuming LMMb and $\bm{b}_{d}^{t}=\bm{k}^{t}_{d}-\bm{o}_{d}^{t}\bm{X}_{d}$ with $\bm{o}^{t}_{d}=\bm{m}^{t}_{d}\bm{G}\bm{Z}^{t}_{d}\bm{V}^{-1}_{d}$, the $\mathrm{MSE}(\tilde{\mu}_{d})$ reduces to

[TABLE]

where $g_{1d}$ accounts for the variability of $\tilde{\mu}_{d}$ once $\bm{\beta}$ is known, and $g_{2d}$ for the estimation of $\tilde{\bm{\beta}}$. The second term in (6) is intractable, but there exists a vast literature which deals with the estimation of it (see Rao and Molina,, 2015, for a review). The third term disappears under normality of errors and random effects; it is rarely considered. Following Chatterjee et al., (2008), we suggest the construction of a SPI using $\bm{g}_{1}(\hat{\bm{\theta}})=(g_{11}(\hat{\bm{\theta}}),\dots,g_{1D}(\hat{\bm{\theta}}))^{t}$, where $g_{1d}(\hat{\bm{\theta}})$ is defined in (7) with $\bm{\theta}$ replaced by a consistent estimator. In fact, simulations studies in Reluga, (2020) indicate that alternative measurements of variability do not improve the performance of SPI based on the max-type statistic.

3 SPI and MT for mixed parameters using max-type statistics

We concentrate on the construction of SPI and MT procedures for the mixed parameter in (4) considering a confidence region $\mathcal{I}_{1-\alpha}=\bigtimes_{d=1}^{D}\mathcal{I}_{d,1-\alpha}$ such that $P(\mu_{d}\in\mathcal{I}_{1-\alpha}\;\forall d\in[D])=1-\alpha$, $[D]=\{1,\dots,D\}$. This is equivalent to finding a critical value $c_{S_{0}}(1-\alpha)$ which satisfies

[TABLE]

where we denote by $\hat{\sigma}(\hat{\mu}_{d})$ the estimated variability of $\hat{\mu}_{d}$ (for example, the estimated square root of $\mathrm{MSE}(\hat{\mu}_{d})$). The critical value $c_{S_{0}}(1-\alpha)$ is in fact the $(1-\alpha)^{th}$-quantile of the studentized statistic

[TABLE]

It follows that with probability $1-\alpha$, a region defined as

[TABLE]

covers all mixed parameters. Since the probability density function (pdf) of $S_{0}$ is right skewed, we suggest to consider its upper quantile and construct symmetric $\mathcal{I}^{S}_{d,1-\alpha}$, $d\in[D]$. This approach can be regarded as a variation of the studentized maximum modulus method of Tukey, (1953). At this point, we formally define CPI to circumvent all possible doubts concerning its relation to SPI. Let $c_{d}(1-\alpha)\coloneqq\inf\{t\in\mathbb{R}:P(S_{0d}\leqslant t)\geqslant 1-\alpha\}$. CPI is defined as

[TABLE]

which covers $\mu_{d}$ with probability $1-\alpha$. Due to the central limit theorem, the most common choice is $c_{d}(1-\alpha)=\Phi^{-1}(1-\alpha)$, that is, a quantile from a normal distribution. Furthermore, due to the correspondence between interval estimation and hypothesis testing, our methodology is applicable for the latter. Consider a following pair of hypotheses:

[TABLE]

where $\bm{A}\in\mathbb{R}^{D^{\prime}\times D}$ with $D^{\prime}\leqslant D$ and $\bm{h}\in\mathbb{R}^{D^{\prime}}$ is a vector of constants. A test based on a max-type statistic $t_{H}$ rejects $H_{0}$ at the $\alpha$-level if $t_{H}\geqslant c_{H_{0}}(1-\alpha)$ with $c_{H_{0}}(1-\alpha)\coloneqq\inf\{t\in\mathbb{R}:P(S_{H_{0}}\leqslant t)\geqslant 1-\alpha\}$,

[TABLE]

where $\bm{\mu}^{H}=(\mu^{H}_{1},\dots,\mu^{H}_{D^{\prime}})^{t}=\bm{A}\bm{\mu}\in\mathbb{R}^{D^{\prime}}$ and $\hat{\bm{\mu}}^{H}$ its estimated counterpart. In other words, $\bm{h}\notin\mathcal{I}_{1-\alpha}^{H_{0}}$ with $\mathcal{I}_{1-\alpha}^{H_{0}}=\bigtimes_{d=1}^{D}\mathcal{I}_{d,1-\alpha}^{H_{0}}$, where $\mathcal{I}_{d,1-\alpha}^{H_{0}}=\left\{\hat{\mu}^{H}_{d}\pm c_{H_{0}}(1-\alpha)\hat{\sigma}(\hat{\mu}^{H}_{d})\right\}$. In practice, a standard problem is to test for statistical differences between various clusters with respect to some characteristic. Our test is based on a single step procedure and exhibits a weak control of a family-wise error (FWER). If one aims at testing multiple hypotheses with a strong control of FWER, the step-down technique of Romano and Wolf, (2005) could be implemented. As this is beyond the scope of this paper, details and related simulation results are deferred to the supplementary material.

4 Bootstrap-based SPI and MT procedure

It is challenging to estimate the distribution of $S_{0}$ in (8) and to recover critical values because, among others, mixed effects $\mu_{d}$, $d=1,...,D$, are unknown.Nevertheless, an almost straightforward way to approximate critical value $c_{S_{0}}(1-\alpha)$ is to use a parametric bootstrap procedure which circumvents a direct application of the normal asymptotic distribution (González-Manteiga et al.,, 2008). It can also provide faster convergence (Hall and Maiti,, 2006; Chatterjee et al.,, 2008). Let $B$ be the number of bootstrap samples $(\bm{y}^{*(b)},\bm{X},\bm{Z})$. The bootstrap analogue of expression in (8) is

[TABLE]

The critical value can be consistently approximated by the $(1-\alpha)^{th}$-quantile of (11), i.e., $c_{BS}(1-\alpha)\coloneqq\inf\{t^{*}\in\mathbb{R}:P(S^{*}_{B}\leqslant t^{*})\geqslant 1-\alpha\}$.Consequently, the bootstrap SPI is defined as

[TABLE]

Our choice of $\hat{\sigma}(\hat{\mu}_{d})=\sqrt{g_{1d}(\hat{\bm{\theta}})}$ is motivated by the asymptotic analysis of Chatterjee et al., (2008). Validity of the above bootstrap method is shown by adapting Theorem 3.1 of these authors (henceforth Theorem CLL, provided in the supplementary material) and combining it with some results from the extreme value theory.

Proposition 1**.**

Suppose that the assumptions of Theorem CLL and regularity conditions R.1-R.7 from Section A.1 hold. Then

[TABLE]

An important implication of Proposition 1 is the coverage probability of $\mathcal{I}_{1-\alpha}^{BS}$.

Corollary 1**.**

Under Proposition 1 it holds that

[TABLE]

According to the theoretical developments for max-type statistics, the Kolmogorov distance defined in Proposition 1 converges to [math] at best at polynomial rate $(\log(\cdot))^{c_{2}}/n^{c_{3}}$, where $(\cdot)$ is the number of parameters for which we wish to obtain the maximum (in our case $D$), and $c_{2},c_{3}$ some constants, see e.g. Chernozhukov et al., (2013). We are not aware of results for max-type statistics that one could employ to obtain second order correctness without such $\log(\cdot)$ term. Observe that in our setting $n\rightarrow\infty$ is equivalent to $D\rightarrow\infty$, because we assumed that $n_{d}$ is bounded (see Appendix A.1). Certainly, our result would still hold if both, $n_{d}$ and $D$ grow. In case of a fixed $D$, we would replace Proposition 1 using explicitly a variation of studentized maximum modulus distribution (cf. Stoline and Ury,, 1979).

Due to the relation between interval estimation and tests, the methodology developed for SPI can be used to find a critical value for MT procedure in (9). In particular, consider slightly modified bootstrap statistics

[TABLE]

with $\bm{\mu}^{*H(b)}=(\mu^{*H(b)}_{1},\dots,\mu^{*H(b)}_{D^{\prime}})^{t}\coloneqq\bm{A}\bm{\mu}^{*(b)}\in\mathbb{R}^{D^{\prime}}$, and its estimated versions

[TABLE]

where $\bm{a}_{d}\in\mathbb{R}^{D}$ are the rows of $\bm{A}$. These are applied to find a bootstrap approximation for the critical value $c_{H_{0}}(1-\alpha)$ of our test, namely $c_{BH_{0}}(1-\alpha)\coloneqq\inf\{t\in\mathbb{R}:P(S^{*}_{BH_{0}}\leqslant t)\geqslant 1-\alpha\}$. It is worth mentioning that we do not need to generate bootstrap samples under $H_{0}$ to obtain the critical values of our test.

In Section 2 we defined NERM and FHM as popular examples of LMM. We describe a parametric bootstrap procedure that yields promising results when constructing SPI under these models. Under NERM and FHM, we use simplified versions of $g_{1d}$ in (7) derived by Prasad and Rao, (1990) as the estimators of $\hat{\sigma}^{2}(\hat{\mu}_{d})$: $g^{N}_{1d}(\hat{\bm{\theta}})=\hat{\sigma}^{2}_{u}/(\hat{\sigma}^{2}_{u}+\hat{\sigma}^{2}_{e}/n_{d})(\hat{\sigma}^{2}_{e}/n_{d})$ for NERM and $g^{F}_{1d}(\hat{\bm{\theta}})=\hat{\sigma}^{2}_{u}\sigma^{2}_{e_{d}}/(\hat{\sigma}^{2}_{u}+\sigma^{2}_{e_{d}})$ for FHM. Under NERM the bootstrap algorithm is

From the original sample, obtain consistent estimators $\hat{\bm{\beta}}$ and $\hat{\bm{\theta}}=(\hat{\sigma}^{2}_{e},\hat{\sigma}^{2}_{u})$. 2. 2.

Generate $D$ independent copies of $W_{1}\sim N(0,1)$. Construct $\bm{u}^{*}=(u^{*}_{1},u^{*}_{2},\dots,u^{*}_{D})$ with $u^{*}_{d}=\hat{\sigma}_{u}W_{1}$, $d=[D]$. 3. 3.

Generate $n$ independent copies of $W_{2}\sim N(0,1)$. Construct $\bm{e}^{*}=(e^{*}_{1},e^{*}_{2},\dots,e^{*}_{n})$ with $e^{*}_{j}=\hat{\sigma}_{e}W_{2}$, $j=[n]$. 4. 4.

Create a bootstrap sample $\bm{y}^{*}=\bm{X}\hat{\bm{\beta}}+\bm{u}^{*}+\bm{e}^{*}$. 5. 5.

Fit the model to the bootstrap sample and obtain bootstrap estimates $\hat{\bm{\beta}}^{*}$, $\hat{\bm{\theta}}^{*}=(\hat{\sigma}^{2*}_{e},\hat{\sigma}^{2*}_{u})$, $\mu_{dj}^{*}$ and $\hat{\mu}_{dj}^{*}$. 6. 6.

Repeat Steps 2-5 $B$ times. Calculate $S^{*(b)}_{B}$, $b=1,\dots,B$, using $g^{N*(b)}_{1d}(\hat{\bm{\theta}}^{*(b)})$ to obtain $c_{BS}(1-\alpha)$ and $\mathcal{I}^{BS}_{1-\alpha}$.

Here, $g^{N*(b)}_{1d}(\hat{\bm{\theta}}^{*(b)})=\hat{\sigma}^{2*(b)}_{u}/(\hat{\sigma}^{2*(b)}_{u}+\hat{\sigma}^{2*(b)}_{e}/n_{d})(\hat{\sigma}^{2*(b)}_{e}/n_{d})$ is the bootstrap equivalent of $g^{N}_{1d}(\hat{\bm{\theta}})$. To implement the analogous bootstrap under FHM, we need to slightly modify step 1 defining $\hat{\bm{\theta}}=\hat{\sigma}^{2}_{u}$ and $g^{F{*(b)}}_{1d}(\hat{\bm{\theta}}^{*(b)})=\hat{\sigma}^{2*(b)}_{u}\sigma^{2}_{e_{d}}/(\hat{\sigma}^{2*(b)}_{u}+\sigma^{2}_{e_{d}})$. Additionally, we need to replace step 3 by

1. 3.’

Generate $D$ independent copies of a variable $W_{2}\sim N(0,1)$. Construct vector $\bm{e}^{*}=(e^{*}_{1},e^{*}_{2},\dots,e^{*}_{D})$ with elements $e^{*}_{d}={\sigma}_{e_{d}}W_{2}$, $d=[D]$.

The parametric bootstrap algorithm can be modified to accommodate more complex models, for example with spatial or temporal correlation, by adapting accordingly the process of generating errors and random effects. An indisputable advantage of the bootstrap approach is its generality. As soon as we can mimic a data generating process for the assumed model, it can be implemented and applied to construct SPI and carry out MT for any kind of estimator. In addition, bootstrap SPI are relatively robust to model misspecifications (see results in Table 3), in particular when the number of units in each cluster grows. This is in alignment with related remarks of Jiang, (1998). On the other hand, it is usually more computer intensive than an analytical derivation.

We conclude this section with a practical extension of our results. In the testing problem (9) we have already allowed for a scenario where only $D^{\prime}<D$ hypotheses were considered, even though all data were used to estimate fixed parameters and predict random effects. Similarly, one might be interested in the construction of SPI with a joint coverage probability for a subset of $D^{\prime}<D$ clusters. Without loss of generality, we assume that our goal is to construct SPI for the first $D^{\prime}$ cluster-level mixed parameters. Then, in the definition of $S_{0}$ in (8) and $S^{*}_{B}$ (11) one replaces $\max_{d=1,\dots,D}$ by $\max_{d=1,\dots,D^{\prime}}$ and proceeds along the same lines as for $D$ areas. If $D^{\prime}=O(D)$, we can evoke the same results from the extreme value theory as in case of Proposition 1 to prove a result similar to Corollary 1, that is:

Corollary 2**.**

Let $D^{\prime}<D$, $d\in[D^{\prime}]$, $D^{\prime}=O(D)$. Consider $S^{*(b)}_{B^{\prime}}\coloneqq\max_{d=1,\dots,D^{\prime}}\left\lvert S^{*(b)}_{Bd}\right\rvert$, $c_{B^{\prime}S}(1-\alpha)\coloneqq\inf\{t^{*}\in\mathbb{R}:P(S^{*}_{B^{\prime}}\leqslant t^{*})\geqslant 1-\alpha\}$ where $S^{*(b)}_{Bd}$ as defined in (11), and $\mathcal{I}^{B^{\prime}S}_{1-\alpha}=\bigtimes_{d=1}^{D^{\prime}}\mathcal{I}^{B^{\prime}S}_{d,1-\alpha}$ with $\mathcal{I}^{B^{\prime}S}_{d,1-\alpha}=\left\{\hat{\mu}_{d}\pm c_{B^{\prime}S}(1-\alpha)\hat{\sigma}(\hat{\mu}_{d})\right\}$. Then, under Proposition 1 it holds that

[TABLE]

5 Alternative methods for SPI and MT

Although, to the best of our knowledge, we are the first who introduce SPI and MT procedures for mixed parameters, various approaches have been put forward to tackle the problem of simultaneous confidence bands for linear regression surfaces. For example, Bonferroni t-statistics are a straightforward tool to compare a set of fixed parameters. Furthermore, other authors, such as Working and Hotelling, (1929), Scheffé, (1953), Sun and Loader, (1994) or Beran, (1988), to mention a few, developed equally important methodologies for the simultaneous inference of fixed parameters. In the rest of this section, we adapt some of the methods from the linear regression and nonparametric curve estimation to our setting. One could thus treat them as alternative approaches to our proposal.

5.1 The volume-of-tube procedure

Consider LMMb defined in (1). One way to obtain BLUP estimates for $\bm{\beta}$ and $\bm{u}$ is to solve the mixed model equations of Henderson, (1950):

[TABLE]

which can be re-expressed in the following simplified form:

[TABLE]

with $\tilde{\bm{\phi}}=\left(\tilde{\bm{\beta}}^{t},\tilde{\bm{u}}^{t}\right)^{t}$, $\bm{C}=\left[\bm{X}\>\bm{Z}\right]$. From (13) and (14) we obtain a straightforward formula for the estimates $\tilde{\bm{\phi}}=\bm{K}^{-1}\bm{C}^{t}\bm{R}^{-1}\bm{y}$ and $\hat{\bm{\phi}}=\tilde{\bm{\phi}}(\hat{\bm{\theta}})$. For some $\bm{x}=\left(1,x_{1},\dots,x_{p}\right)^{t}$ with $x_{1},\dots,x_{p}\in\mathcal{X}\subset\mathbb{R}^{p}$, $\bm{z}=\left(z_{1},\dots,z_{q}\right)^{t}\in\mathcal{Z}\subset\mathbb{R}^{q}$ and $\bm{c}=(\bm{x}^{t},\bm{z}^{t})^{t}\in\mathcal{X}\times\mathcal{Z}\eqqcolon\mathcal{C}$ one has $\bm{x}^{t}\tilde{\bm{\beta}}+\bm{z}^{t}\tilde{\bm{u}}=\bm{c}^{t}\tilde{\bm{\phi}}\equiv\bm{l}(\bm{x},\bm{\theta})^{t}\bm{y}=\sum_{i=1}^{n}l_{i}(\bm{x},\bm{\theta})y_{i}$ where $\bm{l}(\bm{x},\bm{\theta})^{t}=(l_{1}(\bm{x},\bm{\theta}),l_{2}(\bm{x},\bm{\theta}),\dots,l_{n}(\bm{x},\bm{\theta}))=\bm{c}^{t}(\bm{C}^{t}\bm{R}^{-1}\bm{C}+\bm{G}^{+})^{-1}\bm{C}^{t}\bm{R}^{-1}$ i.e., $\bm{l}(\bm{x},\bm{\theta})$ is an $n$-vector. In addition, BLUP fitted values are $\tilde{\bm{y}}=\bm{L}\bm{y}$ where $\bm{L}=\bm{C}(\bm{C}^{t}\bm{R}^{-1}\bm{C}+\bm{G}^{+})^{-1}\bm{C}^{t}\bm{R}^{-1}$ which can be considered as a ridge regression formulation of the BLUP. Having reformulated the LMMb, we extend the approach of Sun et al., (1999) and Krivobokova et al., (2010). To simplify the notation, let $\bm{R}=\sigma^{2}_{e}\bm{I}_{n}$ which leads to $\bm{l}(\bm{x},\bm{\theta})^{t}=\bm{c}^{t}\left(\bm{C}^{t}\bm{C}+\sigma_{e}^{2}\bm{G}^{+}\right)^{-1}\bm{C}^{t}$. Assuming normality for errors and random effects one obtains

[TABLE]

We conclude that $Z$ is a nonsingular Gaussian random variable with mean 0 and variance 1. Consequently, the following expressions can be retrieved from equation (15):

[TABLE]

The problem of finding a $(1-\alpha)$ SPI for $\hat{\mu}_{d}$  $\forall d\in[D]$ boils down to the choice of a critical value $c_{VT}(1-\alpha)$ such that

[TABLE]

Let $\bm{l}_{M}=\bm{l}_{M}(\bm{x},\bm{\theta})$, $\hat{\bm{l}}_{M}=\bm{l}_{M}(\bm{x},\bm{\hat{\theta}})$, $\bm{e}_{M}=\bm{e}_{M}(\bm{x},\bm{\theta})$, $\hat{\bm{e}}_{M}=\bm{e}_{M}(\bm{x},\hat{\bm{\theta}})$ and $\lambda_{M}=\hat{\bm{e}}_{M}-\bm{e}_{M}$. We state Proposition 2 for $p=2$, whereas other cases are considered in the supplementary material.

Proposition 2**.**

Suppose that $\sigma^{2}_{e}$ is estimated by some consistent estimator. Define $\mathcal{Q}=\bm{l}_{M}/\left\lVert\bm{l}_{M}\right\rVert$, $\xi=\inf\limits_{\bm{c}\in\mathcal{C}}||\hat{\bm{l}}_{M}||/||\bm{l}_{M}||$ and $\eta=\sup\limits_{\bm{c}\in\mathcal{C}}|(\hat{\bm{l}}_{M}-\bm{l}_{M})^{t}\bm{e}_{M}+\hat{\bm{l}}^{t}_{M}(\hat{\bm{e}}_{M}-\bm{e}_{M})|/\sigma_{e}||\bm{l}_{M}||$. We assume that $\exists$ $\xi_{0}>0$, $\eta_{0}>0$ such that $P(\xi\leqslant\xi_{0})=o(\alpha)$ and $P(\eta\leqslant\eta_{0})=o(\alpha)$ as $n\rightarrow\infty$ and $\alpha\rightarrow 0$. Thus one can approximate (16) as follows:

[TABLE]

where $t_{\nu}$ is a t-distributed random variable with $\nu$ degrees of freedom, $\kappa_{0}=\int_{\bm{c}\in\mathcal{C}}\left\lVert\mathcal{Q}^{\prime}(\bm{x})\right\rVert d{\bm{x}}$ the volume of the manifold $\mathcal{M}=\left\{\mathcal{Q}(\bm{c}),\bm{c}\in\mathcal{C}\right\}$, and $\zeta_{0}$ the boundary area of $\mathcal{M}$. Finally, $\mathcal{E}$ is the Euler-Poincaré characteristic of $\mathcal{M}$.

Due to its limited practical relevance, the proof of Proposition 2 and a spelled-out expression for function $Q$ in (17) are deferred to our supplementary material. By assuming $\bm{R}=\sigma^{2}_{e}\bm{I}_{n}$ our proposal does not suffer from the loss of generality if we can write $\bm{R}=\sigma^{2}_{e}\bm{R}_{e}$, where $\bm{R}_{e}$ is some positive definite matrix. Then it follows that $\bm{l}(\bm{x},\bm{\theta})=\bm{c}^{t}(\bm{C}^{t}\bm{R}^{-1}\bm{C}+\bm{G}^{+})^{-1}\bm{C}^{t}\bm{R}^{-1}=\bm{c}^{t}(\bm{C}^{t}\bm{R}_{e}^{-1}\bm{C}+\sigma^{2}_{e}\bm{G}^{+})^{-1}\bm{C}^{t}\bm{R}_{e}^{-1}$, and we can still use the distribution of the ratio of $\hat{\sigma}_{e}/\sigma_{e}$ which is essential in the derivation of the confidence bands. In case of FHM, $\bm{R}=\mathrm{diag}(\sigma^{2}_{e_{d}})$ is assumed to be known, but in practice we could estimate it from other sources. Thus the derivation could be applied as well.

Having retrieved the critical value, one could construct the volume-of-tube SPI using $\mathcal{I}^{VT}_{1-\alpha}=\bigtimes_{d=1}^{D}\mathcal{I}^{VT}_{d,1-\alpha}$, where $\mathcal{I}^{VT}_{d,1-\alpha}=\left\{\hat{\mu}_{d}\pm c_{VT}(1-\alpha)\hat{\sigma}_{e}||\hat{\bm{l}}_{M}||\right\}$. The approximation in Proposition 2 is conservative, i.e., the coverage probability is higher than the nominal $1-\alpha$; it approaches $1-\alpha$ as $\alpha\rightarrow 0$ and $n\rightarrow\infty$. Similarly as in in Section 4, the latter assumption is equivalent to $D\rightarrow\infty$ in our setting. The results are valid under different asymptotic regimes too, if we are able to consistently estimate the variance parameters, in particular $\sigma^{2}_{e}$. We immediately see that function $Q$ in (17) contains several constants. Numerical approximation of $\kappa_{0}$ which describes the geometry of the manifold $\mathcal{M}$ may not pose a major problem, but it is not clear how to estimate $\xi_{0}$ and $\eta_{0}$ under LMMb. Some ideas were derived for simpler one-dimensional models. Sun et al., (1999) proposes a derivative and a perturbation method to estimate constant $\xi_{0}$, while Sun and Loader, (1994) suggest estimating $\eta_{0}$ nonparametrically. It is unclear, though, how to extend their implementations to the LMM setting. Bootstrap approximation can be regarded as an alternative. However, in this case it would be easier to use bootstrap directly as described in Section 4. Finally, the application of the volume-of-tube formula results in two sources of errors; from the approximation itself and from the estimation of the constants, making the approximation less reliable.

5.2 Monte Carlo procedure

To deal with the problem of approximating unknown constants in (17), consider mixed model equations in (13) and (14) in Section 5.1. When using LMM for spline regression, Ruppert et al., (2003) proposed a simple numerical approach to construct confidence bands of one-dimensional nonparametric curves by the empirical approximation of (15), that is:

[TABLE]

Likewise, we apply expression (18) to simulate the distribution of $S_{0}$ in (8), and set

[TABLE]

where $\bar{\bm{c}}_{d}=(\bm{k}_{d}^{t},\bm{m}_{d}^{t})^{t}$. Afterwards, we draw $K$ realisations from normal distribution in (18), estimate the critical value $c_{S_{0}}(1-\alpha)$ by the $([(1-\alpha)K]+1)^{th}$ order statistic of $S_{MC}$ and construct MC SPI as follows

[TABLE]

We can similarly obtain a critical value for MT. The consistency of $\mathcal{I}^{MC}_{1-\alpha}$ follows from equation (15) which is a standard result for mixed models. The same results from the extreme value theory as in the proof of Proposition 1 might be invoked to prove the consistency for the maxima. Monte Carlo SPI are easy to implement and less computer intensive than bootstrap. Yet, they are less robust to departures from the normality of errors and random effects (cf. Section 6).

5.3 Bonferroni procedure

Classical simultaneous inference has been considered via Bonferroni procedure. If all statistics $(\tilde{\mu}_{d}-\mu_{d})/\sigma(\tilde{\mu}_{d})$ were independent Gaussian pivots, the critical value to construct SPI or MT could be selected as $c_{BO}(1-\alpha)=\Phi^{-1}(1-\alpha/2D)$. One may use quantiles from the normal instead of the t-distribution, because the number of mixed parameters is allowed to grow to infinity such that the latter distribution converges to the former, cf. the high-dimensional regression setting in Chernozhukov et al., (2013). Having retrieved the value of interest, a Bonferroni SPI is defined as

[TABLE]

While the same critical value might be used in MT procedure (9), it provides a weak control of FWER. Using Bonferroni’s methodology, we do not approximate the true distribution of statistic $S_{0}$ in equation (8). Hence, we work with $\hat{\sigma}(\hat{\mu}_{d})=\sqrt{\mathrm{mse}(\hat{\mu}_{d})}$ which is an estimated version of MSE defined in (6). An application of this procedure is simple and does not require almost any computational effort. It will be our benchmark under asymptotic independence of parameters. However, the results of Romano and Wolf, (2005) confirm that the method of Bonferroni performs poorly for correlated random variables, a problem that is even aggravated when allowing for spatio- and/or temporal dependencies, see our discussion in Section 8. Nevertheless, similarly as bootstrap SPI, Bonferroni bands are fairly robust to the distributional departures from normality of errors and random effect if the number of units in each cluster grows to infinity.

5.4 Beran procedure

Beran, (1988) developed a procedure to obtain balanced simultaneous intervals with an overall coverage probability $1-\alpha$ within the context of models without random effects. His technique is based on so called roots and bootstrapping to approximate their cumulative distribution functions (cdfs). We can follow Beran’s methodology and evaluate its performance under LMM. Suppose that $S_{d}$ is a root and consider $\max_{d=1,\dots,D}S_{d}$ as defined in equation (8). Let $F_{S_{d}}$ and $F_{S}$ be their respective cdfs. Furthermore, we denote with $F^{-1}_{S_{d}}(a)$ and $F^{-1}_{S}(a)$ the largest $a^{th}$ quantiles of $F_{S_{d}}$ and $F_{S}$. Beran, (1988) suggested bootstrap approximations $F^{*}_{S_{d}}$ and $F^{*}_{S}$ to obtain $D$ critical values defined as $c_{BEd}(1-\alpha)=F^{*-1}_{S_{d}}\{F^{*-1}_{S}(1-\alpha)\}$. Further details can be found in his paper. The SPI is then

[TABLE]

The critical values $c_{BEd}(1-\alpha)$ are not directly applicable for MT. Moreover, Beran’s method is as computer intensive as bootstrap SPI, but in comparison to the former it might provide a poorer coverage rate as its convergence in sup-norm is not guaranteed, cf. results in Section 6. Last but not least, it is not necessarily robust to the distributional departures from normality of errors and random effects.

6 Simulation experiments

We carry out simulations to examine finite sample properties of bootstrap (BS), Monte Carlo (MC), Beran (BE) and Bonferroni (BO) SPIs as well as to evaluate the empirical power of MT procedures under various scenarios. In particular, we analysed them under NERM and FHM. As far as the former is concerned, we set $x_{dj1}=1$, $x_{dj2}\sim U(0,1)$ $\forall$ $d\in[D]$ and $j\in[n_{d}]$, whereas under the FHM we set $x_{d1}=1$, $x_{d2}\sim U(0,1)$ $\forall d\in[D]$ with $\bm{\beta}=(1,1)^{t}$ in both models. The number of simulation runs is $I=2500$, each with $B=1000$ bootstrap samples. The covariates are fixed in all simulation runs. We considered small to medium numbers of clusters with $D\in\{15,30,60,90\}$.

When NERM is considered, we first set $n_{d}=5$ $\forall d\in[D]$, $e_{dj}\sim N(0,\sigma^{2}_{e})$, $u_{d}\sim N(0,\sigma^{2}_{u})$ such that the intraclass correlation coefficient $\mathrm{ICC}=\sigma^{2}_{u}/(\sigma^{2}_{u}+\sigma^{2}_{e})$ equals 1/3, 1/2 or 2/3 (see the first column of Table 1). Then we relax the modelling assumptions by allowing $e_{dj}$ and $u_{d}$ to deviate from normality to become heavy-tailed or asymmetric. Namely, we draw them from centred chi-square distribution with 5 degrees of freedom, student-t distribution with 6 degrees of freedom and skewed student-t distribution with 5 degrees of freedom and the skewness parameter equal to 1.25. We then rescale them to variances $\sigma^{2}_{e}$ and $\sigma^{2}_{u}$ indicated in parentheses in Table 3. Furthermore, we allow the number of units to grow with the number of clusters, cf. Jiang, (1998). The unusual choice of the skewed t-distribution is motivated by the data example in Section 7. In particular, it aims to mimic the pdf of estimated errors in our application, see the middle panel of Figure 1o. We consider a scenario with $D=52$ and $n_{d}=100$, i.e., the number of areas in the data example with $n_{d}=100$ being close to the median of the number of units across counties. We also evaluate the performance of our method for a smaller sample size with $D=26$ and $n_{d}=50$. Since the results hardly differ when estimating $\bm{\theta}$ using restricted maximum likelihood (REML) or the method of moments, we skip the latter.

We apply a similar setting as in Datta et al., (2005) in the simulation study with FHM. We suppose that random effects and errors are independent, centred and normally distributed with unknown variance $\sigma_{u}^{2}=1$ and known $\sigma^{2}_{e_{d}}$ as follows. Each fifth part of the total number of clusters is assigned to a different value for $\sigma^{2}_{e_{d}}$; in Scenario 1: $0.7$, $0.6$, $0.5$, $0.4$, $0.3$ and in Scenario 2: $2.0$, $0.6$, $0.5$, $0.4$, $0.2$. That is, we consider the case of known heteroscedasticity for errors. Variance $\sigma_{u}^{2}$ is estimated using REML, Henderson’s method (Prasad and Rao,, 1990) and the method of Fay and Herriot, (1979). We present results only for the former as other methods perform similarly for SPI and MT. All simulated scenarios are almost optimal settings for the Bonferroni procedure as the mixed parameter estimates are asymptotically independent. Therefore we can take it as a benchmark (cf. comments in Section 5.3).

We use three criteria to evaluate the performance of different methods to construct SPI: the empirical coverage probability (ECP), the average width (WS), and the average variance of widths (VS)

[TABLE]

ECP is the percentage of times all cluster-level parameters are inside their SPI. On the other hand, WS is calculated for each cluster over the widths of the intervals from $I$ simulations, and averaged over all clusters to obtain an aggregated indicator. Lower values of WS are preferable. Finally, for assessing their variability, we compute the variance of widths over the simulations, and average them over all clusters (VS). We prefer lower values of VS, as they would indicate that the length of intervals is stable.

Last but not least, in practice, $c_{BS}(1-\alpha)$ and $c_{BE}(1-\alpha)$ are approximated by $[\{(1-\alpha)B\}+1]^{th}$ order statistics of the empirical bootstrap distribution. In addition, to construct $\mathcal{I}^{MC}_{1-\alpha}$ in (19) we can use $\bm{g}_{1}$ or the variance expression from the denominator in (15). Since the numerical differences were negligible, we present results only for the latter.

Table 1 shows the numerical results of our criteria to compare the performance of different methods when errors and random effects are normally distributed. Under these scenarios, BS attains the nominal level of 95% even for a small number of clusters ($D=15$). Yet, due to the overestimation of variability of the cluster parameters, this method suffers from an overcoverage when ICC$=1/3$ for $D=15$ and $D=30$. Furthermore, although our simulations constitute a nearly optimal design for the Bonferroni method, BO exhibits almost always undercoverage. MC has worse performance, the convergence to the nominal level is slower with ECP oscillating around 94% only for $D=60$ and $D=90$ when ICC$=2/3$ or ICC$=1/2$. It does not attain the nominal coverage under the third scenario. Moreover, BE diverges which might by explained by the lack of convergence in sup-norm of the Beran’s procedure. The second part of Table 1 summarizes results for WS and VS. As expected, the width increases with growing $D$. Nonetheless, the speed of this increase is moderate; with a growing number of areas, the SPI has to cover more parameters, but the estimate of variability decreases (for more details, see Reluga et al.,, 2021). When we consider VS, we conclude that BS is more variable than other methods for $D=15$, but this difference decreases for increasing $D$. In SAE, undercoverage is often considered a more severe type of error than overcoverage, partly due to the difficulties to detect and alleviate it (Yoshimori,, 2015). On the other hand, overcoverage is often a result of an excessive variability in small samples which is illustrated in Table 1. Having this in mind, we conclude that BS is the most satisfactory method.

Table 2 shows the performance of SPI constructed for a subset of $D^{\prime}=D/5$ clusters. It illustrates finite sample performance of Corollary 2. Since the simulations under other scenarios led to the same conclusions, we only consider $\sigma^{2}_{u}=\sigma^{2}_{e}=1$. In each simulation run we constructed SPI for $D^{\prime}$ areas, but used all data to compute variances, fixed and random effects. The ECP is similar as in Table 1. In contrast, the widths of SPI are narrower than those in Table 1, because they are constructed to cover simultaneously only $D^{\prime}<D$ mixed parameters. This empirical study confirms a practical relevance of our proposal. In fact, it shows that one can construct reliable SPI or conduct MT for an arbitrary subset of cluster-level parameters.

While the vast majority of the SAE literature relies heavily on the normality of random effects and errors, especially regarding MSE estimation and CPI construction, this assumption may be violated in practice. Thus, we conduct a robustness study regarding departures from normality of errors and/or random effects. The first part of the empirical results of our criteria under this setting is presented in Table 3. Further simulation results are deferred to the supplementary material. First, but not surprisingly, the overall performance of all methods is worse than in Table 1, especially for asymmetric $\chi^{2}$ distributed errors. Second, BS and BO are still superior to all other methods. In addition, in case of chi-square distributed departures, the coverage is higher for ICC$=1/3$ with $D=15$ due to overestimated variability, then it drops for $D=30$ and $D=60$, and increases for $D=90$ in accordance with asymptotic theory. Importantly, under the scenario which mimics the data application, i.e., with skewed t-distributed errors and normal random effects, the ECP is close to the nominal level. However, we must conclude that the considered SPIs do not attain the nominal coverage probability if errors exhibit more severe deviations from normality than we observe in our application, irrespective of the presence of deviations from normality of random effects. The issue of undercoverage might be alleviated by the use of more sophisticated bootstrapping scheme requiring different theoretical derivations which is beyond the scope of this paper (cf. Reluga,, 2020). Simulations therein and in our supplementary material confirm that the deviation from normality of random effects hardly affects the coverage of SPI. A similar conclusion was drawn by McCulloch and Neuhaus, (2011) in the study of the bias of estimated fixed effects and EBLUPs. When it comes to the right hand side of Table 3, WS decreases with a growing sample size due to the increase of $n_{d}$. Even though the critical values increase with the growing number of clusters, $\hat{\sigma}^{2}(\hat{\mu}_{d})$ decreases at a faster rate. For this reason the average width of intervals is decreasing too.

Let us revisit the differences between CPI and SPI. Figure 6 displays 95% bootstrap SPI in light blue and CPIs of Chatterjee et al., (2008) in dark blue. The critical values for CPIs have been calculated using parametric bootstrap, cf. Chatterjee et al., (2008) for details. In comparison to CPI, SPI covers all clusters with a certain probability. Black dots represent the true mixed parameters $\mu_{d}$. Out of thirty, three cluster-level parameters (eighth, twenty-second and twenty-third) are clearly outside of their CPI, and another four (first, seventh, eleventh and twenty-fifth) are on their boundary. It does not happen by chance or by the simulation design, but by the construction of CPIs: for $100(1-\alpha)$% CPI, about $100\alpha$% of the true mixed parameters are not covered by their intervals. Figure 6 illustrates even more severe case with 10% of the true parameters not covered by CPIs. In contrast, SPI contains all of the true mixed parameters. Moreover, SPI is not excessively wide compared to CPI. In fact, SPI seems just as wide as necessary; twenty-third cluster-level mean is right at the boundary. Undoubtedly, CPI and SPI are methodologically different and constructed to cover distinct sets with a certain probability. One can thus argue that their direct comparison is flawed and should not be investigated. We do not claim otherwise; rather, Figure 6 serves as an illustration of a practical relevance of SPI as a valid tool for comparing mixed parameters across clusters. Moreover, Figure 6 demonstrates that the cluster-wise inference can lead to erroneous conclusions once applied to perform joint statements or comparisons.

Regarding MT, Figure 1l displays the empirical power of bootstrap and MC based max-type tests for $H_{0}:\bm{\mu}=\bm{h}$ vs $H_{1}:\bm{\mu}=\bm{h}+\bm{1}_{D}\delta$. In the simulations, we simply set $\bm{h}:=\bm{\mu}$ under $H_{0}$, whereas under under $H_{1}$ we added a constant $\delta\in[-2,2]$ to each element of $\bm{h}$. As expected, ICC influences the Type II error – the curves are the steepest and almost not distinguishable for ICC$=2/3$. The bootstrap test performs significantly better when ICC$=1/2$ and ICC$=1/3$ for small and medium $D$. In contrast, MC based tests do not attain the nominal level under $H_{0}$ for small sample sizes. For larger $D$, the curves almost coincide under all three scenarios.

Let us finally turn to the analysis under FHM. Since the performance of MT under FHM leads to similar conclusions as for NERM, we restrict ourselves to present ECP, WS and VS for different SPIs. Table 4 displays the results. Bootstrap SPI suffers from overcoverage for small $D$, similarly to Bonferroni’s SPI. The overcoverage is probably caused by the same reasons as for NERM. Surprisingly, Bonferroni’s intervals fail to achieve the nominal level for larger numbers of clusters. Beran’s SPI diverge, whereas MC SPI exhibits undercoverage for $D=60$ and $D=90$.

7 Application to the household income data of Galicia

We consider the household income data of the Structural Survey of Homes of Galicia (SSHG) which contains many potentially correlated covariates. It is of great interest for the Galician Institute of Statistics (IGS), and the regional government alike, to study the household income across counties (comarcas), e.g., to adjust regional policies and resource allocations. The IGS provides direct design-based estimates and/or EBLUP of the average household income accompanied by their variability measures or area-specific confidence intervals. However, the joint consideration of county-level parameters or comparisons between them is often important too. We start from the classical design- and model-based area-wise analysis. Afterwards we complete it with the simultaneous inference for counties in Galicia.

The SSHG contains data on 23628 individuals within 9203 households which were collected in 2014 and published in 2015. It comprises information about the total income as well as different characteristics on individual and household level. The variable of interest is the monthly household income. This variable was obtained by taking the twelfth of the total yearly income which consists of paid work, own professional activity and miscellaneous benefits. Following Lombardía et al., (2018), we consider covariates age, education level, type of household, and variables indicating financial difficulties of the household at the end of a month. Galicia is divided into four provinces (A Coruña, Lugo, Ourense and Pontevedra) which are further divided into 53 counties, the small areas that constitute our clusters. There are eighteen counties in A Coruña, thirteen in Lugo, twelve in Ourense and ten in Pontevedra. As the SSGH does not contain data from the county Quiroga in Lugo, we limit the study to the remaining 52 counties. Table 5 displays descriptive statistics of the number of units across the counties of each province. The model based approach is motivated by the scarcity of data (in some counties less than 20 observations were collected).

Even though the SSHG does not produce official estimates of totals $Y^{dir}_{d}$, $X^{dir}_{di}$ and means $\bar{Y}^{dir}_{d}$, $\bar{X}^{dir}_{di}$ at the county level, we calculated them using:

[TABLE]

where $\hat{N}^{dir}_{d}$ stands for the estimate of the county size $N^{dir}_{d}$, $\mathcal{R}_{d}$ is the sample in county $d$ and $w_{j}$ is an official calibrated sample weight. In addition, we have $w_{j}=1/\pi_{j}$ where $\pi_{j}\neq 0$ is the first-order inclusion probability. We used the same design-based direct variance estimator as Lombardía et al., (2018), that is:

[TABLE]

Furthermore, we calculated the coefficient of variation (CV) of direct estimates at the county level. In twelve counties CV $>10\%$, and in three of them CV $>15\%$. A direct estimate is considered official, and thus publishable, if its CV is lower than a certain threshold set by a statistical office. For example, the Office for National Statistics in the UK sets this threshold to 20% for the labour force statistics (Lombardía et al.,, 2018). Although the CV of our estimates does not exceed this threshold, it is still high enough to consider a model-based framework.

We constructed design and model-based point and CPI estimates of monthly incomes. We employed (22) and (23) for the design-based estimation. Within the model-based framework, we considered (4) as a target parameter with $\bm{k}_{d}=\bar{\bm{X}}^{dir}_{d}$, $\bm{m}_{d}=1$ and calculated EBLUP in (5) by $\hat{\mu}_{d}=\hat{\bar{\bm{X}}}^{dir}_{d}\hat{\bm{\beta}}+\hat{\bm{u}}_{d}$ $\forall d\in[D]$. Since SSHG contains information on the household level, we can fit NERM to these data. Figure 7 shows design-and model-based point estimates of monthly household incomes together with 95% CPIs. Model-based CPIs were constructed using the parametric bootstrap (cf. Chatterjee et al.,, 2008). We can use CPIs to compare different methods for the same cluster-level parameter, but not to make comparisons across different counties. For a better presentation, we divided the plot into five panels based on the number of units in each county. First, we can see that the widths of both direct- and model-based intervals decrease with increasing sample size. Second, the widths of direct CPIs are much larger that their model-based counterparts. In fact, direct estimates for certain areas (for example, the fourth and sixth area in the first panel) are too wide to make any informative conclusion. This confirms the necessity of a model-based framework.

Table 6 displays the covariates with their standard deviations as well as the estimated coefficients with standard errors and p-values. We performed a variable selection in two stages. First, we selected a subset of covariates that exhibited the highest Spearman’s rank correlation with the household income. Afterwards, we applied a generalised AIC which uses a quasi-likelihood with generalised degrees of freedom, see Lombardía et al., (2017). The algorithm selected covariates describing characteristics of the household and characteristics of the head of household. The estimates of variance parameters are $(\hat{\sigma}_{e}^{2},\hat{\sigma}_{u}^{2})=(758558.60,19746.24)$.

It is well known that income data are right skewed. Unsurprisingly, our dependent variable exhibits this feature too. It is therefore popular to consider log-income or more sophisticated transformations. Since the naive back-transformation of the dependent variable could cause a serious bias due to the Jensen inequality, different estimation and inference methods were recently suggested by Rojas-Perilla et al., (2020) and references therein. In fact, log transformation of household income in our data example did not help to overcome the problem of skewness. We thus decided to proceed with the household income on the original scale, and assess the sensitivity of our method to the departures from normality of errors and/or random effects after fitting LMM. If the undertaken inference is not compromised be these departures, we shall continue with further analysis. We carried out statistical tests and analysed diagnostic plots. The left panel of Figure 1o displays a diagnostic plot of Lange and Ryan, (1989) using standardised empirical Bayes estimates of the random effects in a weighted normal QQ plot; it supports the adequacy of the normality assumption. Moreover, the p-values of Kolmogorov-Smirnov and Shapiro-Wilk tests, which are $0.997$ and $0.944$, confirm this conclusion. Regarding the normality of errors, the middle and the right panels of Figure 1o present Cholesky residuals (Jacqmin-Gadda et al.,, 2007). The uncorrelated Cholesky residuals are constructed by multiplying $\bm{y}-\bm{X}\hat{\bm{\beta}}$ by the Cholesky square root of the variance matrix. A right tail is visible in both panels. The p-values of Kolmogorov-Smirnov and Shapiro-Wilk tests are $<0.000$. Nevertheless, in the middle panel of Figure 1o we can see that the kernel density estimate of the pdf of the skewed errors has a long, but not a thick tail. Our simulation results in Table 3 indicate that such departure is not problematic for our bootstrap-based SPI. In fact, SPI is quite robust to these departures, and a good coverage probability is still provided in comparison to other methods. Last but not least, we did not assess the robustness of bootstrap CPIs to the departure from normality of errors, because they are not the topic of this article; we present them for illustrative reasons.

Figure 7 displays bootstrap CPI as developed by Chatterjee et al., (2008), together with BS SPI for the county-level averages of monthly household income in Galicia. We can see a lot of variability over the estimates. Evaluating the results of CPI (blue) versus SPI (light blue), it is apparent that the cluster-wise prediction intervals are not adequate to address either a joint consideration or a comparison of the counties. If we consider, for example, the counties of A Fisterra and Noia (7th and 8th regions of the second panel in the black rectangle), the CPIs indicate significantly different incomes, whereas the SPIs do not support this claim. Moreover, there are other counties (practically in each panel) for which CPIs would insinuate significant differences whereas statistically valid SPIs do not confirm this conclusion. Nevertheless, SPIs are not unnecessarily wide for practical use. We detect significant and valid differences between several interval estimates.

Figure 1q presents maps with lower and upper limits of bootstrap SPI. The boundaries are classified into one of five categories which were built using $0.2$, $0.4$, $0.6$ and $0.8$ quantiles of the point estimates. We observe a substantial variation of average household income over the counties. Lower and upper boundaries of the interval estimates for the counties of A Coruña, Lugo, Vigo (with a large number of units) and Lower Miño are classified into the richest category; they are indicated with ellipsoid in the second and the last panel of Figure 7. In contrast, there is a group of eight counties (three in the centre, one in the west and four in the south of Galicia) which are classified to the poorest category in the left panel and the second poorest category in the middle panel. Right panel of Figure 1q presents CPIs and two different types of SPIs for eight poorest counties. In particular, SPI refers to the interval estimate constructed for 52 counties (estimates for all 52 counties are plotted in Figure 7). In contrast, SPI S refers to the interval estimate constructed for a subset of eight poorest counties, but using all data to estimate fixed parameters and predict random effects (cf. Corollary 2 and simulations in Table 2). As expected, SPI S is still wider than CPI, but much narrower than SPI. Moreover, SPI S permits for a valid comparative inference for a subset of these poorest areas. In fact, we can conclude that the differences in the average household income are not statistically significant. Finally, CPI should not be used to make maps in Figure 1q, as this would suggest that we were allowed to compare them.

Finally, it would be interesting to investigate whether the monthly income of the households which reported difficulties in coming to the end of the month is significantly different from the monthly income of the households which did not struggle with this issue. For example, in touristic areas households spend more such that they might face some difficulties without being poorer. The outcome of such test might then be used to develop a more targeted policy. Our MT procedure might be readily applied to support or disprove the hypothesis of no difference in monthly income between two mentioned types of households. To test this hypothesis we take clusters created from a cross-section of counties and difficulty status. Therefore we apply our developed methodology to $2\times 52=104$ counties by difficulty status. More specifically, we consider $\bm{\mu}\in\mathbb{R}^{104}$ and test $H_{0}:\bm{A\mu}=\bm{0}_{104}$ versus $H_{1}:\bm{A\mu}\neq\bm{0}_{104}$, where $\bm{A}\in\mathbb{R}^{52\times 104}$ with rows that are composed of $104$-dimensional vectors $\bm{a}$ with a $1$ on the $2d-1$ place, $-1$ on the $2d$ place but [math] otherwise, where $d$ stands for a particular county. The test statistic is $t_{H}=\max_{d=1,\dots,D}|t_{H_{d}}|=7.569$ whereas the critical value is $c_{BH_{0}}(1-\alpha)=4.220$. That is, we clearly reject $H_{0}$ of no difference. The rejection of the null hypothesis is consistent with the variable selection procedure which has suggested to find this test outcome.

8 Conclusions

We introduce a practical method to construct SPI and MT procedures for mixed parameters under LMM. We illustrate its use and relevance in simulation studies and a data application within the framework of SAE. We theoretically derive two techniques based on bootstrap approximation of the distribution of the max-type statistic and the volume-of-tube formula. However, we proved that the latter is not directly operational. We further discussed various alternatives and assessed their empirical performance in the simulation study. Though slightly conservative for very small samples, the bootstrap-based SPI yield the most satisfactory results in our simulations. Moreover, it is quite robust to certain deviations from the normality assumptions. In addition, our bootstrap based max-type statistic is readily applicable for testing multiple statistical hypotheses which was illustrated by our simulation studies.

Accounting for the joint coverage probability (or the Type I error) for several or all cluster-level parameters makes SPI wider than CPI. However, only SPI are statistically valid for joint statements or comparisons between cluster-level parameters. Moreover, if one conducted studies with several surveys, SPI would contain all true parameters in $100(1-\alpha)\%$ of all studies, whereas CPI would not cover about $D\alpha$ of them in each survey. Our tools are equally applicable to any subset of the clusters while using all data for estimation and prediction.

Our method can be extended to account for more complex data structures such as LMM with spatial and/or temporal dependencies (Pratesi and Salvati,, 2008; Morales and Santamaría,, 2019) or highly skewed response variables (Moura et al.,, 2017). Furthermore, our general idea could be also applied for a comparative analysis of nonlinear indicators of inequality or poverty obtained after the transformation of the dependent variable (Rojas-Perilla et al.,, 2020) or benchmarked estimators under restrictions (Ugarte et al.,, 2009). Yet, for these cases the estimation of the variance components or MSE would have to be adjusted. Furthermore, within more complex modelling frameworks one would need to extend the asymptotic theory. In addition, different bootstrap schemes might be necessary to mimic the data-generation process in a suitable way (Field et al.,, 2008). The above-mentioned extensions are beyond the scope of this paper, but are an open field for future research.

Appendix A Appendix

A.1 Regularity conditions

1. R.1

$\bm{X}_{d}$ and $\bm{Z}_{d}$ are uniformly bounded such that $\sum_{d=1}^{D}\bm{X}^{t}_{d}\bm{V}_{d}^{-1}\bm{X}_{d}=\left\{O(D)\right\}_{(p+1)\times(p+1)}$. 2. R.2

Covariance matrices $\bm{G}_{d}$ and $\bm{R}_{d}$ have a linear structure in $\bm{\theta}$. 3. R.3

Convergence: $D\rightarrow\infty$, $\sup_{d\geqslant 1}n_{d}<<\infty$ and $\sup_{d\geqslant 1}q_{d}<<\infty$. 4. R.4

To ensure the nonsingularity of $\bm{\Sigma}_{\bm{\theta}}$, $0<\inf_{d\leqslant 1}\sigma^{2}_{e_{d}}\leqslant\sup_{d\leqslant 1}\sigma^{2}_{e_{d}}<\infty$ and $\sigma^{2}_{u}\in(0,\infty)$. 5. R.5

$\bm{b}_{d}^{t}=\bm{k}^{t}_{d}-\bm{o}_{d}^{t}\bm{X}_{d}$ with $b_{di}=O(1)$ for $i=1,\dots,p+1$. 6. R.6

$\{\frac{\partial}{\partial\theta_{j}}\bm{o}_{d}^{t}\bm{X}_{d}\}_{i}=O(1)$ for $j=1,\dots,h$ and $i=1,\dots,p+1$. 7. R.7

$\hat{\bm{\theta}}$ satisfies: $(i)$ $\hat{\bm{\theta}}-\bm{\theta}=O_{p}(D^{-1/2})$, $(ii)$ $\hat{\bm{\theta}}(\bm{y})=\hat{\bm{\theta}}(-\bm{y})$ and $(iii)$ $\hat{\bm{\theta}}(\bm{y}+\bm{Xr})=\hat{\bm{\theta}}(\bm{y})$ for any $\bm{r}\in\mathbb{R}^{p+1}$.

Furthermore, we will evoke Assumptions 1-4 from Chatterjee et al., (2008), which are quite technical but largely irrelevant in practice. For the sake of completeness, they are provided in our supplementary material.

A.2 Proof of Proposition 1

We concentrate on the consistency of the bootstrap SPI. The consistency of MT procedure follows straightforwardly with some changes of notation due to the correspondence between tests and interval estimates (for more details, see Corollary 2 in Reluga et al.,, 2021). To demonstrate Proposition 1, we make use of the result in Theorem CLL of Chatterjee et al., (2008). In this section, we use some notation from their paper if it is not in conflict with ours. Consider $S_{0d}$ in (8) and $S^{*}_{Bd}$ in (11) and let $\mathcal{L}_{d}(q)=P(S_{0d}\leqslant q)$ and $\mathcal{L}^{*}_{d}(q)=P^{*}(S^{*}_{Bd}\leqslant q)$, where $P^{*}(\cdot)$ stands for a probability measure induced by a parametric bootstrap. Under suitable regularity conditions Chatterjee et al., (2008) proved that $\mathcal{L}_{d}(q)=\Phi(q)+\gamma(q,\bm{\beta},\bm{\theta},n)+O(n^{-1})$ where $\gamma(\cdot)$ is some smooth function. Furthermore, $\mathcal{L}^{*}_{d}(q)$ admits almost identical, equally short expansion with $\bm{\beta}$, $\bm{\theta}$ replaced by $\hat{\bm{\beta}}$, $\hat{\bm{\theta}}$ . Since $\mathcal{L}_{d}(q)\approx\Phi(q)$, we can follow the same steps as Reluga et al., (2021) to prove the consistency of SPI. In particular, observe that

[TABLE]

The same arguments follow for $S^{*}_{B}$ with $P$ replaced by $P^{*}$. Moreover, the cdf of the standardised maxima of the normal distribution is in the domain of attraction of the Gumbel law. If we notice that the last term in (24) can be approximated by the Gumbel distribution if suitably standardised, the consistency of SPI follows by applying Poyla’s theorem which relates the convergence in law with sup-norm convergence (see Reluga et al.,, 2021, for more details).

Remark 1**.**

Chatterjee et al., (2008)** estimated fixed parameters using an ordinary least squares method, whereas in our paper we used a generalized least squares. As pointed out by the authors, an asymptotic expansion still holds as soon as the weighting matrices are smooth functions of $\bm{\theta}$, which we assume in R.2 above.

Remark 2**.**

In their original proof, Chatterjee et al., (2008) considered a modified version of $S_{0d}$, that is $S_{0d}=\hat{\sigma}_{T}^{-1}(T-\hat{\mu}_{T})$ where $T=\bm{f}^{t}(\bm{X\beta}+\bm{Zu})$ for a given fixed vector $\bm{f}$, and $\hat{\mu}_{T}=\bm{f}^{t}(\bm{X\hat{\beta}}+\bm{Z\hat{u}})$. This modification leads to some smoothing effects and results in a faster convergence rate. Nevertheless, as pointed out in Remark 6 in their paper, the analysis of the area-specific mixed effects leads to a sightly slower convergence, but is equally valid, and can be carried out along the same lines.

Appendix B Supplementary material

This section contains supplementary material to the main article which introduced simultaneous prediction intervals and multiple testing procedures for a joint, comparative analysis of mixed parameters under linear mixed models. We introduce an extension of the multiple testing procedurein order to guarantee a strong control of a family-wise error. We further state the theorem of Chatterjee et al., (2008) together with its assumptions and some additional details. In addition, we present a proof of Proposition 2. Finally, we discuss and analyze the consequences of deviations from standard model assumptions. The extension of the multiple testing procedure and the analysis of the departures from the model assumptions are studied empirically in simulations which complement numerical results presented in the main document.

Appendix C The step-down testing procedure: description and numerical results

The testing procedure described in Section 3 of the main document controls weakly for the family-wise error rate (FWER). If one aims at testing multiple hypotheses like

[TABLE]

then we can use the step-down technique of Romano and Wolf, (2005) which assures a strong control of the FWER. Originally, Romano and Wolf, (2005) developed his multiple testing procedure for fixed effects in a regression context. We adapt it to make it suitable for mixed parameters within the framework of LMM. More specifically, let $\Omega$ be the space for all data-generating processes and $w$ be the true one. Each $H_{0d}$ is equivalent to $\omega\in\Omega_{d}$ for some $\Omega_{d}\subseteq\Omega$. For any $\vartheta\subseteq[D]$, denote $\Omega^{\vartheta}=(\cap_{d\in\vartheta}\Omega_{d})\cap(\cap_{d\notin\vartheta}\Omega^{c}_{d})$ with $\Omega^{c}_{d}=\Omega\setminus\Omega_{d}$. The strong control of the FWER implies that

[TABLE]

We summarise the step-down procedure of Romano and Wolf, (2005) in the following way. Consider $t_{H_{d}}$, $d\in[D]$ defined in Section 3 of the main document and denote by $c_{\vartheta}(1-\alpha)$ an estimate for the $(1-\alpha)$-quantile of $\max_{d\in\vartheta}|t_{H_{d}}|$. Let $\vartheta(1)=[D]$ and reject all hypotheses $H_{0d}$ for which $t_{H_{d}}>c_{\vartheta(1)}(1-\alpha)$. If no hypothesis is rejected, then stop. Otherwise, let $\vartheta(2)$ be the set of null hypotheses not rejected at the first step. At step $k$, let $\vartheta(k)\subseteq[D]$ be the set of hypotheses not rejected at step $k-1$. Reject all hypotheses $H_{0d}$, $d\in\vartheta(k)$ for which $t_{H_{d}}>c_{\vartheta(k)}(1-\alpha)$. If no hypothesis is rejected, stop. Proceed in this way until the algorithm stops. Romano and Wolf, (2005) proved that their procedure controls strongly for FWER, that is

[TABLE]

A direct consequence of Proposition 1 from the main document is the following corollary.

Corollary 3**.**

Under assumptions of Proposition 1, the step-down procedure of Romano and Wolf, (2005) with bootstrap approximations of the critical value $c_{\vartheta}(1-\alpha)$ provides a strong control of the FWER, satisfying condition in (26).

In other words, Corollary 3 says that we can use bootstrap approximations of the critical value defined in Section 4 of the main document in the step-down procedure of Romano and Wolf, (2005) which controls strongly for FWER.

We compare the finite sample performance of the step-down method of Romano and Wolf, (2005) with Bonferroni’s procedure which is treated as a benchmark in the main document. To this end, we consider multiple two-sided testing hypotheses $H_{0,d}:\mu_{d}=h_{d}$ for all $d\in[D]$ for three scenarios under NERM without departures from normality. We assume that $\mu_{d}=h_{d}+1$ for $d=1,\dots,D/5$. The performance criterion is an empirical FWER while the nominal level of FWER is set to $\alpha=0.05$. The analysis of Table 7 confirms that the bootstrap step-down procedure controls strongly for FWER, whereas Bonferroni’s method fails to do it under $\sigma^{2}_{e}=0.5$, $\sigma^{2}_{u}=1$ for $D=60$. Except for the third scenario and $D=15$, bootstrap procedure provides lower or equally good FWER as the Bonferroni’s technique. Recall that the latter works relatively well, because we designed the simulations such that the simultaneous tests are asymptotically independent.

Appendix D Theorem CLL, its assumptions and additional details

To demonstrate Proposition 1, we made use of Theorem CLL of Chatterjee et al., (2008). For the sake of completeness, we state it below with the technical assumptions adapted to our notation and modelling problem. Similarly as in the main document, we use some of the notation of Chatterjee et al., (2008) whenever it is not in conflict with ours.

Theorem 1**.**

(Theorem 3.1 of Chatterjee et al., (2008)). Let $T=\bm{f}^{t}(\bm{X\beta}+\bm{Zu})$, $\bm{f}\in\mathbb{R}^{n}$ a vector and a conditional distribution of $T$ given $\bm{y}$ be $N(\mu_{T},\sigma^{2}_{T})$, where

[TABLE]

Let $T^{*}=\bm{f}^{t}(\bm{X}\hat{\bm{\beta}}+\bm{Z}\bm{u}^{*})$ be a bootstrap equivalent of $T$, $\mathcal{L}_{n}$ the cdf of $\hat{\sigma}_{T}^{-1}(T-\hat{\mu}_{T})$ with $\bm{\theta}$ replaced by $\hat{\bm{\theta}}$, and $\mathcal{L}^{*}_{n}$ the cdf of $\hat{\sigma}^{-1*}(T^{*}-\hat{\mu}^{*}_{T})$. Suppose Assumptions 1-4 hold and $p+1+h=o(n)$. Then, for any continuity point $v$

[TABLE]

Moreover, $\mathcal{L}_{n}(v)$ admits a short asymptotic expansion, i.e., $\mathcal{L}_{n}(v)=\Phi(v)+h^{2}n^{-1}\gamma(q,\bm{\beta},\bm{\theta})+O(h^{3}n^{-3/2})$ where $\Phi$ stands for the cdf of a standard Gaussian random variable.

Before stating the technical assumptions, recall that $\bm{V}=\bm{R}+\bm{ZG}\bm{Z}^{t}$ where $\bm{R}$ and $\bm{G}$ denote the covariance matrices of errors $\bm{e}$ and random effects $\bm{u}$, respectively. In addition, the vector of unknown variance components is $\bm{\theta}\in\mathbb{R}^{h}$. For the readers convenience, we keep the numbering of the assumptions as in Chatterjee et al., (2008).

Assumption 1: Conditions on vector $\bm{f}$.

[TABLE]

Furthermore, $\sigma_{T}^{2}=\bm{f}^{t}\bm{Z}(\bm{G}-\bm{GZ}^{t}\bm{V}^{-1}\bm{ZG})\bm{Z}^{t}\bm{f}>M>0$ for a constant $M>0$.

Assumption 2: Regularity conditions on design matrix $\bm{X}$.

[TABLE]

In addition, the smallest eigenvalue of $\bm{X^{t}X}/n$ is bounded away from zero.

Assumption 3: Conditions on the variance matrices $\bm{G}$, $\bm{R}$ and their estimates $\bm{G}(\hat{\bm{\theta}})$, $\bm{R}(\hat{\bm{\theta}})$.

The eigenvalues of $\bm{G}$, $\bm{R}$ lie in $(L^{-1},L)$ for some $L>1$, and the eigenvalues of $\bm{G}(\hat{\bm{\theta}})$, $\bm{R}(\hat{\bm{\theta}})$ lie in $(L^{-1}/2,2L)$. Furthermore, the eigenvalues of $\bm{V}$ lie in a compact set on the positive half of the real line. Finally, let $\bm{\Lambda}_{G}(\hat{\bm{\theta}})$, $\bm{\Lambda}_{R}(\hat{\bm{\theta}})$ be diagonal matrices implicitly defined via $\bm{G}(\hat{\bm{\theta}})=\bm{G}_{1}^{1/2}\bm{\Lambda}_{G}(\hat{\bm{\theta}})\bm{G}_{2}^{t/2}$ with $\bm{G}=\bm{G}_{1}^{1/2}\bm{G}_{2}^{t/2}$ and $\bm{R}(\hat{\bm{\theta}})=\bm{R}_{1}^{1/2}\bm{\Lambda}_{R}(\hat{\bm{\theta}})\bm{R}_{2}^{t/2}$ with $\bm{R}=\bm{R}_{1}^{1/2}\bm{R}_{2}^{t/2}$ depending only on $\bm{\theta}\in\mathbb{R}^{h}$. It holds that all entries of $\bm{\Lambda}_{G}(\hat{\bm{\theta}})$, $\bm{\Lambda}_{R}(\hat{\bm{\theta}})$ have three bounded continuous derivatives (partial derivatives with respect to the elements of $\hat{\bm{\theta}}$).

Assumption 4: Conditions on the estimator of the variance parameters $\hat{\bm{\theta}}$.

Define $\bm{E}=\sqrt{h/n}(\hat{\bm{\theta}}-\bm{\theta})$, then all moments of $\|\bm{E}\|$ are of order $O(1)$ and for the elements $E_{j}$ of vector $E$, it holds that $\forall\ j,l=1,...,h$

[TABLE]

Appendix E Alternative argument for consistency of bootstrap SPI

The consistency of $\mathcal{I}_{1-\alpha}^{BS}$ and Proposition 1 can be also based on a heuristic argument of Hall and Pittelkow, (1990) combined with the asymptotic expansion in Theorem 1. Ideally, $c_{S_{0}}(1-\alpha)$ would be determined from

[TABLE]

Since $\pi\left\{c_{S_{0}}(1-\alpha)\right\}$ is unknown, we approximate it by bootstrap such that

[TABLE]

with $\mathcal{W}=\{(y_{dj},\bm{x}_{dj},\bm{z}_{dj}),d\in[D],j\in[n_{d}]\}$. If we prove that $\pi$ and $\pi^{*}$ are asymptotically close up to order $O_{P}((h^{2}n^{-3/2})$, then it implies the same order of accuracy for $\mathcal{I}^{BS}_{1-\alpha}$. Define

[TABLE]

which can be represented as a finite number of unions and intersection of convex sets. This number is bounded uniformly for $D\geqslant 2$ and $c_{S_{0}}(1-\alpha)>0$. Observe that $\pi=\int_{\mathbb{R}}d\mathcal{L}_{n}$. Theorem 1 shows that for all continuity points $v$ the cdf’s $\mathcal{L}_{n}$ and $\mathcal{L}^{*}_{n}$ converge to the same limit at the desired speed, and the same speed is maintained in the convergence of $\pi$ and $\pi^{*}$. Since $\pi$ is defined as an integral of $d\mathcal{L}_{n}$ over $\mathbb{R}$, a direct consequence of Lemma 1 is $\sup_{k\in\mathbb{R}}|\pi(k)-\pi^{*}(k)|=O_{P}(h^{2}n^{-3/2})$.

Appendix F Proof of Proposition 2

To derive approximation formulas in Proposition 2, we suppose for the moment that the manifold $\mathcal{M}=\left\{\mathcal{Q}(\bm{c}),\bm{c}\in\mathcal{C}\right\}$ has no boundary, that is, Euler-Poincaré characteristic $\mathcal{E}=0$. In addition, $\bm{l}_{M}$, $\hat{\bm{l}}_{M}$, $\bm{e}_{M}$, $\hat{\bm{e}}_{M}$, $\lambda=\hat{\bm{e}}_{M}-\bm{e}_{M}$ and the other assumptions remain as defined in the main document. For some $\bm{x}\in\mathcal{X}$, $\bm{z}\in\mathcal{Z}$ and $\bm{c}\in\mathcal{C}$ the difference between $\bm{x}^{t}\bm{\beta}+\bm{z}^{t}\bm{u}$ and its estimate can be bounded by

[TABLE]

where $\eta(\bm{x})=\left|\left(\hat{\bm{l}}_{M}-\bm{l}_{M}\right)^{t}\bm{e}_{M}+\hat{\bm{l}}^{t}_{M}\lambda_{M}\right|$. If we normalize the first term on the right hand side, it is straightforward to see that

[TABLE]

which coincides with the formula of a Gaussian random variable in Section 5.1 of the main document. Here, $\bm{e}_{M}/\sigma_{e}$ is an $n$-vector of normally distributed random variables. Following the derivation of Sun and Loader, (1994) and Sun et al., (1999), one needs to choose $c_{VT}(1-\alpha)$ such that

[TABLE]

where $\xi=\inf\limits_{\bm{c}\in\mathcal{C}}\frac{||\hat{\bm{l}}_{M}||}{||\bm{l}_{M}||}$ is the minimum of the ratio between estimated $\hat{\bm{l}}_{M}$ and the true $\bm{l}_{M}$, and $\eta=\sup\limits_{\bm{c}\in\mathcal{C}}\frac{\eta(\bm{x})}{\sigma_{e}||\bm{l}_{M}||}$ accounts for the difference in the estimation of variance parameters. When $\bm{\theta}$ is obtained using some consistent estimator (for example REML), then Jiang, (1998) proved that $\hat{\sigma}_{e}$ is asymptotically independent of $Z$ and

[TABLE]

Therefore, $\xi$ and $\eta$ can be bounded by positive constants $\xi\leqslant\xi_{0}$ and $\eta\leqslant\eta_{0}$ in probability as $n\rightarrow\infty$, and we obtain the approximation

[TABLE]

Under our setting, the variable $\nu^{1/2}\hat{\sigma}/\sigma$ is approximately $\chi$ distributed with $\nu$ degrees of freedom and a pdf

[TABLE]

Since we deal with a Guassian random variable, we can adjust the approximations formulas of Sun, (1993) to retrieve a value of $c_{VT}(1-\alpha)$ for $p=1$, $p=2$ and $p\geqslant 3$. First of all, let us focus on the cases $p=1$ and $p=2$, where we need to find a solution of

[TABLE]

We develop two expressions using Taylor expansions. Let $c^{\prime}_{VT}(1-\alpha)=c_{VT}(1-\alpha)\xi_{0}$ and $\frac{c^{\prime}_{VT}(1-\alpha)w}{\nu^{1/2}}=j(w)$. Then for any $\eta\rightarrow 0$ we have

[TABLE]

Using a $\chi$ distribution, we have for ${\cal A}_{j}:=\int_{0}^{\infty}A_{j}f(w,\nu)dw$, $j=1,2,3$:

[TABLE]

To find an approximation for the model with $p\geqslant 3$, we modify a following expression

[TABLE]

Having calculated all of the necessary terms, we obtain final approximations;

for $p=1$:

[TABLE]

for $p=2$:

[TABLE]

for $p\geqslant 3$:.

[TABLE]

where $t_{\nu}$ is a t-distributed random variable with $\nu$ degrees of freedom, $F_{d_{1},d_{2}}$ an F-distributed random variable with parameters $d_{1}$ and $d_{2}$, $\kappa_{0}=\int_{\bm{c}\in\mathcal{C}}\left\lVert\mathcal{Q}^{\prime}(\bm{x})\right\rVert d{\bm{x}}$ the volume of the manifold $\mathcal{M}=\left\{\mathcal{Q}(\bm{c}),\bm{c}\in\mathcal{C}\right\}$, and $\zeta_{0}$ the boundary area of $\mathcal{M}$. Furthermore, $\kappa_{2}$ and $\zeta_{1}$ measure the curvatures of $\mathcal{M}$ and $\partial\mathcal{M}$ respectively, whereas $m_{0}$ measures the rotation angles of $\partial^{2}\mathcal{M}$. Finally, $\mathcal{E}$ is the Euler-Poincaré characteristic of $\mathcal{M}$. For the manifold $\mathcal{M}$ with the boundary we need to add the boundary’s correction terms $\mathcal{E}P\{|t_{v}|>c_{VT}(1-\alpha)\}$ and $2\mathcal{E}P\{|t_{v}|>c_{VT}(1-\alpha)\}$ for $p=1$ and $p=2$ respectively.

Appendix G Departures from model assumptions

In what follows we first extend simulation results in order to complement the findings regarding departures from normality of the distribution of errors and random effects. Afterwards we discuss potential deviations from the mixed model assumptions, their consequences and remedies.

Various authors (see, among others, a review paper of McCulloch and Neuhaus,, 2011) found that the departures from normality of random effects do not have a strong effect. In contrast, the lack of normality of error terms can have a serious impact, for example on the coverage probability of prediction intervals. Nevertheless, the results of the simulation experiment in Table 3 of the main document indicate that our bootstrap SPI works reasonably well under the departures from normality which we seem to face in our case study. However, the performance of all methods deteriorates if we deal with more sever departures from normality. Therefore we studied this problem further. Tables 8 and 9 present additional simulation results of the evaluating criteria introduced in Section 6 of the main document. In the simulation experiments shown in the main document in Table 3 we discovered the same phenomena. In particular, the coverage probability depends not just on the sample size (both, $D$ and $n_{d}$) but also on ICC and the combination of distributions for $e_{dj}$ and $u_{d}$. Furthermore, in accordance with the remarks of McCulloch and Neuhaus, (2011), it seems that our SPI is more affected and exhibits more severe undercoverage when error terms deviate from normality than when this happens for the random effects $\bm{u}_{d}$. As we mentioned in the main document, Reluga, (2020) studied the performance of alternative bootstrap procedures which seem to achieve a stronger robustness against non-normalities, but they require several changes of the proposed procedure. Moreover, they would require a different asymptotic theory. Furthermore, parametric bootstrap SPI fails when the distribution of errors and random effects is contaminated by outliers belonging to a different distribution with a large variance. In fact, this was confirmed by some simulations not shown in this document. In case of contaminated data, outlier robust estimator should be applied; see Chambers and Tzavidis, (2006); Jacqmin-Gadda et al., (2007); Sinha and Rao, (2009). However, the development of simultaneous inference techniques for robust estimators is beyond the scope of this article.

Many other violations of our model assumptions are possible. The proposed methods for constructing simultaneous prediction intervals and tests are model based. In fact, this characteristic is widely shared by inferential methods in many fields and is not particular to the bootstrap based technique which we developed. Nonetheless it is worth to discuss the consequences of model violations. Apart from the most studied departures from distributional assumptions, the second interesting deviation might be the existence of nonlinearities and interactions.

First, Reluga et al., (2021) consider an extension to a nonlinear link function and generalised linear mixed model. Assuming this modelling framework, the authors did not use EBLUP, but examined the empirical best predictor (Jiang and Lahiri,, 2001). More sophisticated nonlinearities, potentially together with interactions, are in the focus of interest when switching to non- and semiparametric mixed models, see, for example, Lombardía and Sperlich, (2008) and Opsomer et al., (2008) for estimation with kernels and splines, respectively. Moreover, González-Manteiga et al., (2013) discussed in detail the application of bootstrap in samiparametric mixed effects models, whereas Sperlich and Lombardía, (2010) investigated nonparametric specification testing for small area statistics. One could certainly imagine an extension of our methods to those models, but both, the statistical properties as well as the computational challenges would be more involved. Yet, we consider these to be interesting subjects for future research.

We should also mention model specifications which deal with potential dependencies between the random components, or dependencies between covariates and area effects. In the context of the former, we can think about models with temporal and/or spatial correlation, see Signh et al., (2005) and Moura and Migon, (2002), whereas the latter is examined by Lombardía and Sperlich, (2012). They propose a semiparametric device to control for potential correlation between covariates and cluster effects. Parametric counterparts to their semiparametric device fall into the class of models we consider in our main document.

Finally, the question of model (or variable) selection in mixed models is a related aspect to consider once extending the scope of our theory (see Müller et al.,, 2013, for an extensive review of model selection techniques). Nevertheless, due to the additional variability coming from the selection process, one would need to first theoretically investigate a post-selection distribution of mixed effects in order to develop statistically valid simultaneous intervals and tests (see, for example Charkhi and Claeskens,, 2018, for a valid post-selection inference for fixed parameters).