MAP estimation

To compute the MAP estimates of the parameters in a regression model, the command MAP.estimation can be used. To apply this command, the data has to be in a specific form and the inverse covariance matrix of the Gaussian prior needs to be chosen. The analysis below is for the combined data set Nurse. The estimates in the separate hospitals can be obtained with the same commands, but with the local data sets instead.

The command inv.prior.cov creates a diagonal inverse covariance matrix for the prior distribution of the correct dimension. Based on the characteristics of the covariates (continuous or categorical) in M and the number of nuisance parameters, the number of model parameters is computed (the number of regression parameters for a categorical variable equals the number of levels minus one and a linear model (family="gaussian") has one nuisance parameter, the variance of the measurement error). The argument lambda=0.01 means that all elements on the diagonal of $\Lambda $ equal 0.01. The arguments of MAP.estimation are the outcome variable Nurses$stress, the covariate data in M, the type of the model (family="gaussian") and the matrix Lambda. The inference results are stored in the list fit. A summary of it can be found with summary(fit).

When applying the BFI approach, the analyses are performed in every hospital and the results in fit are sent to the central server. There, the results from the different hospitals are combined. This is explained below.

BFI for homogeneous populations

Suppose that all hospitals have sent their output to the central server. For ease of notation, we assume these outputs are stored in fit1, fit2, …, fit25. From every output the relevant elements need to be selected and combined. If the number of hospitals is high, it is easier to work with a for-loop. With the following code, all relevant elements of 25 local centers are created and combined by the main function bfi():

Here Lambda is the inverse covariance matrix of the prior for the (fictive) combined data. The command bfi combines the estimates from the different hospitals into the BFI estimates. The outcome fitbfi_homo is a list with the BFI estimates $\widehat {\mathbf {\theta }}_{\mathrm { {BFI}}}$ and $\widehat {\mathbf {A}}_{\mathrm { {BFI}}}$ . The command summary(fitbfi_homo) gives the BFI estimates (and more information).

BFI for heterogeneous populations

Different types of heterogeneity have been discussed in Section 3. Below we will explain how to do the analyses in R.

Heterogeneity of population characteristics

Heterogeneity across population characteristics in the centers implies that the value of the parameter $\mathbf {\theta }_2$ differs across centers. Because the bfi-command estimates the parameter $\mathbf {\theta }_1$ (and its curvature matrix $\widehat {A}_1$ ), and these estimates are not affected by $\mathbf {\theta }_2$ , the R-code explained in the previous subsection can still be applied.

Heterogeneity across outcome means

Suppose the intercepts differ across hospitals. To take this variation into account we allow a hospital specific intercept in the regression model. Instead of one general intercept there are $L=25$ intercepts; an increase of $L-1$ parameters. The dimension of the inverse covariance matrix Lambda for the fictive combined data set changes as well. For a diagonal matrix with 0.01 at the diagonal, this matrix can be obtained by

These commands replace the two corresponding commands above. The argument L=25 has to be added to indicate the number of centers, and thus the number of location specific intercepts. This matrix should be appended to the list priors instead. The MAP estimates can be obtained with the command bfi, but it needs to be made explicit that the hospitals may have different intercepts:

For this stratified analysis extra arguments have been added: stratified=TRUE and strat_par=1. The first argument indicates that the full model stratifies with respect to the different hospitals. The default is stratified=FALSE. If strat_par=1 there is stratification with respect to the intercept and if strat_par=2 this is the case for the variance of the measurement error in a linear regression model. A summary of the results can be obtained by |summary(fitbfi_hetero)|.This gives a list with estimates, starting with the hospital specific intercepts.

Heterogeneity due to clustering

An example of a cluster variable is hospital size. For all nurses in a hospital this covariate is constant and, as a consequence, the effect of hospital size on stress cannot be estimated within a hospital. However, the model for the (fictive) combined data could include this covariate if there is variation across the hospitals (Section 3.5). Here we explain how to do the analyses in R. In practice, every local hospital sends its size (small, medium, large) to the central server. Then, a vector with all sizes is defined in R. Suppose this vector is named Hsize. After fitting the local models (like explained before), the estimated model for the (fictive) combined data can be obtained with:

The commands return a list with categorical specific intercepts and the estimates of the remaining parameters.

Appendix II.A Homogeneity across centers

In this appendix we derive expressions of the BFI estimators under the assumption that the variables $(\mathbf {X}_{\ell i},Y_{\ell i}), i=1,\ldots ,n_\ell , \ell =1,\ldots , L$ are independent and identically distributed. In equations (1) and (2) in Section 2 we have seen that the log posterior densities for the (fictive) combined data set $\mathbf {D}$ and for the subset $\mathbf {D}_\ell $ equal

(A.1) $$ \begin{align} \log \big\{ p(\boldsymbol{\theta}|\mathbf{D}) \big\} &=\; \log \big\{p(\boldsymbol{\theta}_1)\big\} \;+\; \log \big\{p(\boldsymbol{\theta}_2)\big\} \;+\; \sum_{\ell=1}^L \sum_{i=1}^{n_\ell}\log \big\{p(y_{\ell i}| \mathbf{x}_{\ell i},\boldsymbol{\theta}_1)\big\}\nonumber\\&\quad \;+\; \sum_{\ell=1}^L \sum_{i=1}^{n_\ell}\log \big\{ p(\mathbf{x}_{\ell i}|\boldsymbol{\theta}_2)\big\} \;-\; \log \big\{p(\mathbf{D})\big\} \end{align} $$

(A.2) $$ \begin{align} \log \big\{p_\ell(\boldsymbol{\theta}|\mathbf{D}_\ell)\big\} &=\; \log \big\{p_\ell(\boldsymbol{\theta}_1)\big\} \;+\; \log \big\{p_\ell(\boldsymbol{\theta}_2)\big\} \;+\; \sum_{i=1}^{n_\ell}\log \big\{p(y_{\ell i}| \mathbf{x}_{\ell i},\boldsymbol{\theta}_1)\big\}\nonumber\\&\quad \;+\; \sum_{i=1}^{n_\ell}\log \big\{p(\mathbf{x}_{\ell i}|\boldsymbol{\theta}_2)\big\} \;-\; \log \big\{p_\ell(\mathbf{D}_\ell)\big\}. \end{align} $$

By reordering the terms in equation (A.2), it follows that for every center $\ell $

$$ \begin{align*} &\sum_{i=1}^{n_\ell}\log \big\{p(y_{\ell i}| \mathbf{x}_{\ell i},\boldsymbol{\theta}_1)\big\} \;+\; \sum_{i=1}^{n_\ell}\log \big\{p(\mathbf{x}_{\ell i}|\boldsymbol{\theta}_2)\big\} \\&\qquad=\; \log \big\{p_\ell(\boldsymbol{\theta}|\mathbf{ D}_\ell)\big\} - \log \big\{p_\ell(\boldsymbol{\theta}_1)\big\} \;-\; \log \big\{p_\ell(\boldsymbol{\theta}_2)\big\} + \log \big\{p_\ell(\mathbf{D}_\ell)\big\}. \end{align*} $$

Next, summing over all centers yields

(A.3) $$ \begin{align} &{\sum_{\ell=1}^L\sum_{i=1}^{n_\ell}\log \big\{p(y_{\ell i}| \mathbf{x}_{\ell i},\boldsymbol{\theta}_1)\big\} \;+\; \sum_{\ell=1}^L\sum_{i=1}^{n_\ell}\log \big\{p(\mathbf{x}_{\ell i}|\boldsymbol{\theta}_2)\big\}} \nonumber\\ &\quad =\; \sum_{\ell=1}^L\log \big\{p_\ell(\boldsymbol{\theta}|\mathbf{D}_\ell)\big\} - \log \Big\{\prod_{\ell=1}^L p_\ell(\boldsymbol{\theta}_1)\Big\} \;-\; \log \Big\{\prod_{\ell=1}^L p_\ell(\boldsymbol{\theta}_2)\Big\} \;+\; \log \Big\{\prod_{\ell=1}^L p_\ell(\mathbf{D}_\ell)\Big\}. \end{align} $$

By inserting the right hand side of equation (A.3) into the right hand side of equation (A.1) this yields

(A.4) $$ \begin{align} \log \big\{p(\boldsymbol{\theta}|\mathbf{D})\big\} &\;=\; \sum_{\ell=1}^L \log \big\{p_\ell(\boldsymbol{\theta}|\mathbf{D}_\ell)\big\}+\log\Bigg\{\frac{p(\boldsymbol{\theta}_1)}{\prod_{\ell=1}^L p_\ell(\boldsymbol{\theta}_1)}\Bigg\}\nonumber\\&\quad+\log\Bigg\{\frac{p(\boldsymbol{\theta}_2)}{\prod_{\ell=1}^L p_\ell(\boldsymbol{\theta}_2)}\Bigg\}-\log\Bigg\{\frac{ p(\mathbf{D})}{\prod_{\ell=1}^L p_\ell(\mathbf{D}_\ell)}\Bigg\}. \end{align} $$

We expressed the log posterior densities of the combined data, $\log \{p(\mathbf {\theta }|\mathbf {D})\}$ , in terms of the log posterior densities of the local data sets, $\log \{p_\ell (\mathbf {\theta }|\mathbf {D}_\ell )\}$ . However, the final aim is to express the MAP estimator $\widehat {\mathbf {\theta }}$ based on the (fictive) combined data set $\mathbf {D}$ in terms of the MAP estimators based on the local data sets $\mathbf {D}_\ell $ . This will be done next. We approximate the log posterior densities for the data set $\mathbf {D}_\ell $ by a Taylor expansion up to the quadratic order in $\mathbf {\theta }$ around the MAP estimator $\widehat {\mathbf {\theta }}_\ell $ :

$$ \begin{align*} \log \big\{p_\ell\big(\boldsymbol{\theta}|\mathbf{D}_\ell\big)\big\} \;=\; \log \big\{p_\ell\big(\widehat{\boldsymbol{\theta}}_\ell|\mathbf{D}_\ell\big)\big\} -\tfrac{1}{2}\big(\boldsymbol{\theta}-\widehat{\boldsymbol{\theta}}_\ell\big)^t \widehat{\mathbf{A}}_\ell(\boldsymbol{\theta}-\widehat{\boldsymbol{\theta}}_\ell\big) + O_p\big(\|\widehat{\boldsymbol{\theta}}_\ell-\boldsymbol{\theta}\|^3\big), \end{align*} $$

with $\widehat {\mathbf {A}}_\ell $ equal to minus the second derivative of $\log \{p_\ell (\mathbf {\theta }|\mathbf {D}_\ell )\}$ with respect to $\mathbf {\theta }$ , evaluated at $\widehat {\mathbf {\theta }}_\ell $ . The linear term in the Taylor expansion is equal to zero and therefore missing in the expansion; the MAP estimator maximizes the log posterior density and the first derivative evaluated at $\widehat {\mathbf {\theta }}_\ell $ is therefore equal to zero. The last term in the Taylor expansion is equal to $O_p(\|\widehat {\mathbf {\theta }}_\ell -\mathbf {\theta }\|^3)=\|\widehat {\mathbf {\theta }}_\ell -\mathbf {\theta }\|^3 O_p(1)$ , where $O_p(1)$ represents a term that is bounded in probability for the sample size going to infinity.Reference van der Vaart²⁸ For $\mathbf {\theta }$ in a small neighborhood of $\widehat {\mathbf {\theta }}_\ell $ , the term $\|\widehat {\mathbf {\theta }}_\ell -\mathbf {\theta }\|^3$ will be close to zero (in probability), and the remainder term $O_p(\|\widehat {\mathbf {\theta }}_\ell -\mathbf {\theta }\|^3)$ is small compared to the other terms in the Taylor expansion which are of an order of at most $\|\widehat {\mathbf {\theta }}_\ell -\mathbf {\theta }\|^2$ .

