# Bayesian Federated Inference
for Statistical Models

Marianne A Jonker
Department for Health Evidence, Section Biostatistics,
Radboudumc, Nijmegen, The Netherlands
and
Hassan Pazira
Department for Health Evidence, Section Biostatistics,
Radboudumc, Nijmegen, The Netherlands
and
Anthony CC Coolen
Donders Institute, Faculty of Science, Radboud University,
Nijmegen, The Netherlands
Saddle Point Science Europe, Mercator Science Park,
Nijmegen, The Netherlands corresponding author<EMAIL_ADDRESS>
###### Abstract

Identifying predictive factors via multivariable statistical analysis is for
rare diseases often impossible because the data sets available are too small.
Combining data from different medical centers into a single (larger) database
would alleviate this problem, but is in practice challenging due to regulatory
and logistic problems. Federated Learning (FL) is a machine learning approach
that aims to construct from local inferences in separate data centers what
would have been inferred had the data sets been merged. It seeks to harvest
the statistical power of larger data sets without actually creating them. The
FL strategy is not always feasible for small data sets. Therefore, in this
paper we refine and implement an alternative Bayesian Federated Inference
(BFI) framework for multi center data with the same aim as FL. The BFI
framework is designed to cope with small data sets by inferring locally not
only the optimal parameter values, but also additional features of the
posterior parameter distribution, capturing information beyond that is used in
FL. BFI has the additional benefit that a single inference cycle across the
centers is sufficient, whereas FL needs multiple cycles. We quantify the
performance of the proposed methodology on simulated and real life data.

Keywords: Federated Learning, data integration, multi center data, MAP
estimator, small data sets.

##  1 Introduction

Medical data are not always easily shared for research by those who hold them
(companies, hospitals, research centers, universities, consortia). This
hampers their integrated use in estimation, inference and prediction. It
reduces our ability to identify prognostic factors reliably and to predict the
clinical outcomes for future patients, due to possible mis-identification of
prognostic factors and overfitting of the statistical model. This is
especially a problem for rare diseases, for which only small data sets will
typically be available for research. Not surprisingly, therefore, the
prognosis of patients with rare diseases and in particular rare cancers is
worse than that of patients with more common ones. It follows that we must
either create more effective mechanisms and incentives for data sharing
between institutions, or focus on technology for the integration of individual
analysis outcomes obtained on localized data sets. If we go for the second
route, the scientific problem is how to extract from the results of analyses
on data subsets the result that would have been found had the subsets been
combined into a single data set.

Federated learning (FL) is an on-line (i.e. non-batch) machine learning
approach, developed several years ago in the context of interacting mobile
devices (McMahan et al., 2017), in which local centers use local data for
training machine learning systems on site by optimizing a pre-specified loss
function, and only the estimated parameters of the trained systems are sent
out for integration at a central server. The data, in contrast, stay at their
owners’ institutions. This procedure is ‘cycled’ around the centers
iteratively; all centers update their estimates in an on-line fashion until a
convergence criterion is met (see Figure 1, left plot). Typically, deep
learning neural networks are employed, where the circulated parameters are
strengths of interactions between nodes and activation thresholds of nodes.
While FL performs excellently in image analysis (see e.g. Rieke et al., 2020;
Sheller et al., 2020; Gafni et al., 2021), it has limitations in application
to rare disease data. First, in its standard form FL still requires relatively
large data sets (given the large numbers of deep learning parameters). Second,
FL does not generate rigorous error bars, since only the most probable
parameters are iterated (as opposed to posterior distributions). Third, it is
not clear that this process will always converge to a satisfactory final
state, since the latent heterogeneity across centers may generate
prohibitively conflicting gradients for the model parameters in the
maximization procedure, leading to an end result that works for none of the
local data sets.

In recent years much progress has been made in FL (e.g. improved optimization
in each center, aggregation of the local models at the central server and
dealing with heterogeneity of the local populations (Li et al., 2020; Chen and
Chao, 2020; Zhu et al., 2021). An overview of the recent developments
including references is given in Liu et al. (2021). Multiple researchers have
proposed FL in a Bayesian setting for deep learning models (e.g. Maddox et
al., 2019; Al-Shedivat et al., 2020): In each local center the posterior
distribution is estimated and communicated to the central server, where
aggregation takes place. However, this Bayesian procedure turns out to be
challenging, especially for deep learning models in which the parameter
dimensionality is high. Proposals to address this include approximating the
estimated local and global posterior distribution with MCMC (Zhang et al.,
2019; Izmailov et al., 2021), variational inference (Zhang et al., 2018)), or
the use of multivariate Gaussian distributions with a Laplace approximation
(Maddox et al. 2019; Al-Shedivat et al. 2020). For a linear model (as an
example) Al-Shedivat et al. (2020) approximated the local posterior by a
multivariate Gaussian. Under the assumption of a flat prior for the model
parameters in the local centers, the global posterior distribution in the
central server is multivariate Gaussian as well, and estimates of its
parameters can be computed from the local sample means and covariance
matrices. Although this seems to be an interesting idea on paper, the authors
directly note that this way of estimating the global posterior distribution is
not feasible in a general setting where models are neural networks with
millions of parameters, as it requires estimating and computing the inverse of
an $d\times d$ matrix, where $d$ is the number of parameters.

In the field of statistical modelling (e.g. via generalized linear models) the
challenges in FL are not about the computational burden due to the high
dimensionality of the model parameters, but rather about the complexities of
(medical) data, such as heterogeneity of populations, random or structural
missingness of covariates, presence of confounding factors, overfitting, and
the interpretation of the parameters in the estimated models. In this paper a
Bayesian Federated Inference (BFI) framework is proposed and developed for
arbitrary generalized linear regression models, in which the challenges
described above are addressed. More specifically, the posterior distributions
in the local centers and for the fictive combined data are approximated by
multivariate Gaussian distributions around the maximum a posteriori (MAP)
estimate. By choosing either a multivariate Gaussian or uninformative
distribution for the prior distribution, the outcome of the inference on the
combined set can be expressed directly in terms of the outcomes of the
separate inferences on the local data sets. With the proposed methodology,
there is no need to do inference on the full data set, as its inference
results can be computed a posteriori from the inference results in the subsets
in only one round, in contrast to traditional FL where very many iterative
inference cycles across the centers are needed (see plot on the right in
Figure 1). The proposed BFI framework is also accurate for small data sets by
choosing the prior distribution wisely.

This paper is organized as follows. In Section 2 the framework for the BFI
strategy is explained in a general setting. Next, in Section 3 we focus on BFI
for generalized linear models. Here we also address population heterogeneity
and missing covariates. Simulation studies based on real life data have been
performed to quantify the efficiency in statistical inference if the data are
only locally available for analysis, and the impact of complexities related to
heterogeneity and small sample sizes in the local centers (Section 4). The
paper ends with a discussion in Section 5.

Figure 1: Left: Visualization of an iterative FL procedure. Right:
Visualization of the BFI procedure proposed in this paper. One cycle is
sufficient for this aggregation.

##  2 Bayesian Federated Inference (BFI)

###  2.1 Problem definition and setting

Suppose that the stochastic variable $Y$ is distributed with parametric
density $y\to p(y|\mbox{\boldmath$\theta$})$ that is known up to the parameter
$\theta$, which itself has a density $\mbox{\boldmath$\theta$}\to
p(\mbox{\boldmath$\theta$})$. Data for statistical inference is available from
$L$ non-overlapping data sets from different centers (e.g. hospitals). These
data sets, $D_{1},\ldots,D_{L}$, and their union are defined as:

$\displaystyle D_{\ell}=\\{y_{\ell 1},\ldots,y_{\ell
n_{\ell}}\\},~{}~{}~{}~{}~{}~{}~{}~{}n_{\ell}=|D_{\ell}|,~{}~{}~{}~{}~{}~{}~{}~{}D=\bigcup_{\ell=1}^{L}D_{\ell},$

where $y_{\ell i}$ is the $i^{th}$ observation in the $\ell^{th}$ data subset
$D_{\ell}$. The observations within the local data sets, and also across the
$L$ data sets, are assumed to be independent. We imagine the scenario where
the data from the $L$ subsets are prohibited from being combined into one data
set $D$; inference can only be performed for the $L$ subsets separately, and
only the latter inference results can be combined. In the next subsection we
aim to express the outcome of inference on the combined set $D$ in terms of
the outcomes of the $L$ separate inferences on the constituent sets
$D_{\ell}$.

###  2.2 Formulae for Bayesian subset inference

In a Bayesian analysis based on any data set $D=\\{y_{1},\ldots,y_{n}\\}$ from
a density $y\to p(y|\mbox{\boldmath$\theta$})$ with prior
$\mbox{\boldmath$\theta$}\to p(\mbox{\boldmath$\theta$})$, the posterior
density equals 111We use the letter $p$ for any density. From the arguments it
follows which density is meant.:

$\displaystyle
p(\mbox{\boldmath$\theta$}|D)\;=\;\frac{p(D|\mbox{\boldmath$\theta$})p(\mbox{\boldmath$\theta$})}{Z(D)},~{}~{}~{}~{}~{}Z(D)\;=\;\int\\!p(D|\mbox{\boldmath$\theta$})p(\mbox{\boldmath$\theta$})\mathrm{d}\mbox{\boldmath$\theta$},~{}~{}~{}~{}~{}p(D|\mbox{\boldmath$\theta$})=\prod_{i=1}^{n}p(y_{i}|\mbox{\boldmath$\theta$}).$
(1)

The maximum a posteriori (MAP) estimator for $\theta$ is obtained by
maximizing $\mbox{\boldmath$\theta$}\to p(\mbox{\boldmath$\theta$}|D)$. Since
the complete data set $D$ is the union of $L$ data sets from different
centers, we can write the posterior density as

$\displaystyle
p(\mbox{\boldmath$\theta$}|D)\;=\;\frac{p(\mbox{\boldmath$\theta$})}{Z(D)}\prod_{\ell=1}^{L}\prod_{i=1}^{n_{\ell}}p(y_{\ell
i}|\mbox{\boldmath$\theta$}).$ (2)

The expressions in (1) hold for any data set, so they apply also to the sets
$D_{\ell},\ell=1,\ldots,L$. By rewriting the expression in (2) with $D$
replaced by $D_{\ell}$, and with the subset priors denoted by
$p_{\ell}(\mbox{\boldmath$\theta$})$, we find that

$\displaystyle\prod_{i=1}^{n_{\ell}}p(y_{\ell i}|\mbox{\boldmath$\theta$})$
$\displaystyle\\!=\\!$
$\displaystyle\frac{Z_{\ell}(D_{\ell})}{p_{\ell}(\mbox{\boldmath$\theta$})}\;p_{\ell}(\mbox{\boldmath$\theta$}|D_{\ell}),~{}~{}~{}~{}~{}~{}~{}~{}Z_{\ell}(D_{\ell})\;=\;\int\\!p(D_{\ell}|\mbox{\boldmath$\theta$})p_{\ell}(\mbox{\boldmath$\theta$})\mathrm{d}\mbox{\boldmath$\theta$}.$

By substituting this into (2), we obtain:

$\displaystyle p(\mbox{\boldmath$\theta$}|D)$ $\displaystyle\\!=\\!$
$\displaystyle\frac{1}{C}\;\bigg{(}\frac{p(\mbox{\boldmath$\theta$})}{\prod_{\ell=1}^{L}p_{\ell}(\mbox{\boldmath$\theta$})}\bigg{)}\prod_{\ell=1}^{L}p_{\ell}(\mbox{\boldmath$\theta$}|D_{\ell}),~{}~{}~{}~{}\mbox{
with }~{}~{}~{}~{}C=\frac{Z(D)}{\prod_{\ell=1}^{L}Z_{\ell}(D_{\ell})}.$ (3)

The constant $C$ can always be recovered via normalization. Hence we can in a
relatively simple way express the posterior parameter density
$p(\mbox{\boldmath$\theta$}|D)$ in terms of the $L$ local densities
$p_{\ell}(\mbox{\boldmath$\theta$}|D_{\ell})$. The next question is under
which conditions we can obtain accurate (possibly approximate) expressions for
the MAP estimator $\hat{}\mbox{\boldmath$\theta$}$ of $\theta$ (and its
accuracy), which is based on data set $D$, from the $L$ MAP estimators
$\hat{}\mbox{\boldmath$\theta$}_{\ell},\ell=1,\ldots,L$ (and their accuracy),
which are solely inferred from the data subsets $D_{\ell}$. This is the topic
for the next subsection.

###  2.3 Bayesian Federated Inference

If the distribution of the data comes from an exponential family, BFI is
straightforward; the MAP estimator based on the data $D$ can be obtained
directly from summary statistics in the $L$ subsets. Once models become more
complex approximating the posterior density $p(\mbox{\boldmath$\theta$}|D)$ by
a multivariate Gaussian density may be needed. The Bernstein-von Mises theorem
states that under certain regularity conditions and for sufficiently large
sample size $n$, the posterior distribution can be approximated well by a
multivariate Gaussian distribution centered at the maximum likelihood
estimator of $\theta$ (see e.g. van der Vaart, 1998). In practice, the sample
sizes in the centers are finite and may even be small compared to the number
of covariates (especially in the case of rare diseases). We, therefore, center
around the MAP estimator. We approximate the logarithm of the posterior
density (given the combined data set $D$)

$\displaystyle\log p(\mbox{\boldmath$\theta$}|D)\;=\;\log
p(\mbox{\boldmath$\theta$})+\sum_{\ell=1}^{L}\sum_{i=1}^{n_{\ell}}\log
p(y_{\ell i}|\mbox{\boldmath$\theta$})-\log Z(D)$

by a Taylor expansion up to a quadratic order in $\theta$ around its MAP
parameter estimator $\hat{}\mbox{\boldmath$\theta$}$ which is based on the
combined data set $D$:

$\displaystyle\log p(\mbox{\boldmath$\theta$}|D)\;=\;\log
p(\hat{}\mbox{\boldmath$\theta$}|D)-\tfrac{1}{2}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})^{t}\hat{}\mbox{\boldmath$A$}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})+O_{p}((\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})^{3}),$
(4)

