3.1. Formulation

This section introduces the area-level measurement error Poisson mixed (MEPM) model. Let U denote a finite population that is segmented into D domains indexed by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1, \ldots , D$$\end{document} . We assume a situation where (i) survey data for a variable of interest y, and (ii) either administrative records or survey data regarding auxiliary variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}= \{x_1, \ldots , x_p\}$$\end{document} are available. The area-level MEPM model is composed of three hierarchical stages. In the first stage a model, called sampling model, gives the distribution of the sample counts \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_d$$\end{document} . Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_d$$\end{document} be the domain probability of an event of interest and let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _d=n_d$$\end{document} be the domain sample size (or a size parameter in other counting setups). The sampling model indicates that the distribution of the target variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{d}$$\end{document} is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} y_{d}\sim \text{ Poiss }(\nu _{d}p_{d}),\quad d=1,\ldots ,D. \end{aligned}$$\end{document}

That is to say, the target variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_d$$\end{document} gives the observed domain counts in the sample for the event of interest. Against the background of Sect. 2, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_d$$\end{document} refers to the number of individuals in the sample that have a MDD in domain d. In the second stage, the logarithm of the mean parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _d=\nu _dp_d$$\end{document} (the natural parameter in Poisson regression models) is assumed to vary linearly with the p area-level auxiliary variables contained in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}$$\end{document} . In the second stage, we consider a set of random effects \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{v_{d}:\,\,d=1,\ldots ,D\}$$\end{document} that are independent and identically distributed (i.i.d.) according to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_d \sim N(0,1)$$\end{document} . In matrix notation, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{v}=\underset{1\le d \le D}{\text{ col }}(v_{d})\sim N_D(\varvec{0},\varvec{I}_D), \quad f_v(\varvec{v})=(2\pi )^{-D/2}\exp \left\{ -\frac{1}{2}\,\varvec{v}^\prime \varvec{v}\right\} , \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_v(\varvec{v})$$\end{document} gives the probability density function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{v}$$\end{document} . The linking model is specified as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \log \mu _{d}=\log \nu _{d}+\log p_{d}=\log \nu _{d}+ {{\tilde{\varvec{x}}}}_{0,d}\varvec{\beta }_0+{{\tilde{\varvec{x}}}}_{1,d}\varvec{\beta }_1+\phi v_{d},\quad d=1,\ldots ,D, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\tilde{\varvec{x}}}}_{0,d}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\tilde{\varvec{x}}}}_{1,d}$$\end{document} are a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1\times p_1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1\times p_2$$\end{document} row vectors, respectively, containing the true aggregated values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=p_1+p_2$$\end{document} auxiliary variables for domain d. The terms

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{\beta }_0=\underset{1\le k \le p_1}{\hbox {col}}(\beta _{k}), \quad \varvec{\beta }_1=\underset{p_1+1\le k \le p}{\hbox {col}}(\beta _{k}), \end{aligned}$$\end{document}

are the corresponding column vectors of regression parameters. By defining \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\tilde{\varvec{x}}}}_{d}=({{\tilde{\varvec{x}}}}_{0,d},{{\tilde{\varvec{x}}}}_{1,d})$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }=(\varvec{\beta }_0^\prime ,\varvec{\beta }_1^\prime )^\prime $$\end{document} , the linking model can be written in the simpler form.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \log \mu _{d}=\log \nu _{d}+\log p_{d}=\log \nu _{d}+ {{\tilde{\varvec{x}}}}_{d}\varvec{\beta }+\phi v_{d},\quad d=1,\ldots ,D. \end{aligned}$$\end{document}

We assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\tilde{\varvec{x}}}}_{0,d}$$\end{document} is known and equal to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_{0,d}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1,\ldots ,D$$\end{document} . With respect to practice, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_{0,d}$$\end{document} may correspond to administrative records that are taken from a register. Further, we assume that the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\tilde{\varvec{x}}}}_{1,d}$$\end{document} ’s are unknown random vectors with predictors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_{1,d}$$\end{document} ’s calculated from independent data sources. These data sources might be data sets obtained from surveys with larger sample sizes than the survey where the dependent variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_d$$\end{document} ’s are recorded. In the third stage, we consider the measurement error model

(1) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{\tilde{\varvec{x}}}}_{1,d}^\prime =\varvec{x}_{1,d}^\prime +\varvec{u}_d,\quad \varvec{u}_d\sim N_{p_1}(\varvec{0}, \varvec{\Sigma }_d),\quad d=1,\ldots ,D, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_{1,d}$$\end{document} is a row vector containing unbiased predictors of the components of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\tilde{\varvec{x}}}}_{1,d}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{u}_d$$\end{document} is a column vector of random measurement errors, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_d$$\end{document} is the known covariance matrix of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{u}_d$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{u}_d$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_{1,d}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1,\ldots ,D$$\end{document} , are independent. In practice, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_d$$\end{document} is typically unknown and subsequently estimated from the same unit-level survey data as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_{1,d}$$\end{document} . We finally assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_{1,d}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{u}_d$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_d$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1,\ldots ,D$$\end{document} , are independent, but we only introduce inference procedures conditionally on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_{d}=(\varvec{x}_{0,d},\varvec{x}_{1,d})$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1,\ldots ,D$$\end{document} . In matrix notation, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{y}=\underset{1\le d \le D}{\text{ col }}(y_{d}),\quad \varvec{X}=\underset{1\le d \le D}{\hbox {col}}(\varvec{x}_{d}),\quad \varvec{u}=\underset{1\le d \le D}{\text{ col }}(\varvec{u}_{d}), \quad \varvec{u}_d=\underset{1\le k \le p_1}{\text{ col }}(u_{dk}). \end{aligned}$$\end{document}