Since the log posterior density in equation (A.2) is decomposed in terms that depend on either $\mathbf {\theta }_1$ or $\mathbf {\theta }_2$ , but never on both, the matrices $\widehat {\mathbf {A}}_\ell , \ell =1,\ldots ,L$ are diagonal block matrices:

$$ \begin{align*} \widehat{\mathbf{A}}_\ell=\left(\!\begin{array}{cc} \widehat{\mathbf{A}}_{1,\ell} & \mathbf{0}\\ \mathbf{0} & \widehat{\mathbf{A}}_{2,\ell}\end{array}\! \right), \nonumber \end{align*} $$

with the blocks $\widehat {\mathbf {A}}_{1,\ell }$ and $\widehat {\mathbf {A}}_{2,\ell }$ equal to minus the second derivative matrices for $\mathbf {\theta }_1$ and $\mathbf {\theta }_2$ , respectively, and the log posterior densities can be approximated by

$$ \begin{align*} \log \big\{p_\ell\big(\boldsymbol{\theta}|\mathbf{D}_\ell\big)\big\} &=\; \log \big\{p_\ell\big(\widehat{\boldsymbol{\theta}}_\ell|\mathbf{D}_\ell\big)\big\} \;-\; \tfrac{1}{2}\big(\boldsymbol{\theta}_1-\widehat{\boldsymbol{\theta}}_{1,\ell}\big)^t\widehat{\mathbf{A}}_{1,\ell}\big(\boldsymbol{\theta}_1-\widehat{\boldsymbol{\theta}}_{1,\ell}\big)\\&\quad \;-\; \tfrac{1}{2}\big(\boldsymbol{\theta}_2-\widehat{\boldsymbol{\theta}}_{2,\ell}\big)^t\widehat{\mathbf{A}}_{2,\ell}\big(\boldsymbol{\theta}_2-\widehat{\boldsymbol{\theta}}_{2,\ell}\big) + O_p\big(\|\widehat{\boldsymbol{\theta}}_\ell-\boldsymbol{\theta}\|^3\big). \end{align*} $$

By substituting this expansion for $\log \{ p_\ell (\mathbf {\theta }|\mathbf {D}_\ell )\},\ell =1,\ldots ,L$ , into the relation (A.4), we obtain:

(A.5) $$ \begin{align} \log \big\{p(\boldsymbol{\theta}|\mathbf{D})\big\} &=\; -\tfrac{1}{2}\sum_{\ell=1}^L \big(\boldsymbol{\theta}_1-\widehat{\boldsymbol{\theta}}_{1,\ell}\big)^t\widehat{\mathbf{ A}}_{1,\ell}\big(\boldsymbol{\theta}_1-\widehat{\boldsymbol{\theta}}_{1,\ell}\big) -\; \tfrac{1}{2}\sum_{\ell=1}^L\big(\boldsymbol{\theta}_2-\widehat{\boldsymbol{\theta}}_{2,\ell}\big)^t\widehat{\mathbf{ A}}_{2,\ell}\big(\boldsymbol{\theta}_2-\widehat{\boldsymbol{\theta}}_{2,\ell}\big) \nonumber\\ &\quad \;+\;\log\Bigg\{\frac{p(\boldsymbol{\theta}_1)}{\prod_{\ell=1}^L p_\ell(\boldsymbol{\theta}_1)}\Bigg\} \;+\;\log\Bigg\{\frac{p(\boldsymbol{\theta}_2)}{\prod_{\ell=1}^L p_\ell(\boldsymbol{\theta}_2)}\Bigg\} \;+\; B \;+\; O_p\Big(\sum_{\ell=1}^L\|\widehat{\boldsymbol{\theta}}_\ell-\boldsymbol{\theta}\|^3\Big). \end{align} $$

where B is a term that depends on the data, but is not a function of $\mathbf {\theta }=(\mathbf {\theta }_1,\mathbf {\theta }_2)$ . Now choose the prior densities $\mathbf {\theta }_1 \rightarrow p(\mathbf {\theta }_1)$ and $ \mathbf {\theta }_2 \rightarrow p(\mathbf {\theta }_2)$ in the combined data set and $\mathbf {\theta }_1 \rightarrow p_\ell (\mathbf {\theta }_1)$ and $\mathbf {\theta }_2 \rightarrow p_\ell (\mathbf {\theta }_2)$ in center $\ell $ to be Gaussian with mean zero and inverse covariance matrices $\boldsymbol {\Lambda }_1$ and $\boldsymbol {\Lambda }_2$ in the combined data set, and $\boldsymbol {\Lambda }_{1,\ell }$ and $\boldsymbol {\Lambda }_{2,\ell }$ in center $\ell $ : e.g., $p(\mathbf {\theta }_1) = (\det \boldsymbol {\Lambda }_1 / (2\pi )^d)^{1/2}\exp (-\tfrac {1}{2}\mathbf {\theta }_1^t\boldsymbol {\Lambda }_1\mathbf {\theta }_1)$ . Inserting the expressions of the densities into (A.5) yields

(A.6) $$ \begin{align} \log \big\{p(\boldsymbol{\theta}|\mathbf{D})\big\} &=\; -\tfrac{1}{2}\sum_{\ell=1}^L \big(\boldsymbol{\theta}_1-\widehat{\boldsymbol{\theta}}_{1,\ell}\big)^t\widehat{\mathbf{ A}}_{1,\ell}\big(\boldsymbol{\theta}_1-\widehat{\boldsymbol{\theta}}_{1,\ell}\big) \;-\; \tfrac{1}{2} \sum_{\ell=1}^L \big(\boldsymbol{\theta}_2-\widehat{\boldsymbol{\theta}}_{2,\ell}\big)^t\widehat{\mathbf{ A}}_{2,\ell}\big(\boldsymbol{\theta}_2-\widehat{\boldsymbol{\theta}}_{2,\ell}\big)\nonumber\\ &\quad-\tfrac{1}{2} \boldsymbol{\theta}_1^t \Big(\boldsymbol{\Lambda}_1-\sum_{\ell=1}^L \boldsymbol{\Lambda}_{1,\ell}\Big)\boldsymbol{\theta}_1 \;-\tfrac{1}{2} \boldsymbol{\theta}_2^t \Big(\boldsymbol{\Lambda}_2-\sum_{\ell=1}^L \boldsymbol{\Lambda}_{2,\ell}\Big)\boldsymbol{\theta}_2 \;+\; B'~\;+\; O_p\Big(\sum_{\ell=1}^L\|\widehat{\boldsymbol{\theta}}_\ell-\boldsymbol{\theta}\|^3\Big) \nonumber\\ &=: \; \Omega_{{\mathrm{{BFI}}}}(\theta) \;+\; O_p\Big(\sum_{\ell=1}^L\|\widehat{\boldsymbol{\theta}}_\ell-\boldsymbol{\theta}\|^3\Big), \end{align} $$

for $B'$ a term that depends on the data, but not of $\mathbf {\theta }_1$ and $\mathbf {\theta }_2$ . The function $\mathbf {\theta } \rightarrow \Omega _{{\mathrm { {BFI}}}}(\mathbf {\theta })$ in equation (A.6) is quadratic function of $\mathbf {\theta }_1$ and $\mathbf {\theta }_2$ . Maximizing $\mathbf {\theta } \rightarrow \Omega _{{\mathrm { {BFI}}}}(\mathbf {\theta })$ with respect to $\mathbf {\theta }=(\mathbf {\theta }_1,\mathbf {\theta }_2)$ yields the BFI estimators

$$ \begin{align*} \hspace{-5mm} & \widehat{\boldsymbol{\theta}}_{1,{\mathrm{{BFI}}}} \;:=\; \big(\widehat{\mathbf{A}}_{1,{\mathrm{{BFI}}}}\big)^{-1}\sum_{\ell=1}^L \widehat{\mathbf{ A}}_{1,\ell}\widehat{\boldsymbol{\theta}}_{1,\ell}, \quad \widehat{\mathbf{A}}_{1,{\mathrm{{BFI}}}} := \sum_{\ell=1}^L \widehat{\mathbf{A}}_{1,\ell}+\boldsymbol{\Lambda}_1-\sum_{\ell=1}^L \boldsymbol{\Lambda}_{1,\ell}, \\ \hspace{-5mm} & \widehat{\boldsymbol{\theta}}_{2,{\mathrm{{BFI}}}} \;:=\; \big(\widehat{\mathbf{A}}_{2,{\mathrm{{BFI}}}}\big)^{-1}\sum_{\ell=1}^L \widehat{\mathbf{ A}}_{2,\ell}\widehat{\boldsymbol{\theta}}_{2,\ell}, \quad\widehat{\mathbf{A}}_{2,{\mathrm{{BFI}}}} \;:=\; \sum_{\ell=1}^L \widehat{\mathbf{A}}_{2,\ell}+\boldsymbol{\Lambda}_2-\sum_{\ell=1}^L \boldsymbol{\Lambda}_{2,\ell}, \end{align*} $$

where $\widehat {\mathbf {A}}_{1,{\mathrm { {BFI}}}}$ and $\widehat {\mathbf {A}}_{2,{\mathrm { {BFI}}}}$ equal minus the second derivative of $\Omega _{{\mathrm { {BFI}}}}$ with respect to $\mathbf {\theta }_1$ and $\mathbf {\theta }_2$ . In Appendix III.B the asymptotic distribution of the BFI estimators are derived.

Appendix II.B Heterogeneity across centers, center-specific parameter

Suppose that the vector of regression parameters can be subdivided into two parts. One part is equal across the centers and the other part may vary. A special case is the situation in which the intercepts vary. In the calculations of the BFI estimator, we assume that the covariates are statistically independent between the individuals within and across the centers. We, moreover, assume that the outcome variables given the covariates and the center are independent.

Suppose that the vector $\mathbf {\theta }$ can be decomposed as $\mathbf {\theta }=(\mathbf {\theta }_1,\mathbf {\theta }_2)=(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b},\mathbf {\theta }_2)$ , where, as before, $\mathbf {\theta }_2$ is the vector of parameters that specifies the distribution of the covariates. The parameter $\mathbf {\theta }_1=(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b})$ is decomposed so that $\mathbf {\theta }_{1a}$ is the vector of (regression) parameters that is assumed to be equal across the centers, and $\mathbf {\theta }_{1b}$ is the vector of (regression) parameters that may vary. In this appendix we consider the situation in which every center has its own specific vector of parameters: $\mathbf {\theta }_{1b,1}, \ldots , \mathbf {\theta }_{1b,L}$ for the L centers. The vector of parameters in the combined data set is equal to $\mathbf {\theta }=(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b,1},\ldots ,\mathbf {\theta }_{1b,L},\mathbf {\theta }_{2})$ , where $\mathbf {\theta }_{1b,\ell }$ is the parameter vector in center $\ell $ . If only the intercepts vary across the centers, $\mathbf {\theta }_{1b,\ell }$ is one-dimensional, but for now we allow $\mathbf {\theta }_{1b,\ell }$ to be a vector.

For simplicity of notation we assume that $\mathbf {\theta }_{1a}, \mathbf {\theta }_{1b}$ and $\mathbf {\theta }_2$ are independent: $p(\mathbf {\theta })=p(\mathbf {\theta }_{1a})p(\mathbf {\theta }_{2})\prod _{\ell =1}^L p(\mathbf {\theta }_{1b,\ell })$ for the combined data set, and in center $\ell $ : $p_\ell (\mathbf {\theta }_{1a},\mathbf {\theta }_{1b,\ell },\mathbf {\theta }_{2})=p_\ell (\mathbf {\theta }_{1a})p_\ell (\mathbf {\theta }_{1b,\ell })p_\ell (\mathbf {\theta }_{2})$ . As before, the log posterior densities can be written as