where $\hat{}\mbox{\boldmath$A$}$ is equal to minus the curvature matrix of
$\log p(\mbox{\boldmath$\theta$}|D)$, i.e. the matrix of minus second order
derivatives with respect to $\theta$, evaluated at
$\hat{}\mbox{\boldmath$\theta$}$. Note that there is no linear term in the
Taylor expansion (4), because the gradient of $\log
p(\mbox{\boldmath$\theta$}|D)$ at $\hat{}\mbox{\boldmath$\theta$}$ equals zero
by definition (as $\hat{}\mbox{\boldmath$\theta$}$ is the MAP estimator).
Hence, under certain regularity assumptions, for $\theta$ in a neighbourhood
of $\hat{}\mbox{\boldmath$\theta$}$,

$\displaystyle
p(\mbox{\boldmath$\theta$}|D)\;\propto\;\exp(-\tfrac{1}{2}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})^{t}\hat{}\mbox{\boldmath$A$}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$}))+O_{p}((\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})^{3}),$

and thus near $\hat{\mbox{\boldmath$\theta$}}$ one can approximate the
posterior distribution by a Gaussian density:

$\displaystyle
p(\mbox{\boldmath$\theta$}|\mbox{\boldmath$D$})\;\approx\;\Big{(}\frac{{\rm
det}\;\hat{}\mbox{\boldmath$A$}}{(2\pi)^{d}}\Big{)}^{\frac{1}{2}}\exp\big{(}-\tfrac{1}{2}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})^{t}\hat{}\mbox{\boldmath$A$}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})\big{)},$
(5)

where $d={\rm dim}(\mbox{\boldmath$\theta$})$ and ${\rm
det}\hat{}\mbox{\boldmath$A$}$ denotes the determinant of the matrix
$\hat{}\mbox{\boldmath$A$}$. In a similar way, each of the $L$ posterior
distributions for the subsets $\ell=1\ldots L$ is approximately Gaussian
around $\hat{}\mbox{\boldmath$\theta$}_{\ell}$ (the MAP estimator based on
data set $D_{\ell}$) and with covariance matrix equal to the inverse of
$\hat{}\mbox{\boldmath$A$}_{\ell}$. By substituting expression (5) for
$p(\mbox{\boldmath$\theta$}|D)$, and the equivalent expressions for
$p_{\ell}(\mbox{\boldmath$\theta$}|D_{\ell}),\ell=1,\ldots,L$, into relation
(3), we obtain:

$\displaystyle\Big{(}\frac{{\rm
det}\;\hat{}\mbox{\boldmath$A$}}{(2\pi)^{d}}\Big{)}^{\frac{1}{2}}\exp(-\tfrac{1}{2}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})^{t}\hat{}\mbox{\boldmath$A$}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$}))$
$\displaystyle=$
$\displaystyle\frac{\prod_{\ell=1}^{L}Z_{\ell}(D_{\ell})}{Z(D)}~{}\frac{p(\mbox{\boldmath$\theta$})}{\prod_{\ell=1}^{L}p_{\ell}(\mbox{\boldmath$\theta$})}\prod_{\ell=1}^{L}\Big{(}\frac{{\rm
det}\;\hat{}\mbox{\boldmath$A$}_{\ell}}{(2\pi)^{d}}\Big{)}^{\frac{1}{2}}\exp\big{(}-\tfrac{1}{2}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$}_{\ell})^{t}\hat{}\mbox{\boldmath$A$}_{\ell}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$}_{\ell})\big{)}.$

If the prior densities $\mbox{\boldmath$\theta$}\to
p(\mbox{\boldmath$\theta$})$ and $\mbox{\boldmath$\theta$}\to
p_{\ell}(\mbox{\boldmath$\theta$})$ are chosen to be Gaussian, with mean zero
and covariance matrices $\mbox{\boldmath$\Lambda$}^{-1}$ and
$\mbox{\boldmath$\Lambda$}_{\ell}^{-1}$, i.e.
$p(\mbox{\boldmath$\theta$})=(\det\mbox{\boldmath$\Lambda$}/(2\pi)^{d})^{\frac{1}{2}}\exp(-\frac{1}{2}\mbox{\boldmath$\theta$}^{t}\mbox{\boldmath$\Lambda$}\mbox{\boldmath$\theta$})$
and
$p_{\ell}(\mbox{\boldmath$\theta$})=(\det\mbox{\boldmath$\Lambda$}_{\ell}/(2\pi)^{d})^{\frac{1}{2}}\exp(-\frac{1}{2}\mbox{\boldmath$\theta$}^{t}\mbox{\boldmath$\Lambda$}_{\ell}\mbox{\boldmath$\theta$})$,
the above equation becomes

$\displaystyle\Big{(}\frac{{\rm
det}\;\hat{}\mbox{\boldmath$A$}}{(2\pi)^{d}}\Big{)}^{\frac{1}{2}}\exp(-\tfrac{1}{2}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})^{t}\hat{}\mbox{\boldmath$A$}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$}))\;=\;\;\frac{\prod_{\ell=1}^{L}Z_{\ell}(D_{\ell})}{(2\pi)^{\frac{1}{2}d}Z(D)}\Bigg{(}\frac{(\det\mbox{\boldmath$\Lambda$})\prod_{\ell=1}^{L}({\rm
det}\;\hat{}\mbox{\boldmath$A$}_{\ell})}{\prod_{\ell=1}^{L}(\det\mbox{\boldmath$\Lambda$}_{\ell})}\Bigg{)}^{\\!\frac{1}{2}}\;\;$
(7)
$\displaystyle\times\exp\Big{(}-\tfrac{1}{2}\mbox{\boldmath$\theta$}^{t}\big{(}\mbox{\boldmath$\Lambda$}-\sum_{\ell=1}^{L}\mbox{\boldmath$\Lambda$}_{\ell}\big{)}\mbox{\boldmath$\theta$}\Big{)}\prod_{\ell=1}^{L}\exp\big{(}-\tfrac{1}{2}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$}_{\ell})^{t}\hat{}\mbox{\boldmath$A$}_{\ell}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$}_{\ell})\big{)}.$

Taking the logarithm on both sides yields, when viewed as a function of
$\theta$:

$\displaystyle(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})^{t}\hat{}\mbox{\boldmath$A$}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$})$
$\displaystyle=$
$\displaystyle\mbox{\boldmath$\theta$}^{t}\Big{(}\mbox{\boldmath$\Lambda$}-\sum_{\ell=1}^{L}\mbox{\boldmath$\Lambda$}_{\ell}\Big{)}\mbox{\boldmath$\theta$}+\sum_{\ell=1}^{L}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$}_{\ell})^{t}\hat{}\mbox{\boldmath$A$}_{\ell}(\mbox{\boldmath$\theta$}-\hat{}\mbox{\boldmath$\theta$}_{\ell})+B,$
(8)

with $B$ representing a term that does not depend on $\theta$. Both sides are
quadratic functions of $\theta$, and must have the same expansion coefficients
for the linear and quadratic terms in $\theta$. This implies that

$\displaystyle{\rm quadratic~{}terms:}$
$\displaystyle\sum_{j,k=1}^{d}\theta_{j}\theta_{k}\Big{(}\hat{A}_{jk}-\Lambda_{jk}+\sum_{\ell=1}^{L}\Lambda_{\ell,jk}-\sum_{\ell=1}^{L}\hat{A}_{\ell,jk}\Big{)}=0~{}~{}~{}~{}~{}\forall\mbox{\boldmath$\theta$}$
$\displaystyle{\rm linear~{}terms:}$
$\displaystyle\sum_{j=1}^{d}\theta_{j}\Big{(}\hat{}\mbox{\boldmath$A$}\hat{}\mbox{\boldmath$\theta$}-\sum_{\ell=1}^{L}\hat{}\mbox{\boldmath$A$}_{\ell}\hat{}\mbox{\boldmath$\theta$}_{\ell}\Big{)}_{j}=0~{}~{}~{}~{}~{}\forall\mbox{\boldmath$\theta$}.$

The solutions of these equations are:

$\displaystyle\hat{\mbox{\boldmath$A$}}=\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}_{\ell}+\mbox{\boldmath$\Lambda$}-\sum_{\ell=1}^{L}\mbox{\boldmath$\Lambda$}_{\ell},~{}~{}~{}~{}~{}~{}~{}~{}~{}\hat{\mbox{\boldmath$\theta$}}=\Big{(}\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}_{\ell}+\mbox{\boldmath$\Lambda$}-\sum_{\ell=1}^{L}\mbox{\boldmath$\Lambda$}_{\ell}\Big{)}^{-1}\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}_{\ell}\hat{\mbox{\boldmath$\theta$}}_{\ell}$
(9)

Again there is no need to do inference on the combined data set $D$ to find
the estimator $(\hat{}\mbox{\boldmath$\theta$},\hat{}\mbox{\boldmath$A$})$
from the combined set; we can compute approximations a posteriori from the
inference results on the subsets.

Even if one only wishes to do MAP inference, one would still need the matrix
$\hat{}\mbox{\boldmath$A$}$ to compute error bars or regions for the MAP
estimator $\hat{}\mbox{\boldmath$\theta$}$. More specifically, for the
$k^{th}$ element of $\theta$ its approximate $(1-2\alpha)100\%$ confidence
interval equals
$\hat{}\mbox{\boldmath$\theta$}_{k}\pm\xi_{\alpha}\;\sqrt{(\hat{}\mbox{\boldmath$A$}^{-1})_{k,k}},$
for $\xi_{\alpha}$ the upper $\alpha$-quantile of the standard Gaussian
distribution. Hypothesis testing is also straightforward by the asymptotic
normality of $\hat{}\mbox{\boldmath$\theta$}$.

The formulae in (9) do not say anything about the plausibility of the subsets
describing similar subpopulations. However, once we have computed the
estimators $(\hat{}\mbox{\boldmath$\theta$},\hat{}\mbox{\boldmath$A$})$ we
should find that $\hat{}\mbox{\boldmath$\theta$}$ is compatible with each
‘local’ estimator $\hat{}\mbox{\boldmath$\theta$}_{\ell}$, given the error
bars coded in the matrices $\hat{}\mbox{\boldmath$A$}$ and
$\hat{}\mbox{\boldmath$A$}_{\ell}$. One way to address this is to consider the
coordinate wise approximate $(1-2\alpha)100\%$ confidence intervals for the
difference between the true $\theta$-value in all subsets except subset $\ell$
and the true parameter value in subset $\ell$:

$\displaystyle(\hat{}\mbox{\boldmath$\theta$}_{-\ell}-\hat{}\mbox{\boldmath$\theta$}_{\ell})_{k}\;\pm\;\xi_{\alpha}\;\sqrt{(\hat{}\mbox{\boldmath$A$}^{-1}_{-\ell}+\hat{}\mbox{\boldmath$A$}_{\ell}^{-1})_{kk}},$

where subscript $-\ell$ indicates that the estimator excluded data from subset
$\ell$.

###  2.4 Nuisance parameters differ across centers

Suppose that the vector $\theta$ can be decomposed as
$\mbox{\boldmath$\theta$}^{t}=(\mbox{\boldmath$\theta$}_{a}^{t},\mbox{\boldmath$\theta$}_{b}^{t})$,
where $\mbox{\boldmath$\theta$}_{a}$ denotes the core parameters of interest
of dimension $d_{1}$, and $\mbox{\boldmath$\theta$}_{b}$ the nuisance
parameters of dimension $d_{2}$, with $d_{1}+d_{2}=d$. The parameter vector of
interest $\mbox{\boldmath$\theta$}_{a}$ is assumed to be equal in the $L$
different sub-populations, but $\mbox{\boldmath$\theta$}_{b}$ may vary across
the sub-populations. Let $\mbox{\boldmath$\theta$}_{b,\ell}$ be the vector of
nuisance parameters in population $\ell$. Then, the parameter vector in the
model for the combined data set is equal to
$\mbox{\boldmath$\theta$}^{t}=(\mbox{\boldmath$\theta$}_{a}^{t};\mbox{\boldmath$\theta$}_{b,1}^{t},\ldots,\mbox{\boldmath$\theta$}_{b,L}^{t})$,
and is of dimension $d^{\prime}=d_{1}+Ld_{2}$. Let
$\hat{\mbox{\boldmath$\theta$}}_{b,\ell}$ be the estimate of
$\mbox{\boldmath$\theta$}_{b,\ell}$ based solely on data in the $\ell^{th}$
subset, and $\hat{\mbox{\boldmath$\theta$}}^{\sim}_{b,\ell}$ be the estimate
based on the combined data set. We choose simple priors of the form
$p(\mbox{\boldmath$\theta$})=p(\mbox{\boldmath$\theta$}_{a})\prod_{\ell=1}^{L}p(\mbox{\boldmath$\theta$}_{b,\ell})$
(for the combined data set), and
$p(\mbox{\boldmath$\theta$}_{a},\mbox{\boldmath$\theta$}_{b,\ell})=p_{\ell}(\mbox{\boldmath$\theta$}_{a})p_{\ell}(\mbox{\boldmath$\theta$}_{b,\ell})$
in data subset $\ell$. We then have

$\displaystyle\log p(\mbox{\boldmath$\theta$}|D)\\!$ $\displaystyle=$
$\displaystyle\\!\\!\log
p(\mbox{\boldmath$\theta$}_{a})+\sum_{\ell=1}^{L}\sum_{i=1}^{n_{\ell}}\log
p(y_{\ell
i}|\mbox{\boldmath$\theta$}_{a},\mbox{\boldmath$\theta$}_{b,\ell})+\sum_{\ell=1}^{L}\log
p(\mbox{\boldmath$\theta$}_{b,\ell})-\log Z(D),$ $\displaystyle\log
p_{\ell}(\mbox{\boldmath$\theta$}_{a},\mbox{\boldmath$\theta$}_{b,\ell}|D_{\ell})\\!$
$\displaystyle=$ $\displaystyle\\!\\!\log
p_{\ell}(\mbox{\boldmath$\theta$}_{a})+\sum_{i=1}^{n_{\ell}}\log p(y_{\ell
i}|\mbox{\boldmath$\theta$}_{a},\mbox{\boldmath$\theta$}_{b,\ell})+\log
p_{\ell}(\mbox{\boldmath$\theta$}_{b,\ell})-\log Z_{\ell}(D_{\ell}),$