Based on these definitions, the probability density function of the measurement errors is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f_u(\varvec{u})=(2\pi )^{-Dp_1/2}\bigg (\prod _{d=1}^D|\varvec{\Sigma }_{d}|\bigg )^{-1/2} \exp \bigg \{-\frac{1}{2}\,\sum _{d=1}^D\varvec{u}_d^\prime \varvec{\Sigma }_d^{-1}\varvec{u}_d\bigg \}. \end{aligned}$$\end{document}

The area-level MEPM model can be expressed as a single model in the form

(2) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} y_{d}|_{v_{d},\varvec{u}_d,\varvec{x}_d}\sim \text{ Poiss }(\nu _{d}p_{d}),\,\, \log \mu _{d}=\log \nu _{d}+\log p_{d}=\log \nu _{d}+\varvec{x}_d\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_{d}, \end{aligned}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{d}=\exp \left\{ \varvec{x}_{d}\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_{d}\right\} $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1,\ldots ,D$$\end{document} . From the model assumptions, it holds that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&P(y_{d}|\varvec{v},\varvec{u},\varvec{X})=P(y_{d}|v_d,\varvec{u}_d,\varvec{x}_d)=\frac{1}{y_{d}!}\exp \{-\nu _{d}p_{d}\} \nu _{d}^{y_{d}} p_{d}^{y_{d}}, \\&P(\varvec{y}|\varvec{v},\varvec{u},\varvec{X})=\prod _{d=1}^DP(y_{d}|v_d,\varvec{u}_d,\varvec{x}_d), \\&P(\varvec{y}|\varvec{X})=\int _{R^{D}}\int _{R^{Dp_1}} P(\varvec{y}|v_d,\varvec{u}_d,\varvec{x}_d) f_{v}(\varvec{v})f_u(\varvec{u})\,d\varvec{v}\,d\varvec{u}=\int _{R^{D(p_1+1)}} \psi (\varvec{y},\varvec{v},\varvec{u})\,d\varvec{v}\,d\varvec{u}, \end{aligned}$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \psi (\varvec{y},\varvec{v},\varvec{u})= & {} (2\pi )^{-\frac{D}{2}}\exp \left\{ \frac{-\varvec{v}^\prime \varvec{v}}{2}\right\} (2\pi )^{-Dp_1/2}\bigg (\prod _{d=1}^D|\varvec{\Sigma }_{d}|\bigg )^{-1/2}\exp \bigg \{-\frac{1}{2}\,\sum _{d=1}^D\varvec{u}_d^\prime \varvec{\Sigma }_d^{-1}\varvec{u}_d\bigg \} \\&\cdot&\prod _{d=1}^D\frac{\exp \{-\nu _{d}p_{d}\} \nu _{d}^{y_{d}} \exp \left\{ y_{d}(\varvec{x}_{d}\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_{d})\right\} }{y_{d}!} \\= & {} (2\pi )^{-\frac{D}{2}}\exp \left\{ \frac{-\varvec{v}^\prime \varvec{v}}{2}\right\} (2\pi )^{-Dp_1/2}\bigg (\prod _{d=1}^D|\varvec{\Sigma }_{d}|\bigg )^{-1/2}\exp \bigg \{-\frac{1}{2}\,\sum _{d=1}^D\varvec{u}_d^\prime \varvec{\Sigma }_d^{-1}\varvec{u}_d\bigg \} \\&\cdot&\big (\prod _{d=1}^Dy_{d}!\big )^{-1} \exp \left\{ \sum _{d=1}^D\big \{-\nu _{d}\exp \{\varvec{x}_{d}\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_{d}\}+y_{d}\log \nu _{d}\big \}\right\} \\&\cdot&\exp \left\{ \sum _{k=1}^p\big (\sum _{d=1}^Dy_{d}x_{dk}\big )\beta _k+\sum _{k=1}^{p_1}\big (\sum _{d=1}^Dy_{d}u_{dk}\big )\beta _{ k} +\phi \sum _{d=1}^Dy_{d}v_{d}\right\} . \end{aligned}$$\end{document}

3.2. Empirical Best Prediction