$$ \begin{align*} {\log \big\{p(\boldsymbol{\theta}|\mathbf{D})\big\}}&\;=\; \log \big\{p(\boldsymbol{\theta}_{1a})\big\} + \sum_{\ell=1}^L \log \big\{p(\boldsymbol{\theta}_{1b,\ell})\big\} + \log \big\{p(\boldsymbol{\theta}_{2})\big\} + \sum_{\ell=1}^L\sum_{i=1}^{n_\ell} \log \big\{p(y_{\ell i}|\boldsymbol{\theta}_{1a},\boldsymbol{\theta}_{1b,\ell},\mathbf{ x}_{\ell i})\big\}\\&\quad + \sum_{\ell=1}^L\sum_{i=1}^{n_\ell} \log \big\{p(\mathbf{x}_{\ell i}|\boldsymbol{\theta}_{2})\big\} -\log \big\{p(\mathbf{D})\big\}, \end{align*} $$

and

$$ \begin{align*} {\log \big\{p_\ell(\boldsymbol{\theta}_{1a},\boldsymbol{\theta}_{1b,\ell},\boldsymbol{\theta}_2|\mathbf{D}_\ell)\big\}} &=\; \log \big\{p_\ell(\boldsymbol{\theta}_{1a})\big\} + \log \big\{p_\ell(\boldsymbol{\theta}_{1b,\ell})\big\} + \log \big\{p_\ell(\boldsymbol{\theta}_{2})\big\}\\&\quad + \sum_{i=1}^{n_\ell} \log \big\{p(y_{\ell i}|\boldsymbol{\theta}_{1a},\boldsymbol{\theta}_{1b,\ell},\mathbf{x}_{\ell i})\big\} + \sum_{i=1}^{n_\ell} \log \big\{p(\mathbf{x}_{\ell i}|\boldsymbol{\theta}_{2})\big\} -\log \big\{p_\ell(\mathbf{D}_\ell)\big\}. \end{align*} $$

Previously, and in the formulas above, we see that the log posterior density is decomposed into terms that depend on $\mathbf {\theta }_1$ or $\mathbf {\theta }_2$ , but never on both. That means that the BFI estimator for $\mathbf {\theta }_1$ is not affected by the estimator $\mathbf {\theta }_2$ and vice versa. Therefore, in this setting, the BFI estimator $\widehat {\mathbf {\theta }}_2$ can be expressed in terms of the local MAP estimators $\widehat {\mathbf {\theta }}_{2,\ell }$ and $\widehat {\mathbf {A}}_{2,\ell }$ as in (4). In the remainder of the derivation we focus on $\mathbf {\theta }_1$ only and leave out the terms with $\mathbf {\theta }_2$ from the expressions.

Like in the homogeneous setting, the log posterior density in the full data set can be written in terms of the local log posterior densities:

(A.7) $$ \begin{align} \hspace{-3mm} \log \big\{p(\boldsymbol{\theta}|\mathbf{D})\big\} \;=\; \sum_{\ell=1}^L \log \big\{p_\ell(\boldsymbol{\theta}_{1a},\boldsymbol{\theta}_{1b,\ell}|\mathbf{D}_\ell)\big\}+\log\Bigg\{\frac{p(\boldsymbol{\theta}_{1a})}{\prod_{\ell=1}^L p_\ell(\boldsymbol{\theta}_{1a})}\Bigg\} \;+ \; \log\Bigg\{\frac{\prod_{\ell=1}^L p(\boldsymbol{\theta}_{1b,\ell})}{\prod_{\ell=1}^L p_\ell(\boldsymbol{\theta}_{1b,\ell})}\Bigg\}+B \end{align} $$

with B a term that depends on the data and on $\mathbf {\theta }_2$ , but is not a function of $\mathbf {\theta }_1$ .

Let $\widehat {\mathbf {\theta }}_{1a,\ell }$ and $\widehat {\mathbf {\theta }}_{1b,\ell }$ be the MAP estimators of $\mathbf {\theta }_{1a}$ and $\mathbf {\theta }_{1b,\ell }$ based on the data set $\mathbf {D}_\ell $ . Moreover, let $\widehat {\mathbf {A}}_{1a,\ell }$ and $\widehat {\mathbf {A}}_{1b,\ell }$ be minus the second derivative of $\log \{p_\ell (\mathbf {\theta }|\mathbf {D}_\ell )\}$ with respect to $\mathbf {\theta }_{1a}$ and $\mathbf {\theta }_{1b,\ell }$ respectively, and let $\widehat {\mathbf {A}}_{1ab,\ell }$ be minus the second derivative with respect to both $\mathbf {\theta }_{1a}$ and $\mathbf {\theta }_{1b,\ell }$ all evaluated at $\widehat {\mathbf {\theta }}_{1,\ell } =(\widehat {\mathbf {\theta }}_{1a,\ell },\widehat {\mathbf {\theta }}_{1b,\ell })$ .

The Taylor expansions up to the quadratic term of $\log \{p_\ell (\mathbf {\theta }_{1a},\mathbf {\theta }_{1b,\ell }|\mathbf {D}_\ell )\}$ around $\widehat {\mathbf {\theta }}_{1,\ell }$ is given by

$$ \begin{align*} \log \big\{p_\ell\big(\boldsymbol{\theta}_{1a},\boldsymbol{\theta}_{1b,\ell}|\mathbf{D}_\ell\big)\big\} &=\; \log \big\{p_\ell\big(\widehat{\boldsymbol{\theta}}_{1a,\ell},\widehat{\boldsymbol{\theta}}_{1b,\ell}|\mathbf{D}_\ell\big)\big\} -\tfrac{1}{2}\big(\boldsymbol{\theta}_{1a}-\widehat{\boldsymbol{\theta}}_{1a, \ell}\big)^t \widehat{\mathbf{A}}_{1a, \ell}\big(\boldsymbol{\theta}_{1a}-\widehat{\boldsymbol{\theta}}_{1a, \ell}\big) \nonumber\\[5pt] &\qquad -\;\tfrac{1}{2} \big(\boldsymbol{\theta}_{1b,\ell}-\widehat{\boldsymbol{\theta}}_{1b,\ell}\big)^t \widehat{\mathbf{ A}}_{1b,\ell}\big(\boldsymbol{\theta}_{1b,\ell}-\widehat{\boldsymbol{\theta}}_{1b,\ell}\big)\\&\qquad-\big(\boldsymbol{\theta}_{1a}-\widehat{\boldsymbol{\theta}}_{1a, \ell}\big)^t \widehat{\mathbf{ A}}_{1ab,\ell}\big(\boldsymbol{\theta}_{1b,\ell}-\widehat{\boldsymbol{\theta}}_{1b,\ell}\big) \;+\; O_p\big(\|\widehat{\boldsymbol{\theta}}_{1,\ell} - \boldsymbol{\theta}_1\|^3\big). \end{align*} $$

Next, we insert the Taylor expressions into (A.7). For the combined data $\mathbf {D}$ we assume a Gaussian prior with mean zero and inverse covariance matrix $\boldsymbol {\Lambda }_{1a}$ for $\mathbf {\theta }_{1a}$ , and a zero mean Gaussian prior with inverse covariance matrix $\boldsymbol {\Lambda }_{1b\ell }$ for $\mathbf {\theta }_{1b,\ell }, \ell =1,\ldots ,L$ . For center $\ell $ , also zero mean Gaussian priors are chosen, but with inverse covariance matrices $\boldsymbol {\Lambda }_{1a,\ell }$ and $\boldsymbol {\Lambda }_{1b,\ell }$ . The dimension of $\boldsymbol {\Lambda }_{1b\ell }$ and $\boldsymbol {\Lambda }_{1b,\ell }$ depends on the number of parameters that may vary across the centers. If only the intercepts vary, the matrices are scalars. After inserting these densities in expression (A.7) as well, we obtain

(A.8) $$ \begin{align} \log \big\{p(\boldsymbol{\theta}|\mathbf{D})\big\} &=\; -\tfrac{1}{2}\sum_{\ell=1}^L \big(\boldsymbol{\theta}_{1a}-\widehat{\boldsymbol{\theta}}_{1a,\ell}\big)^t \widehat{\mathbf{A}}_{1a,\ell}\big(\boldsymbol{\theta}_{1a}-\widehat{\boldsymbol{\theta}}_{1a,\ell}\big) -\tfrac{1}{2}\boldsymbol{\theta}_{1a}^t\Big(\boldsymbol{\Lambda}_{1a}-\sum_{\ell=1}^L\boldsymbol{\Lambda}_{1a,\ell}\Big)\boldsymbol{\theta}_{1a}\nonumber\\&\quad \;-\tfrac{1}{2}\sum_{\ell=1}^L \big(\boldsymbol{\theta}_{1b,\ell}-\widehat{\boldsymbol{\theta}}_{1b,\ell}\big)^t \widehat{\mathbf{A}}_{1b,\ell}\big(\boldsymbol{\theta}_{1b,\ell}-\widehat{\boldsymbol{\theta}}_{1b,\ell}\big) \nonumber\\ &\quad -\tfrac{1}{2} \sum_{\ell=1}^L \boldsymbol{\theta}_{1b,\ell}^t\big(\boldsymbol{\Lambda}_{1b\ell}-\boldsymbol{\Lambda}_{1b,\ell}\big)\boldsymbol{\theta}_{1b,\ell} \;-\; \sum_{\ell=1}^L \big(\boldsymbol{\theta}_{1a}-\widehat{\boldsymbol{\theta}}_{1a,\ell}\big)^t \widehat{\mathbf{A}}_{1ab,\ell}\big(\boldsymbol{\theta}_{1b,\ell}-\widehat{\boldsymbol{\theta}}_{1b,\ell}\big)\nonumber\\&\quad \;+\; B' \;+\; O_p\Big(\sum_{\ell=1}^L \|\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\|^3\Big) \nonumber\\ &=:\; \Omega_{{\mathrm{{BFI}}}}(\boldsymbol{\theta}_1) + O_p\Big(\sum_{\ell=1}^L \|\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\|^3\Big) \end{align} $$

with $B'$ representing a term that does not depend on $\mathbf {\theta }_1$ . The function $\mathbf {\theta }_1\rightarrow \Omega _{{\mathrm { {BFI}}}}(\mathbf {\theta }_1)$ is a quadratic function of $\mathbf {\theta }_1$ . Maximization of this function with respect to $\mathbf {\theta }_1$ by setting its derivative equal to zero, yields the BFI estimators

(A.9) $$ \begin{align} &\widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}}:=\; \Big( \widehat{\mathbf{A}}_{1a,{\mathrm{{BFI}}}} -\sum_{\ell=1}^L \widehat{\mathbf{A}}_{1ab,\ell}(\widehat{\mathbf{A}}_{1b,\ell,{\mathrm{{BFI}}}})^{-1}(\widehat{\mathbf{A}}_{1ab,\ell,{\mathrm{{BFI}}}})^t \Big)^{\!-1}\times \nonumber\\ & \quad \sum_{\ell=1}^L \Big[ \Big( \widehat{\mathbf{A}}_{1a,\ell} - \widehat{\mathbf{A}}_{1ab,\ell}(\widehat{\mathbf{A}}_{1b,\ell,{\mathrm{{BFI}}}})^{-1}(\widehat{\mathbf{A}}_{1ab,\ell})^t\Big)\widehat{\boldsymbol{\theta}}_{1a,\ell} + \widehat{\mathbf{A}}_{1ab,\ell}\Big(\mathbf{1} - (\widehat{\mathbf{A}}_{1b,\ell,{\mathrm{{BFI}}}})^{-1}\widehat{\mathbf{A}}_{1b,\ell}\Big)\widehat{\boldsymbol{\theta}}_{1b,\ell} \Big] \end{align} $$