and hence, like in (3),

$\displaystyle\log p(\mbox{\boldmath$\theta$}|D)$ $\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\log
p_{\ell}(\mbox{\boldmath$\theta$}_{a},\mbox{\boldmath$\theta$}_{b,\ell}|D_{\ell})+\log\Big{(}\frac{p(\mbox{\boldmath$\theta$}_{a})}{\prod_{\ell=1}^{L}p_{\ell}(\mbox{\boldmath$\theta$}_{a})}\Big{)}$
(10)
$\displaystyle\quad\qquad+\;\log\Big{(}\frac{\prod_{\ell=1}^{L}p(\mbox{\boldmath$\theta$}_{b,\ell})}{\prod_{\ell=1}^{L}p_{\ell}(\mbox{\boldmath$\theta$}_{b,\ell})}\Big{)}-\log\Big{(}\frac{Z(D)}{\prod_{\ell=1}^{L}Z_{\ell}(D_{\ell})}\Big{)}.$

The quadratic approximations of the a posteriori densities are now

$\displaystyle\log p(\mbox{\boldmath$\theta$}|D)\;\approx\;\log
p(\hat{}\mbox{\boldmath$\theta$}|D)-\tfrac{1}{2}(\mbox{\boldmath$\theta$}_{a}-\hat{}\mbox{\boldmath$\theta$}_{a})^{t}\hat{}\mbox{\boldmath$A$}_{a}(\mbox{\boldmath$\theta$}_{a}-\hat{}\mbox{\boldmath$\theta$}_{a})$
$\displaystyle\hskip
28.45274pt-\;\tfrac{1}{2}\sum_{\ell=1}^{L}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{}\mbox{\boldmath$\theta$}^{\sim}_{b,\ell})^{t}\hat{\mbox{\boldmath$A$}}^{\sim}_{b,\ell}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{}\mbox{\boldmath$\theta$}^{\sim}_{b,\ell})-(\mbox{\boldmath$\theta$}_{a}-\hat{\mbox{\boldmath$\theta$}}_{a})^{t}\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}^{\sim}_{ab,\ell}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{\mbox{\boldmath$\theta$}}^{\sim}_{b,\ell}),$
$\displaystyle\log
p_{\ell}(\mbox{\boldmath$\theta$}_{a},\mbox{\boldmath$\theta$}_{b,\ell}|D_{\ell})\;\approx\;\log
p_{\ell}(\hat{\mbox{\boldmath$\theta$}}_{a},\hat{\mbox{\boldmath$\theta$}}_{b,\ell}|D_{\ell})-\tfrac{1}{2}(\mbox{\boldmath$\theta$}_{a}-\hat{\mbox{\boldmath$\theta$}}_{a,\ell})^{t}\hat{\mbox{\boldmath$A$}}_{a,\ell}(\mbox{\boldmath$\theta$}_{a}-\hat{\mbox{\boldmath$\theta$}}_{a,\ell})$
$\displaystyle\hskip
28.45274pt-\;\tfrac{1}{2}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{\mbox{\boldmath$\theta$}}_{b,\ell})^{t}\hat{}\mbox{\boldmath$A$}_{b,\ell}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{\mbox{\boldmath$\theta$}}_{b,\ell})-(\mbox{\boldmath$\theta$}_{a}-\hat{\mbox{\boldmath$\theta$}}_{a,\ell})^{t}\hat{\mbox{\boldmath$A$}}_{ab,\ell}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{\mbox{\boldmath$\theta$}}_{b,\ell})$

where $\hat{\mbox{\boldmath$\theta$}}_{b,\ell}^{\sim}$ is the MAP estimator of
$\mbox{\boldmath$\theta$}_{b}$ found via inference on the full set $D$, and
$\hat{\mbox{\boldmath$A$}}^{\sim}_{b,\ell}$ and
$\hat{\mbox{\boldmath$A$}}_{ab,\ell}^{\sim}$ denote the matrices of minus the
second derivatives of $\log p(\mbox{\boldmath$\theta$}|D)$ with respect to the
components of $\mbox{\boldmath$\theta$}_{b,\ell}$ and with respect to both
$\mbox{\boldmath$\theta$}_{a}$ and $\mbox{\boldmath$\theta$}_{b,\ell}$,
respectively. Insertion into (10), assuming the usual Gaussian zero-average
priors for $\mbox{\boldmath$\theta$}_{a}$, and Gaussian zero-average priors
for $\mbox{\boldmath$\theta$}_{b,\ell}$ with inverse covariance matrices
$\mbox{\boldmath$\Lambda$}_{b}$ in $p(\mbox{\boldmath$\theta$}_{b,\ell})$ and
$\mbox{\boldmath$\Lambda$}_{b,\ell}$ in
$p_{\ell}(\mbox{\boldmath$\theta$}_{b,\ell})$ gives

$\displaystyle(\mbox{\boldmath$\theta$}_{a}-\hat{}\mbox{\boldmath$\theta$}_{a})^{t}\hat{}\mbox{\boldmath$A$}_{a}(\mbox{\boldmath$\theta$}_{a}-\hat{}\mbox{\boldmath$\theta$}_{a})+\sum_{\ell=1}^{L}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{}\mbox{\boldmath$\theta$}^{\sim}_{b,\ell})^{t}\hat{\mbox{\boldmath$A$}}^{\sim}_{b,\ell}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{}\mbox{\boldmath$\theta$}^{\sim}_{b,\ell})$
$\displaystyle\hskip
170.71652pt+\;2(\mbox{\boldmath$\theta$}_{a}-\hat{\mbox{\boldmath$\theta$}}_{a})^{t}\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}^{\sim}_{ab,\ell}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{\mbox{\boldmath$\theta$}}^{\sim}_{b,\ell})$
$\displaystyle=\sum_{\ell=1}^{L}(\mbox{\boldmath$\theta$}_{a}-\hat{\mbox{\boldmath$\theta$}}_{a,\ell})^{t}\hat{\mbox{\boldmath$A$}}_{a,\ell}(\mbox{\boldmath$\theta$}_{a}-\hat{\mbox{\boldmath$\theta$}}_{a,\ell})+\mbox{\boldmath$\theta$}_{a}^{t}\Big{(}\mbox{\boldmath$\Lambda$}_{a}-\sum_{\ell=1}^{L}\mbox{\boldmath$\Lambda$}_{a,\ell}\Big{)}\mbox{\boldmath$\theta$}_{a}$
(11) $\displaystyle\hskip
28.45274pt+\;\sum_{\ell=1}^{L}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{\mbox{\boldmath$\theta$}}_{b,\ell})^{t}\hat{}\mbox{\boldmath$A$}_{b,\ell}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{\mbox{\boldmath$\theta$}}_{b,\ell})\;+\;\sum_{\ell=1}^{L}\mbox{\boldmath$\theta$}_{b,\ell}^{t}(\mbox{\boldmath$\Lambda$}_{b}-\mbox{\boldmath$\Lambda$}_{b,\ell})\mbox{\boldmath$\theta$}_{b,\ell}$
$\displaystyle\hskip
28.45274pt+\;2\sum_{\ell=1}^{L}(\mbox{\boldmath$\theta$}_{a}-\hat{\mbox{\boldmath$\theta$}}_{a,\ell})^{t}\hat{\mbox{\boldmath$A$}}_{ab,\ell}(\mbox{\boldmath$\theta$}_{b,\ell}-\hat{\mbox{\boldmath$\theta$}}_{b,\ell})+B$

with $B$ representing a term that does not depend on $\theta$. Both sides in
the above equation are quadratic functions of $\theta$, and this equation must
hold for all $\theta$. We must therefore equate the coefficients on either
side of the specific linear and quadratic terms. Identification of the
quadratic terms in $\theta$ leads to the following equations:

$\displaystyle\hat{\mbox{\boldmath$A$}}_{a}=\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}_{a,\ell}+\mbox{\boldmath$\Lambda$}_{a}-\sum_{\ell=1}^{L}\mbox{\boldmath$\Lambda$}_{a,\ell},~{}~{}~{}~{}~{}\hat{\mbox{\boldmath$A$}}_{b,\ell}^{\sim}=\hat{\mbox{\boldmath$A$}}_{b,\ell}+\mbox{\boldmath$\Lambda$}_{b}-\mbox{\boldmath$\Lambda$}_{b,\ell},~{}~{}~{}~{}~{}\hat{\mbox{\boldmath$A$}}_{ab,\ell}^{\sim}=\hat{\mbox{\boldmath$A$}}_{ab,\ell},$
(12)

whereas identification of the linear terms in $\theta$ gives the two equations

$\displaystyle\hat{}\mbox{\boldmath$A$}_{a}\hat{}\mbox{\boldmath$\theta$}_{a}+\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}^{\sim}_{ab,\ell}\hat{\mbox{\boldmath$\theta$}}^{\sim}_{b,\ell}$
$\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}_{a,\ell}\hat{\mbox{\boldmath$\theta$}}_{a,\ell}+\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}_{ab,\ell}\hat{\mbox{\boldmath$\theta$}}_{b,\ell},$
$\displaystyle\hat{\mbox{\boldmath$A$}}^{\sim}_{b,\ell}\hat{}\mbox{\boldmath$\theta$}^{\sim}_{b,\ell}+(\hat{\mbox{\boldmath$A$}}^{\sim}_{ab,\ell})^{t}\hat{\mbox{\boldmath$\theta$}}_{a}$
$\displaystyle=$
$\displaystyle\hat{}\mbox{\boldmath$A$}_{b,\ell}\hat{\mbox{\boldmath$\theta$}}_{b,\ell}+(\hat{\mbox{\boldmath$A$}}_{ab,\ell})^{t}\hat{\mbox{\boldmath$\theta$}}_{a,\ell}.$

These equations can again be solved, giving formulae for the estimators
$\hat{\mbox{\boldmath$\theta$}}_{a}$ and
$\hat{\mbox{\boldmath$\theta$}}_{b}^{\sim}$ that would be obtained based on
the combined data in terms of estimators in the separate centers. For the core
parameters one finds:

$\displaystyle\hat{\mbox{\boldmath$\theta$}}_{a}=\Big{(}\hat{}\mbox{\boldmath$A$}_{a}-\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}_{ab,\ell}(\hat{\mbox{\boldmath$A$}}^{\sim}_{b,\ell})^{-1}(\hat{\mbox{\boldmath$A$}}_{ab,\ell})^{t}\Big{)}^{\\!-1}\;\times\;$
$\displaystyle\sum_{\ell=1}^{L}\Big{[}\Big{(}\hat{\mbox{\boldmath$A$}}_{a,\ell}-\hat{\mbox{\boldmath$A$}}_{ab,\ell}(\hat{\mbox{\boldmath$A$}}^{\sim}_{b,\ell})^{-1}(\hat{\mbox{\boldmath$A$}}_{ab,\ell})^{t}\Big{)}\hat{\mbox{\boldmath$\theta$}}_{a,\ell}+\hat{\mbox{\boldmath$A$}}_{ab,\ell}\Big{(}1\\!\\!{\rm
I}-(\hat{\mbox{\boldmath$A$}}^{\sim}_{b,\ell})^{-1}\hat{}\mbox{\boldmath$A$}_{b,\ell}\Big{)}\hat{\mbox{\boldmath$\theta$}}_{b,\ell}\Big{]}~{}~{}$

with $1\\!\\!{\rm I}$ the $d_{2}\times d_{2}$ unit matrix and the expressions
for the matrices $\hat{\mbox{\boldmath$A$}}_{a}$ and
$\hat{\mbox{\boldmath$A$}}^{\sim}_{b,\ell}$ as given in (12). If required, the
full set estimators for the nuisance parameters can subsequently be calculated
from

$\displaystyle\hat{}\mbox{\boldmath$\theta$}^{\sim}_{b,\ell}=(\hat{\mbox{\boldmath$A$}}^{\sim}_{b,\ell})^{-1}\Big{[}\hat{}\mbox{\boldmath$A$}_{b,\ell}\hat{\mbox{\boldmath$\theta$}}_{b,\ell}+(\hat{\mbox{\boldmath$A$}}_{ab,\ell})^{t}(\hat{\mbox{\boldmath$\theta$}}_{a,\ell}-\hat{\mbox{\boldmath$\theta$}}_{a})\Big{]}.$

We conclude that also if nuisance parameters are different across centers,
there is no need to combine the data subsets; we can compute all relevant full
set estimators a posteriori from the estimators obtained from the subsets.

The above equations are seen to simplify if we choose identical nuisance
parameter priors throughout, i.e.
$\mbox{\boldmath$\Lambda$}_{b,\ell}=\mbox{\boldmath$\Lambda$}_{b}$ for all
$\ell$. In that case
$\hat{\mbox{\boldmath$A$}}^{\sim}_{b,\ell}=\hat{\mbox{\boldmath$A$}}_{b,\ell}$,
and hence

$\displaystyle\hat{\mbox{\boldmath$\theta$}}_{a}=\Big{(}\hat{}\mbox{\boldmath$A$}_{a}-\sum_{\ell=1}^{L}\hat{\mbox{\boldmath$A$}}_{ab,\ell}(\hat{\mbox{\boldmath$A$}}_{b,\ell})^{-1}(\hat{\mbox{\boldmath$A$}}_{ab,\ell})^{t}\Big{)}^{\\!-1}\\!\sum_{\ell=1}^{L}\Big{(}\hat{\mbox{\boldmath$A$}}_{a,\ell}-\hat{\mbox{\boldmath$A$}}_{ab,\ell}(\hat{\mbox{\boldmath$A$}}_{b,\ell})^{-1}(\hat{\mbox{\boldmath$A$}}_{ab,\ell})^{t}\Big{)}\hat{\mbox{\boldmath$\theta$}}_{a,\ell}.$