This section derives the EBPs of the domain probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_d$$\end{document} and the mean parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _d = \nu _d p_d$$\end{document} . We proceed according to the following basic strategy. First, the best predictors (BPs) are calculated under the model under the preliminary assumption that the model parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }= (\varvec{\beta }', \phi )'$$\end{document} are known. Afterward, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} is replaced by a consistent estimator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\varvec{\theta }}} = ({\hat{\varvec{\beta }}}', {\hat{\phi }})'$$\end{document} to obtain the EBP. We show how to perform consistent model parameter estimation in Sect. 4. We start with the EBP of the domain probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_d$$\end{document} . Recall that the conditional distribution of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{y}$$\end{document} , given \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{v}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{u}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{X}$$\end{document} , is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(\varvec{y}|\varvec{v},\varvec{u},\varvec{X})=\prod _{d=1}^D P(y_{d}|v_{d},\varvec{u}_d,\varvec{x}_d)$$\end{document} , where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P(y_{d}|v_{d},\varvec{u}_d,\varvec{x}_d)= & {} \frac{1}{y_{d}!}\exp \{-\nu _{d}p_{d}\} \nu _{d}^{y_{d}} p_{d}^{y_{d}} \\= & {} \frac{\nu _d^{y_d}}{y_d!} \exp \left\{ y_{d}(\varvec{x}_{d}\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_d)-\nu _{d}\exp \{\varvec{x}_{d}\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_{d}\}\right\} . \end{aligned}$$\end{document}

The probability density functions of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_d$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{u}_d$$\end{document} are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(v_d)=(2\pi )^{-1/2}\exp \big \{-\frac{1}{2}\, v_{d}^2\big \},\quad f(\varvec{u}_d)=(2\pi )^{-p_1/2}|\varvec{\Sigma }_{d}|^{-1/2} \exp \Big \{-\frac{1}{2}\,\varvec{u}_d^\prime \varvec{\Sigma }_d^{-1}\varvec{u}_d\Big \}. \end{aligned}$$\end{document}

The BP of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{d}$$\end{document} is its conditional expectation under the model, that is, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{p}}_{d}(\varvec{\theta })=E_\theta [p_{d}|\varvec{y},\varvec{X}]$$\end{document} . In this case, we have that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\theta [p_{d}|\varvec{y},\varvec{X}]=E_\theta [p_{d}|y_{d},\varvec{x}_d]$$\end{document} and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} E_\theta [p_{d}|y_{d},\varvec{x}_d]= & {} \frac{ \int _{R^{p_1+1}}\exp \{\varvec{x}_{d}\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_{d}\} P(y_{d}|v_{d},\varvec{u}_d,\varvec{x}_d)f(v_{d})f(\varvec{u}_d)\,\mathrm{d}v_{d} \mathrm{d}\varvec{u}_{d}}{\int _{R^{p_1+1}}P(y_{d}|v_{d},\varvec{u}_d,\varvec{x}_d)f(v_{d})f(\varvec{u}_d)\,\mathrm{d}v_{d}\mathrm{d}\varvec{u}_{d}}\\= & {} \frac{N_{d}(y_{d},\varvec{\theta })}{D_{d}(y_{d},\varvec{\theta })}. \end{aligned}$$\end{document}

The numerator and denominator of the fraction are given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} N_{d}(y_{d},\varvec{\theta })= & {} \int _{R^{p_1+1}}\exp \Big \{(y_{d}+1)(\varvec{x}_{d}\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_{d}) - \nu _{d}\exp \big \{\varvec{x}_{d}\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_{d}\big \}\Big \}\\&\quad \times f(v_{d})f(\varvec{u}_d)\,\mathrm{d}v_{d}\mathrm{d}\varvec{u}_{d}, \\ D_{d}(y_{d},\varvec{\theta })= & {} \int _{R^{p_1+1}}\exp \Big \{y_{d}(\varvec{x}_{d}\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_{d}) - \nu _{d}\exp \{\varvec{x}_{d}\varvec{\beta }+\varvec{u}_{d}^\prime \varvec{\beta }_1+\phi v_{d}\}\Big \}\\&\quad \times f(v_{d})f(\varvec{u}_d)\,\mathrm{d}v_{d}\mathrm{d}\varvec{u}_{d}. \end{aligned}$$\end{document}

In principle, the EBP is then calculated according to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{p}}_d = N_d(y_d, {\hat{\varvec{\theta }}})/ D_d(y_d, {\hat{\varvec{\theta }}})$$\end{document} . However, both integrals contained in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\theta [p_{d}|y_{d},\varvec{x}_d]$$\end{document} do not have a closed-form solution. As a consequence, the EBP cannot be calculated exactly, but has to be approximated instead. This can be done by using Monte Carlo integration. That is to say, we choose a number of 2L symmetric locations at which numerator and denominator are evaluated within the support of the multivariate probability density of both random effects and measurement errors. Afterward, we average the obtained functional values and obtain an approximation to the corresponding integrals. This is summarized in the subsequent procedure.