with $\mathbf {1}$ the unit matrix and the matrices $\widehat {\mathbf {A}}_{1a,{\mathrm { {BFI}}}}$ and $\widehat {\mathbf {A}}_{1b,\ell ,{\mathrm { {BFI}}}}$ as given in (A.11) below and

(A.10) $$ \begin{align} \widehat{\boldsymbol{\theta}}_{1b,\ell,{\mathrm{{BFI}}}} \; :=\;\big(\widehat{\mathbf{A}}_{1b,\ell,{\mathrm{{BFI}}}}\big)^{-1}\Big[\widehat{\mathbf{A}}_{1b,\ell}\widehat{\boldsymbol{\theta}}_{1b,\ell} + \big(\widehat{\mathbf{A}}_{1ab,\ell}\big)^t\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}}\big)\Big] \end{align} $$

with

(A.11) $$ \begin{align} \widehat{\mathbf{A}}_{1a,{\mathrm{{BFI}}}}:=\sum_{\ell=1}^L \widehat{\mathbf{A}}_{1a,\ell}+\boldsymbol{\Lambda}_{1a}-\sum_{\ell=1}^L\boldsymbol{\Lambda}_{1a,\ell}, \quad \widehat{\mathbf{A}}_{1b,\ell, {\mathrm{{BFI}}}}:=\widehat{\mathbf{A}}_{1b,\ell}+\boldsymbol{\Lambda}_{1b\ell}-\boldsymbol{\Lambda}_{1b,\ell}, \quad \widehat{\mathbf{A}}_{1ab,\ell, {\mathrm{{BFI}}}}:=\widehat{\mathbf{A}}_{1ab,\ell}, \end{align} $$

where $\widehat {\mathbf {A}}_{1a,{\mathrm { {BFI}}}}, \widehat {\mathbf {A}}_{1b,\ell , {\mathrm { {BFI}}}}$ and $\widehat {\mathbf {A}}_{1ab,\ell , {\mathrm { {BFI}}}}$ equal minus the second derivatives of $\Omega _{{\mathrm { {BFI}}}}$ with respect to $\mathbf {\theta }_{1a}, \mathbf {\theta }_{1b}$ and the mix $\mathbf {\theta }_{1a}$ and $\mathbf {\theta }_{1b}$ .

Appendix II.C Heterogeneity across centers due to clustering

In this appendix we consider the situation in which the centers are clustered by, for example, due to location or type (academic / non-academic medical center). Another example is the covariate hospital size in the nurse-data set, where the clusters are: small, medium, large. Within a hospital/center, all nurses are in the same cluster and have the same covariate value, so that the covariate can not be included in the local regression model as it is collinear with the intercept.

In the calculations of the BFI estimators, we assume that the covariates are independent between the individuals within and across the centers. We, moreover, assume that the outcome variables given the covariates and the cluster level are independent. Suppose that the vector of model parameters in center $\ell $ is equal to $\mathbf {\theta }_\ell = (\mathbf {\theta }_{1a},\theta _{1b,\ell },\mathbf {\theta }_2)$ , where, as before, $\mathbf {\theta }_2$ is the parameter vector that specifies the distribution of the covariates. The parameter $\mathbf {\theta }_{1a}$ is the vector of regression parameters which are assumed to be equal in all centers, but excluding the intercept which may vary across the centers. The parameter $\theta _{1b,\ell } \in \{\theta _{1b1},\ldots ,\theta _{1bK}\}$ for $\ell =1,\ldots ,L$ and with $K\leq L$ is the intercept of the model in center $\ell $ . So, $\theta _{1b,\ell }$ (with a comma in the subscript) is the parameter in center $\ell $ , whereas $\theta _{1bk}$ (without a comma in the subscript) is the parameter for the kth category of the center-specific covariate. If $K=L$ , $\theta _{1b,\ell }\neq \theta _{1b,\ell '}$ for $\ell \neq \ell '$ and we are in the situation of Appendix II.B, where every center has its own specific intercept value. If $K<L$ , there are centers $\ell $ and $\ell '$ with $\ell \neq \ell '$ with $\theta _{1b,\ell }= \theta _{1b,\ell '}$ . In the example, the covariate “hospital size” has three levels (small, medium, large). That means that $K=3$ and the three parameters represent the three intercepts for the three classes of centers. The parameter vector in the (fictive) combined data set $\mathbf {D}$ is defined as $\mathbf {\theta }=(\mathbf {\theta }_{1a}, \theta _{1b1}, \ldots , \theta _{1bK}, \mathbf {\theta }_2)$ .

For simplicity of notation we assume (again) that $\mathbf {\theta }_{1a}, \theta _{1b}$ and $\mathbf {\theta }_2$ are independent: $p(\mathbf {\theta })=p(\mathbf {\theta }_{1a})p(\mathbf {\theta }_{2})\prod _{k=1}^K p(\theta _{1bk})$ for the combined data set, and locally $p_\ell (\mathbf {\theta }_{1a},\theta _{1b,\ell },\mathbf {\theta }_{2})=p_\ell (\mathbf {\theta }_{1a})p_\ell (\theta _{1b,\ell })p_\ell (\mathbf {\theta }_{2})$ in data subset $\ell $ . For the combined data $\mathbf {D}$ we assume a Gaussian prior with mean zero and inverse covariance matrices $\boldsymbol {\Lambda }_{1a}$ for $\mathbf {\theta }_{1a}$ , and a zero mean Gaussian prior with inverse variance $\boldsymbol {\Lambda }_{1bk}$ for $\mathbf {\theta }_{1bk}, k=1,\ldots ,K$ . Also for center $\ell $ zero mean Gaussian priors are chosen, but with inverse covariance matrix $\boldsymbol {\Lambda }_{1a,\ell }$ and inverse variance $\boldsymbol {\Lambda }_{1b,\ell }$ . Similar notation and calculations as in Appendix II.B lead to the equation below, instead of the equation (A.8):

$$ \begin{align*} \log \big\{p(\boldsymbol{\theta}|\mathbf{D})\big\} &=\; -\tfrac{1}{2}\sum_{\ell=1}^L \big(\boldsymbol{\theta}_{1a}-\widehat{\boldsymbol{\theta}}_{1a,\ell}\big)^t \widehat{\mathbf{A}}_{1a,\ell}\big(\boldsymbol{\theta}_{1a}-\widehat{\boldsymbol{\theta}}_{1a,\ell}\big) \;-\;\tfrac{1}{2}\boldsymbol{\theta}_{1a}^t\Big(\boldsymbol{\Lambda}_{1a}-\sum_{\ell=1}^L\boldsymbol{\Lambda}_{1a,\ell}\Big)\boldsymbol{\theta}_{1a} \\&\quad\;-\;\tfrac{1}{2}\sum_{\ell=1}^L \big(\theta_{1b}-\widehat{\theta}_{1b,\ell}\big) \widehat{A}_{1b,\ell}\big(\theta_{1b,\ell}-\widehat{\theta}_{1b,\ell}\big) \nonumber \\ & \hspace{5mm} \;-\;\tfrac{1}{2} \sum_{k=1}^K \theta_{1bk}\Lambda_{1bk}\theta_{1bk} \;+\;\tfrac{1}{2} \sum_{\ell=1}^L \theta_{1b,\ell}\Lambda_{1b,\ell}\theta_{1b,\ell}\\&\quad \;-\; \sum_{\ell=1}^L \big(\boldsymbol{\theta}_{1a}-\widehat{\boldsymbol{\theta}}_{1a,\ell}\big)^t \widehat{\mathbf{A}}_{1ab,\ell}\big(\boldsymbol{\theta}_{1b,\ell}-\widehat{\boldsymbol{\theta}}_{1b,\ell}\big) \;+\; B' \;+\; O_p\Big(\sum_{\ell=1}^L \|\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\|^3\Big)\\ &=:\; \Omega_{{\mathrm{{BFI}}}}(\boldsymbol{\theta}_1) \;+\; O_p\Big(\sum_{\ell=1}^L \|\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\|^3\Big)\\[-18pt] \end{align*} $$

with $B'$ representing a term that does not depend on $\mathbf {\theta }_1$ . Let $z_\ell $ denote the category of center $\ell $ for the center-specific covariate. So $z_\ell \in \{1,\ldots ,K\}$ . Differentiating $\Omega _{{\mathrm { {BFI}}}}(\mathbf {\theta }_1)$ with respect to $\mathbf {\theta }_{1a}$ and $\theta _{1b}$ and setting the derivatives equal to zero, yields the BFI estimators:

$$ \begin{align*} \widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}} &:=\; \Big(\widehat{\mathbf{A}}_{1a,{\mathrm{{BFI}}}}-\sum_{k=1}^K\widehat{\mathbf{A}}_{1abk,{\mathrm{{BFI}}}}\big(\widehat{A}_{1bk,{\mathrm{{BFI}}}}\big)^{-1}\;(\widehat{\mathbf{ A}}_{1abk,{\mathrm{{BFI}}}})^t\Big)^{-1} \\ &\quad\times \Bigg( \sum_{\ell=1}^L \widehat{\mathbf{A}}_{1a,\ell}\widehat{\boldsymbol{\theta}}_{1a,\ell} +\sum_{\ell=1}^L \widehat{\mathbf{A}}_{1ab,\ell}\widehat{\theta}_{1b,\ell}\\&\qquad-\sum_{k=1}^K \widehat{\mathbf{A}}_{1abk,{\mathrm{{BFI}}}}\big(\widehat{A}_{1bk,{\mathrm{{BFI}}}}\big)^{-1}\Big[\sum_{\ell=1: z_\ell=k}^L\widehat{A}_{1b,\ell}\widehat{\theta}_{1b,\ell} +\sum_{\ell=1:z_\ell=k}^L(\widehat{\mathbf{A}}_{1ab,\ell})^t\widehat{\boldsymbol{\theta}}_{1a,\ell} \Big]\Bigg)\\[-18pt] \end{align*} $$

with $\widehat {\mathbf {A}}_{1a,{\mathrm { {BFI}}}}$ , $\widehat {A}_{1bk,{\mathrm { {BFI}}}}$ and $\widehat {\mathbf {A}}_{1abk,{\mathrm { {BFI}}}}$ as given in (A.12) below and the estimator $\widehat {\theta }_{1bk,{\mathrm { {BFI}}}}, k=1,\ldots ,K$ is given by

$$ \begin{align*} &\widehat{\theta}_{1bk,{\mathrm{{BFI}}}} \;:=\;\\ &\quad \big(\widehat{A}_{1bk,{\mathrm{{BFI}}}}\big)^{-1}\Big[\sum_{\ell=1: z_\ell=k}^L\widehat{A}_{1b,\ell}\widehat{\theta}_{1b,\ell} \;+\;\sum_{\ell=1: z_\ell=k}^L(\widehat{\mathbf{A}}_{1ab,\ell})^t\widehat{\boldsymbol{\theta}}_{1a,\ell} \;-\;(\widehat{\mathbf{A}}_{1abk,{\mathrm{{BFI}}}})^t\widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}} \Big]. \end{align*} $$

with minus the second derivatives of $\Omega _{{\mathrm { {BFI}}}}$ equal to