In Section 4 an example is given in which the data in three centers follow a
logistic regression model in which the center specific intercepts are the
nuisance parameters and the remaining regression parameters are the parameters
of interest. This results in an aggregated model that contains an intercept
for each of the three populations.

##  3 Bayesian regression and subset inference

###  3.1 Bayesian regression

The most relevant problem areas in medicine are those where each data sample
is a pair $(\mbox{\boldmath$X$},Y)$ of input vectors $X$ (the covariates) and
output values $Y$ (e.g. clinical outcome or treatment response). A
parametrized model relates $X$ to $Y$, and the available data are used to
infer the parameters for this model. In this section the BFI methodology
described in Section 2 is applied to parametric regression models, and
associated challenges such as heterogeneity across the populations are
addressed. Suppose that $Y|\mbox{\boldmath$X$}=\mbox{\boldmath$x$}$ and $X$
have densities $y\to p(y|\mbox{\boldmath$x$},\mbox{\boldmath$\theta$}_{1})$
and $\mbox{\boldmath$x$}\to
p(\mbox{\boldmath$x$}|\mbox{\boldmath$\theta$}_{2})$, respectively, so that
for
$\mbox{\boldmath$\theta$}^{t}=(\mbox{\boldmath$\theta$}_{1}^{t},\mbox{\boldmath$\theta$}_{2}^{t})$
with unknown parameter vectors
$\mbox{\boldmath$\theta$}_{1}\in\mathbb{R}^{d_{1}}$ and
$\mbox{\boldmath$\theta$}_{2}\in\mathbb{R}^{d_{1}}$ we seek to infer:
$(y,\mbox{\boldmath$x$})\to
p(y,\mbox{\boldmath$x$}|\mbox{\boldmath$\theta$})\;=\;p(y|\mbox{\boldmath$x$},\mbox{\boldmath$\theta$}_{1})p(\mbox{\boldmath$x$}|\mbox{\boldmath$\theta$}_{2}).$
An example is the class of generalized linear models (GLMs). In GLMs the
outcome variable $Y$ is related to the covariates $x$ only via the linear
predictor $\mbox{\boldmath$\beta$}^{t}\mbox{\boldmath$x$}$. By setting the
first covariate equal to 1, an intercept can be included trivially in the
model.

Let $D_{\ell}$ be the $\ell^{th}$ data subset, and $D$ be the union of the $L$
subsets:

$\displaystyle D_{\ell}=\\{(\mbox{\boldmath$x$}_{\ell 1},y_{\ell
1}),\ldots,(\mbox{\boldmath$x$}_{\ell n_{\ell}},y_{\ell
n_{\ell}})\\},~{}~{}~{}~{}~{}~{}~{}~{}n_{\ell}=|D_{\ell}|,~{}~{}~{}~{}~{}~{}~{}~{}D=\bigcup_{\ell=1}^{L}D_{\ell}.$

As in the previous section, we aim to express the outcome of inference on the
parameters in the combined set $D$ in terms of the outcomes of $L$ separate
inferences on the constituent sets $D_{\ell},\ell=1,\ldots,L$.

For simplicity we assume statistically independent
$\mbox{\boldmath$\theta$}_{1}$ and $\mbox{\boldmath$\theta$}_{2}$, i.e.
$p(\mbox{\boldmath$\theta$}_{1},\mbox{\boldmath$\theta$}_{2})=p(\mbox{\boldmath$\theta$}_{1})p(\mbox{\boldmath$\theta$}_{2})$
and
$p_{\ell}(\mbox{\boldmath$\theta$}_{1},\mbox{\boldmath$\theta$}_{2})=p_{\ell}(\mbox{\boldmath$\theta$}_{1})p_{\ell}(\mbox{\boldmath$\theta$}_{2})$
for all $\ell$. For the combined data set $D$ and for the subsets $D_{\ell}$
the log posteriori distribution is then given by, respectively,

$\displaystyle\log p(\mbox{\boldmath$\theta$}|D)\;=\;\log
p(\mbox{\boldmath$\theta$})+\sum_{\ell=1}^{L}\sum_{i=1}^{n_{\ell}}\log
p(y_{\ell i},\mbox{\boldmath$x$}_{\ell i}|\mbox{\boldmath$\theta$})-\log Z(D)$
$\displaystyle=$ $\displaystyle\log
p(\mbox{\boldmath$\theta$}_{1})+\sum_{\ell=1}^{L}\sum_{i=1}^{n_{\ell}}\log
p(y_{\ell i}|\mbox{\boldmath$x$}_{\ell i},\mbox{\boldmath$\theta$}_{1})+\log
p(\mbox{\boldmath$\theta$}_{2})+\sum_{\ell=1}^{L}\sum_{i=1}^{n_{\ell}}\log
p(\mbox{\boldmath$x$}_{\ell i}|\mbox{\boldmath$\theta$}_{2})-\log Z(D)$
$\displaystyle\log p_{\ell}(\mbox{\boldmath$\theta$}|D_{\ell})\;=\;\log
p_{\ell}(\mbox{\boldmath$\theta$})+\sum_{i=1}^{n_{\ell}}\log p(y_{\ell
i},\mbox{\boldmath$x$}_{\ell i}|\mbox{\boldmath$\theta$})-\log
Z_{\ell}(D_{\ell})$ $\displaystyle=$ $\displaystyle\log
p_{\ell}(\mbox{\boldmath$\theta$}_{1})+\sum_{i=1}^{n_{\ell}}\log p(y_{\ell
i}|\mbox{\boldmath$x$}_{\ell i},\mbox{\boldmath$\theta$}_{1})+\log
p_{\ell}(\mbox{\boldmath$\theta$}_{2})+\sum_{i=1}^{n_{\ell}}\log
p(\mbox{\boldmath$x$}_{\ell i}|\mbox{\boldmath$\theta$}_{2})-\log
Z_{\ell}(D_{\ell}).$

So, the log posterior densities are decomposed into terms that depend on
either $\mbox{\boldmath$\theta$}_{1}$, or on $\mbox{\boldmath$\theta$}_{2}$
(or neither), and maximization with respect to $\mbox{\boldmath$\theta$}_{1}$
and $\mbox{\boldmath$\theta$}_{2}$ to obtain their MAP estimators can be
performed independently. Similarly as in (3), the link between the log
posterior density for the combined data set and for the subsets takes the form

$\displaystyle\log p(\mbox{\boldmath$\theta$}|D)$ $\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\log
p_{\ell}(\mbox{\boldmath$\theta$}|D_{\ell})+\log\Big{(}\frac{p(\mbox{\boldmath$\theta$}_{1})}{\prod_{\ell=1}^{L}p_{\ell}(\mbox{\boldmath$\theta$}_{1})}\Big{)}+\log\Big{(}\frac{p(\mbox{\boldmath$\theta$}_{2})}{\prod_{\ell=1}^{L}p_{\ell}(\mbox{\boldmath$\theta$}_{2})}\Big{)}-\log\Big{(}\frac{Z(D)}{\prod_{\ell=1}^{L}Z_{\ell}(D_{\ell})}\Big{)}.$

Again the log posterior densities $\log p(\mbox{\boldmath$\theta$}|D)$ and
$\log p_{\ell}(\mbox{\boldmath$\theta$}|D_{\ell})$ can for sufficiently large
$n$ be approximated by quadratic expansions for $\theta$ near the MAP
estimators $\hat{}\mbox{\boldmath$\theta$}$ and
$\hat{\mbox{\boldmath$\theta$}}_{\ell}$, leading to Gaussian approximations
for the posterior densities themselves. Since the log posterior densities can
both be decomposed in terms dependent on either $\mbox{\boldmath$\theta$}_{1}$
or $\mbox{\boldmath$\theta$}_{2}$ (there are no mixture terms), the matrices
$\hat{\mbox{\boldmath$A$}}$ and $\hat{\mbox{\boldmath$A$}}_{\ell}$ are
diagonal block matrices of the form

$\displaystyle\hat{\mbox{\boldmath$A$}}=\left(\\!\begin{array}[]{cc}\hat{\mbox{\boldmath$A$}}_{1}&{\bf
0}\\\ {\bf
0}&\hat{\mbox{\boldmath$A$}}_{2}\end{array}\\!\right),~{}~{}~{}~{}~{}~{}~{}~{}\hat{\mbox{\boldmath$A$}}_{\ell}=\left(\\!\begin{array}[]{cc}\hat{\mbox{\boldmath$A$}}_{1,\ell}&{\bf
0}\\\ {\bf 0}&\hat{\mbox{\boldmath$A$}}_{2,\ell}\end{array}\\!\right),$ (20)

in which the blocks
$\\{\hat{\mbox{\boldmath$A$}}_{1},\hat{\mbox{\boldmath$A$}}_{1,\ell}\\}$ and
$\\{\hat{\mbox{\boldmath$A$}}_{2},\hat{\mbox{\boldmath$A$}}_{2,\ell}\\}$
represent the minus curvature matrices for $\mbox{\boldmath$\theta$}_{1}$ and
$\mbox{\boldmath$\theta$}_{2}$, respectively. The quadratic log-likelihood
approximation for $p(\mbox{\boldmath$\theta$}|D)$ is

$\displaystyle\log p(\mbox{\boldmath$\theta$}|D)\\!$
$\displaystyle\\!\approx\\!$ $\displaystyle\\!\log
p(\hat{}\mbox{\boldmath$\theta$}|D)-\tfrac{1}{2}(\mbox{\boldmath$\theta$}_{1}-\hat{}\mbox{\boldmath$\theta$}_{1})^{t}\hat{}\mbox{\boldmath$A$}_{1}(\mbox{\boldmath$\theta$}_{1}-\hat{}\mbox{\boldmath$\theta$}_{1})-\tfrac{1}{2}(\mbox{\boldmath$\theta$}_{2}-\hat{}\mbox{\boldmath$\theta$}_{2})^{t}\hat{}\mbox{\boldmath$A$}_{2}(\mbox{\boldmath$\theta$}_{2}-\hat{}\mbox{\boldmath$\theta$}_{2}),$

and a similar expression holds for
$p_{\ell}(\mbox{\boldmath$\theta$}|D_{\ell})$. If all prior parameter
densities are chosen to be Gaussian with mean zero and inverse covariance
matrices $\mbox{\boldmath$\Lambda$}_{1}$ and
$\mbox{\boldmath$\Lambda$}_{1,\ell}$ (for $\mbox{\boldmath$\theta$}_{1}$) and
$\mbox{\boldmath$\Lambda$}_{2}$ and $\mbox{\boldmath$\Lambda$}_{2,\ell}$ (for
$\mbox{\boldmath$\theta$}_{2}$), then insertion of these quadratic
approximations into (3.1) gives, for $k=1,2$

$\displaystyle(\mbox{\boldmath$\theta$}_{k}-\hat{}\mbox{\boldmath$\theta$}_{k})^{t}\hat{}\mbox{\boldmath$A$}_{k}(\mbox{\boldmath$\theta$}_{k}-\hat{}\mbox{\boldmath$\theta$}_{k})=\mbox{\boldmath$\theta$}_{k}^{t}\Big{(}\mbox{\boldmath$\Lambda$}_{k}-\sum_{\ell=1}^{L}\mbox{\boldmath$\Lambda$}_{k,\ell}\Big{)}\mbox{\boldmath$\theta$}_{k}+\sum_{\ell=1}^{L}(\mbox{\boldmath$\theta$}_{k}-\hat{\mbox{\boldmath$\theta$}}_{k,\ell})^{t}\hat{\mbox{\boldmath$A$}}_{k,\ell}(\mbox{\boldmath$\theta$}_{k}-\hat{\mbox{\boldmath$\theta$}}_{k,\ell})+B_{k},$

with $B_{1}$ and $B_{2}$ not dependent on $\theta$. Expressions for
$\hat{}\mbox{\boldmath$\theta$}_{1},\hat{}\mbox{\boldmath$\theta$}_{2},\hat{}\mbox{\boldmath$A$}_{1}$
and $\hat{}\mbox{\boldmath$A$}_{2}$ in terms of their subset analogons can be
now be extracted similarly as in (9), giving

$\displaystyle\hat{}\mbox{\boldmath$A$}_{1}=\sum_{\ell=1}^{L}\hat{}\mbox{\boldmath$A$}_{1,\ell}+\mbox{\boldmath$\Lambda$}_{1}-\sum_{\ell=1}^{L}\mbox{\boldmath$\Lambda$}_{1,\ell},~{}~{}~{}~{}~{}~{}~{}~{}~{}\hat{}\mbox{\boldmath$\theta$}_{1}=\hat{}\mbox{\boldmath$A$}_{1}^{-1}\sum_{\ell=1}^{L}\hat{}\mbox{\boldmath$A$}_{1,\ell}\hat{}\mbox{\boldmath$\theta$}_{1,\ell},$
(21)
$\displaystyle\hat{}\mbox{\boldmath$A$}_{2}=\sum_{\ell=1}^{L}\hat{}\mbox{\boldmath$A$}_{2,\ell}+\mbox{\boldmath$\Lambda$}_{2}-\sum_{\ell=1}^{L}\mbox{\boldmath$\Lambda$}_{2,\ell},~{}~{}~{}~{}~{}~{}~{}~{}~{}\hat{}\mbox{\boldmath$\theta$}_{2}=\hat{}\mbox{\boldmath$A$}_{2}^{-1}\sum_{\ell=1}^{L}\hat{}\mbox{\boldmath$A$}_{2,\ell}\hat{}\mbox{\boldmath$\theta$}_{2,\ell}.$
(22)

From these formulae it is also clear that the distribution of the covariates
does not affect the estimator
$(\hat{}\mbox{\boldmath$\theta$}_{1},\hat{A}_{1})$.