1. 1. Estimate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\varvec{\theta }}}=({{\hat{\varvec{\beta }}}}',{\hat{\phi }})$$\end{document} .

2. 2. Generate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{d}^{(\ell )}$$\end{document} i.i.d. N(0, 1), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{u}_{d}^{(\ell )}$$\end{document} i.i.d. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{p_1}(\varvec{0},\varvec{\Sigma }_{d})$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{d}^{(L+\ell )}=-v_{d}^{(\ell )}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{u}_{d}^{(L+\ell )}=-\varvec{u}_{d}^{(\ell )}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1,\ldots ,D$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell =1,\ldots ,L$$\end{document} .

3. 3. Calculate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{p}}_{d}({{\hat{\varvec{\theta }}}})={\hat{N}}_{d}/{\hat{D}}_{d}$$\end{document} , where

   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\hat{N}}_{d}= & {} \frac{1}{2L}\sum _{\ell =1}^{2L}\exp \left\{ (y_{d}+1) (\varvec{x}_{d}{{\hat{\varvec{\beta }}}}+\varvec{u}_{d}^{\prime (\ell )}{{\hat{\varvec{\beta }}}}+{{\hat{\phi }}} v_{d}^{(\ell )}) - \nu _{d}\exp \{\varvec{x}_{d}{{\hat{\varvec{\beta }}}}+\varvec{u}_{d}^{\prime (\ell )}{{\hat{\varvec{\beta }}}}+{{\hat{\phi }}} v_{d}^{(\ell )}\}\right\} , \\ {\hat{D}}_{d}= & {} \frac{1}{2L}\sum _{\ell =1}^{2L}\exp \left\{ y_{d}(\varvec{x}_{d}{{\hat{\varvec{\beta }}}}+\varvec{u}_{d}^{\prime (\ell )}{{\hat{\varvec{\beta }}}}+{{\hat{\phi }}} v_{d}^{(\ell )}) - \nu _{d}\exp \{\varvec{x}_{d}{{\hat{\varvec{\beta }}}}+\varvec{u}_{d}^{\prime (\ell )}{{\hat{\varvec{\beta }}}}+{{\hat{\phi }}} v_{d}^{(\ell )}\}\right\} . \end{aligned}$$\end{document}

It can be shown that for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L \rightarrow \infty $$\end{document} , it holds that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{N}}_{d} / {\hat{D}}_{d} \overset{a.s.}{\rightarrow } N_d(y_d, {\hat{\varvec{\theta }}})/ D_d(y_d, {\hat{\varvec{\theta }}})$$\end{document} (Calfisch, Reference Calfisch1998). That is to say, for L chosen sufficiently large, the upper algorithm can approximate the EBP up to arbitrary order. Based on these developments, we can conclude that the EBP of the mean parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _d=\nu _d p_d$$\end{document} is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\mu }}_{d}({{\hat{\varvec{\theta }}}})=\nu _d{\hat{p}}_{d}({{\hat{\varvec{\theta }}}})$$\end{document} . Note that there may be applications where the researcher wants to avoid Monte Carlo integration, for instance when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_1$$\end{document} is large and heavy computational burden shall be avoided. For such cases, we state synthetic predictors as an alternative that can be calculated very efficiently. The synthetic predictor of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_d$$\end{document} is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{p}}_d^{syn}=\exp \{ \varvec{x}_{d}{\hat{\varvec{\beta }\}}}$$\end{document} . Likewise, the synthetic predictor of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _d=\nu _d p_d$$\end{document} is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\mu }}_d^{syn} = \nu _d \exp \{ \varvec{x}_{d}{\hat{\varvec{\beta }\}}}$$\end{document} . These predictors have acceptable accuracy when both measurement errors and random effects are negligible in the distribution of the target variable.

3.3. Mean Squared Error Estimation

This section presents a parametric bootstrap algorithm that estimates the mean squared error (MSE) of the EBP \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\mu }}_d={\hat{\mu }}_d({\hat{\varvec{\theta }}})$$\end{document} , which is characterized by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} MSE({\hat{\mu }}_d) = E \big [ ({\hat{\mu }}_d - \mu _d)^2 \big ], \quad d=1, \ldots , D. \end{aligned}$$\end{document}

Let B be the number of bootstrap replicates. The procedure is given as follows.

1. 1. Obtain estimates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\varvec{\theta }}}=({\hat{\varvec{\beta }}}, {\hat{\phi }})'$$\end{document} based on the observations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y_1, \varvec{x}_1), \ldots , (y_D, \varvec{x}_D)$$\end{document}

2. 2. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=1, \ldots , B$$\end{document} , do

   • (a)Generate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{d}^{(b)} \overset{i.i.d.}{\sim } N(0,1)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{u}_{d}^{(b)}\overset{i.i.d.}{\sim }N_p(\varvec{0},\varvec{\Sigma }_{d})$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1, \ldots , D$$\end{document}

   • (b)Calculate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _d^{(b)} = \nu _d p_d^{(b)}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_d^{(b)} = \exp \{\varvec{x}_d{\hat{\varvec{\beta }}} + \varvec{u}_d^{(b)\prime }{\hat{\varvec{\beta }}}+{\hat{\phi }}v_d^{(b)}\}$$\end{document} , and draw \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_d^{(b)} \sim \text {Poiss}(\mu _d^{(b)})$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1, \ldots , D$$\end{document}

   • (c)Obtain estimates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\varvec{\theta }}}^{(b)}=({\hat{\varvec{\beta }}}^{(b)}, {\hat{\phi }}^{(b)})$$\end{document} based on the observations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y_1^{(b)}, \varvec{x}_1), \ldots , (y_D^{(b)}, \varvec{x}_D)$$\end{document}

   • (d)Calculate the EBP \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\mu }}_d^{(b)} = {\hat{\mu }}_d({\hat{\varvec{\theta }}}^{(b)})$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1, \ldots , D$$\end{document}