(A.12) $$ \begin{align} \widehat{\mathbf{A}}_{1a,{\mathrm{{BFI}}}} &:=\; \sum_{\ell=1}^L \widehat{\mathbf{A}}_{1a,\ell}\;+\;\boldsymbol{\Lambda}_{1a}\;-\;\sum_{\ell=1}^L\boldsymbol{\Lambda}_{1a,\ell}, \nonumber\\ \widehat{A}_{1bk,{\mathrm{{BFI}}}} &:=\;\sum_{\ell=1: z_\ell=k}^L \widehat{A}_{1b,\ell}\;+\;\Lambda_{1bk}\;-\;\sum_{\ell=1: z_\ell=k}^L\Lambda_{1b,\ell}, \\ \widehat{\mathbf{A}}_{1abk,{\mathrm{{BFI}}}} &:=\; \sum_{\ell=1: z_\ell=k}^L\widehat{\mathbf{A}}_{1ab,\ell}, \nonumber \end{align} $$

and 0 for the remaining terms.

Appendix III.A Asymptotic distribution of the MAP estimator based on the combined data

In this section we study the asymptotic distribution of the MAP estimator $\widehat {\mathbf {\theta }}_1=(\widehat {\mathbf {\theta }}_{1a},\widehat {\mathbf {\theta }}_{1b})$ for $\mathbf {\theta }_1=(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b})$ in the combined data set. The asymptotic distribution for the MAP estimator for $\mathbf {\theta }_2$ can be derived similarly.

From literature (Bernstein–Von Mises TheoremReference van der Vaart²⁸) it is known that the MAP estimator is asymptotically Gaussian:

$$ \begin{align*} \sqrt{n} \begin{pmatrix} \begin{pmatrix} \widehat{\boldsymbol{\theta}}_{1a}\\ \widehat{\boldsymbol{\theta}}_{1b} \end{pmatrix} - \begin{pmatrix} \boldsymbol{\theta}_{1a}\\ \boldsymbol{\theta}_{1b} \end{pmatrix} \end{pmatrix} \;\leadsto \; {\mathcal N}\big(\mathbf{0}, J_{1}^{-1} \big), \end{align*} $$

where “ $\leadsto $ ” means convergence in distribution for the sample size to infinity. Further, $\mathbf {0}$ is a vector of zeroes, and $J_{1}^{-1}$ is the inverse Fisher information matrix for the combined populations from all centers. The Fisher information matrix and its inverse have the form

$$ \begin{align*} J_{1} &=\; \begin{pmatrix}J_{1a} & J_{1ab} \\(J_{1ab})^t & J_{1b}\end{pmatrix}\\[10pt](J_{1})^{-1} &=\; \begin{pmatrix}H & ~~~ -H J_{1ab} (J_{1b})^{-1} \\-(J_{1b})^{-1} (J_{1ab})^t H & ~~~~~~~~ (J_{1b})^{-1} + (J_{1b})^{-1} (J_{1ab})^t H J_{1ab} (J_{1b})^{-1}\end{pmatrix}\end{align*} $$

with $H=(J_{1a}-J_{1ab} (J_{1b})^{-1} (J_{1ab})^t)^{-1}$ and where $J_{1a}$ and $J_{1b}$ equal the expectation of the second derivatives of $-\log \{p(\mathbf {D}|\mathbf {\theta })\}$ with respect to $\mathbf {\theta }_{1a}$ and $\mathbf {\theta }_{1b}$ , evaluated at the true value of $\mathbf {\theta }_1$ . The matrix $J_{1ab}$ equals the expectation of the derivative of $-\log \{p(\mathbf {D}|\mathbf {\theta })\}$ with respect to $\mathbf {\theta }_{1a}$ and $\mathbf {\theta }_{1b}$ , evaluated at the true value of $\mathbf {\theta }_1$ .

In the following we rewrite the four submatrices of $(J_1)^{-1}$ in terms of the Fisher information matrices in the local centers. In the next appendices it will be proven that the asymptotic covariance matrices of the BFI-estimators equal these expressions and the BFI estimators are therefore asymptotically efficient.

The data from the different centers are independent. Therefore, the Fisher information matrix can be written as a weighted sum of the Fisher information matrices in the different centers. By the law of large numbers

$$ \begin{align*} J_1 = \lim_{n_1,\ldots,n_L\rightarrow \infty} \frac{\partial^2}{\partial \boldsymbol{\theta}_1^2} \Big(- \; \frac{1}{n}\log \big\{p(\mathbf{D} | \boldsymbol{\theta}_1)\big\}\Big) &= \lim_{n_1,\ldots,n_L\rightarrow \infty} \frac{\partial^2}{\partial \boldsymbol{\theta}_1^2} \Big(- \; \sum_{\ell=1}^L \frac{1}{n}\log \big\{p(\mathbf{D}_\ell | \boldsymbol{\theta}_1)\big\}\Big)\\[5pt] &= \lim_{n_1,\ldots,n_L\rightarrow \infty} \; \sum_{\ell=1}^L \frac{n_\ell}{n} \frac{\partial^2}{\partial \boldsymbol{\theta}_1^2} \Big(-\frac{1}{n_\ell}\log \big\{p(\mathbf{D}_\ell | \boldsymbol{\theta}_1)\big\}\Big)\\& = \sum_{\ell=1}^L w_\ell J_{1,\ell} \end{align*} $$

with $J_{1,\ell }$ the Fisher information matrix for $\mathbf {\theta }_1$ in center $\ell $ and $n_\ell /n \rightarrow w_\ell $ if $n_\ell , n \rightarrow \infty $ . So $J_1 = \sum _{\ell =1}^L w_\ell J_{1,\ell }$ .

In the homogeneous setting all parameters are included in $\mathbf {\theta }_{1a}$ (there is no parameter $\mathbf {\theta }_{1b}$ ). Then, for $I_{1,\ell }$ the Fisher information matrix for $\mathbf {\theta }_{1a}$ in center $\ell $ , it follows that $J_{1,\ell }=I_{1,\ell }=I_1, \ell =1,2,\ldots ,L$ and

(A.13) $$ \begin{align} J_1 \;=\; \sum_{\ell=1}^L w_\ell J_{1,\ell} \;=\; \sum_{\ell=1}^L w_\ell I_{1,\ell} \;=\; I_1, \quad (J_1)^{-1} \;=\; \Big(\sum_{\ell=1}^L w_\ell I_{1,\ell}\Big)^{-1} \;=\; (I_1)^{-1}. \end{align} $$

The heterogeneous setting is more complex. Suppose the parameter $\mathbf {\theta }_{1a}$ is assumed to be same across all centers, but $\mathbf {\theta }_{1b}=(\mathbf {\theta }_{1b,1},\ldots ,\mathbf {\theta }_{1b,L})$ is a vector with center-specific parameters (the index refers to the center). The log likelihood function for center $\ell $ is a function of $\mathbf {\theta }_{1b,\ell }$ , but not of $\mathbf {\theta }_{1b,k}$ with $k\neq \ell $ . Therefore, for the Fisher information matrix $J_{1,\ell }$ for $(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b})=(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b,1},\ldots ,\mathbf {\theta }_{1b,L})$ in center $\ell $ , the columns and rows that are related to $\mathbf {\theta }_{1b,k}, k\neq \ell $ contain zeroes only. The matrix $I_{1,\ell }$ is the Fisher information matrix for $(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b,\ell })$ in center $\ell $ (so not of $(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b})$ like $J_{1,\ell }$ ), with the blocks $I_{1a,\ell }, I_{1b,\ell }$ and $I_{1ab,\ell }$ , defined in a similar way as in $J_1$ . Since $\mathbf {\theta }_{1a}$ is the same across the centers, $J_{1a,\ell }=I_{1a,\ell }$ . However, $J_{1b,\ell }\neq I_{1b,\ell }$ , since $J_{1b,\ell }$ is the Fisher information matrix for $\mathbf {\theta }_{1b}=(\mathbf {\theta }_{1b,1}, \ldots , \mathbf {\theta }_{1b,L})$ in center $\ell $ , whereas $I_{1b,\ell }$ is the Fisher information matrix for $\mathbf {\theta }_{1b,\ell }$ in center $\ell $ ; the dimensions of the matrices are different.

Since the parameter $\mathbf {\theta }_{1b}$ is a vector with center-specific parameters, the matrix $J_{1b}$ has a block diagonal matrix with center-specific blocks. Because of this form, it follows that

$$ \begin{align*} J_{1ab} (J_{1b})^{-1} (J_{1ab})^t = \sum_{\ell=1}^L w_{\ell} I_{1ab,\ell} (I_{1b,\ell})^{-1} (I_{1ab,\ell})^t. \end{align*} $$

Since $J_{1a}=\sum _{\ell =1}^L w_\ell J_{1a,\ell }=\sum _{\ell =1}^L w_\ell I_{1a,\ell }$ (the parameter vector $\mathbf {\theta }_{1a}$ is shared across all centers),

$$ \begin{align*} J_{1a}-J_{1ab} (J_{1b})^{-1} (J_{1ab})^t = \sum_{\ell=1}^L w_\ell \big(I_{1a,\ell}-I_{1ab,\ell} (I_{1b,\ell})^{-1} (I_{1ab,\ell})^t \big) \end{align*} $$

and the asymptotic covariance matrix for $\widehat {\mathbf {\theta }}_{1a}$ is equal to

(A.14) $$ \begin{align} \Big(J_{1a}-J_{1ab} (J_{1b})^{-1} (J_{1ab})^t\Big)^{-1} = \Big(\sum_{\ell=1}^L w_\ell \big(I_{1a,\ell}-I_{1ab,\ell} (I_{1b,\ell})^{-1} (I_{1ab,\ell})^t \big)\Big)^{-1}. \end{align} $$

The asymptotic covariance matrix for the MAP estimator $\widehat {\mathbf {\theta }}_{1b}$ equals:

(A.15) $$ \begin{align} (J_{1b})^{-1} + (J_{1b})^{-1} (J_{1ab})^t \big(J_{1a}-J_{1ab} (J_{1b})^{-1} (J_{1ab})^t\big)^{-1} J_{1ab} (J_{1b})^{-1}. \end{align} $$

For parameter $\mathbf {\theta }_{1b,\ell }$ the asymptotic covariance matrix equals the $\ell $ th diagonal block of this matrix. The corresponding block of the matrix $J_{1b}$ equals $w_\ell I_{1b,\ell }$ . By the structure of $J_{1b}$ and the equation (A.14) it follows that $\ell ^{th}$ diagonal block of the matrix in (A.15) is given by

(A.16) $$ \begin{align} \big(w_\ell I_{1b,\ell}\big)^{-1}+\big(I_{1b,\ell}\big)^{-1} \big(I_{1ab,\ell}\big)^{t} \Big(\sum_{k=1}^L w_k \big(I_{1a,k}-I_{1ab,k}(I_{1b,k})^{-1}(I_{1ab,k}\big))^t\Big)^{-1} I_{1ab,\ell} \big(I_{1b,\ell}\big)^{-1}. \end{align} $$

Appendix III.B Asymptotic distribution in homogeneous setting

Asymptotic distribution BFI estimator

In Appendix II.A we have derived an expression for the BFI estimator in the homogeneous setting. In this setting, all parameters are included in the vector $\mathbf {\theta }_{1a}$ ; there is no vector $\mathbf {\theta }_{1b}$ . The BFI estimator $\widehat {\mathbf {\theta }}_{1,{\mathrm { {BFI}}}}$ is defined as

$$ \begin{align*} \widehat{\boldsymbol{\theta}}_{1,{\mathrm{{BFI}}}} = \big(\widehat{\mathbf{A}}_{1,{\mathrm{{BFI}}}}\big)^{-1}\sum_{\ell=1}^L \widehat{\mathbf{A}}_{1,\ell}\widehat{\boldsymbol{\theta}}_{1,\ell} \quad\mbox{with} \quad\widehat{\mathbf{A}}_{1,{\mathrm{{BFI}}}} = \sum_{\ell=1}^L \widehat{\mathbf{A}}_{1,\ell}+\boldsymbol{\Lambda}_1-\sum_{\ell=1}^L \boldsymbol{\Lambda}_{1,\ell}. \end{align*} $$