It is well known that there is a direct link between Bayesian inference for
regression models and penalized regression. If one chooses a Gaussian prior
with a diagonal inverse covariance matrix $\Lambda$ with all diagonal elements
equal to $\lambda$, the MAP estimator for $\theta$ is identical to the ridge
penalized ML estimator for $\theta$ with regularizer $\lambda$ (Wu et al.
2018). The higher the value of $\lambda$, the stronger the penalization, or
the smaller the variance of the prior distribution. When the sample size in
one or more local subsets is small compared to the number of covariates in the
model, the model in that particular center may be overfitted. Increasing the
value of the regularizer $\lambda_{\ell}$ (i.e. of the inverse variance of the
prior in subset $\ell$) can prevent overfitting in the local models and
increase the inference accuracy. This is illustrated in Section 4. Note that
in most ridge regression analyses only the regression (or association)
parameters are penalized, not the intercept.

###  3.2 Heterogeneity across populations

So far we have assumed that the distributions of the covariates and the
conditional distributions of the outcome given the covariates are the same in
all centers. In practice this means that we assumed that the subsets were
sampled from statistically identical populations. This assumption might be
violated, for instance if the different centers represent different
nationalities or hospitals. In that case, the $L$ subsets are sampled from
structurally different populations, and the distributions of the covariates
and/or of the (conditional) outcomes may well differ across centers. If we
assume that the regression parameters $\mbox{\boldmath$\theta$}_{1}$ do not
vary across centers but that the parameters of the covariate distributions,
$\mbox{\boldmath$\theta$}_{2}$, may do so, and choose simple priors of the
form
$p(\mbox{\boldmath$\theta$})=p(\mbox{\boldmath$\theta$}_{1})\prod_{\ell=1}^{L}p(\mbox{\boldmath$\theta$}_{2,\ell})$,
then the calculations as before yields equation (21) derived earlier for all
estimators related to $\mbox{\boldmath$\theta$}_{1}$. Estimation of
$\mbox{\boldmath$\theta$}_{1}$ is hence not affected by having subset-specific
parameters $\mbox{\boldmath$\theta$}_{2}$. Equation (22), in contrast, is now
replaced by

$\displaystyle\hat{\mbox{\boldmath$A$}}^{\sim}_{2,\ell}=\hat{\mbox{\boldmath$A$}}_{2,\ell}+\mbox{\boldmath$\Lambda$}_{2}-\mbox{\boldmath$\Lambda$}_{2,\ell},~{}~{}~{}~{}~{}~{}~{}~{}~{}\hat{\mbox{\boldmath$\theta$}}^{\sim}_{2,\ell}=\big{(}\hat{\mbox{\boldmath$A$}}_{2,\ell}+\mbox{\boldmath$\Lambda$}_{2}-\mbox{\boldmath$\Lambda$}_{2,\ell}\big{)}^{-1}\hat{\mbox{\boldmath$A$}}_{2,\ell}\hat{\mbox{\boldmath$\theta$}}_{2,\ell},$

where $\hat{\mbox{\boldmath$\theta$}}_{2,\ell}^{\sim}$ is the MAP estimator of
$\mbox{\boldmath$\theta$}_{2}$ found via inference on the full set $D$, and
$\hat{\mbox{\boldmath$A$}}^{\sim}_{2,\ell}$ denotes the matrix of minus second
derivatives of $\log p(\mbox{\boldmath$\theta$}|D)$ with respect to
$\mbox{\boldmath$\theta$}_{2,\ell}$.

###  3.3 Sets of covariates differ across centers

Not only the covariate distributions, but also the actual set of available
covariates may differ across centers. For instance, it might be that in some
hospitals the body mass index of patients is measured, whereas in other
hospitals it is not. If this covariate is not of interest for the study, it
can be excluded from the analysis. Including the covariate in the regression
models in some centers and not in others might give non-interpretable results,
as covariates may be correlated.

In general, if an individual covariate is missing, imputation methods can be
applied. Common methods are single and multiple regression imputation (van
Buuren, 2018): a value for a missing covariate is predicted by a regression
model that was fitted based on the observed covariates and outcome values. In
our present context a similar strategy can be applied. A prediction model for
the unobserved covariate is fitted in the centers in which the covariate is
measured, combined by means of the BFI strategy in order to obtain one
regression model (according to the formulae above) and subsequently used to
predict and impute the covariate values in the centers where the covariate was
not measured. After imputation, the final BFI model can then be constructed as
before. Of course, this strategy only works if the simultaneous distributions
of the covariates are identical across centers.

##  4 Simulation Studies and Application

###  4.1 Study aims

In this section we describe several simulation studies aimed at quantifying
the performance of the BFI methodology. More specifically, the agreement
between the inference results of the BFI strategy and of those obtained when
actually combining all the data of $L$ subsets are studied.

Regression models are employed for multiple purposes, e.g. for predicting
outcomes of new patients or subjects, or for studying association between a
factor and the outcome. Therefore, we focus on quantifying the agreement of
parameter estimates and on outcome predictions in the settings where there is
homogeneity or heterogeneity across local centers, low sample sizes and
unmeasured covariates. We use existing real life data for our simulation
studies.

###  4.2 Real life data

Our data set consists of data of trauma patients from different hospitals. The
outcome variable of interest is mortality (binary) and the covariates are age,
sex, the Injury Severity Score (ISS, ranging from 1 (low) to 75 (high)), the
Glasgow Coma Scale (GCS, which expresses the level of consciousness, ranging
from 3 (low) to 15 (high)). The data originate from multiple hospitals which
can be categorised in three groups as: peripheral hospital without a neuro-
surgical unit, peripheral hospital with a neuro-surgical unit, and academic
medical center. They were collected between October 1984 and October 1985, in
12 hospitals from one of the three categories above. The original aim in
collecting these data was to compare hospitals from the three categories, in
terms of differences in rates of management error (Draaisma et al., 1989). In
the data in their present form (that have since been made available for
educational purposes) there is no longer any reference to individual
hospitals; there is only a variable indicating in which of the three hospital
categories the data were recorded. In the absence of more detailed information
we assume, for simplicity and in the terminology of the present study, that
there are only three data subsets or centers, which correspond to three
categories. Because the data in each subset come from a distinct hospital
type, the patient populations and characteristics of the patient traumas may
be different, i.e. the distributions of the covariates may vary across
subsets. Nevertheless, the regression parameters of the corresponding
covariates may be the same for the three subsets. Some data characteristics
are given in Table 1.

Data subsets | number | mortality | age | sex | ISS | GCS
---|---|---|---|---|---|---
| $n_{\ell}$ | % | median | % females | median | median
peripheral no NSU | 49 | 43 | 42 | 22 | 41 | 11
peripheral with NSU | 106 | 40 | 34 | 24 | 33 | 14
academic hospitals | 216 | 22 | 35 | 30 | 29 | 11
combined data | 371 | 30 | 36 | 27 | 30 | 12
Table 1: Data characteristics in the three subsets (here indicating different
hospital types) and in the combined data set. NSU stands for ‘neuro-surgical
unit’.

We combined the data, normalized every covariates (subtracted the sample mean
and divided by the sample standard deviation of the full data set), and fitted
a Bayesian logistic regression model in the combined data set. The prior
distribution for the vector of regression parameters was chosen to be Gaussian
with mean zero and a diagonal covariance matrix with the diagonal elements
equal to $10$ (the diagonal elements of
$\mbox{\boldmath$\Lambda$}=\mbox{\boldmath$\Lambda$}_{\ell}$ equal
$\lambda=0.1$). The MAP-estimated regression parameters for the combined set,
together with their estimated standard errors, were
$\hat{\beta}_{0}=-1.70\;(0.22)$ (intercept), $\hat{\beta}_{1}=-0.15\;(0.18)$
(sex), $\hat{\beta}_{2}=1.36\;(0.20)$ (age), $\hat{\beta}_{3}=0.55\;(0.19)$
(ISS), and $\hat{\beta}_{4}=-1.98\;(0.24)$ (GCS).

###  4.3 Homogeneity across population

Simulation procedure
In this subsection we study the performance of the BFI methodology for the
case where the data in the subsets come from populations with identical
characteristics, i.e. with identical covariate distributions and identical
true parameter values (including those of the regression parameters). Since
this property does not seem to hold in our data set (see Table 1), we first
randomly assign each patient to one of the three subsets, keeping the sample
sizes for the three subsets fixed. After this randomization we can indeed
assume that the three patient groups have the same covariate distributions and
regression parameters. We then perform regression analyses at subset level and
at full data set level, and quantify the agreement between the results
obtained on the full combined data set versus those obtained via the BFI
protocol for integrating the outcomes of the subset regressions. This
procedure is repeated multiple times. Below we describe the quantities used
for measuring agreement. For the prediction comparisons some extra steps are
taken to ensure that estimation and prediction are not based on the same data.
In every cycle and every subset, 10% of the patients are randomly left out of
the estimation procedure. After obtaining the estimates for the regression
coefficients with the BFI strategy, mortality probabilities are predicted for
the 10% patients that were left out (based on their covariates and estimated
regression parameters). More details are given below.

Quantities to measure agreement
The outcome variable is binary (mortality) and, therefore, a logistic
regression model is fitted. The parameter of interest is the vector of
regression parameters $\beta$ (including the intercept). For simplicity of
interpretation of the regression parameters and the choice of the prior
distribution, we first normalize every covariate by subtracting from each
covariate value in the data set the sample mean of the covariate values in the
combined data set and dividing by the sample standard deviation (computed in
the combined data set). This is done for each covariate. Note that this sample
mean and standard deviation in the combined data set can be computed without
combining all local data, because they are functions of the first and second
sample moments in every local subset only (and of the local sample size).
Therefore, the local data themselves do not need to be combined and
standardization is also possible in practice.

Let $(\hat{}\mbox{\boldmath$\beta$}_{b},\hat{}\mbox{\boldmath$A$}_{b})$ and
$(\hat{}\mbox{\boldmath$\beta$}_{c},\hat{}\mbox{\boldmath$A$}_{c})$ denote the
estimates of $(\mbox{\boldmath$\beta$},\mbox{\boldmath$A$})$ based on the BFI
strategy and found after combining all data (the subscript $b$ stands for BFI,
$c$ for combined data). To quantify the agreement between the regression
estimates, the mean squared error (MSE) is computed for every regression
coefficient with the formula:

$\displaystyle
MSE_{\beta_{k}}=\frac{1}{M}\sum_{m=1}^{M}\Big{(}\hat{\beta}_{b,k}^{(m)}-\hat{\beta}_{c,k}\Big{)}^{2}.$
(23)

Here $k=1,\ldots,d$ labels the $d$ components of the vector $\beta$, and the
superscript $m$ labels the $M$ independent randomization cycles of the
patients over the three subsets. Note that the values $\hat{\beta}_{c,k}$ are
the same in every cycle (i.e. is independent of $m$), since this estimate is
based on the combined data. A small value of $MSE_{\beta_{k}}$ means that
there is hardly any loss when computing the MAP estimates with the BFI method
(i.e. without combining all the data). Furthermore, the MSE is calculated for
the square root of every diagonal element of minus the inverse of the
curvature matrix (representing estimators of the expected errors in the
regression coefficients):

$\displaystyle
MSE_{A_{k}}=\frac{1}{M}\sum_{m=1}^{M}\Big{(}\sqrt{((\hat{}\mbox{\boldmath$A$}_{b}^{(m)})^{-1})_{kk}}-\sqrt{(\hat{}\mbox{\boldmath$A$}_{c}^{-1})_{kk}}\Big{)}^{2}.$
(24)

For testing the accuracy of BFI’s patient outcome predictions we use similar
quantities. For the 10% individuals that were left out from the estimation
part, the prediction probabilities are estimated. In each cycle $m$, for an
individual $i$ from center $\ell$, these probabilities are computed as

$\displaystyle\hat{p}_{b,\ell,i}^{(m)}=\frac{\exp((\hat{}\mbox{\boldmath$\beta$}^{(m)}_{b})^{t}\mbox{\boldmath$x$}_{\ell,i})}{1+\exp((\hat{}\mbox{\boldmath$\beta$}^{(m)}_{b})^{t}\mbox{\boldmath$x$}_{\ell,i})},~{}~{}~{}~{}~{}~{}~{}~{}\hat{p}_{c,\ell,i}^{(m)}=\frac{\exp((\hat{}\mbox{\boldmath$\beta$}^{(m)}_{c})^{t}\mbox{\boldmath$x$}_{\ell,i})}{1+\exp((\hat{}\mbox{\boldmath$\beta$}^{(m)}_{c})^{t}\mbox{\boldmath$x$}_{\ell,i})}$
(25)

The estimators $\hat{\mbox{\boldmath$\beta$}}^{(m)}_{b}$ and
$\hat{\mbox{\boldmath$\beta$}}^{(m)}_{c}$ are now both obtained from 90% of
the patients (a selection that is different for each cycle $m$, so both depend
on $m$). To compare the results from the two regression routes (regression
after data integration versus BFI recombination of subset results), we compute

$\displaystyle MSE_{p}=\frac{1}{37M}\sum_{m=1}^{M}\sum_{\ell=1}^{L}\sum_{i\in
S_{\ell}^{(m)}}\Big{(}\hat{p}_{b,\ell,i}^{(m)}-\hat{p}_{c,\ell,i}^{(m)}\Big{)}^{2},$
(26)

where $S_{\ell}^{(m)}$ is the set of indices of individuals selected for
prediction in subset $\ell$ in cycle $m$. The value of 37 is equal to 10% of
the total sample size, and hence gives the number of individuals selected for
prediction in every cycle.

Results: Parameter estimation
The prior distribution for the vector of regression parameters $\beta$ is
chosen to be multivariate Gaussian with mean zero and a diagonal covariance
matrices
$\mbox{\boldmath$\Lambda$}^{-1}=\mbox{\boldmath$\Lambda$}_{\ell}^{-1}$ of
which every element equals $\lambda$. In every subset $D_{\ell}$ we fit a
logistic regression model with the same prior (i.e.
$\mbox{\boldmath$\Lambda$}_{\ell}=\mbox{\boldmath$\Lambda$}$ for all $\ell$),
and we aggregate the results by applying the formulae in (9). The
$MSE_{\beta_{k}}$ are in the upper part of Table 2.