3. 3. Output: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$mse^*({\hat{\mu }}_d) = \frac{1}{B} \sum _{b=1}^B ({\hat{\mu }}_d^{(b)} - \mu _d^{(b)})^2$$\end{document}

The MSE estimator based on this naive parametric bootstrap approach has generally \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(D^{-1})$$\end{document} bias, which is not always acceptable in practice. More refinement can be achieved by double bootstrapping, as proposed Hall and Maiti (Reference Hall and Maiti2006). We have not implemented the double bootstrap approach because it is computationally very intensive and has a great impact on simulation times. Nevertheless, in some real data setups with small D and high and accurate computing capacity, the double bootstrap approach is recommendable.

4.1. Method of Moments

A natural set of equations for the MM approach is

(3) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} 0= & {} f_k(\varvec{\theta })=M_k(\varvec{\theta })-{\hat{M}}_k=\sum _{d=1}^DE_\theta [y_{d}]x_{dk} -\sum _{d=1}^Dy_{d}x_{dk},\quad k=1,\ldots p, \end{aligned}$$\end{document}

(4) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} 0= & {} f_{p+1}(\varvec{\theta })=M_{p+1}(\varvec{\theta })-{\hat{M}}_{p+1}=\sum _{d=1}^DE_\theta [y_{d}^2]-\sum _{d=1}^Dy_{d}^2, \end{aligned}$$\end{document}

where the model-based expectations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{\theta }[y_d]$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{\theta }[y_d^2]$$\end{document} are given in Appendix A. For solving the system of nonlinear MM equations, we use a Newton–Raphson algorithm. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=0,1,2, \ldots $$\end{document} denote the index of iterations. Further, let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{f}(\varvec{\theta }^{(i)})$$\end{document} be the set of MM equations based on the model parameter vector in the i-th iteration. Likewise, let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{H}(\varvec{\theta }^{(i)})$$\end{document} denote the corresponding Jacobian matrix. The Newton–Raphson algorithm is performed as follows:

1. 1. Set the initial values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${i}=0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }^{(0)}=(\varvec{\beta }^{(0)},\phi ^{(0)})$$\end{document} .

2. 2. Repeat the following steps till convergence

   • (a)Update \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }^{({i})}$$\end{document} by using the equation

     \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{\theta }^{({i}+1)}=\varvec{\theta }^{({i})}-\varvec{H}^{-1}(\varvec{\theta }^{({i})})\varvec{f}(\varvec{\theta }^{({i})}), \end{aligned}$$\end{document}

   • (b)Update the iteration index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${i}\leftarrow {i}+1$$\end{document} .

3. 3. Output: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\varvec{\theta }}} = \varvec{\theta }^{({i}+1)}$$\end{document} .

Appendix A demonstrates how to compute the components of the Jacobian matrix. Note that the MM estimator is consistent in probability, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\varvec{\theta }}} \overset{P}{\rightarrow } \varvec{\theta }$$\end{document} as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D \rightarrow \infty $$\end{document} . We proof this result in Appendix B based on the developments of Jiang (Reference Jiang1998).

4.2. Variance Estimation and Significance Testing

Hereafter, we show how to estimate the variance of the model parameter estimates obtained from the MM method. This step is important with respect to practice when significance tests for covariates that are measured with error shall be conducted. Let us define

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{M}(\varvec{\theta })=\underset{1\le k \le p+1}{\hbox {col}}(M_k(\varvec{\theta })), \quad {{\hat{\varvec{M}}}}=\underset{1\le k \le p+1}{\hbox {col}}({\hat{M}}_k), \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_k$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{M}}_k$$\end{document} are defined as in Sect. 4.1. The asymptotic variance of the MM estimators can be approximated by a Taylor expansion of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{M}(\varvec{\theta })$$\end{document} around \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} . This is to say,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\hat{\varvec{M}}}=\varvec{M}({{\hat{\varvec{\theta }}}})\approx \varvec{M}(\varvec{\theta })+\varvec{H}(\varvec{\theta })({{\hat{\varvec{\theta }}}}-\varvec{\theta }),\quad {{\hat{\varvec{\theta }}}}-\varvec{\theta }\approx \varvec{H}^{-1}(\varvec{\theta })({{\hat{\varvec{M}}}}-\varvec{M}(\varvec{\theta })). \end{aligned}$$\end{document}

Under some regularity conditions, it holds that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text{ var }({{\hat{\varvec{\theta }}}})=E[({{\hat{\varvec{\theta }}}}-\varvec{\theta })({{\hat{\varvec{\theta }}}}-\varvec{\theta })^\prime ]\approx \varvec{H}^{-1}(\varvec{\theta })\text{ var }({{\hat{\varvec{M}}}})\varvec{H}^{-1}(\varvec{\theta }). \end{aligned}$$\end{document}