Below, we derive the asymptotic distribution of $\sqrt {n}\big (\widehat {\mathbf {\theta }}_{1,{\mathrm { {BFI}}}} - \mathbf {\theta }_1\big )$ :

$$ \begin{align*} \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1,{\mathrm{{BFI}}}} - \boldsymbol{\theta}_1\big) &= \sqrt{n} \Big\{ \Big(\sum_{\ell=1}^L \widehat{\mathbf{A}}_{1,\ell}+\boldsymbol{\Lambda}_1-\sum_{\ell=1}^L \boldsymbol{\Lambda}_{1,\ell}\Big)^{-1}\Big(\sum_{\ell=1}^L \widehat{\mathbf{A}}_{1,\ell}\widehat{\boldsymbol{\theta}}_{1,\ell}\Big)-\boldsymbol{\theta}_1\Big\} \\[4pt] &= \sqrt{n} \Big\{ \Big(\sum_{\ell=1}^L \frac{1}{n}\widehat{\mathbf{A}}_{1,\ell}+\frac{1}{n}\boldsymbol{\Lambda}_1-\sum_{\ell=1}^L \frac{1}{n}\boldsymbol{\Lambda}_{1,\ell}\Big)^{-1}\Big(\sum_{\ell=1}^L \frac{1}{n}\widehat{\mathbf{ A}}_{1,\ell}\widehat{\boldsymbol{\theta}}_{1,\ell}\Big)-\boldsymbol{\theta}_1\Big\}. \end{align*} $$

Because the term

$$ \begin{align*} \Big(\sum_{\ell=1}^L \frac{1}{n}\widehat{\mathbf{A}}_{1,\ell}+\frac{1}{n}\boldsymbol{\Lambda}_1-\sum_{\ell=1}^L \frac{1}{n}\boldsymbol{\Lambda}_{1,\ell}\Big)^{-1} = \Big(\sum_{\ell=1}^L \frac{1}{n}\widehat{\mathbf{ A}}_{1,\ell}\Big)^{-1} + O_p\Big(\frac{1}{n}\Big), \end{align*} $$

it holds that

$$ \begin{align*} \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1,{\mathrm{{BFI}}}} - \boldsymbol{\theta}_1\big) &= \sqrt{n} \Big\{ \Big(\sum_{\ell=1}^L \frac{1}{n}\widehat{\mathbf{A}}_{1,\ell}\Big)^{-1}\Big(\sum_{\ell=1}^L \frac{1}{n}\widehat{\mathbf{A}}_{1,\ell}\widehat{\boldsymbol{\theta}}_{1,\ell}\Big)-\boldsymbol{\theta}_1\Big\} + O_p\Big(\frac{1}{\sqrt{n}}\Big)\\[4pt] &= \Big(\sum_{\ell=1}^L \frac{1}{n}\widehat{\mathbf{A}}_{1,\ell}\Big)^{-1}\Big(\sum_{\ell=1}^L \frac{1}{n}\widehat{\mathbf{A}}_{1,\ell} \;\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\big)\Big) + O_p\Big(\frac{1}{\sqrt{n}}\Big)\\[4pt] &= \Big(\sum_{\ell=1}^L \frac{n_\ell}{n} \frac{1}{n_\ell}\widehat{\mathbf{A}}_{1,\ell}\Big)^{-1}\Big(\sum_{\ell=1}^L \frac{\sqrt{n_\ell}}{\sqrt{n}} \; \frac{1}{n_\ell}\widehat{\mathbf{A}}_{1,\ell} \;\sqrt{n_\ell}\big(\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\big)\Big) + O_p\Big(\frac{1}{\sqrt{n}}\Big). \end{align*} $$

Asymptotically, the MAP estimator and the maximum likelihood estimator are equivalent. It follows that the MAP estimator in center $\ell $ , $\widehat {\mathbf {\theta }}_{1,\ell }$ , is asymptotically normalReference van der Vaart^28, Reference Bijma, Jonker and van der Vaart²⁹: $\sqrt {n_\ell }(\widehat {\mathbf {\theta }}_{1,\ell }-\mathbf {\theta }_{1,\ell }) \leadsto {\mathcal N}(0,(I_{1,\ell })^{-1})$ for $I_{1,\ell }$ the Fisher information matrix in center $\ell $ . Remember that $\widehat {\mathbf { A}}_{1,\ell }$ is defined as the second derivative of $-\log \{p(\mathbf {D}_\ell |\mathbf {\theta }_1)\}$ evaluated at $\widehat {\mathbf {\theta }}_{1,\ell }$ . If this second derivative is sufficiently smooth near $\mathbf {\theta }_{1,\ell }$ , it follow by the law of large numbers, that $n_\ell ^{-1}\widehat {\mathbf {A}}_{1,\ell }$ converges in probability to $I_{1,\ell }$ .Reference van der Vaart^28, Reference Bijma, Jonker and van der Vaart²⁹ By Slutsky’s lemma, it follows that, for every center $\ell $

$$ \begin{align*} \frac{1}{n_\ell}\widehat{\mathbf{A}}_{1,\ell} \; \sqrt{n_\ell}\big(\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\big) \;=\; I_{1,\ell} \; \sqrt{n_\ell}\big(\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\big) + o_P(1) \;\leadsto\; {\mathcal N}\big(0, I_{1,\ell}\big). \end{align*} $$

Since the data across the L centers are assumed to be independent, it follows that

$$ \begin{align*} \sum_{\ell=1}^L \frac{\sqrt{n_\ell}}{\sqrt{n}} \; \frac{1}{n_\ell}\widehat{\mathbf{A}}_{1,\ell} \; \sqrt{n_\ell}\big(\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\big) \;=\; \sum_{\ell=1}^L \sqrt{w_\ell} \; I_{1,\ell} \; \sqrt{n_\ell}\big(\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\big) + o_P(1) \;\leadsto\; {\mathcal N}\Big(0,\sum_{\ell=1}^L w_\ell I_{1,\ell}\Big). \end{align*} $$

Further, the term

$$ \begin{align*} \Big(\sum_{\ell=1}^L \frac{n_\ell}{n} \frac{1}{n_\ell}\widehat{\mathbf{A}}_{1,\ell}\Big)^{-1} \;=\; \Big(\sum_{\ell=1}^L w_\ell I_{1,\ell}\Big)^{-1} \;+\; o_p(1). \end{align*} $$

Combining the results, yields

$$ \begin{align*} \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1,{\mathrm{{BFI}}}} - \boldsymbol{\theta}_1\big) \leadsto {\mathcal N}\Big(0,\Big(\sum_{\ell=1}^L w_\ell I_{1,\ell}\Big)^{-1}\Big). \end{align*} $$

The asymptotic covariance matrix equals $J_1^{-1}$ as defined in (A.13), which equals the asymptotic covariance matrix of the MAP estimator based on the combined data. The BFI estimator is asymptotically efficient; no information is lost if the data from the centers can not be combined. Under homogeneity the matrices $I_{1,\ell }=I_1=J_1, \ell =1,\ldots ,L$ , and because $\sum _{\ell =1}^L w_\ell =1$ ,

$$ \begin{align*} \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1,{\mathrm{{BFI}}}} - \boldsymbol{\theta}_1\big) \leadsto {\mathcal N}\Big(0,\big(I_{1}\big)^{-1}\Big). \end{align*} $$

Further, we have seen that the BFI estimator

$$ \begin{align*} \frac{1}{n}\widehat{\mathbf{A}}_{1,{\mathrm{{BFI}}}} &=\; \frac{1}{n}\sum_{\ell=1}^L \widehat{\mathbf{A}}_{1,\ell}+\frac{1}{n}\boldsymbol{\Lambda}_1-\frac{1}{n}\sum_{\ell=1}^L \boldsymbol{\Lambda}_{1,\ell} \\ &=\; \sum_{\ell=1}^L \frac{n_\ell}{n}\; \frac{1}{n_\ell}\widehat{\mathbf{A}}_{1,\ell}\;+\;O_p\Big(\frac{1}{n}\Big) \;=\; \sum_{\ell=1}^L w_\ell I_{1,\ell} \;+\;o_p(1), \end{align*} $$

converges in probability to $\sum _{\ell =1}^L w_\ell I_{1,\ell }$ .

If the number of centers L increases to infinity as well, but $L=o_p(n)$ (i.e., the number of centers is smaller in rate than the total sample size), the asymptotic results remain valid, but the derivation needs to be adjusted slightly.

Asymptotic distribution of the WAV and Single center estimators

Suppose the MAP estimator $\widehat {\mathbf {\theta }}_{1,\ell }$ in center $\ell $ is used for estimating the parameter $\mathbf {\theta }_1$ , then we obtain

$$ \begin{align*} \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1,\ell} - \boldsymbol{\theta}_1\big) \;=\; \frac{\sqrt{n}}{\sqrt{n_\ell}} \; \sqrt{n_\ell}\big(\widehat{\boldsymbol{\theta}}_{1,\ell} - \boldsymbol{\theta}_1\big) \leadsto {\mathcal N}\Big(0,\big(w_\ell I_{1,\ell}\big)^{-1}\Big), \end{align*} $$

if $w_\ell> 0$ . That means that the single-center estimator is not optimal for estimating $\mathbf {\theta }_1$ unless $w_\ell =1$ (there is only one center).

For the weighted average estimator $\sum _{\ell =1}^L \frac {n_\ell }{n} \widehat {\mathbf {\theta }}_{1,\ell }$ for estimating $\mathbf {\theta }_1$ , we find

$$ \begin{align*} \sqrt{n}\Big(\sum_{\ell=1}^L \frac{n_\ell}{n}\; \widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\Big) \;=\; \sum_{\ell=1}^L \frac{n_\ell}{n}\; \sqrt{n} \big(\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\big) \;=\; \sum_{\ell=1}^L \frac{\sqrt{n_\ell}}{\sqrt{n}} \sqrt{n_\ell} \big(\widehat{\boldsymbol{\theta}}_{1,\ell}-\boldsymbol{\theta}_1\big) \leadsto {\mathcal N}\Big(0,\sum_{\ell=1}^Lw_\ell \big( I_{1,\ell}\big)^{-1}\Big), \end{align*} $$

with, in the homogeneous setting $I_{1,\ell }=I_1, \ell =1,\ldots ,L$

$$ \begin{align*} \sum_{\ell=1}^Lw_\ell \big( I_{1,\ell}\big)^{-1}=(I_1)^{-1} = (J_1)^{-1} \end{align*} $$

as defined in (A.13). In the homogeneous setting, the weighted average estimator is asymptotically efficient as well.

Appendix III.C Asymptotic distribution of the BFI estimator in heterogeneous setting

Asymptotic distribution BFI estimator

Let $\widehat {\mathbf {A}}_{1a,\ell }$ and $\widehat {\mathbf {A}}_{1b,\ell }$ be the second derivatives of $-\log \{p_\ell (\mathbf {\theta }_1 | \mathbf {D}_\ell )\}$ with respect to $\mathbf {\theta }_{1a}$ and $\mathbf {\theta }_{1b}$ , respectively, and let $\widehat {\mathbf {A}}_{1ab,\ell }$ be the second derivative with respect to both $\mathbf {\theta }_{1a}$ and $\mathbf {\theta }_{1b}$ , all evaluated at the MAP estimator $\widehat {\mathbf {\theta }}_{1,\ell }$ . If these second derivatives are sufficiently smooth in the neighborhood of $\mathbf {\theta }_{1}$ , it follows by the law of large numbers that