|  | intercept | sex | age | ISS | GCS
---|---|---|---|---|---|---
$\lambda=0.1$ | $\hat{\beta}_{c}~{}~{}(SE)$ | -1.70 (0.22) | -0.15 (0.18) | 1.36 (0.20) | 0.55 (0.19) | -1.98 (0.24)
| $MSE_{\beta_{k}}$ |  1.11$\times 10^{-2}$ |  0.06$\times 10^{-2}$ |  0.890$\times 10^{-2}$ |  0.10$\times 10^{-2}$ |  2.00$\times 10^{-2}$
| $MSE_{A_{k}}$ |  0.13$\times 10^{-4}$ |  0.25$\times 10^{-4}$ |  0.13$\times 10^{-4}$ |  0.11$\times 10^{-4}$ |  0.25$\times 10^{-4}$
$\lambda=0.01$ | $\hat{\beta}_{c}~{}~{}(SE)$ | -1.72 (0.22) | -0.15 (0.18) | 1.37 (0.21) | 0.55 (0.19) | -2.00 (0.24)
| $MSE_{\beta_{k}}$ |  1.57$\times 10^{-2}$ |  0.10$\times 10^{-2}$ |  1.26$\times 10^{-2}$ |  0.17$\times 10^{-2}$ |  2.77$\times 10^{-2}$
| $MSE_{A_{k}}$ |  0.43 $\times 10^{-4}$ |  0.49$\times 10^{-4}$ |  0.29$\times 10^{-4}$ |  0.20$\times 10^{-4}$ |  0.28$\times 10^{-4}$
$\lambda_{1}=1,\lambda_{2}=0.1$ | $MSE_{\beta_{k}}$ | $1.64\times 10^{-2}$ | $0.05\times 10^{-2}$ | $1.34\times 10^{-2}$ | $0.11\times 10^{-2}$ | $2.75\times 10^{-2}$
$\lambda_{3}=\lambda=0.01$ | $MSE_{A_{k}}$ | $1.57\times 10^{-4}$ | $0.12\times 10^{-4}$ | $1.16\times 10^{-4}$ | $0.28\times 10^{-4}$ | $2.85\times 10^{-4}$
$\lambda=1$ | $\hat{\beta}_{c}~{}~{}(SE)$ | -1.58 (0.20) | -0.13 (0.17) | 1.25 (0.19) | 0.54 (0.18) | -1.84 (0.21)
| $MSE_{\beta_{k}}$ |  3.01$\times 10^{-2}$ |  0.04$\times 10^{-2}$ |  2.69$\times 10^{-2}$ |  0.19$\times 10^{-2}$ |  5.32 $\times 10^{-2}$
| $MSE_{A_{k}}$ |  15.87$\times 10^{-4}$ |  2.79$\times 10^{-4}$ |  10.91$\times 10^{-4}$ |  4.02$\times 10^{-4}$ |  20.37$\times 10^{-4}$
$\lambda=0.1$ | $MSE_{\beta_{k}}$ |  7.42 $\times 10^{-2}$ |  0.20 $\times 10^{-2}$ |  6.54 $\times 10^{-2}$ |  0.50 $\times 10^{-2}$ |  14.50 $\times 10^{-2}$
| $MSE_{A_{k}}$ |  0.35 $\times 10^{-4}$ |  2.23 $\times 10^{-4}$ |  0.35$\times 10^{-4}$ |  0.41 $\times 10^{-4}$ |  1.25$\times 10^{-4}$
$\lambda=0.01$ | $MSE_{\beta_{k}}$ |  18.34$\times 10^{-2}$ |  0.45$\times 10^{-2}$ |  14.20$\times 10^{-2}$ |  1.10$\times 10^{-2}$ |  31.02$\times 10^{-2}$
| $MSE_{A_{k}}$ |  6.37$\times 10^{-4}$ |  9.67$\times 10^{-4}$ |  3.93$\times 10^{-4}$ |  3.46$\times 10^{-4}$ |  2.01$\times 10^{-4}$
Table 2:  Results of the simulation study 1) for homogeneous populations (rows
1 and 2) as described in Subsection 4.3; 2) for different priors in the
subsets (row 3) with $\lambda_{i}$ the value in subset $i$ and $\lambda$ in
the combined data set (Subsection 4.3); 3) with small subsets (rows 4-6) as
described in Subsection 4.4. In all cases $M=1000$.

In this table we see that the estimates based on the BFI strategy are
accurate. As expected, for higher values of $\lambda$ the estimates of the
regression coefficients are closer to zero, and the estimated standard
deviations are smaller. Furthermore, the $MSE$ values are also closer to zero
for higher values of $\lambda$. However, better agreement does not necessarily
mean that the estimates are closer to the true values.

Results: Prediction
The BFI-estimated logistic regression model is used to predict mortality for
new patients. To quantify its prediction performance we carried out simulation
studies as described above: in each of the $M=1000$ simulation cycles 10% of
the patients were left out from estimation, and used instead to predict the
mortality probability. The quantity $MSE_{p}$ in (26) was calculated,
resulting in the values $MSE_{p}=0.24\times 10^{-3}$ for $\lambda=0.01$, and
$MSE_{p}=0.16\times 10^{-3}$ for $\lambda=0.1$. In addition, in Figure 2(a) we
show a scatter plot with the predicted probabilities based on the BFI model
(vertical axis) and the combined data model (horizontal axis), for
$\lambda=0.1$. This plot shows excellent agreement. As expected for the
present model definitions, the data points in Figure 2(a) are not distributed
homogeneously along the diagonal, since Gaussian distributed uncertainty in
parameter estimates is mapped nonlinearly to probabilities via the expressions
in (25). This effect is more clearly present when the local data sets are
smaller (next subsection), where the estimate uncertainties are larger.

(a)

(b)

(c)

Figure 2: Prediction probabilities based on the estimates
$\hat{}\mbox{\boldmath$\beta$}_{b}$ (BFI, vertical axis) versus those based on
$\hat{}\mbox{\boldmath$\beta$}_{c}$ (regression on combined data set,
horizontal axis), for $\lambda=0.1$. Perfect agreement corresponds to all
points being located on the diagonal. (a) Homogeneous populations. (b)
Heterogenous populations. (c) Missing covariate.

Results: parameter estimation with different priors
We also considered the situation in which the priors for the vector of
regression parameters $\beta$ in the subsets are different. We took
multivariate Gaussian distributions with mean zero again, but with inverse
covariance matrices with on the diagonal $\lambda_{1}=1,\lambda_{2}=0.1$ and
$\lambda_{3}=0.01$ for the three subsets and $\lambda=0.01$ for the combined
data set (so the lower the sample size the stronger the regularisation). The
$MSE_{\beta_{k}}$ are in row 3 of Table 2. The estimates based on the BFI
strategy are accurate. We also performed a simulation study for mortality
prediction. The results are very similar to those in Figure 2(a).

###  4.4 Small subsets

If the sample sizes of the subsets are small compared to the number of
covariates in the model, there is a risk of overfitting. Our present data set
consists of three data subsets, ranging in size from 49 to 216 patients.
Especially the subset with 49 patients is small, and overfitting is here
addressed by taking values for $\lambda$ that are not (too) close to zero. In
this subsection we consider a more extreme situation in which multiple subsets
are small. We create such data subsets from our original hospital data by
randomly dividing all patients over eight subsets of size 40 each, plus one
additional subset of 51 patients. The combined data set $D$ remains the same,
hence the same is true for all estimates of regression parameters based on
this combined data set. We now perform the simulation study described in
Subsection 4.3 to study the effect of modest subset sample sizes on the
estimation and prediction agreement. We do this for different values of
$\lambda$. The values of $MSE_{\beta_{k}}$ are given in Table 2. For the
patient outcome predictions we found the following MSE values:
$MSE_{p}=3.04\times 10^{-3}$ for $\lambda=0.01$, $MSE_{p}=1.08\times 10^{-3}$
for $\lambda=0.1$ and $MSE_{p}=0.43\times 10^{-3}$ for $\lambda=1$. Upon
increasing the value of $\lambda$, the estimates of the regression parameters
are shrunk towards zero and agreement increases. This can be seen by the
decreasing values of $MSE_{p}$, and also from the scatter plots in Figure 3.
For $\lambda=0.01$ (left plot in Figure 3) the agreement is weaker, and one
observes what appears to be overfitting (rotation of the data cloud relative
to the diagonal). For larger values of $\lambda$ the agreement between BFI
recombination and regression results on the full data set improves, as
expected. As mentioned before, better agreement does not necessarily mean that
the estimates are closer to the true values; we just conclude that the BFI
protocol reliably predicts the results of regression on the integrated data
set.

We repeated the simulation study for $\lambda=0.01$ in the combined data set
and $\lambda=0.1$ or $\lambda=1$ in the subsets. The results are very similar
as before.

Figure 3: Prediction probabilities based on the estimates
$\hat{}\mbox{\boldmath$\beta$}_{b}$ (BFI protocol, vertical axis) versus those
based on $\hat{}\mbox{\boldmath$\beta$}_{c}$ (regression on combined data set,
horizontal axis), in the setting with small sample sizes in nine subsets. Here
$n_{\ell}=40$ for $\ell=1\ldots 8$, and $n_{9}=51$. From left to right:
$\lambda=0.01$, $\lambda=0.1$, and $\lambda=1$.

###  4.5 Heterogeneity across populations

Different covariate distributions and effects
In this subsection we no longer use randomizations, but return to the original
data to study agreement between BFI and direct regression on the combined set
in terms of parameter estimates and outcome predictions, in the setting where
(as is the case in our original data) the populations from the different
hospital types are expected to differ. The combined set priors and subset
priors were chosen to be identical, Gaussian with zero average and
$\mbox{\boldmath$\Lambda$}=\mbox{\boldmath$\Lambda$}_{\ell}$ for all $\ell$,
with $\Lambda$ a diagonal matrix with $\lambda$ on the diagonal. The parameter
estimates obtained via the BFI protocol and those found via regression on the
combined data set are given in the Table 3. From this table we conclude that,
even if there are differences in the covariate effects between the populations
where the subsets where sampled from, the BFI-estimated regression parameters
are still very similar to those that would have been found upon first
combining all data in a single set. For the comparisons of the predictions, as
before in every cycle 10% of the individuals are selected randomly for
prediction, and not used for estimation. We quantify the agreement in patient
outcome probabilities again by (26). For $\lambda=0.1$ we found
$MSE_{p}=0.36\times 10^{-3}$. A scatter plot for the comparison of the
prediction probabilities is shown in Figure 2(b). It can be seen that the
agreement between the predictions is high.

| intercept | sex | age | ISS | GCS
---|---|---|---|---|---
subset 1 | -2.82 (0.94) | -1.74 (0.90) | 1.54 (0.56) | 1.90 (0.70) | -2.06 (0.65)
subset 2 | -0.97 (0.33) | -0.34 (0.31) | 1.90 (0.45) | 0.52 (0.36) | -1.82 (0.40)
subset 3 | -2.22 (0.36) | 0.07 (0.25) | 1.27 (0.29) | 0.30 (0.26) | -2.36 (0.40)
BFI ($\mbox{\boldmath$\beta$}_{b}$) | -1.51 (0.22) | -0.12 (0.18) | 1.23 (0.21) | 0.51 (0.19) | -1.74 (0.24)
Combined ($\mbox{\boldmath$\beta$}_{c}$) | -1.70 (0.22) | -0.15 (0.18) | 1.36 (0.20) | 0.55 (0.19) | -1.98 (0.24)
Table 3: Estimates of $\beta$ in each subset, based on the BFI strategy, and
obtained after first combining the data. Here $\lambda=0.1$. The subsets
represent: peripheral hospital without NSU (subset 1), peripheral hospital
with NSU (subset 2) and academic medical centers (subset 3).

Different prevalence
In Table 1 we see that the mortalities in the three subsets (representing
hospital types) differ. This may not be explained completely by the covariates
in the model. The prevalence of mortality is reflected in the intercepts of
the logistic regression models for the three subsets. To address this observed
mortality difference we allow the intercept to vary across hospital types, and
we aggregate the estimates from the different subsets as explained in
Subsection 2.4. This yields a single logistic regression model that includes a
hospital type specific intercept, but no general intercept. For comparison, we
also fitted a logistic regression model with hospital type intercepts based on
the combined data (instead of a model with a general intercept and a reference
hospital type). The combined set priors and subset priors were chosen to be
identical, Gaussian with zero mean and
$\mbox{\boldmath$\Lambda$}=\mbox{\boldmath$\Lambda$}_{\ell}$ for all $\ell$,
and $\Lambda$ a diagonal matrix with $\lambda$ on the diagonal. For
$\lambda=0.1$ and $\lambda=1$ the results are shown in Table 4. We see very
satisfactory agreement between the estimates. We also see, as expected, that
the estimates $\hat{}\mbox{\boldmath$\beta$}_{c}$ are typically shrunk to zero
for $\lambda=1$ more than for $\lambda=0.1$.

$\lambda$ |  | subset 1 | subset 2 | subset 3 | sex | age | ISS | GCS
---|---|---|---|---|---|---|---|---
$\lambda=0.1$ |  BFI ($\mbox{\boldmath$\beta$}_{b}$) | -1.45 | -1.07 | -1.94 | -0.12 | 1.30 | 0.53 | -1.82
|  Combined ($\mbox{\boldmath$\beta$}_{c}$) | -1.61 | -1.15 | -2.09 | -0.16 | 1.40 | 0.57 | -1.96
$\lambda=1$ |  BFI ($\mbox{\boldmath$\beta$}_{b}$) | -1.17 | -0.94 | -1.79 | -0.13 | 1.20 | 0.54 | -1.69
|  Combined ($\mbox{\boldmath$\beta$}_{c}$) | -1.34 | -1.10 | -1.95 | -0.14 | 1.28 | 0.55 | -1.82
Table 4: Estimates of the regression coefficients computed with the BFI
strategy, and those obtained after first combining the data, when taking into
account the differences in death prevalence in the three hospital types by
including subset specific intercepts.