An estimator of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{ var }({{\hat{\varvec{\theta }}}})$$\end{document} is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \widehat{\text{ var }}({{\hat{\varvec{\theta }}}})=\varvec{H}^{-1}({{\hat{\varvec{\theta }}}})\widehat{\text{ var }}({{\hat{\varvec{M}}}})\varvec{H}^{-1}({{\hat{\varvec{\theta }}}}), \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\hat{\varvec{\theta }}}}=\varvec{\theta }^{(r+1)}$$\end{document} is taken from the output of the MM algorithm and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{\text{ var }}({{\hat{\varvec{M}}}})$$\end{document} is an estimator of the covariance matrix of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\hat{\varvec{M}}}}$$\end{document} . An algorithm to estimate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{ var }({{\hat{\varvec{\theta }}}})$$\end{document} is

1. 1. Fit the model to the sample and calculate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\hat{\varvec{\theta }}}}$$\end{document} .

2. 2. Generate bootstrap samples \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{y_{d}^{(b)}: d=1,\ldots ,D\}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=1,\ldots ,B$$\end{document} , from the fitted model.

3. 3. Calculate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\varvec{M}}}^{(b)}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=1,\ldots ,B$$\end{document} , and

   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\overline{\varvec{M}}}=\frac{1}{B}\sum _{b=1}^B{\hat{\varvec{M}}}^{(b)},\quad {\widehat{var}}_B({{\hat{\varvec{M}}}})=\frac{1}{B}\sum _{b=1}^B({\hat{\varvec{M}}}^{(b)}-{\overline{\varvec{M}}})({\hat{\varvec{M}}}^{(b)}-{\overline{\varvec{M}}})^\prime . \end{aligned}$$\end{document}

4. 4. Output: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{\text{ var }}_A({{\hat{\varvec{\theta }}}})=\varvec{H}^{-1}({{\hat{\varvec{\theta }}}})\widehat{\text{ var }}_B({{\hat{\varvec{M}}}})\varvec{H}^{-1}({{\hat{\varvec{\theta }}}})$$\end{document} .

Alternatively, steps 3 and 4 can be replaced by

1. 1. Fit the model to the bootstrap samples and calculate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\varvec{\theta }}}^{(b)}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=1,\ldots ,B$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{\varvec{\theta }}}=\frac{1}{B}\sum _{b=1}^B{\hat{\varvec{\theta }}}^{(b)}$$\end{document} .

2. 2. Output: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{\text{ var }}_B({{\hat{\varvec{\theta }}}})=\frac{1}{B}\sum _{b=1}^B ({\hat{\varvec{\theta }}}^{(b)}-{\overline{\varvec{\theta }}})({\hat{\varvec{\theta }}}^{(b)}-{\overline{\varvec{\theta }}})^\prime $$\end{document} .

From this procedure, for a given significance level \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \in (0, 1)$$\end{document} and the corresponding t-value, a model parameter estimate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\theta }}_k$$\end{document} is considered to be significantly different from zero if

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\hat{\theta }}_k \notin \left( {\hat{\theta }}_k - t_{D-p, 1-\frac{\alpha }{2}} \sqrt{\widehat{\text{ var }}({\hat{\theta }}_k)},\,\, {\hat{\theta }}_k + t_{D-p, 1-\frac{\alpha }{2}} \sqrt{\widehat{\text{ var }}({\hat{\theta }}_k)}\right) . \end{aligned}$$\end{document}

5.1. Set Up

A Monte Carlo simulation with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I=500$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1, \ldots , I$$\end{document} , iterations is conducted. We generate a population of D domains, where D varies over scenarios. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1, \ldots , D$$\end{document} , we define

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mu _d= \nu _d p_d, \quad y_d \sim \text {Poiss}(\mu _d), \quad p_d = \exp \{\beta _0+\varvec{x}_{1,d} \varvec{\beta }_1 + \varvec{x}_{2,d} \varvec{\beta }_2 + \varvec{u}_{1,d}' \varvec{\beta }_1 + \varvec{u}_{2,d}' \varvec{\beta }_2 + \phi v_d\}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _d = 300$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _0 = -4$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }_1 = (0.5, 0.5)'$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }_2 = - \varvec{\beta }_1$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi = 0.3$$\end{document} , as well as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_d \sim N(0,1)$$\end{document} . Accordingly, we have an intercept and four covariates that are measured with error. Note that the random effects are generated in every iteration individually. The unbiased covariate predictors are drawn from uniform distributions according to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{dj} \sim U(1.0,1.4)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1, \ldots , 4$$\end{document} , and held fixed over all Monte Carlo iterations. The covariate measurement errors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{u}_d = (\varvec{u}_{1,d}', \varvec{u}_{2,d}')'$$\end{document} are drawn in each iteration individually according to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{u}_d \sim N_{4}({\varvec{0}}, \varvec{\Sigma }_d), \quad \varvec{\Sigma }_d = \begin{pmatrix} \sigma _{1,d}^2 &{} \quad \sigma _{12,d} &{} \quad \sigma _{13,d} &{} \quad \sigma _{14,d}\\ \sigma _{21,d} &{} \quad \sigma _{2,d}^2 &{} \quad \sigma _{23,d} &{} \quad \sigma _{24,d}\\ \sigma _{31,d} &{} \quad \sigma _{32,d} &{} \quad \sigma _{3,d}^2 &{} \quad \sigma _{34,d}\\ \sigma _{41,d} &{} \quad \sigma _{42,d} &{} \quad \sigma _{34,d} &{} \quad \sigma _{4,d}^2\\ \end{pmatrix} , \quad d=1, \ldots , D, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{j,d}^2 \sim U(0.05, 0.15)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{jk,d} = \rho _{jk} \sigma _{j,d}^2 \sigma _{k,d}^2$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{jk} =0.5 $$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=2,3,4$$\end{document} , as well as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{jk}=-0.3$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=2,3$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=3,4$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j \ne k$$\end{document} . Overall, we consider four simulation scenarios arising from the four different values for D. The scenarios are defined as follows.