$$ \begin{align*} \frac{1}{n_\ell}\widehat{\mathbf{A}}_{1a,\ell} \rightarrow I_{1a,\ell} \quad \frac{1}{n_\ell}\widehat{\mathbf{A}}_{1b,\ell} \rightarrow I_{1b,\ell} \quad \frac{1}{n_\ell}\widehat{\mathbf{A}}_{1ab,\ell} \rightarrow I_{1ab,\ell}, \end{align*} $$

with the matrices $I_{1a,\ell }, I_{1b,\ell }$ , and $I_{1ab,\ell }$ for center $\ell $ .

Asymptotic behaviour of the BFI estimator for $\mathbf {\theta _{1a}}$

In Appendix II.B we computed that, under the assumption that $\Lambda _{1b\ell }=\Lambda _{1b,\ell }$ , the BFI estimator for $\mathbf {\theta }_{1a}$ is equal to:

(A.17) $$ \begin{align} \widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}}&:=\Big( \sum_{\ell=1}^L \Big(\widehat{\mathbf{A}}_{1a,\ell}- \widehat{\mathbf{A}}_{1ab,\ell}(\widehat{\mathbf{A}}_{1b,\ell})^{-1}(\widehat{\mathbf{ A}}_{1ab,\ell})^t\Big) +\boldsymbol{\Lambda}_{1a}-\sum_{\ell=1}^L\boldsymbol{\Lambda}_{1a,\ell} \Big)^{\!-1} \times \nonumber\\ &\quad\sum_{\ell=1}^L \Big(\widehat{\mathbf{A}}_{1a,\ell} - \widehat{\mathbf{A}}_{1ab,\ell}(\widehat{\mathbf{A}}_{1b,\ell})^{-1}(\widehat{\mathbf{A}}_{1ab,\ell})^t\Big)\widehat{\boldsymbol{\theta}}_{1a,\ell}. \end{align} $$

Define $\widehat {\Gamma }_\ell := n_\ell ^{-1}\widehat {\mathbf {A}}_{1a,\ell }- n_\ell ^{-1}\widehat {\mathbf {A}}_{1ab,\ell }(n_\ell ^{-1}\widehat {\mathbf {A}}_{1b,\ell })^{-1}(n_\ell ^{-1}\widehat {\mathbf {A}}_{1ab,\ell })^t$ . Then,

(A.18) $$ \begin{align} \widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}}=\Big(\sum_{\ell=1}^L \frac{n_\ell}{n}\;\widehat{\Gamma}_\ell + \frac{1}{n}\;\boldsymbol{\Lambda}_{1a}-\sum_{\ell=1}^L\frac{1}{n}\;\boldsymbol{\Lambda}_{1a,\ell} \Big)^{\!-1} \sum_{\ell=1}^L \frac{n_\ell}{n}\;\widehat{\Gamma}_\ell \;\widehat{\boldsymbol{\theta}}_{1a,\ell}. \end{align} $$

If the derivatives are sufficiently smooth, it follows by the law of large numbers, continuous mapping theorem, and Slutsky’s lemma,Reference van der Vaart²⁸ that $\widehat {\Gamma }_\ell ^{-1}$ converges in probability to $\big (I_{1a,\ell }-I_{1ab,\ell }(I_{1b,\ell })^{-1}(I_{1ab,\ell })^t\big )^{-1}$ , which equals the left-upper block of the inverse of the Fisher information matrix $I_{1,\ell }$ in center $\ell $ . Now, based on similar calculations as in the homogeneous setting

$$ \begin{align*} \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}} - \boldsymbol{\theta}_{1a}\big) \;\leadsto\; {\mathcal N}\Big(0,\Big(\sum_{\ell=1}^L w_\ell \big(I_{1a,\ell}-I_{1ab,\ell}(I_{1b,\ell})^{-1}(I_{1ab,\ell})^t\big)\Big)^{-1}\Big), \end{align*} $$

with the asymptotic covariance matrix equal to the asymptotic covariance matrix of the MAP estimator based on all data, given in Equation (A.14). The BFI estimator $\widehat {\mathbf {\theta }}_{1a,{\mathrm { {BFI}}}}$ is asymptotically efficient for estimating $\mathbf {\theta }_{1a}$ .

If $I_{1a,\ell }, I_{1b,\ell }$ and $I_{1ab,\ell }$ equal across the centers:

$$ \begin{align*} \Big(\sum_{\ell=1}^L w_\ell \big\{I_{1a,\ell}-I_{1ab,\ell}(I_{1b,\ell})^{-1}(I_{1ab,\ell})^t\big\}\Big)^{-1} = \Big(I_{1a,\ell}-I_{1ab,\ell}(I_{1b,\ell})^{-1}(I_{1ab,\ell})^t\Big)^{-1}. \end{align*} $$

Asymptotic behaviour of the BFI estimator of $\mathbf {\theta _{1b,\ell }}$

Under the assumption that $\boldsymbol {\Lambda }_{1b\ell }=\boldsymbol {\Lambda }_{1b,\ell }$ and $w_\ell>0$ , the BFI estimator for $\mathbf {\theta }_{1b,\ell }$ is given by

$$ \begin{align*} \widehat{\boldsymbol{\theta}}_{1b,\ell,{\mathrm{{BFI}}}}& \;=\; \big(\widehat{\mathbf{A}}_{1b,\ell}\big)^{-1}\Big[\widehat{\mathbf{A}}_{1b,\ell}\widehat{\boldsymbol{\theta}}_{1b,\ell} + \big(\widehat{\mathbf{A}}_{1ab,\ell}\big)^t\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}}\big)\Big]\\& \;=\; \widehat{\boldsymbol{\theta}}_{1b,\ell} + \big(\widehat{\mathbf{A}}_{1b,\ell}\big)^{-1}\big(\widehat{\mathbf{A}}_{1ab,\ell}\big)^t\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}}\big). \end{align*} $$

Then,

(A.19) $$ \begin{align} \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1b,\ell,{\mathrm{{BFI}}}}-\boldsymbol{\theta}_{1b,\ell}\big) &=\; \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1b,\ell}-\boldsymbol{\theta}_{1b,\ell}\big) + \big(\widehat{\mathbf{ A}}_{1b,\ell}\big)^{-1}\big(\widehat{\mathbf{A}}_{1ab,\ell}\big)^t\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}}\big)\nonumber \\[4pt] &=\; \sqrt{n}(\widehat{\boldsymbol{\theta}}_{1b,\ell}-\boldsymbol{\theta}_{1b,\ell}) + \big(\widehat{\mathbf{A}}_{1b,\ell}\big)^{-1}\big(\widehat{\mathbf{A}}_{1ab,\ell}\big)^t\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \boldsymbol{\theta}_{1a}\big)\nonumber\\&\quad \;+\; \big(\widehat{\mathbf{A}}_{1b,\ell}\big)^{-1}\big(\widehat{\mathbf{A}}_{1ab,\ell}\big)^t\sqrt{n}\big(\boldsymbol{\theta}_{1a} - \widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}}\big). \end{align} $$

We first leave out the last term and consider the asymptotic behaviour of

$$ \begin{align*} \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1b,\ell}-\boldsymbol{\theta}_{1b,\ell}\big) + \big(\widehat{\mathbf{A}}_{1b,\ell}\big)^{-1}\big(\widehat{\mathbf{A}}_{1ab,\ell}\big)^t\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \boldsymbol{\theta}_{1a}\big). \end{align*} $$

As explained before, this term equals

$$ \begin{align*} \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1b,\ell}-\boldsymbol{\theta}_{1b,\ell}\big) &+ \big(I_{1b,\ell}\big)^{-1}\big(I_{1ab,\ell}\big)^t\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \boldsymbol{\theta}_{1a}\big) \;+\; o_p(1)\\[4pt] &=\; w_\ell^{-1/2} \sqrt{n_\ell}\big(\widehat{\boldsymbol{\theta}}_{1b,\ell}-\boldsymbol{\theta}_{1b,\ell}\big) + w_\ell^{-1/2}\big(I_{1b,\ell}\big)^{-1}\big(I_{1ab,\ell}\big)^t\sqrt{n_\ell}\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \boldsymbol{\theta}_{1a}\big) \;+\; o_p(1). \end{align*} $$

In center $\ell $ the MAP estimator $\widehat {\mathbf {\theta }}_{1,\ell }$ is asymptotically normal (Bernstein–Von Mises TheoremReference van der Vaart²⁸):

$$ \begin{align*} \sqrt{n_\ell} \begin{pmatrix} \begin{pmatrix} \widehat{\boldsymbol{\theta}}_{1a,\ell}\\ \widehat{\boldsymbol{\theta}}_{1b,\ell} \end{pmatrix} - \begin{pmatrix} \boldsymbol{\theta}_{1a}\\ \boldsymbol{\theta}_{1b,\ell} \end{pmatrix} \end{pmatrix} \;\leadsto \; {\mathcal N}\big(\mathbf{0}, (I_{1,\ell})^{-1} \big), \end{align*} $$

Define the function $g(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b,\ell }) = \mathbf {\theta }_{1b,\ell } + \big (I_{1b,\ell }\big )^{-1}\big (I_{1ab,\ell }\big )^t \mathbf {\theta }_{1a}$ . Then, by the continuous mapping theorem

$$ \begin{align*} \sqrt{n_\ell}(\widehat{\boldsymbol{\theta}}_{1b,\ell}-\boldsymbol{\theta}_{1b,\ell}) &+ \big(I_{1b,\ell}\big)^{-1}\big(I_{1ab,\ell}\big)^t\sqrt{n_\ell}\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \boldsymbol{\theta}_{1a}\big) =\; \sqrt{n}_\ell \begin{pmatrix} g\begin{pmatrix} \widehat{\boldsymbol{\theta}}_{1a,\ell}\\ \widehat{\boldsymbol{\theta}}_{1b,\ell} \end{pmatrix} - g\begin{pmatrix} \boldsymbol{\theta}_{1a}\\ \boldsymbol{\theta}_{1b,\ell} \end{pmatrix} \end{pmatrix}\\[8pt] &\leadsto \; {\mathcal N}\big(\mathbf{0}, g'(\boldsymbol{\theta}_{1a},\boldsymbol{\theta}_{1b,\ell})^t \big(I_{1,\ell}\big)^{-1} g'(\boldsymbol{\theta}_{1a},\boldsymbol{\theta}_{1b,\ell}) \big) \end{align*} $$

with $g'(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b,\ell })$ the derivative of the function g in $(\mathbf {\theta }_{1a},\mathbf {\theta }_{1b,\ell })$ . Straightforward calculations show that

$$ \begin{align*}g'(\boldsymbol{\theta}_{1a},\boldsymbol{\theta}_{1b,\ell})^t \big(I_{1,\ell}\big)^{-1} g'(\boldsymbol{\theta}_{1a},\boldsymbol{\theta}_{1b,\ell}) = (I_{1b,\ell})^{-1},\end{align*} $$

which equals the asymptotic covariance matrix of the Gaussian limit distribution for $\widehat {\mathbf {\theta }}_{1b,\ell }$ if the parameter $\mathbf {\theta }_{1a}$ is known. Apparently, leaving out the last term in equation (A.19), means that we assume that the BFI estimator for $\mathbf {\theta }_{1a}$ is (almost) equal to the true value $\mathbf {\theta }_{1a}$ , and thus that $\mathbf {\theta }_{1a}$ is (almost) known. The result is interesting; the asymptotic accuracy of the BFI estimator for $\mathbf {\theta }_{1b,\ell }$ is increased, because the parameter $\mathbf {\theta }_{1a}$ can be estimated more accurately with the BFI estimator (and thus using information from the other centers) compared to the situation in which $\mathbf {\theta }_{1a}$ is estimated based on data from center $\ell $ only.

We go back to the expression in (A.19):