###  4.6 Sets of covariates differ across centers

In this subsection we study the performance of the BFI methodology in the
situation in which one of the covariates is missing in one of the subsets, and
is imputed by the data from the other subsets. We do this for the setting in
which we assume that the population is homogeneous, so we first randomly
assign the patients to the three subsets. Next, we delete all patient’s
covariate values for ”GCS” in subset 1, and estimate a linear regression model
with the dependent variable ”GCS” and the independent variables ”mortality”,
”sex”, ”age”, and ”ISS”, based on the data from subsets 2 and 3 with the BFI
strategy with identical priors with $\lambda=0.1$. This estimated model is
used to predict the GCS-scores in subset 1. Next, the BFI protocol with
identical priors (with $\lambda=0.1$) can be applied as discussed before. This
is repeated $M=1000$ times. The outcome predictions of the BFI protocol and of
regression performed on the combined data set are visualized in Figure 2(c).

##  5 Discussion

In this paper we have discussed a one-shot BFI strategy for estimating
parameters in parametric models, where the data in multiple centers can not be
combined for analysis. This is for instance of interest in a medical setting,
where the data sets in individual medical centers or hospitals are often
small, and can not be combined because of privacy legislation. A major
advantage of the proposed BFI strategy is that, in contrast to most FL
strategies, not only parameter estimates but also the uncertainties of these
estimates are computed systematically. These uncertainties reflect the
differences and uncertainties in the individual centers. In medicine,
complementing parameter estimates with error bars is essential in order to
interpret correctly the estimation results.

In this study we have performed numerical simulations in multiple settings,
with real medical data in which the patient outcomes were binary, and with
realistic sample sizes. They reveal a very good agreement between the
parameter estimates and patient outcome predictions obtained with the BFI
procedure, compared to those found after first combining all data in a single
integrated set. We conclude that for the given data and the chosen regression
model (logistic) the inference results based on the combined data set can be
computed reliably a posteriori from the inference results on the local
centers; hardly any information is lost by not being able to combine the data
sets.

In the proposed BFI strategy we approximated the posterior parameter
distribution around the MAP estimator $\hat{}\mbox{\boldmath$\theta$}$ by a
Gaussian distribution with mean $\hat{}\mbox{\boldmath$\theta$}$ and a
covariance matrix equal to minus the inverse curvature matrix in
$\hat{}\mbox{\boldmath$\theta$}$. For large sample sizes, the MAP estimator
will approach the maximum likelihood estimator and minus the inverse curvature
matrix will approach the Fisher information matrix (Bernstein-von Mises
theorem (van der Vaart, 1989)). If the sample sizes in the centers are small
(compared to the number of covariates) these identification may no longer
hold. In particular, if the sample size is close to or smaller than the number
of covariates, the ML estimator will be inaccurate or not even defined,
whereas the MAP estimator is still well defined by the presence of the ridge
penalty in the likelihood (for a Gaussian prior).

While in most of our simulations we have for simplicity used the same inverse
covariance matrices for the data subsets and the combined data set, the
derivations of sections 2 and 3 are in fact more general and allow these
matrices to be distinct. This freedom offered by the theory to vary the
covariances of the priors could for instance be employed to
regularize/penalize smaller data subsets more than larger ones as was shown in
subsections 4.3 and 4.4.

In this paper we have developed and applied the BFI methodology for parametric
(regression) models, and in particular for GLMs. This is not a limitation of
the BFI methodology, which is also applicable in more complex models and
settings, such as models for survival data that are subject to censoring, or
random effects models. In particular, the development of BFI protocols for
semi-parametric models, like the Cox proportional hazards model (Cox, 1972)
which is widely used in medical data analysis, will be our next project.

In conclusion, we have shown that by harnessing systematically and accurately
the cumulative power of multiple disjunct data sets without actually combining
these data sets, the BFI methodology can reduce significantly the sizes of
data sets required for extracting reliably statistically significant
predictive or prognostic patterns (if such patterns are present). While with
conventional methods only very strong regularities and associations can be
identified in small data sets, with the new approaches more subtle ones may
become detectable. It may thereby bring rare diseases, and in particular rare
cancers, closer to more common disease types in terms of our ability to
predict outcomes.

Acknowledgements
We thank J. Draaisma for making his data available.

Funding
This research was supported by an unrestricted grant of Stichting Hanarth
Fonds, The Netherlands.

## References

  * [1] Al-Shedivat, M., Gillenwater, J., Xing, E., Rostamizadeh, A., (2020), ”Federated learning via posterior averaging: A new perspective and practical algorithms,” In International Conference on Learning Representations.. url: https://openreview.net/forum?id=GFsU8a0sGB
  * [2] van Buuren, S., (2018), Flexible Imputation of Missing Data, Chapman and Hall/CRC
  * [3] Chen, H.Y., Chao, W.L., (2020), ”Fedbe: Making Bayesian model ensemble applicable to federated learning,” In International Conference on Learning Representations. url: https://arxiv.org/abs/2009.01974
  * [4] Cox, D.R., (1972), ”Regression Models and Life-Tables,” Journal of the Royal Statistical Society, Series B. 34 (2): 187–220.
  * [5] Draaisma, J.M.Th., de Haan, A.F.J., Goris, R.J.A., (1989),  ”Preventable Trauma Deaths in the Netherlands - A prospective Multicenter Study,”  The journal of Trauma 29(11): 1552-1557.
  * [6] Gafni, T., Shlezinger, N., Cohen, K., Eldar, Y.C., Poor, H.V., (2021), ”Federated Learning: A Signal Processing Perspective,” url: https://arxiv.org/abs/2103.17150
  * [7] Izmailov, P., Vikram, S., Hoffman, M.D., Wilson, A.G.G., (2021), ”What are bayesian neural network posteriors really like?” url: https://arxiv.org/abs/2104.14421.
  * [8] Kairouz, P., McMahan, H.B., Avent, B., Bellet, A., Bennis, M., Bhagoji, A.N., et al., (2021), ”Advances and open problems in federated learning,” Foundations and Trends in Machine Learning, 14(1–2):1–210.
  * [9] Li, T., Sahu, A.K., Zaheer, M., Sanjabi, M., Talwalkar, A., Smith, V., (2020), ”Federated optimization in heterogeneous networks,” Proceedings of Machine Learning and Systems, 2:429–450.
  * [10] Liu, L., Zheng, F., Chen, H., Qi, G-J., Huang, H., Shao, L., (2021),  ”A Bayesian Federated Learning Framework with Online Laplace Approximation” url: https://arxiv.org/abs/2102.01936
  * [11] Maddox, W.J., Izmailov, P., Garipov, T., Vetrov, D.P., Wilson, A.G., (2019), ”A simple baseline for Bayesian uncertainty in deep learning,” Advances in Neural Information Processing Systems, 32.
  * [12] McMahan, H.B.., Moore, E., Ramage, D., Hampson, S., Arcas, B.A., (2017), ”Communication-efficient learning of deep networks from decentralized data,” In Artificial Intelligence and Statistics, pp. 1273–1282.
  * [13] Rieke, N., Hancox, J., Li, W., Milletari, F., Roth, H.R., Albarqouni, S., et al., (2020), ”The future of digital health with federated learning,” NPJ Digital Medicine, 3, 119.
  * [14] Sheller, M.J., Edwards, B., Reina, A.G., Martin, J., Pati, S., Kotrotsou, A., et al., (2020), ”Federated learning in medicine: facilitating multi-institutional collaborations without sharing patient data,” Scientific Reports Nature research, 10:12598
  * [15] van der Vaart, A.W., (1998),  Asymptotic Statistics. Cambridge University Press. ISBN 0-521-49603-9.
  * [16] Zhang, C., Bütepage, J., Kjellström, H., Mandt, S., (2018), ”Advances in variational inference,” IEEE transactions on pattern analysis and machine intelligence, 41(8):2008–2026.
  * [17] Zhang, R., Li, C., Zhang, J., Chen, C., Wilson, A.G.. (2019), ”Cyclical stochastic gradient MCMC for Bayesian deep learning,” In International Conference on Learning Representations
  * [18] Zhu, Z., Hong, J., Zhou, J.. (2021), ”Data-free knowledge distillation for heterogeneous federated learning,” In International Conference on Machine Learning, pp. 12878–12889.


# An overview of the quantitative causality analysis and causal graph
reconstruction based on a rigorous formalism of information flow

X. San Liang<EMAIL_ADDRESS>Fudan University, Shanghai 200438,
China  Shanghai Qi Zhi Institute
(Andrew C. Yao Institute for Artificial Intelligence),
Shanghai 200232, China

(September 24, 2021)

###### Abstract

Inference of causal relations from data now has become an important field in
artificial intelligence. During the past 16 years, causality analysis has been
developed independently in physics from first principles, and, moreover, in a
quantitative sense. This short note presents a brief summary of this line of
work, including part of the theory and several representative applications.

Quantitative causality, information flow, causal graph, correlation, dynamical
system, normalization, self loop

Presented at the 1st Int’l AIxIA Workshop on Causality, Causal-ITALY
Italian Conference on Artificial Intelligence, November 30, 2021

##  I Introduction

The recent rush in artificial intelligence has stimulated enormous interest in
causal inference, particularly after the connection of independent causal
mechanisms to semi-supervised learning, thanks to Schölkopf et al. (2012).
Historically causal inference has been formulated as a problem of statistical
testing; see Granger (1969) and Pearl (2009) for the classics. In parallel, it
has also been investigated as a physical problem. Researches along this line
include Schreiber (2000) or Palu$\rm\check{s}$ et al. (2001), and Liang and
Kleeman (2005). Particularly, the latter is the first one formulated on a
rigorous footing within the framework of dynamical systems, which yields an
explicit solution in closed form, allowing for quantifying and normalizing
with ease the causality between dynamical events. This short note presents a
summary of this line of work, including the theory and several representative
applications.

##  II A gentle stroll through part of the theory

Although causality has long been studied ever since Granger’s seminal work,
its “mathematization is a relatively recent development,” said Peters et al.
(2017). On the other hand, Liang (2016) argued that causality is actually “a
real physical notion that can be rigorously derived ab initio.” This line of
work begins with Liang and Kleeman (2005), where a discovery about the
information flow with two-dimensional deterministic systems was presented. A
comprehensive study with generic systems has been fulfilled recently, with
explicit formulas attained in closed forms; see Liang (2008) and Liang (2016).
These formulas have been validated with benchmark systems such as baker
transformation, Hénon map, Kaplan-Yorke map, Rössler system, to name a few.
They have also been applied to real world problems in the diverse disciplines
such as climate science, meteorology, turbulence, neuroscience, financial
economics, etc. The following is a brief introduction of the theory.

Consider a dynamical system, which may be either a discrete-time mapping or a
continuous-time system. For easy presentation, hereafter only the results with
the latter are shown. (The former requires the aid of the Frobenius-Perron
operator; see Liang (2016).) Let

$\displaystyle d{{\bf X}}={\bf F}({\bf X},t)dt+{\bf B}({\bf X},t)d{{\bf W}},$
(1)

be a $d$-dimensional continuous-time stochastic system for ${\bf
X}=(X_{1},...,X_{d})$, where ${\bf F}=(F_{1},...,F_{d})$ may be arbitrary
nonlinear functions of ${\bf X}$ and $t$, ${{\bf W}}$ is a vector of standard
Wiener processes, and ${\bf B}=(b_{ij})$ is the matrix of perturbation
amplitudes which may also be any functions of ${\bf X}$ and $t$. Assume that
${\bf F}$ and ${\bf B}$ are both differentiable with respect to ${\bf X}$ and
$t$. Now define the rate of information flow, or simply information flow, from
a component $X_{j}$ to another component $X_{i}$ as the contribution of
entropy from $X_{j}$ per unit time in increasing the marginal entropy of
$X_{i}$. We then have the following theorem (Liang, 2016):

###### Theorem 1

For the system (1), the rate of information flowing from $X_{j}$ to $X_{i}$
(in nats per unit time) is

$\displaystyle T_{j\to i}$ $\displaystyle=$
$\displaystyle-E\left[\frac{1}{\rho_{i}}\int_{{}^{d-2}}\frac{\partial(F_{i}\rho_{{\backslash
j}})}{\partial x_{i}}d{\bf x}_{{\backslash i}{\backslash
j}}\right]+\frac{1}{2}E\left[\frac{1}{\rho_{i}}\int_{{}^{d-2}}\frac{\partial^{2}(g_{ii}\rho_{{\backslash
j}})}{\partial x_{i}^{2}}d{\bf x}_{{\backslash i}{\backslash j}}\right],$ (2)
$\displaystyle=$
$\displaystyle-\int_{{}^{d}}\rho_{j|i}(x_{j}|x_{i})\frac{\partial(F_{i}\rho_{{\backslash
j}})}{\partial x_{i}}d{\bf
x}+\frac{1}{2}\int_{{}^{d}}\rho_{j|i}(x_{j}|x_{i})\frac{\partial^{2}(g_{ii}\rho_{\backslash
j})}{\partial x_{i}^{2}}d{\bf x},$ (3)

where $d{\bf x}_{{\backslash i}{\backslash j}}$ signifies
$dx_{1}...dx_{i-1}dx_{i+1}...dx_{j-1}dx_{j+1}...dx_{n}$, $E$ stands for
mathematical expectation, $g_{ii}=\sum_{k=1}^{n}b_{ik}b_{ik}$,
$\rho_{i}=\rho_{i}(x_{i})$ is the marginal probability density function (pdf)
of $X_{i}$, $\rho_{j|i}$ is the pdf of $X_{j}$ conditioned on $X_{i}$, and
$\rho_{{\backslash j}}=\int\rho({\bf x})dx_{j}$.

Equation (2) has a nice property, which forms the basis of the information
flow-based causality analysis (Liang, 2008). Without loss of generality, here
the subscripts 1 and 2 are substituted by $i$ and $j$; same below.

###### Theorem 2

If in (1) neither $F_{1}$ nor $g_{11}$ depends on $X_{2}$, then $T_{2\to
1}=0$.