Table 1

Overview of simulation scenarios

The simulation is conducted with the statistic software package R. The objective is to estimate the mean parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _d$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1, \ldots , D$$\end{document} . We compare the following fitting methods:

• • PM: maximum likelihood Laplace method for the Poisson mixed model, with the R function glmer.

• • MEPM: methods of moments for the measurement error Poisson mixed model, with own programming.

5.2. Results of Model Parameter Estimation

This section studies the behavior of the MM algorithm in Sect. 4.1. We consider MSE and bias as performance measures. For parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta \in \{\beta _0, \varvec{\beta }_1, \varvec{\beta }_2, \phi \}$$\end{document} , they are given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} MSE(\theta ) = \frac{1}{I} \sum _{i=1}^I ({\hat{\theta }}^{(i)} - \theta )^2, \quad BIAS(\theta ) = \frac{1}{I} \sum _{i=1}^I ({\hat{\theta }}^{(i)} - \theta ). \end{aligned}$$\end{document}

Note that since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }_1 = (\beta _1, \beta _2)'$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _1 = \beta _2$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }_2 = (\beta _3, \beta _4)'$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _3=\beta _4$$\end{document} , we average the performance measures for each of for the sake of a compact presentation. Table 2 presents the MSE results, and Table 3 presents the bias results. We start with the regression parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }= (\beta _0, \varvec{\beta }_1', \varvec{\beta }_2')'$$\end{document} . It can be seen that the MEPM model outperforms the PM model in every scenario and for all elements of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} . The MSE of the estimates obtained from the MM approach are always smaller than those obtained from the maximum likelihood Laplace method under the PM model. The additional information obtained from the anticipation of the covariate measurement errors allows for efficiency gains in these settings. The effect is particular evident when looking at the intercept \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _0$$\end{document} . Here, the standard PM model cannot clearly identify the contribution of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _0$$\end{document} to the functional description of the target variable, which causes excessive variation in the respective estimates. The MEPM model, on the other hand, shows good performance on that regard. For the bias of the regression parameter estimates, we have mixed results. This is due to the fact that measurement error distributions are symmetric around zero. There is no clear tendency which of the two approaches obtains better results. However, when looking at the overall efficiency (in terms of the MSE), the MEPM model is better.

Table 2

Mean squared error of model parameter estimation

Table 3

Bias of model parameter estimation

We now consider the standard deviation parameter estimates. Here, the PM model obtains more efficient results. By looking at the bias, we see that the MEPM approach underestimates the true random effect standard deviation. This is likely due to the additional anticipation of the covariate measurement errors. The method seems to attribute some of the random effect variance to the uncertainty stemming from the auxiliary data, which causes the underestimation tendency and the overall MSE inefficiency. However, we will see in the next section that the predictions obtained under the MEPM model are superior to those under the PM model regardless of this effect.

5.3. Results of Mean Parameter Prediction

This section studies the behavior of the EBP under the area-level MEPM model presented in Sect. 3.2 against the EBP derived by Boubeta et al. (Reference Boubeta, Lombardía and Morales2016) under the area-level Poisson mixed model. We consider relative mean squared error (RMSE), relative root mean squared error (RRMSE), absolute bias (ABIAS), as well as relative absolute bias (RABIAS) as performance measures. They are given as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&RMSE_{d} = \left( \frac{1}{I} \sum _{i=1}^I (\mu _d^{(i)} - {\hat{\mu }}_d^{(i)})^2 \right) ^{1/2}, \quad RRMSE_d = \frac{RMSE_d}{{\bar{\mu }}_d}, \quad {\bar{\mu }}_d = \frac{1}{I} \sum _{i=1}^I \mu _d^{(i)},\\&ABIAS_d = \frac{1}{I} \sum _{i=1}^I |\mu _d^{(i)} - {\hat{\mu }}_d^{(i)}|, \quad RABIAS_d = \frac{ABIAS_d}{{\bar{\mu }}_d}. \end{aligned}$$\end{document}

We further define the proportion of efficient predictions (PoEP)

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} PoEP = \frac{1}{D} \sum _{d=1}^D \mathbbm {1}_{(RMSE_d < RMSE_d^*)} \cdot 100\%, \end{aligned}$$\end{document}

as well as the subsequent aggregated measures

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&RMSE = \frac{1}{D} \sum _{d=1}^D RMSE_d, \quad ABIAS = \frac{1}{D} \sum _{d=1}^D ABIAS_d,\\&RRMSE = \frac{1}{D} \sum _{d=1}^R RRMSE_d \cdot 100, \quad RABIAS = \frac{1}{D} \sum _{d=1}^D RABIAS_d \cdot 100\%, \end{aligned}$$\end{document}