$$ \begin{align*} \sqrt{n}(\widehat{\boldsymbol{\theta}}_{1b,\ell,{\mathrm{{BFI}}}}&-\boldsymbol{\theta}_{1b,\ell}) \;=\; \sqrt{n}(\widehat{\boldsymbol{\theta}}_{1b,\ell}-\boldsymbol{\theta}_{1b,\ell}) + \big(\widehat{\mathbf{A}}_{1b,\ell}\big)^{-1}\big(\widehat{\mathbf{A}}_{1ab,\ell}\big)^t\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \widehat{\boldsymbol{\theta}}_{1a,{\mathrm{{BFI}}}}\big)\\[4pt] &=\; \sqrt{n}(\widehat{\boldsymbol{\theta}}_{1b,\ell}-\boldsymbol{\theta}_{1b,\ell}) - \big(I_{1b,\ell}\big)^{-1}\big(I_{1ab,\ell}\big)^t \times\\ &\quad\Bigg(\Big(\sum_{k=1}^L \frac{n_k}{n}\;\widehat{\Gamma}_k \Big)^{\!-1} \sum_{k=1}^L \frac{n_k}{n}\;\widehat{\Gamma}_k \;\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,k} - \widehat{\boldsymbol{\theta}}_{1a,\ell}\big)\Bigg) + o_P(1)\\[8pt] &=\; \sqrt{n}(\widehat{\boldsymbol{\theta}}_{1b,\ell}-\boldsymbol{\theta}_{1b,\ell}) - \big(I_{1b,\ell}\big)^{-1}\big(I_{1ab,\ell}\big)^t \times\\ &\quad\Bigg(\Big(\sum_{k=1}^L w_k\;\Gamma_k \Big)^{\!-1} \sum_{k=1,k\neq \ell}^L w_k\;\Gamma_k \;\big(\sqrt{n}(\widehat{\boldsymbol{\theta}}_{1a,k} - \boldsymbol{\theta}_{1a}) - \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \boldsymbol{\theta}_{1a}) \big)\Bigg) + o_P(1)\\ &=\; \sqrt{n}(\widehat{\boldsymbol{\theta}}_{1b,\ell}-\boldsymbol{\theta}_{1b,\ell}) + \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \boldsymbol{\theta}_{1a}\big) \big(I_{1b,\ell}\big)^{-1}\big(I_{1ab,\ell}\big)^t\Big(\sum_{k=1}^L w_k\;\Gamma_k \Big)^{\!-1} \sum_{k=1,k\neq \ell}^L w_k\;\Gamma_k\\ & \quad -\big(I_{1b,\ell}\big)^{-1}\big(I_{1ab,\ell}\big)^t\Big(\sum_{k=1}^L w_k\;\Gamma_k \Big)^{\!-1} \sum_{k=1,k\neq \ell}^L w_k\;\Gamma_k \;\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,k} - \boldsymbol{\theta}_{1a}\big) \; + o_P(1). \end{align*} $$

The asymptotic distribution can be obtained with the continuous mapping theorem and Slutky’s lemma again. The third and last term depends on data from all centers except from center $\ell $ . Since the data from the different centers are assumed to be independent, it is sufficient to show that the asymptotic distributions of the sum of the first and second terms and of the third term are asymptotically normal and next add the asymptotic mean (which are zero) and the variances. The asymptotic distribution of the sum of the first and second term can obtained with the continuous mapping theorem, like before. The third term is asymptotically normal:

$$ \begin{align*} \big(I_{1b,\ell}\big)^{-1}\big(I_{1ab,\ell}\big)^t\Big(\sum_{k=1}^L w_k\;\Gamma_k \Big)^{\!-1} \sum_{k=1,k\neq \ell}^L \sqrt{w_k}\;\Gamma_k \;\sqrt{n_k}\big(\widehat{\boldsymbol{\theta}}_{1a,k} - \boldsymbol{\theta}_{1a}\big) \leadsto {\mathcal N}(\mathbf{0}, \Sigma) \end{align*} $$

with

$$ \begin{align*} \Sigma &= \big(I_{1b,\ell}\big)^{-1} \big(I_{1ab,\ell}\big)^{t} \Big(\sum_{k=1}^L w_k \Gamma_k\Big)^{-1}\Big(\sum_{k\neq l}w_k\Gamma_k\Big)\Big(\sum_{k=1}^L w_k \Gamma_k\Big)^{-1} I_{1ab,\ell} \big(I_{1b,\ell}\big)^{-1} \\ &= \big(I_{1b,\ell}\big)^{-1} \big(I_{1ab,\ell}\big)^{t} \Big(\sum_{k=1}^L w_k \Gamma_k\Big)^{-1}\Big(\sum_{k=1}^Lw_k\Gamma_k - w_\ell\Gamma_\ell\Big)\Big(\sum_{k=1}^L w_k \Gamma_k\Big)^{-1} I_{1ab,\ell} \big(I_{1b,\ell}\big)^{-1}\\ &= \big(I_{1b,\ell}\big)^{-1} \big(I_{1ab,\ell}\big)^{t} \Big(I - \Big(\sum_{k=1}^L w_k \Gamma_k\Big)^{-1}w_\ell\Gamma_\ell\Big)\Big(\sum_{k=1}^L w_k \Gamma_k\Big)^{-1} I_{1ab,\ell} \big(I_{1b,\ell}\big)^{-1}\\ &= \big(I_{1b,\ell}\big)^{-1} \big(I_{1ab,\ell}\big)^{t} \Big(\sum_{k=1}^L w_k \Gamma_k\Big)^{-1} I_{1ab,\ell} \big(I_{1b,\ell}\big)^{-1} \\&\quad- \big(I_{1b,\ell}\big)^{-1} \big(I_{1ab,\ell}\big)^{t} \Big(\sum_{k=1}^L w_k \Gamma_k\Big)^{-1}w_\ell\Gamma_\ell\Big(\sum_{k=1}^L w_k \Gamma_k\Big)^{-1} I_{1ab,\ell} \big(I_{1b,\ell}\big)^{-1}. \end{align*} $$

If we sum up the variances of the asymptotic zero mean Gaussian distributions of all terms, we obtain

$$ \begin{align*} \big(I_{1b,\ell}\big)^{-1} \big(I_{1ab,\ell}\big)^{t} \Big(\sum_{k=1}^L w_k \Gamma_k\Big)^{-1} I_{1ab,\ell} \big(I_{1b,\ell}\big)^{-1} + \big(w_\ell I_{1b,\ell}\big)^{-1}. \end{align*} $$

So,

$$ \begin{align*} &\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1b,\ell,{\mathrm{{BFI}}}} - \boldsymbol{\theta}_{1b,\ell}\big) \leadsto \\&\quad{\mathcal N}\Big(\mathbf{0}, \big(w_\ell I_{1b,\ell}\big)^{-1}+\big(I_{1b,\ell}\big)^{-1} \big(I_{1ab,\ell}\big)^{t} \Big(\sum_{k=1}^L w_k \big(I_{1a,k}-I_{1ab,k}(I_{1b,k})^{-1}(I_{1ab,k})^t\big)\Big)^{-1} I_{1ab,\ell} \big(I_{1b,\ell}\big)^{-1}\Big). \end{align*} $$

The BFI estimator follows, asymptotically, the same distribution as the MAP estimator based on the combined data, see Equation (A.16), and is asymptotically efficient. If $\boldsymbol {\Lambda }_{1b\ell } \neq \boldsymbol {\Lambda }_{1b,\ell }$ the asymptotic distribution of $\widehat {\mathbf {\theta }}_{1b,\ell ,{\mathrm { {BFI}}}}$ will not change, because $n^{-1}\boldsymbol {\Lambda }_{1b\ell }$ and $n^{-1}\boldsymbol {\Lambda }_{1b,\ell }$ converge to zero.

Like in the homogeneous setting, it can be directly seen that the BFI estimator $n^{-1} \widehat {\mathbf {A}}_{1a,{\mathrm { {BFI}}}}$ converges in probability to $\sum _{\ell =1}^L w_\ell I_{1a,\ell }$ . Similarly, the BFI estimator $n^{-1} \widehat {\mathbf {A}}_{1b,\ell ,{\mathrm { {BFI}}}}$ converges in probability to $w_\ell I_{1b,\ell }$ , and $n^{-1} \widehat {\mathbf {A}}_{1ab,\ell ,{\mathrm { {BFI}}}}$ to $w_\ell I_{1ab,\ell }$ .

Asymptotic distribution of the WAV and single center estimators

If the parameter $\mathbf {\theta }_{1a}$ is estimated by the MAP estimator from center $\ell $ , $\widehat {\mathbf {\theta }}_{1a,\ell }$ , the asymptotic distribution equals

$$ \begin{align*} \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1a,\ell} - \boldsymbol{\theta}_{1a}\big) \;\leadsto\; {\mathcal N}\Big(0,w_\ell^{-1}\Big(I_{1a,\ell}-I_{1ab,\ell}(I_{1b,\ell})^{-1}(I_{1ab,\ell})^t\Big)^{-1}\Big). \end{align*} $$

If the weighted average estimator $\sum _{\ell =1}^L \frac {n_\ell }{n} \widehat {\mathbf {\theta }}_{1a,\ell }$ is used, the asymptotic distribution equals, like before

$$ \begin{align*} \sqrt{n}\Big(\sum_{\ell=1}^L \frac{n_\ell}{n}\; \widehat{\boldsymbol{\theta}}_{1a,\ell} - \boldsymbol{\theta}_{1a}\Big) \;\leadsto\; {\mathcal N}\Big(0,\sum_{\ell=1}^Lw_\ell\Big( I_{1a,\ell}-I_{1ab,\ell}(I_{1b,\ell})^{-1}(I_{1ab,\ell})^t\Big)^{-1}\Big). \end{align*} $$

Asymptotically, both estimators are not efficient, unless $I_{1,\ell }=I_1, \ell =1,\ldots ,L$ . In that case the weighted average estimator is asymptotically efficient as well.

If we estimate $\mathbf {\theta }_{1b,\ell }$ by its MAP estimator in the corresponding center, the asymptotic distribution equals:

$$ \begin{align*}& \sqrt{n}\big(\widehat{\boldsymbol{\theta}}_{1b,\ell} - \boldsymbol{\theta}_{1b,\ell}\big) \leadsto \\&\quad{\mathcal N}\Bigg(\mathbf{0}, w_\ell^{-1}\Big( I_{1b,\ell}^{-1} + \big(I_{1b,\ell}\big)^{-1} (I_{1ab,\ell})^t \big(I_{1a,\ell}-I_{1ab,\ell} (I_{1b,\ell})^{-1} (I_{1ab,\ell})^t\big)^{-1} I_{1ab,\ell} (I_{1b,\ell})^{-1}\Big)\Bigg). \end{align*} $$

Since $\mathbf {\theta }_{1b,\ell }$ is center-specific, the WAV estimator equals the single center estimator. The estimator $(\widehat {\mathbf {\theta }}_{1a,\ell },\widehat {\mathbf {\theta }}_{1b,\ell })$ is based on data from center $\ell $ only, and does not use any information from the other centers for estimating the parameter $\mathbf {\theta }_{1a}$ . Using information from the other centers, $\mathbf {\theta }_{1a}$ is estimated more accurately and improves the estimation of $\mathbf {\theta }_{1b}$ . That is why the BFI estimator $\widehat {\mathbf {\theta }}_{1b,\ell ,{\mathrm { {BFI}}}}$ is more accurate than the weighted average estimator.

The asymptotic distribution of the BFI estimators for clustered data (Appendix II.C) can be derived in a similar way, but this is more complicated due to the complex expressions of these estimators. However, since the BFI estimators in the homogeneous and the heterogeneous case with center-specific parameters have shown to be asymptotically efficient and these settings are special cases of the one with clustered data, it is expected that the BFI estimators for the clustered data are asymptotically efficient as well.