The algorithm for the information flow-based causal inference is as follows:
If $T_{j\to i}=0$, then $X_{j}$ is not causal to $X_{i}$; otherwise it is
causal, and the absolute value measures the magnitude of the causality from
$X_{j}$ to $X_{i}$.

Another property regards the invariance upon coordinate transformation,
indicating that the obtained information flow is an intrinsic property in
nature (Liang, 2018).

###### Theorem 3

$T_{2\to 1}$ is invariant under arbitrary nonlinear transformation of
$(X_{3},X_{4},...,X_{d})$.

As shown in Liang (2021) (and other publications), this is very important in
causal graph reconstruction. It together with Theorem 2 makes it promising
toward a solution of the problem of latent confounding.

For linear systems, the formula in (2) can be simplified.

###### Theorem 4

In (1), if ${\bf F}({\bf X})={\bf f}+{\bf A}{\bf X}$, and the matrix ${\bf B}$
is constant, then

$\displaystyle T_{j\to i}=a_{ij}\frac{\sigma_{ij}}{\sigma_{ii}},$ (4)

where $a_{ij}$ is the $(i,j)^{th}$ entry of ${\bf A}$ and $\sigma{ij}$ the
population covariance between $X_{i}$ and $X_{j}$.

Notice if $X_{i}$ and $X_{j}$ are not correlated, then $\sigma_{ij}=0$, which
yields a zero causality: $T_{j\to i}=0$. But conversely it is not true. We
hence have the following corollary:

###### Corollary 1

In the linear sense, causation implies correlation, but not vice versa.

In an explicit expression, this corollary fixes the debate on causation vs.
correlation ever since George Berkeley (1710).

In the case with only time series (no dynamical system is given), we have the
following result (Liang, 2021):

###### Theorem 5

Given $d$ time series $X_{1},X_{2},...,X_{d}$, under the assumption of a
linear model with additive noise, the maximum likelihood estimator (mle) of
(2) for $T_{2\to 1}$ is

$\displaystyle\hat{T}_{2\to 1}=\frac{1}{\det{\bf
C}}\cdot\sum_{j=1}^{d}\Delta_{2j}C_{j,d1}\cdot\frac{C_{12}}{C_{11}},$ (5)

where $C_{ij}$ is the sample covariance between $X_{i}$ and $X_{j}$,
$\Delta_{ij}$ the cofactors of the matrix ${\bf C}=(C_{ij})$, and $C_{i,dj}$
the sample covariance between $X_{i}$ and a series derived from $X_{j}$ using
the Euler forward differencing scheme:
$\dot{X}_{j,n}=(X_{j,n+k}-X_{j,n})/(k{\Delta t})$, with $k\geq 1$ some
integer.

Eq. (5) is rather concise in form, involving only the common statistics, i.e.,
sample covariances. The transparent formula makes causality analysis, which
otherwise would be complicated, very easy and computationally efficient. Note,
however, that Eq. (5) cannot replace (2); it is just the maximum likelihood
estimator (mle) of the latter. One needs to test the statistical significance
before making a causal inference based on the estimator $\hat{T}_{2\to 1}$.

If what are given are not time series, but independent, identically
distributed (i.i.d.) panel data, it has been shown that $\hat{T}_{2\to 1}$ has
the same form as (5); see Rong and Liang (2021).

Besides the information flow between two components, say $X_{1}$ and $X_{2}$,
it is also possible to estimate the influence of one component, say $X_{1}$,
on itself. Following the convention since Liang and Kleeman (2005), write it
as $dH^{*}_{1}/dt$.

###### Theorem 6

Under a linear assumption, the mle of $dH^{*}_{1}/dt$ is

$\displaystyle\widehat{\left(\frac{dH^{*}_{1}}{dt}\right)}=\frac{1}{\det{\bf
C}}\cdot\sum_{j=1}^{d}\Delta_{1j}C_{j,d1}.$ (6)

This result, first obtained in Liang (2014), provides an efficient approach to
identifying self loops in a causal graph (cf. Hyttinen et al., 2012).

Statistical significance tests can be performed for the estimators. This is
done with the aid of a Fisher information matrix. See Liang (2014) and Liang
(2021) for details.

Causality in this sense can be normalized in order to reveal the relative
importance of a causal relation. See Liang (2015) for details. Note that
recently there has some similar developments along this line, e.g.,
R$\phi$ysland (2012), Moij et al. (2013), and Mogensen et al. (2018). We want
to mention that the above formalism appears 7-13 years earlier, and, to our
best knowledge, it is the first one studying causality within the framework of
dynamical systems.

##  III Some representative applications

The above rigorous formalism has been successfully put to application to many
real world problems such as tropical cyclone genesis prediction (Bai et al.,
2018), near-wall turbulence (Liang and Lozano-Durán, 2016), global climate
change (Stips et al., 2016), financial analysis (Lu et al., 2020; Liang,
2015), soil moisture-precipitation interaction (Hagan et al., 2018),
neuroscience problems (Hristopulos et al., 2019), El Niño (Liang et al.,
2021), to name a few. Among these we want to particularly mention the study by
Stips et al. (2016) on CO2 emission vs. global warming. They found that CO2
emission does drive the recent global warming during the past century, and the
causal relation is one-way. However, on a time scale of 1000 years or up, this
one-way causality is is completely reversed, becoming a causality from air
temperature to carbon dioxide. In other words, on the paleoclimate scale, it
is global warming that drives the CO2 emission! This remarkable result is
consistent with that inferred from the recent ice-core data from Antarctica.

Another interesting application (Liang, 2015) regards the relation between the
two corporations IBM and GE, using the time series of US stocks downloaded
from ${YAHOO!}\atop{\rm{finance}}$. Overall the causality between the two is
insignificant, but if a running time causality analysis is performed, there
appears a strong, almost one-way causality from IBM to GE in 70’s, starting
from 1971. This abrupt one-way causality out of blue reveals to us an old
story about “Seven Dwarfs and a Giant”: In 50-60’s, GE was believed to be the
biggest computer user outside the U.S. Federal Government; to avoid relying on
IBM the computer “Giant”, it together with six other companies (“Seven
Dwarfs”) began to build mainframes. GE itself of course is a giant, but in the
computer market, it is just a “Dwarf”. Since it could not beat IBM, in 1970,
it sold its computer division. As a result, starting from 1971, it had to rely
on IBM again. This is the story behind the jump in $T_{IBM\to GE}$ from 1970
to 1971. While the story has almost gone to oblivion, this finding, which is
solely based on a simple causality analysis of two time series with Eq. (5),
is really remarkable.

Figure 1: El Niño prediction has become a benchmark problem for the testing of
machine learning algorithms. The present wisdom is that El Niño may be
predicted at a lead time of 1-2 years. Shown here are 1000 predictions (pink)
of the El Niño Modoki index (EMI) as described in the text. Overlaid are the
observed EMI (blue), the mean of the realizations (cyan). The light shading
marks the period for validation, while the darker shading marks the prediction
period. (From Liang et al., 2021.)

The latest application on the prediction of El Niño Modoki, or Central-Pacific
type El Niño, is also a remarkable one. El Niño is a climate mode that has
been linked to many hazards globewide, e.g., flooding, drought, wild fires,
heat waves, etc. Its accurate forecasting is of great importance to many
sectors of our society such as agriculture, energy, hydrology, to name
several. With its societal importance and the elegant setting, El Niño
prediction has become a testbed for AI algorithms.

Currently the widsom for El Niño prediction is that it may be predictable at a
lead time of 1-2 years. However, there still exists much uncertainty; an
example is the 2014-16 “Monster El Niño,” almost all projections fell off the
mark. Among the El Niño varieties, it is believed that El Niño Modoki is
particularly difficult to predict.

A striking breakthrough has just been made. Liang et al. (2021) took advantage
of the quantitative nature of the above information flow-based causality
analysis, and identified a delayed causal pattern, i.e., the structure of the
information flow from the solar activity to the sea surface temperature, very
similar to the El Niño Modoki mode. They then conjectured that, based on the
series of sunspot numbers, El Niño Modoki should be predictable. This is
indeed the case, and, remarkably, the prediction can be at a lead time of as
long as 10 years or up! This remarkable progress, among others, is a result of
the rigorously formulated quantitative causality analysis.

## References

  * [1] Bai, C., R. Zhang, and Coauthors (2018). Forecasting the tropical cyclone genesis over the Northwest Pacific through identifying the causal factors in cyclone-climate interactions. J. Atmos. Ocean. Tech., 35, 247-259.
  * [2] Berkeley, G. A Treatise Concerning the Principles of Human Knowledge; Aaron Rhames: Dublin, Ireland, 1710.
  * [3] Hyttinen, A., F. Eberhardt, P.O. Hoyer (2012). Learning linear cyclic causal models with latent variables. Journal of Machine Learning Research, 13, 3387-3439.
  * [4] Granger, C.W.J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica, _37_ , 424–438.
  * [5] Hagan, D.F.T.; Wang, G.; Liang, X.S.; Dolman, H.A.J. (2019) A time-varying causality formalism based on the Liang-Kleeman information flow for analyzing directed interactions in nonstationary climate systems. J. Clim., _32_ , 7521–7537.
  * [6] Hristopulos, D.T.; Babul, A.; Babul, S.; Brucar, L.R.; Virji-Babul, N. (2019) Disrupted information flow in resting-state in adolescents with sports related concussion. Front. Hum. Neurosci., _13_ , 419. doi:10.3389/fnhum.2019.00419.
  * [7] Janzing, D., B. Schölkopf (2017). Detecting confounding in multivariate linear models via spectral analysis. J. Causal Inference, 6(1), DOI:10.1515/jci-2017-0013.
  * [8] Liang, X.S., R. Kleeman (2005). Information transfer between dynamical system components. Phys. Rev. Lett., 95(24), 244101.
  * [9] Liang, X.S. Information flow within stochastic dynamical systems. Phys. Rev. E 2008, _78_ , 031113.
  * [10] Liang, X.S. (2014). Unraveling the cause-effect relation between time series. Phys. Rev. E, 90, 052150.
  * [11] Liang, X.S. (2015). Normalizing the causality between time series. Phys. Rev. E 2015, _92_ , 022126.
  * [12] Liang, X.S. (2016). Information flow and causality as rigorous notions ab initio. Phys. Rev. E, 94, 052201.
  * [13] Liang, X.S., A. Lozano-Durán (2016). A preliminary study of the causal structure in fully developed near-wall turbulence. Proceedings of the Summer Program 2016, 233-242. Center for Turbulence Research, Stanford University, CA, USA.
  * [14] Liang, X.S. (2018). Information flow with respect to relative entropy. Chaos, _28_ , 075311.
  * [15] Liang, X.S. (2019). A study of the cross-scale causation and information flow in a stormy model mid-latitude atmosphere. Entropy, 21, 149.
  * [16] Liang, X.S. (2021). Normalized multivariate time series causality analysis and causal graph reconstruction. Entropy, 23, 679.
  * [17] Liang, X.S., F. Xu, Y. Rong, R. Zhang, X. Tang, and F. Zhang (2021). El Niño Modoki can be mostly predicted more than 10 years ahead of time. Sci. Rep., 11, 17860.
  * [18] Lu, X., K. Liu, et al. (2020). The break point-dependent causality between the cryptocurrency and emerging stock markets. Economic Computation and Economic Cybernetics Studies and Research, 54(4), 203-216.
  * [19] Mooij, J.M.; Janzing, D.; Heskes, T.; Schölkopf, B. (2013) From ordinary differential equations to structural causal models: The deterministic case. In Proceedings of the 29th Annual Conference on Uncertainty in Artificial Intelligence, Bellevue, WA, USA, 11–15, July, pp. 440–448.
  * [20] Mogensen, S.W.; Malinksky, D.; Hansen, N.R. (2018) Causal learning for partially observed stochastic dynamical systems. In Proceedings of the 34th Conference on Uncertainty in Artificial Intelligence (UAI), Monterey, CA, USA, 6–10 August 2018.
  * [21] Palu$\rm\check{s}$, M.; Komárek, V.; Hrn$\rm\check{c}$i$\rm\check{r}$, Z.; $\rm\check{S}$t$\rm\check{e}$rbová, K. (2001) Synchronization as adjustment of information rates: Detection from bivariate time series. Phys. Rev. E, _63_ , 046211.
  * [22] Pearl, J. Causality: Models, Reasoning, and Inference, 2nd ed; Cambridge University Press: New York, NY, USA, 2009.
  * [23] Hahs, D.W.; Pethel, S.D. (2011) Distinguishing anticipation from causality: Anticipatory bias in the estimation of information flow. Phys. Rev. Lett., _107_ , 12870.
  * [24] Peters, J., Janzing, D., Schölkopf, B. (2017) Elements of Causal Inference: Foundations and Learning Algorithms. The MIT Press, Cambridge, MA, USA.
  * [25] Rong, Y., X.S. Liang (2021). Panel data causal inference using a rigorous information flow analysis for homogeneous, independent and identically distributed datasets. IEEE Access, 9, 47266-47274.
  * [26] Reichenbach, H. (1956). The direction of time. University of California Press, Berkeley.
  * [27] R${/o}$ysland, K. (2012). Counterfactual analyses with graphical models based on local independence. The Annals of Statistics, 40(4):2162-2194.
  * [28] Schölkopf, B.; Janzing, D.; Peters, J.; Sgouritsa, E.; Zhang, K.; Mooij, J.M. On causal and anticausal learning. In Proceedings of the 29th International Conference on Machine Learning (ICML), Edinburgh, Scotland, June 26-July 1, 2012; pp.1255-1262.
  * [29] Schreiber, T. Measuring information transfer. Phys. Rev. Lett. 2000, _85_ , 461.
  * [30] Spirtes, P., K. Zhang (2016). Causal discovery and inference: concepts and recent methodological advances. Appl. Informatics, _3_ , 3.
  * [31] Stips, A.; Macias, D.; Coughlan, C.; Garcia-Gorriz, E.; Liang, X.S. On the causal structure between CO2 and global temperature. Sci. Rep. 2016, _6_ , 21691.