to allow for a more compact presentation. Here, the term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbbm {1}_{(RMSE_d < RMSE_d^*)}$$\end{document} yields the value 1 when the respective method has achieved a smaller RMSE than the alternative, and 0 else. Table 4 presents the results of mean parameter prediction.

Table 4

Performance of mean parameter prediction

We see that the EBP under the MEPM model outperforms the predictor under the PM model in all scenarios and for all considered performance measures. Further, by looking at the measure PoEP, we see that the MEPM predictions are more efficient than the PM predictions in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$100\%$$\end{document} of the domains for Scenario 1 and 2, and in over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90\%$$\end{document} for Scenario 3 and 4. The anticipation of the covariate measurement errors leads to a significant improvement of the domain count predictions, which is in line with the theoretical developments of Sect. 3. This is further visualized in Fig. 1. It contains boxplots of the measure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$RRMSE_d$$\end{document} over all domains per scenario. The results of the PM model are displayed in blue, while the results of the MEPM model are depicted in red. It becomes evident that both the boxes and the median of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$RRMSE_d$$\end{document} -values are visibly lower under the MEPM model than under the PM model. The strongest efficiency gains are obtained under the scenario with a smaller number of domains. This is important with respect to the empirical application in Sect. 6, where the number of domains is small as well.

Figure. 1

RRMSE of mean parameter prediction

5.4. Results of Mean Squared Error Estimation

This section studies the performance of the parametric bootstrap algorithm for MSE estimation in Sect. 3.3. For this, we use the Monte Carlo MSE obtained measured in the simulation experiment as approximation to the real MSE. That is to say, we define

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} MSE_{d} = \frac{1}{I} \sum _{i=1}^I (\mu _d^{(i)} - {\hat{\mu }}_d^{(i)})^2, \quad mse_d^* = \frac{1}{I} \sum _{i=1}^{I} mse_d^{*(i)}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$mse_d^{*(i)}$$\end{document} is an MSE estimate obtained for the d-th domain in the i-th Monte Carlo iteration. Further, we consider absolute bias and relative absolute bias as performance measures, which are given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} ABIAS_d = \frac{1}{I} \sum _{i=1}^{I} |mse_d^{*(i)} - MSE_d|, \quad RABIAS_d = \frac{ABIAS_d}{MSE_d}\cdot 100\%. \end{aligned}$$\end{document}

As in the previous two simulation experiments, we use aggregated performance measures

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} mse^* = \frac{1}{D} \sum _{d=1}^D mse_d^*, \quad ABIAS = \frac{1}{D} \sum _{d=1}^D ABIAS_d, \quad RABIAS = \frac{1}{D} \sum _{d=1}^D RABIAS_d \end{aligned}$$\end{document}

for a compact presentation. We choose \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=300$$\end{document} for the number of bootstrap replicates used for MSE estimation. Table 5 presents the corresponding simulation results and further visualized in Fig. 2. It can be seen that the parametric bootstrap yields plausible results. There is a slight tendency for underestimation for scenarios where the number of domains is small. However, this tendency decreases when increasing the number of domains, although at a rather slow rate. This is because increasing D also increases the number of unknown MSE values with to be estimated under new domain-specific measurement error distributions. Furthermore, recall the comment on the MSE estimator’s bias at the end of Sect. 3.3. Depending on the application, researcher may consider using double bootstrap approaches in order to improve MSE estimates. However, with a relative absolute bias of roughly \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10\%$$\end{document} , the MSE estimates are relatively accurate given the fact that a GLMM with measurement errors is used for prediction.

Table 5

Simulation results for MSE estimation

Figure 2 studies the convergence behavior of the parametric bootstrap for MSE estimation. It shows the distribution of the difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$mse_d^*(B) - mse_d^*(1000)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=1, \ldots ,D$$\end{document} , when the parametric bootstrap estimator is calculated based on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B \in \{1, \ldots , 1000\}$$\end{document} replicates.

Figure. 2

Convergence behavior of MSE estimators per domain

The light gray band displays the inner \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90\%$$\end{document} of differences over all domains, while the gray band marks the inner \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$70\%$$\end{document} . The mean is plotted by the black line. The blue dashed lines mark \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10\%$$\end{document} differences, and the red dashed line displays \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\%$$\end{document} difference. We see that the mean stabilizes at about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=250$$\end{document} iterations. Accordingly, with the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=300$$\end{document} we employed in the simulation, the results are generally stable and allow for an evaluation of the overall MSE level in practice. However, in case regional analysis shall be conducted where MSE estimates are compared between domains, a higher number of iterations is required. This is due to the additional uncertainty resulting from the covariate measurement errors. The inner \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$70\%$$\end{document} of MSE estimates are within the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10\%$$\end{document} difference range for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=400$$\end{document} , and the inner \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90\%$$\end{document} only for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=600$$\end{document} or higher. That is to say, if regional analysis between domains in terms of MSE estimates shall be conducted, we recommend at least \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=600$$\end{document} .