# Neuroimaging Meta Regression for Coordinate Based Meta Analysis Data with a
Spatial Model

Yifan Yu
Oxford Big Data Institute
University of Oxford
Oxford, UK
& Rosario Pintos Lobo
Department of Psychology
Florida International University
Miami, FL, USA
&Michael Cody Riedel
Department of Physics
Florida International University
Miami, FL, USA
& Katherine Bottenhorn
Department of Population and Public Health Sciences
University of Southern California
Los Angeles, CA, USA
& Angela R. Laird
Center for Imaging Science
Florida International University
Miami, FL, USA
& Thomas E. Nichols
Oxford Big Data Institute
University of Oxford
Oxford, UK
<EMAIL_ADDRESS>
To whom correspondence should be addressed.

###### Abstract

Coordinate-based meta-analysis combines evidence from a collection of
Neuroimaging studies to estimate brain activation. In such analyses, a key
practical challenge is to find a computationally efficient approach with good
statistical interpretability to model the locations of activation foci. In
this article, we propose a generative coordinate-based meta-regression (CBMR)
framework to approximate smooth activation intensity function and investigate
the effect of study-level covariates (e.g., year of publication, sample size).
We employ spline parameterization to model spatial structure of brain
activation and consider four stochastic models for modelling the random
variation in foci. To examine the validity of CBMR, we estimate brain
activation on $20$ meta-analytic datasets, conduct spatial homogeneity tests
at voxel level, and compare to results generated by existing kernel-based
approaches.

##  1 Introduction

Functional neuroimaging includes a number of techniques to image brain
activity, including Positron emission tomography (PET) and functional magnetic
resonance imaging (fMRI). Starting three decades ago, PET studies were used to
compare brain activity between rest and experimental conditions, producing
maps of “activation", images of statistics measuring the strength of the
experimental effect. Especially in the last two decades, the literature of
fMRI activations has grown rapidly, which motivates a need to integrate
findings, establish consistency and heterogeneity across independent but
related studies. However in both PET and fMRI studies, the validity is
challenged by common drawbacks of small sample size, high prevalence of false
positives (approximately $10-20\%$ of reported foci in publications are false
positives (Wager et al., 2007), as well as significant heterogeneity among
studies and unreliable inference due to their diversity in measurements and
types of analysis (Samartsidis et al., 2017b). Meta-analysis is an essential
tool to address these limitations and improve statistical power by pooling
evidence from multiple studies and providing insight into consistent results.

There are also applications of neuroimaging meta-analysis to resting-state
fMRI and structural analysis using voxel-based morphometry. Going forward we
will only reference fMRI, but note that the application extends to other types
of data. Meta-analysis is classified into two categories in neuroimaging
research: image-based meta-analysis (IBMA) which uses the 3D statistic maps of
original studies and coordinate-based meta-analysis (CBMA) which uses the
reported spatial coordinates of activation foci in a standard MNI or Talairach
space. Ideally, only IBMA would be used, as there is substantial information
loss by only using activation foci as compared to full statistics maps, and
further accuracy loss occurs when deactivation foci are ignored (Salimi-
Khorshidi et al., 2009). However, while it is now more common to share entire
statistical maps in published studies, historically this has not been the case
(Salimi-Khorshidi et al., 2009) and there exists large-scale coordinate
databases (e.g., BrainMap (Laird et al., 2005), Neurosynth (Yarkoni et al.,
2011)). Hence, CBMA is still the predominant approach for neuroimaging meta-
analysis.

To identify consensus in brain regions with consistent activation across
studies, researchers have developed a variety of CBMA methods, which are
either kernel-based or model-based. Among those kernel-based methods,
activation likelihood estimation (ALE, with a Gaussian kernel), multilevel
kernel density analysis (MKDA, with a uniform sphere) and signed differential
mapping (SDM, with a Gaussian kernel scaled by effect size) are commonly used
(Turkeltaub et al., 2002; Eickhoff et al., 2012; Wager et al., 2007; Radua et
al., 2012). None of the three methods is based on a formal statistical model,
however, all are able to obtain statistical inferences by reference to a null
hypothesis of total random arrangement of the foci (Samartsidis et al.,
2017b). Voxels with significant p-values are considered to be the regions of
consistent activation. Multiple testing corrected inferences are made by
controlling the family wise error rate using the null maximum distribution
(Westfall and Young, 1993) or the false discovery rate (FDR) (Benjamini-
Hochberg procedure). However, kernel-based methods lack interpretability as
mass-univariate approaches rather than explicit probabilistic models, they do
not generally allow group comparison, do not model the spatial dependence of
activation foci, as well as cannot accommodate study-level covariates to
conduct a meta-regression (Samartsidis et al., 2019).

Bayesian model-based methods address these limitations, and are categorised
into parametric spatial point process models (Kang et al., 2011; Montagna et
al., 2018; Samartsidis et al., 2019) and non-parametric Bayesian models (Yue
et al., 2012; Kang et al., 2014). They use explicit generative models for the
data with testable assumptions. Although they generally provide advances in
interpretability and accuracy over kernel-based methods, they are
computationally intensive approaches and generally require parallel computing
on GPUs (Samartsidis et al., 2019), and only some approaches can conduct meta-
regression to estimate the effect of study-level covariates. Further, it can
be more challenging for practitioners to interpret the spatial posterior
intensity functions and utilise spatial Bayesian models in practice.

In this work, we investigate classical frequentist models that explicitly
account for the spatial structure of the distribution of activation foci.
Specifically, we focus on developing a spatial model that takes the form of a
generalised linear model, where we make use of a spline parameterization to
induce a smooth response and model the entire image jointly, allow for image-
wise study-level regressors and consider different stochastic models to find
the most accurate but parsimonious fit. Although Poisson is the classic
distribution describing independent foci counts, we have previously found
evidence of over-dispersion (Samartsidis et al., 2017a), and thus we further
explore a Negative Binomial model, a Clustered Negative Binomial model and a
Quasi-Poisson model to allow excess variation in counts data.

Our work draws on the existing methods for CBMA, while introducing key
innovations. From the Bayesian work, we take the idea of explicit spatial
models; from the kernel methods, we take the idea of fixing the degree of
spatial smoothness. The contribution of this meta-regression model is both
methodological and practical, it provides a generative regression model that
estimates a smooth intensity function and can have study-level regressors.
Meanwhile, using a crucial memory-saving model factorisation, it is also a
computationally efficient alternative to the existing Bayesian spatial
regression models and provides an accurate estimation of the intensity
function. While our method is suitable for any CBMA data, we are particularly
motivated by studies of cognition. Cognition encompasses various mental
processes, including perception, intelligence, problem solving, social
interactions, and can be affected by substance use. We demonstrate this meta-
regression framework on previously published meta-analyses of $20$ cognitive
and psychological tasks, allowing generalised linear hypothesis testing on
spatial effect, as well as inference on the effect of study-level covariates.

In the reminder of this work, we present the proposed meta-regression
framework, discuss model factorisation and optimisation procedures, as well as
inferences on meta-regression outcomes via statistical tests in Section 2.
Then we explain experiment settings in Section 3 and explore different
variants of stochastic models on the $20$ meta-analytic datasets and describe
multiple goodness-of-fit statistics to identify the most accurate model,
establish valid FPR control via Monte Carlo simulation under the null
hypothesis of spatial homogeneity, followed by a comparison of homogeneity
test with kernel methods in Section 4. Finally, Section 5 summarises our
findings and potential extension of this meta-regression framework in the
future.

##  2 Methods

Generalised linear models are described in terms of their stochastic and
deterministic components. Our deterministic model has a regression structure
with a spatial component utilising a spline parameterization and study-level
covariate component. For the stochastic model, we consider multiple models
motivated by CBMA data characteristics. We then propose a model factorisation
approach to make our methods scalable, before outlining a general inference
framework.

###  2.1 Deterministic model

####  2.1.1 Generic regression structure

Assume there are $N$ voxels in each of $M$ studies, and then our CBMA data at
voxel $j$ for study $i$ is the voxelwise count of foci $Y_{ij}$, collected for
study $i$ as the N-vector
$Y_{i}=\left[Y_{i1},Y_{i2},\cdots,Y_{iN}\right]^{\top}$. We generate a spatial
design matrix $X(N\times P)$ with $P$ cubic B-spline bases (more details to
follow in Section 2.1.2) and construct study-level covariates matrix
$Z(M\times R)$ by extracting $R$ study-level covariates from each of $M$
studies. For the CBMA framework, the central object of interest is the
voxelwise intensity function for study $i$, which considers both effects of
smooth spatial bases and study-level covariates. In this setting, it is most
concise to write the model for study $i$ as

$\log(\mu_{i})=\log\left[\mathbb{E}(Y_{i})\right]=X\beta+(Z_{i}\gamma)\mathbf{1}_{N}$
(1)

where $\beta(P\times 1)$ and $\gamma(R\times 1)$ are regression coefficients
for spatial bases $X$ and study-level covariates $Z$ respectively, $Z_{i}$ is
the $i^{th}$ row of study-level regressors $Z$, $\mathbf{1}_{N}$ is a
$N$-vector of $1$’s; and the estimated intensity is captured via $\mu_{ij}$
for studies $i=1,...,M$ and voxels $j=1,...,N$, collected for study $i$ as the
N-vector $\mu_{i}=\left[\mu_{i1},\mu_{i2},\cdots\mu_{iN}\right]^{\top}$. This
model is identifiable as long as we ensure each covariate variable is mean
zero, letting $X$ capture the overall mean. The GLM for all voxels in all $M$
studies is then

$\log\left[\mathbb{E}(Y)\right]=(\mathbf{1}_{M}\otimes
X)\beta+(Z\otimes\mathbf{1}_{N})\gamma$ (2)

where $Y=[Y_{1},Y_{2},\cdots,Y_{M}]^{\top}$ is a $(M\times N)-$vector,
containing voxelwise foci count for all of $M$ studies, and $\otimes$ is the
Kronecker product. This formulation has millions of rows ($MN$) and the
spatial design matrix has billions of entries ($MN\times P$). In consideration
of implementation complexity and memory requirement, we will propose a
simplified reformulation of this GLM in Section 2.3.

####  2.1.2 Spline parameterization

Previous work on spatial point process modelling of CBMA data have treated
each study’s foci as realisation of a doubly-stochastic Poisson process, also
known as a Cox process. In some of that work, the log intensity function is
parameterised by superimposed Gaussian kernel basis functions (Montagna et
al., 2018), while in others, the log intensity is a Gaussian process
(Samartsidis et al., 2019). Here, we propose tensor product of cubic B-spline
basis for modelling the spatial intensity, as its smoothness, stability and
local support make it an ideal spatial basis for CBMA application. A
1-dimensional cubic B-spline is a piecewise polynomial of order $k=3$, where
pre-specified knots $T=(t_{0},t_{1},\cdots,t_{n})$ determines the
parameterization of basis functions as the intersections of polynomial
sections. The order $k$ B-spline basis functions $B_{ik}$ are defined by
recurrence relations, for $i=0,1,\cdots,n$,

$\begin{split}B_{i1}&=\begin{cases}1\;\text{for }t_{i}\leq t<t_{i+1}\\\
0\;\text{otherwise}\end{cases}\;\text{for}\;k=1\,,\\\
B_{i,k}(t)&=\frac{t-t_{i}}{t_{i+k-1}-t_{i}}B_{i,k-1}(t)+\frac{t_{i+k}-t}{t_{i+k}-t_{i+1}}B_{i+1,k-1}(t)\;\text{for}\;k>1\,.\end{split}$
(3)

The B-spline curve is a linear combination of the B-spline basis function
$B_{ik}$. For our $3D$ lattice, assume there are $v_{x}$ voxels along x
direction, the coefficients of $v_{x}$ voxels evaluated at each of $n_{x}$
B-spline bases construct a coefficient matrix $C_{x}$ (size $v_{x}\times
n_{x}$). Similarly, there exist another two coefficient matrices $C_{y}$ and
$C_{z}$ (size $v_{y}\times n_{y}$ and $v_{z}\times n_{z}$) along y and z
direction. The whole coefficient matrix $C$ of 3-dimensional B-spline bases is
constructed by taking tensor product of the $3$ coefficient matrices (see
Figure 1 for a 2D illustration),

$C=C_{x}\otimes C_{y}\otimes C_{z}$ (4)

The matrix of $C$ is $(v_{x}v_{y}v_{z})\times(n_{x}n_{y}n_{z})$, and is the
basis for the entire $3D$ volume, while the analysis is based on a brain mask
of $N$ voxels. The design matrix $X$ is obtained from $C$ after a three-step
process: First, rows corresponding to voxels outside the brain mask are
removed; then, columns are removed if they correspond to weakly supported
B-spline bases (a B-spline basis is regarded as “weakly supported" if its
maximum value of coefficients evaluated at each voxel is below $0.1$).
Finally, the rows are re-normalised (sum to $1$) to preserve the property of
“partition of unity" of spline bases.

Figure 1: Illustration of tensor product of $2D$ spline bases (with equal
knots spacing)

We define our cubic B-spline bases with equally spaced knots on $x,y,z$
dimension, and thus we parameterise the level of spatial smoothness by the
knot spacing. Larger knots spacing, smaller basis, and greater smoothness;
conversely, closer knots, larger basis, and greater ability to represent fine
details. Conceptually, more flexible parameterizations would allow arbitrary
knots locations, but with the consideration of minimising computational
complexity, we fix the design matrix $X$ based on pre-specified knots spacing
according to prior knowledge. While other spline applications use a dense
array of knots and then control smoothness with a roughness penalty, the
computational and memory requirements of our spatial model demand that we
judiciously select the coarsest spline spacing consistent with our
application.

###  2.2 Stochastic model

Assumed stochastic behaviours of CBMA foci data determine the form of
statistical likelihood we use. We consider a set of four stochastic models for
the distribution of foci counts at voxel level. All of our models take the
form of generalised linear models, where inhomogeneous intensity at each voxel
is captured by the spline bases and any study-level covariate (as per Equation
2). We fit our model either by maximising log-likelihood function iteratively
via L-BFGS algorithm or iteratively re-weighted least squares (IRLS) for
Quasi-likelihood models. To identify the most accurate but parsimonious model,
we will elaborate our meta-regression framework for these models and
illustrate their strengths and limitations.

####  2.2.1 Poisson model

In practice, the count of foci $Y_{ij}$ (for studies $i=1,\cdots,M$, voxels
$j=1,\cdots,N$) is only ever $0$ or $1$, which strictly indicates a Binomial
model. However, inspired by previous success with Poisson point process, and
accuracy of Poisson approximation for low-rate Binomial data (Eisenberg et
al., 1966), we consider a Poisson model.

If foci arise from a realisation of a (continuous) inhomogeneous Poisson
process, the (discrete) voxel-wise counts will be independently distributed as
Poisson random variables, with rate equal to the integral of the (true,
unobserved, continuous) intensity function over each voxel. As the sum of
multiple independent Poisson random variables is also Poisson, also gives rise
to a practical consequence that it is equivalent to either model the set of
$M$ study-level counts or the summed counts at each voxel. Following the
deterministic structure outlined in Equation 1, the intensity for voxel $j$ in
study $i$ is

$\begin{split}\mathbb{E}[Y_{ij}]&=\mu_{ij}\\\
\log(\mu_{ij})&=\eta_{ij}=x^{\top}_{j}\beta+Z_{i}\gamma\end{split}$ (5)

where $Y_{ij}\sim\mathrm{Poisson}(\mu_{ij})$, $x^{\top}_{j}$ is the $j^{th}$
row of spatial design matrix $X(N\times P)$, and $\beta$ is regression
coefficient of spline bases. The data vector $Y$ has length-$(MN)$, which is
impractical to represent explicitly. Under the assumption of independence of
counts across studies, the likelihood function is exactly same if we model the
voxelwise total foci count over studies instead (more details to follow in
S1.1 in supplementary material), which gives rise to the modified Poisson
model on summed data at voxel $j$ over all studies,
$Y_{\centerdot,j}=\sum\limits_{i=1}^{M}Y_{ij}$,

$\begin{split}\mathbb{E}[Y_{\centerdot,j}]&=\mu_{\centerdot,j}\\\
\mu_{\centerdot,j}&=\sum\limits_{i=1}^{M}\mu_{ij}=\sum\limits_{i=1}^{M}\exp\left(x_{j}^{\top}\beta+Z_{i}\gamma\right)=\exp(x_{j}^{\top}\beta)\left(\sum\limits_{i=1}^{M}\exp(Z_{i}\gamma)\right)\end{split}$
(6)

where $\mu_{\centerdot,j}=\sum\limits_{i=1}^{M}\mu_{ij}$ is the expected sum
of intensity at voxel $j$ over studies. Under this formulation, the likelihood
to be optimised is,

$l(\theta)=l(\beta,\gamma)=\sum\limits_{j=1}^{N}\left[Y_{\centerdot,j}\log(\mu_{\centerdot,j})-\mu_{\centerdot,j}-\log(Y_{\centerdot,j}!)\right]$
(7)

####  2.2.2 Negative Binomial model

While Poisson model is widely used in the regression of count data, it is
recognised that counts often display over-dispersion (the variance of response
variable substantially exceeds the mean). Imposition of Poisson model
according to an unrealistic assumption (variance equals mean) may
underestimate the standard error, and give rise to biased estimation of the
regression coefficients. While Barndorff-Nielsen and Yeo (1969) proposed a
formal definition of spatial Negative Binomial model, it involves Gaussian
processes and complexities we sought to avoid. Hence, here we do not pose a
formal point process model, but rather simply assert that the count data at
each voxel follows a Negative binomial (NB) distribution independently, thus
allowing for anticipated excess variance relative to Poisson (Lawless, 1987).

Our NB model uses a single parameter $\alpha$ shared over all studies and all
voxels to index variance in excess of Poisson model. For each study $i$ and
voxel $j$, let $\lambda_{ij}$ follows a Gamma distribution with mean
$\mu_{ij}$ and variance $\alpha\mu_{ij}^{2}$; then conditioned on
$\lambda_{ij}$, let $Y_{ij}$ be Poisson with mean $\lambda_{ij}$. Then it can
be shown that the marginal distribution of $Y_{ij}$ follows a NB distribution
with probability mass function,

$\mathbb{P}(Y_{ij}=y_{ij})=\frac{\Gamma(y_{ij}+\alpha^{-1})}{\Gamma(y_{ij}+1)\Gamma(\alpha^{-1})}(\frac{1}{1+\alpha\mu_{ij}})^{\alpha^{-1}}(\frac{\alpha\mu_{ij}}{1+\alpha\mu_{ij}})^{y_{ij}}.$
(8)

In terms of the success count and probability parameterization,
$\text{NB}(r,p)$, we have
$Y_{ij}\sim\text{NB}(\alpha^{-1},\frac{\mu_{ij}}{\alpha^{-1}+\mu_{ij}})$, with
mean $\mathbb{E}(Y_{ij})=\mu_{ij}$ and variance
$\mathbb{V}(Y_{ij})=\mu_{ij}+\alpha\mu_{ij}^{2}$. Details on derivation of
probability density function of NB model can be found in S1.2 of the
Supplementary material. When $\alpha>0$, we have Poisson-excess variance of
$\alpha\mu_{ij}^{2}$; or analogous to the coefficient of variation, the
coefficient of excess variation is
$\sqrt{\alpha\mu_{ij}^{2}}/\mu_{ij}=\sqrt{\alpha}$, which can be interpreted
roughly as the relative excess standard deviation, relative to a Poisson
model.

Again, the data vector is impractical to represent explicitly, but unlike
Poisson, the sum of multiple independent NB random variables doesn’t follow a
NB distribution. Thus, we propose moment matching approach to approximate the
mean (first moment) and variance (second moment) of this convolution of NB
distributions, which significantly facilitates the simplification of log-
likelihood function. Matching the first two moments, the approximate NB
distribution of total count of foci over all studies at voxel $j$ is given by
$Y_{\centerdot,j}=\sum\limits_{i=1}^{M}Y_{ij}\sim\mathrm{NB}(r^{\prime}_{j},p^{\prime}_{j})$,
where

$r^{\prime}_{j}=\frac{\mu_{\centerdot,j}^{2}}{\alpha\sum\limits_{i=1}^{M}\mu_{ij}^{2}},\;\;p^{\prime}_{j}=\frac{\sum\limits_{i=1}^{M}\mu_{ij}^{2}}{\alpha^{-1}\mu_{\centerdot,j}+\sum\limits_{i=1}^{M}\mu_{ij}^{2}}$

with corresponding excess variance

$\alpha^{\prime}=\alpha\frac{\sum\limits_{i=1}^{M}\mu_{ij}^{2}}{\mu_{\centerdot,j}^{2}},$

which gives rise to the simplified NB log-likelihood function,

$l(\theta)\approx
l(\beta,\alpha^{\prime})=\sum\limits_{j=1}^{N}\left[\log\Gamma(Y_{\centerdot,j}+r^{\prime}_{j})-\log\Gamma(Y_{\centerdot,j}+1)-\log\Gamma(r^{\prime}_{j})+r^{\prime}_{j}\log{(1-p^{\prime}_{j})+Y_{\centerdot,j}\log{p^{\prime}_{j}}}\right]$
(9)

Details on derivations of moment matching approach can be found in S1.3 in the
Supplementary material.

####  2.2.3 Clustered Negative Binomial model

While NB model can be regarded as a kind of “random effects" Poisson model, as
developed above, the latent Gamma random variable introduces independent
variation at each voxel. We could instead assert that the random (Gamma-
distributed) effects are not independent voxelwise effects, but rather latent
characteristics of each study, and represent a shared effect over the entire
brain for a given study. This is, in fact, the approach used by a Bayesian
CBMA method (Samartsidis et al., 2019), and in a non-imaging setting, a
Poisson-Gamma model for two-stage cluster sampling (Geoffroy and Weerakkody,
2001). Therefore, we now consider a third GLM, where at the first stage, we
assume each individual study $i$ is sampled with a global latent value
$\lambda_{i}$ from a Gamma distribution with mean $1$ and variance $\alpha$,
which accommodates excess variance by dispersion parameter $\alpha$
($\lambda_{i}\sim Gamma(\alpha^{-1},\alpha^{-1})$). At the second stage,
conditioned on the global variable $\lambda_{i}$, $Y_{ij}$ are drawn from a
Poisson distribution with mean $\lambda_{i}\mu_{ij}$
($Y_{ij}|\lambda_{i}\sim\mathrm{Poisson}(\lambda_{i}\mu_{ij})$), where
$\mu_{ij}$ is the expected intensity parameterised by spatial regression
parameter $\beta$ and covariates regression parameter $\gamma$. The marginal
distribution of $Y_{ij}$ also follows a NB distribution,

$\mathbb{P}(Y_{ij}=y_{ij})=\frac{\Gamma(y_{ij}+\alpha^{-1})}{\Gamma(y_{ij}+1)\Gamma(\alpha^{-1})}(\frac{\alpha^{-1}}{\mu_{ij}+\alpha^{-1}})^{\alpha^{-1}}(\frac{\mu_{ij}}{\mu_{ij}+\alpha^{-1}})^{y_{ij}}$
(10)

where $Y_{ij}\sim NB(\alpha^{-1},\frac{\mu_{ij}}{\alpha^{-1}+\mu_{ij}})$ with
mean $\mathbb{E}(Y_{ij})=\mu_{ij}$ and variance
$\mathbb{V}(Y_{ij})=\mu_{ij}+\alpha\mu_{ij}^{2}$. Details on derivation of
probability density function of clustered NB model can be found in S1.4 in
Supplementary material. This two-stage hierarchical Clustered NB model also
introduces a covariance structure between foci within a study, which is
determined by the expected intensity of the observations as well as the
dispersion parameter $\alpha$ (see S1.5 in Supplementary material). The
covariance for studies $i$ and $i^{\prime}$, and distinct voxel $j$ and
$j^{\prime}$ is,

$\begin{cases}\mathbb{C}(Y_{ij},Y_{i^{\prime},j^{\prime}})=\alpha\mu_{ij}\mu_{ij^{\prime}},\;\text{if}\;i=i^{\prime}\\\
\mathbb{C}(Y_{ij},Y_{i^{\prime},j^{\prime}})=0,\;\text{if}\;i\neq
i^{\prime}\end{cases}$ (11)

The log-likelihood is the sum of terms over independent studies,

$\begin{split}l(\beta,\alpha,\gamma)&=\sum\limits_{i=1}^{M}\log[f(Y_{i1},Y_{i2},\cdots,Y_{iN})]\\\
&=M\alpha^{-1}\log(\alpha^{-1})-M\log\Gamma(\alpha^{-1})+\sum\limits_{i=1}^{M}\log\Gamma(Y_{i,\centerdot}+\alpha^{-1})\\\
&-\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}\log{Y_{ij}!}-\sum\limits_{i=1}^{M}(Y_{i,\centerdot}+\alpha^{-1})\log(\mu_{i,\centerdot}+\alpha^{-1})+\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}Y_{ij}\log\mu_{ij}\end{split}$
(12)

where $Y_{i,\centerdot}=\sum\limits_{j=1}^{N}Y_{ij}$ is the sum of foci within
study $i$. One limitation of this model, though, is that it doesn’t admit a
factorisation and depends on the length-(MN) data vector (see S1.5 in
Supplementary material).

Despite a good motivation to induce intra-study dependence, the Clustered NB
model depends on the strong assumption that excess variance is captured by the
global dispersion $\lambda_{i}$. If there is voxelwise independent excess
variance, the previous NB model will be preferred; we assess this issue below
with real data evaluations.

####  2.2.4 Quasi-Poisson model

As an alternative to NB model, Quasi-Poisson model also allows over-dispersed
count data, and is a straightforward elaboration of the GLM. Instead of
specifying a well-defined probability distribution for count data, Quasi-
Poisson model only needs a mean model and a variance function,
$\mathbb{V}(Y_{ij})=\theta\mu_{ij}$ (with $\theta\geq 1$). While the variance-
mean relationship is linear for the Quasi-Poisson model, the relationship is
quadratic in NB model. This results in small foci counts being weighted more
and can have greater adjustment effect in Quasi-Poisson model, which
theoretically might be a perfect fit to our scenario that most brain regions
have zero or low foci counts (Ver Hoef and Boveng, 2007).

Quasi-Poisson model can be framed as GLM, with mean and variance for voxel $j$
in study $i$ is given by,

$\begin{split}E[Y_{ij}]&=\mu_{ij}\\\
\mathrm{Var}(Y_{ij})&=\theta\mu_{ij}.\end{split}$ (13)

Without a likelihood function, we instead use ILRS algorithm, with the
$(j+1)^{th}$ iteration given by,

$\begin{split}\hat{\beta}^{[j+1]}&=\beta^{[j]}+(X^{*\top}W^{[j]}X^{*})^{-1}X^{*\top}(Y-\mu^{[j]})\\\
\hat{\gamma}^{[j+1]}&=\hat{\gamma}^{[j]}+(Z^{*\top}W^{[j]}Z^{*})^{-1}Z^{*\top}(Y-\mu^{[j]})\end{split}$
(14)

where
$W=\mathrm{diag}(\frac{\mu_{11}}{\theta},\cdots,\frac{\mu_{1N}}{\theta},\cdots,\frac{\mu_{M1}}{\theta},\cdots,\frac{\mu_{MN}}{\theta})$,
and $X^{*}=\mathbb{1}_{M}\otimes X$, $Z^{*}=\mathbb{1}_{N}\otimes Z$. This
model can be simplified as well, though we again defer that to the next
Section 2.3.

###  2.3 Model factorisation

Having derived the explicit log-likelihood functions for meta-regression with
three stochastic likelihood-based models, as well as the updating equation for
a quasi-likelihood based model, we now consider model factorisation to replace
the full $(MN)$-vector of foci counts by sufficient statistics. Following the
generic formulation of GLM proposed in Section 2.1.1,

$\eta_{ij}=\log(\mu_{ij})=\sum\limits_{k=1}^{P}X_{jk}\beta_{k}+\sum\limits_{s=1}^{R}Z_{is}\gamma_{s}.$
(15)

$\eta_{ij}$ is the estimated linear response from GLM, specific to each voxel
$j$ in each individual study $i$. In this application, there are always at
least $220,000$ voxels ($N$), hundreds or thousands of studies $M$, and as
many as $P=456$ basis elements (with $20mm$ knots spacing), giving rise to
millions of rows ($MN$) and billions of entries ($MN\times(P+R)$) in a GLM
formulation. Thus, we propose a reformulation of this model into a series of
sufficient statistics that are never larger than $M$ or $N$ in dimension.
First, note that the localised spatial effect $\mu^{X}$ and global effect of
study-level covariates $\mu_{i}^{Z}$ for study $i$ factorise $\mu_{ij}$ as

$\mu_{ij}=\exp\left(\sum\limits_{k=1}^{P}X_{jk}\beta_{k}+\sum\limits_{s=1}^{R}Z_{is}\gamma_{s}\right)=\exp\left(\sum\limits_{k=1}^{P}X_{jk}\beta_{k}\right)\exp\left(\sum\limits_{s=1}^{R}Z_{is}\gamma_{s}\right)=\mu^{X}_{j}\,\mu^{Z}_{i}$
(16)

This model is identifiable as long as we ensure each covariate $Z_{s}$ is mean
zero, letting $X$ capture the overall mean. To further simplify the total log-
likelihood function, we also use the fact that $Y_{ij}\leq 1$ (either $0$ or
$1$), as there will never be more than $1$ foci at the same location in a
given study. Define the following notation:

  * •

Let N-vector $\mu^{X}=\exp(X\beta)$ be the vector of localised spatial effect
of studies;

  * •

let $M$-vector $\mu^{Z}=\exp(Z\gamma)$ be the vector of global study-level
covariates effect of studies;

  * •

Let $Y_{\centerdot,j}=\sum\limits_{i=1}^{M}Y_{ij}$ be the sum of foci counts
at voxel $j$ across all studies, and the $N-$vector
$Y_{\centerdot,}=[Y_{\centerdot,1},\cdots,Y_{\centerdot,N}]^{\top}$;

  * •

$Y_{i,\centerdot}=\sum\limits_{j=1}^{N}Y_{ij}$ be the sum of foci counts for
study $i$ across all voxels, and the $M-$vector
$Y_{,\centerdot}=[Y_{1,\centerdot},\cdots,Y_{M,\centerdot}]^{\top}$.

The simplified factorisation of total log-likelihood functions or IRLS
updating equation are specific to each stochastic model (see S2 in
Supplementary material),

  * •

Poisson model:

$l(\beta,\alpha)=Y_{\centerdot,}^{\top}\log(\mu^{X})+Y_{,\centerdot}^{\top}\log(\mu^{Z})-\left[\boldsymbol{1}^{\top}\mu^{X}\right]\left[\boldsymbol{1}^{\top}\mu^{Z}\right],$
(17)

  * •

NB model: As described in Section 2.2.2, we approximate a sum of independent
NB variables again as a NB:

$Y_{\centerdot,j}=\sum\limits_{i=1}^{M}Y_{ij}\sim\mathrm{NB}(r^{\prime}_{j},p^{\prime}_{j})=\mathrm{NB}(\frac{(\mu^{X}_{j})^{2}[\mathbf{1}^{\top}\mu^{Z}]^{2}}{\alpha^{\prime}\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}},\frac{\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}}{{(\alpha^{\prime})}^{-1}\mu^{X}_{j}[\mathbf{1}^{\top}\mu^{Z}]^{\top}+\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}})$
(18)

with dispersion parameter
$\alpha^{\prime}=\frac{\alpha\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}}{(\mu^{X}_{j})^{2}[\mathbf{1}^{\top}\mu^{Z}]^{2}}$.
The log-likelihood function is given by,

$l(\alpha^{\prime},\beta,\gamma)=\sum\limits_{j=1}^{N}\left[\log\Gamma(Y_{\centerdot,j}+r^{\prime}_{j})-\log\Gamma(Y_{\centerdot,j}+1)-\log\Gamma(r^{\prime}_{j})+r^{\prime}_{j}\log{(1-p^{\prime}_{j})+Y_{\centerdot,j}\log{p^{\prime}_{j}}}\right],$
(19)

  * •

Clustered NB model:

$\begin{split}l(\alpha,\beta,\gamma)&=M\alpha^{-1}\log(\alpha^{-1})-M\log\Gamma(\alpha^{-1})+\sum\limits_{i=1}^{M}\log\Gamma(Y_{i,\centerdot}+\alpha^{-1})\\\
&-\sum\limits_{i=1}^{M}(Y_{i,\centerdot}+\alpha^{-1})\log(\alpha^{-1}+\mu_{i,\centerdot})+Y_{\centerdot,}^{\top}\log(\mu^{X})+Y_{,\centerdot}^{\top}\log(\mu^{Z})\end{split}$
(20)

where dispersion parameter $\alpha$ measures the excess variance across all
studies and all voxels,

  * •

Quasi-Poisson model:

$\begin{split}\hat{\beta}^{[j+1]}&=\beta^{[j]}+(X^{\top}W^{[j]}X)^{-1}X^{\top}(Y_{\centerdot,}-(\mu^{X})^{[j]})\\\
\hat{\gamma}^{[j+1]}&=\hat{\gamma}^{[j]}+(Z^{\top}V^{[j]}Z)^{-1}Z^{\top}(Y_{,\centerdot}-(\mu^{Z})^{[j]})\end{split}$
(21)

where
$W=\mathrm{diag}(\frac{\mu_{1}^{X}}{\theta},\cdots,\frac{\mu_{N}^{X}}{\theta})$
and
$V=\mathrm{diag}(\frac{\mu_{1}^{Z}}{\theta},\frac{\mu_{2}^{Z}}{\theta},\cdots,\frac{\mu_{M}^{Z}}{\theta})$.

###  2.4 Fisher Scoring and optimisation via L-BFGS

Based on explicit log-likelihood functions associated with three stochastic
likelihood-based models (Poisson, NB and clustered NB model) in Section
2.2.1-Section 2.2.3, we employ Fisher scoring for iterative optimisation of
parameters in GLMs. Fisher scoring replaces the gradient and Hessian of
Newton’s method with the score and observed Fisher’s information, respectively
(Longford, 1987). Writing $\theta$ for all parameters, the updating equation
at $(j+1)^{th}$ iteration is,

$\theta^{[k+1]}=\theta^{[k]}+I(\theta^{[k])})^{-1}\frac{\partial}{\partial\theta^{[k]}}l(\theta^{[k])})$
(22)

where the observed Fisher information is
$I(\theta^{[k]})=\mathbb{E}\left[-\frac{\partial^{2}l(\theta)}{\partial\theta\partial\theta^{\top}}\right]_{\theta=\theta^{[k]}}$.

For the Poisson model, $\theta=[\beta,\gamma]$, the Fisher information is
given by,

$I(\theta)=I(\beta,\gamma)=\begin{bmatrix}-\frac{\partial^{2}l}{\partial\beta\partial\beta^{\top}}&-\frac{\partial^{2}l}{\partial\beta\partial\gamma^{\top}}\\\
-\frac{\partial^{2}l}{\partial\gamma\partial\beta^{\top}}&-\frac{\partial^{2}l}{\partial\gamma\partial\gamma^{\top}}\end{bmatrix}$
(23)

with negative Hessian matrix of $\beta$,
$\left(-\frac{\partial^{2}l}{\partial\beta\partial\beta^{\top}}\right)_{P\times
P}=X^{\top}\mathrm{diag}(\mu^{X})X^{\top}$; the negative cross term
$\left(-\frac{\partial^{2}l}{\partial\beta\partial\gamma^{\top}}\right)_{P\times
R}=\left(-\frac{\partial^{2}l}{\partial\gamma\partial\beta^{\top}}\right)^{\top}_{R\times
P}=[X^{\top}\mu^{X}][(\mu^{Z})^{\top}Z]$; and negative Hessian matrix of
$\gamma$,
$\left(-\frac{\partial^{2}l}{\partial\gamma\partial\gamma^{\top}}\right)=Z^{\top}\mathrm{diag}(\mu^{Z})Z$.

Our other stochastic models with study-level covariates lead to more
complicated derivations of updating equations via Fisher scoring. Instead, we
use a more efficient quasi-Newton algorithm (L-BFGS algorithm, Shanno, 1970),
which approximates the observed Fisher information from gradient evaluations
instead, and only requires limited memory and reduces the computation
complexity.

###  2.5 Statistical inference

####  2.5.1 Global test of model fitness

Among the proposed stochastic models in Section 2.2, Poisson, NB and clustered
NB model are likelihood-based, while Quasi-Poisson model is Quasi-likelihood
based (its exact likelihood is computationally infeasible). To compare the
goodness of fit from a global perspective, we will utilise likelihood-based
comparison criteria (e.g., LRT, Akaike information criterion (AIC), Bayesian
information criterion (BIC)) with likelihood-based models, as well as other
global model fitness criteria across all stochastic models within this meta-
regression framework.

Likelihood-based model selection criteria LRT uses the difference in log-
likelihoods to test the null hypothesis that true model is the smaller nested
model. As Poisson model is nested in both NB model and clustered NB model with
dispersion parameter $\alpha=0$, for the null hypothesis $H_{0}$: dispersion
parameter $\alpha=0$, the likelihood-ratio test statistic is given by,

$\lambda_{LR}=-2\left[l(\hat{\theta}_{0})-l(\hat{\theta})\right]$

where $l(\hat{\theta})=l(\hat{\alpha},\hat{\beta},\hat{\gamma})$ is maximum
log-likelihood of NB model or clustered NB model without any constraint on
parameters, and
$l(\hat{\theta}_{0})=l(\hat{\alpha}=0,\hat{\beta},\hat{\gamma})$ is maximum
log-likelihood of NB model or clustered NB model with dispersion parameter
$\alpha$ constrained at $0$ (i.e. Poisson model). The test statistic is Chi-
square distributed, with degree of freedom equals to $1$.

AIC and BIC are two alternatives of LRT which also deal with the trade-off
between the goodness of fit and simplicity of the model, and resolve
overfitting problem by penalising the number of parameters in the model. To
measure the goodness of fit of a model $M$ on dataset $D$,

$AIC=2k-2l(\hat{\theta}),\;\;BIC=k\ln(n)-2l(\hat{\theta})$ (24)

where $l(\hat{\theta})$ is the maximised log-likelihood function of the model
$M$, $k$ is the number of parameters in model $M$ and $n$ is the number of
data points in the dataset $D$. The model with smaller AIC or BIC value is
believed to be a better fit to the dataset.

Bias and variance of estimation For the purpose of selecting the best model in
terms of goodness of fit on a variety of datasets, we extend the model
comparisons to all stochastic models proposed in Section 2.2, including Quasi-
Poisson model. As the central outcome of this meta-regression framework is
voxelwise intensity estimation for each study, with the effect of study-level
covariates being considered, it’s natural to utilise bias and variance of
intensity estimation as new criteria stated below,

  * •

Relative bias of estimated total sum of intensity (per study), comparing with
the averaged sum of foci count (per study) across multiple datasets;

  * •

Relative bias of standard deviation (Std) in each of $x,y,z$ dimension,
comparing with the actual Standard deviation in foci count (per study) across
multiple datasets;

  * •

Relative bias of voxel-wise variance between actual foci count (per study) and
intensity estimation (per study).

Here, relative bias is evaluated instead of bias, especially when applied to a
variety of datasets with diverse foci counts.

####  2.5.2 Localised inference with Wald tests on $\mu^{X}_{ij}$ and
$\eta^{X}_{ij}$

While our model is parameterised by $P$ basis elements, users want to make
inference at each of the $N$ voxels. Hence, we will also explore localised
inference on estimated spatial intensity $\mu^{X}_{ij}$ (or
$\eta^{X}_{ij}=\log(\mu^{X}_{ij})$) and regression coefficient of study-level
covariates ($\gamma$) via Wald tests.

Test of spatial homogeneity:  In the CBMA context, the most basic inference is
a test of homogeneity to identify regions where more foci arise than would be
expected if there were no spatial structure. Precisely, we use the null
hypothesis on voxelwise intensity estimation or estimated linear response,
$H_{0}:\mu^{X}_{ij}=\mu_{0}=\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}Y_{ij}/(MN)$
or $\eta^{X}_{ij}=\eta_{0}=\log(\mu_{0})$ at voxel $j$, for study $i$. The
standard error for $\beta$ can be asymptotically estimated from the inverse of
observed Fisher Information matrix, which gives rise to the standard error for
the linear response $\eta^{X}_{ij}$, and thus standard error for
$\mu^{X}_{ij}$ is obtained via delta method. It allows inference via Wald
tests by examining voxelwise intensity estimation against null hypothesis of
homogeneity over space. The signed Wald statistic for $\mu^{X}_{ij}$ or
$\eta^{X}_{ij}$ takes the form:

$Z_{\mu^{X}}=\frac{\mu^{X}_{ij}-\mu_{0}}{\mathrm{SE}(\mu^{X}_{ij})},\;\;Z_{\eta^{X}}=\frac{\eta^{X}_{ij}-\eta_{0}}{\mathrm{SE}(\eta^{X}_{ij})}$
(25)

where $SE(\mu^{X}_{ij})$ is the standard error of estimated spatial intensity
$\mu^{X}_{ij}$, and $SE(\eta^{X}_{ij})$ is the standard error of estimated
linear response $\eta^{X}_{ij}$, and the statistics are Gaussian
asymptotically. Finally, we can create p-value maps that are thresholded to
control FDR at $5\%$ (Benjamini and Hochberg, 1995).

####  2.5.3 Inference on study-level covariates

For regression coefficient $\gamma$ ($s\times 1$) of study-level covariates,
we consider general linear hypothesis (GLH) tests through a contrast matrix
$C_{\gamma}$ ($m\times s$). Under the null hypothesis,

$H_{0}:C_{\gamma}\gamma=\mathbb{0}_{m\times 1}$ (26)

The test statistic follows a $\chi^{2}$ distribution with $m$ degree of
freedom asymtotically,

$(C_{\gamma}\hat{\gamma})^{T}(C_{\gamma}\mathrm{Cov}(\hat{\gamma})C^{T}_{\gamma})^{-1}(C_{\gamma}\hat{\gamma})\xrightarrow{D}\chi^{2}_{m}$
(27)

and in the case of a single contrast $(m=1)$, a signed Z test can be computed.
Details of GLH on study-level covariates can be found in S3.1 in Supplementary
material.

##  3 Experiments

###  3.1 Simulation settings

The statistical analyses of model estimation with CBMA data are conducted at
voxel level: voxelwise test statistics are evaluated to examine the
significance of experimental effect. Therefore, before investigating the model
fitness, we evaluate our models’ false positive rates (FPR) under null
settings. Due to the computationally intensive nature of these evaluations, we
only evaluated the two models that showed promise in other evaluations,
Poisson and NB. Under the null hypothesis of spatial homogeneity, we use Monte
Carlo (MC) simulation to establish the validity of FPR control for the test of
spatial intensity ($\mu^{X}$). Specifically, we will explore meta-regression
with Poisson or NB model, either with non-null study-level covariates or
without study-level covariates. To ensure the validity of FPR control is
applicable to all CBMA data, the sampling mechanism is either model-based or
empirical, with simulated foci count always analogous to the foci count within
a real dataset. Specifically, in model-based sampling, data generating
mechanism matches with regression model, with number of studies and average
foci per study identical to the original dataset; while in empirical sampling,
real data foci locations are randomly shuffled to guarantee the spatial
homogeneity of foci distribution.

###  3.2 Applications to $20$ meta-analytic datasets

Cognition concerns psychological and cognitive processes that focus on
learning people’s perception, interpretation and response to information and
stimuli. It refers to both conscious procedure and unconscious, automatic
mechanisms in the brain that occur as a response to stimuli, and is highly
variable across individuals (Gallagher et al., 2019). Cognition has been
studied intensively to identify brain regions that are involved in cognition
tasks, conducted in an MRI scanner. For the purpose of evaluating the accuracy
and sensitivity of this meta-regression framework, as well as analysing
goodness of fit of stochastic models with respect to different CBMA datasets,
$20$ previously published meta-analytic datasets are used in this article,
which involves multiple aspects of cognition research, as well as other
stimulus-based and diagnosis-based research (as displayed in Table 1).

[]

Table 1: Number of contrasts and foci counts of $20$ meta-analytic datasets

Dataset | number of contrasts | total count of foci | max foci count | average foci count
---|---|---|---|---
1\. Social Processing | $599$ | $4934$ | $47$ | $8.24$
2\. PTSD | $22$ | $154$ | $26$ | $7.00$
3\. Substance Use | $89$ | $657$ | $110$ | $7.38$
4\. Dementia | $28$ | $1194$ | $548$ | $42.64$
5\. Cue Reactivity | $275$ | $3197$ | $58$ | $11.63$
6\. Emotion Regulation | $338$ | $3543$ | $87$ | $10.48$
7\. Decision Making | $145$ | $1225$ | $49$ | $8.45$
8\. Reward | $850$ | $6791$ | $59$ | $7.99$
9\. Sleep Deprivation | $44$ | $454$ | $59$ | $10.32$
10\. Naturalistic | $122$ | $1220$ | $59$ | $10.00$
11\. Problem Solving | $282$ | $3043$ | $44$ | $10.79$
12\. Emotion | $1738$ | $22038$ | $203$ | $12.68$
13\. Cannabis Use | $81$ | $314$ | $16$ | $3.88$
14\. Nicotine Use | $13$ | $77$ | $23$ | $5.92$
15\. Frontal Pole CBP | $795$ | $9525$ | $57$ | $11.98$
16\. Face Perception | $385$ | $2920$ | $50$ | $7.58$
17\. Nicotine Administration | $75$ | $349$ | $24$ | $4.65$
18\. Executive Function | $243$ | $2629$ | $54$ | $10.82$
19\. Finger Tapping | $76$ | $696$ | $27$ | $9.16$
20\. n-Back | $29$ | $640$ | $69$ | $22.07$

Figure 2: Preprocessing pipeline of $20$ meta-analytic datasets before fitting
CBMR framework. Note that panel A and B are applicable to all datasets, which
generate a spatial design matrix $X$, total foci count per voxel
$Y_{\centerdot,}(N\times 1)$ and total foci count per study
$Y_{,\centerdot}(M\times 1)$. While covariates matrix $Z(M\times R)$ in panel
C is integrated into CBMR only if the effect of study-level covariates is
considered.

Preprocessing steps are summarised in Figure 2. The discrete sampling space of
our analysis is the 2$mm^{3}$ MNI atlas, with dimensions $91\times 109\times
91$, and $N=228483$ brain voxels. We first apply this brain mask to remove
foci outside the brain and remove any multiple-foci (while original data peaks
are always distinct, a foci count in excess of 1 cannot occur after rounding
to the 2mm space). We then extract all the sufficient statistics after model
factorisation in Section 2.3, including spatial design matrix $X(N\times P)$
generated from B-spline bases, total foci count per voxel
$Y_{\centerdot,}(N\times 1)$ and total foci count per study
$Y_{,\centerdot}(M\times 1)$ and study-level covariates $Z(M\times R)$ if
considered.

##  4 Results

###  4.1 Simulation results

For each of the $20$ meta-analytic datasets, we simulate foci distribution
under a null hypothesis of spatial homogeneity, estimate spatial intensity and
investigate the distribution of voxel-wise p-values for the eight different
scenarios: fitting Poisson or NB model, use of a model-based or empirical
(random shuffling) data sampling mechanism, and use or omission of study-level
covariates; for all settings we use B-spline knot of spacing $20mm$ in $x,y,z$
direction, producing $P=456$ basis elements. The computation of test
statistics depends on the covariance of regression coefficients, which is
approximated by the inverse of Fisher Information matrix of optimised
parameters at maximised log-likelihood (see Section 2.4). Empirically, we
sometimes observe the Fisher Information matrix is ill-conditioned and
numerically singular, which is associated with datasets having low foci count.
Especially in datasets with few studies, some regions have essentially no
foci, leading some of the $\beta$ coefficients to be driven to negative
infinity and produce an estimated rate of zero. In the experiments, we found
that datasets with a total foci count of at least 200 generally avoided these
singularity problems and produced accurate standard errors.

To establish the validity of spatial homogeneity tests ($\mu_{j}^{X}=\mu_{0}$,
$\forall j=1,\cdots,N$) for each $20$ meta-analytic datasets, we compute
p-values and create P-P plots. We compute 100 null realisations, each
producing N p-values (one for each voxel), with the null expected $-\log_{10}$
p-values ranging from $-\log_{10}(N/(N+1))\approx 0$ to
$-\log_{10}(1/(N+1))=5.359$. For each expected null p-value, we plot the mean
and 95$\%$ prediction interval via a normal approximation (mean $\pm$ 1.96
$\times$ standard deviation, computed over 100 realisations). Since the P-P
plots are very similar for each of the eight scenarios, we only display
results for the setting of CBMR with Poisson model without study-level
covariates, sampled with a model-based approach. Figure 3 shows the four
representative $-\log_{10}$ P-P plots (results for all 20 studies shown in
Figure S7 in Supplementary material), with identity (dashed diagonal line),
$5\%$ significance (dashed horizontal line) and the FDR $5\%$ boundary (solid
diagonal line); gray shaded areas plot the point-wise $95\%$ prediction
intervals. It shows that p-values $<0.05\approx 10^{-1.3}$ are valid, and
extreme p-values can skew liberal; the worst affected cases are datasets with
very few foci (e.g. analysis 14). In general, datasets with total foci counts
less than 200 show poor behaviour.

Figure 3: P-P plot of $p$-value (under $-\log_{10}$ scale) with four
representative meta-analytic datasets (Social Processing, Substance Use,
Cannabis Use and PTSD datasets), estimated by CBMR with Poisson model without
study-level covariates, sampled with model-based approach.

Since multiple testing correction requires valid p-values far smaller than
$0.05$, we focus on control of FDR in these null simulations. None of the 20
datasets have valid FDR control (PP-plots or prediction intervals fall above
the $5\%$ Benjamini-Hochberg threshold). However, the PP plots generally show
valid p-values $<10^{-3}$, and if we truncate p-values by replacing any
p-value smaller than $10^{-3}$ with that value, we obtain valid (if
conservative) FDR control (Table 2). This pragmatic approach could impact
power, but empirical results (Section 4.2) suggest the inferences based on
truncated p-values remain sensitive.

[]

Table 2: The percentage of invalid FDR control (before/after p-value truncated
at $10^{-3}$) in $20$ meta-analytic datasets over $100$ realisations

Dataset | Before | After | Dataset | Before | After
---|---|---|---|---|---
1\. Social Processing | $44\%$ | $0\%$ | 2\. PTSD | $100\%$ | $0\%$
3\. Substance Use | $26\%$ | $0\%$ | 4\. Dementia | $16\%$ | $0\%$
5\. Cue Reactivity | $28\%$ | $0\%$ | 6\. Emotion Regulation | $23\%$ | $0\%$
7\. Decision Making | $18\%$ | $0\%$ | 8\. Reward | $43\%$ | $0\%$
9\. Sleep Deprivation | $30\%$ | $0\%$ | 10\. Naturalistic | $22\%$ | $0\%$
11\. Problem Solving | $26\%$ | $0\%$ | 12\. Emotion | $100\%$ | $0\%$
13\. Cannabis Use | $63\%$ | $0\%$ | 14\. Nicotine Use | $94\%$ | $0\%$
15\. Frontal Pole CBP | $90\%$ | $0\%$ | 16\. Face Perception | $19\%$ | $0\%$
17\. Nicotine Administration | $54\%$ | $0\%$ | 18\. Executive Function | $22\%$ | $0\%$
19\. Finger Tapping | $22\%$ | $0\%$ | 20\. n-Back | $27\%$ | $0\%$

###  4.2 Results from $20$ meta-analytic datasets

We first evaluate the goodness of fit among likelihood-based stochastic models
(Poisson, NB and clustered NB model) via comparisons of maximised log-
likelihood, AIC and BIC. As shown in Figure S8 in S3.3 of Supplementary
material, CBMR with NB model outperforms two other likelihood-based stochastic
models in every dataset. Not surprisingly, as NB model is the only stochastic
approach that allows for anticipated excess variance relative to Poisson at
voxel level; clustered NB is better than Poisson for the majority of these
$20$ meta-analytic datasets, but only by a small margin. It is conceivable
that although a study-wise global dispersion parameter exists in clustered NB
model, CBMA data is completely specified by a Poisson model at voxel level. As
there exists nested relationships between Poisson and both NB and clustered NB
model (with dispersion parameter $\alpha=0$), we also conduct LRT to evaluate
the trade-off between sufficiency and complexity of model. We found the null
hypothesis that the nested model (Poisson) is better than full model (NB) is
rejected for all datasets, with $p$-value less than $10^{-8}$. The clustered
NB model is preferred over Poisson for the majority of the $20$ meta-analytic
datasets (with $p$-value less than $10^{-8}$) (more details to follow in Table
S5 in Section S3.3 of Supplementary material).

(a)

(b)

(c)

Figure 4: Results from bias-related model comparison criteria, fitted with
four stochastic models on each of $20$ meta-analytic datasets: (a) Boxplots of
relative bias of sum of intensity estimation (per study); (b) Boxplots (each
of $x,y,z$ dimensions) of relative bias of standard deviation of intensity
estimation (per study); (c) Boxplots of relative bias of voxelwise variance of
intensity estimation (per study).

Apart from model comparisons via likelihood-based criteria (LRT, AIC and BIC),
we also integrate Quasi-Poisson model into comparisons of model fitness, based
on bias and variance criteria which are applicable to Quasi-likelihood models.
The results in Figure 4(a) suggest that four stochastic models (Poisson, NB,
clustered NB and Quasi-Poisson model) all give rise to accurate intensity
estimation with relative bias of sum of intensity estimation (per study) less
than $0.5\%$, among which Poisson model has the lowest median relative bias
($0.05\%$) across $20$ meta-analytic datasets. Quasi-Poisson and NB model
display both underestimation and overestimation of study-wise sum of intensity
across $20$ meta-analytic datasets, while sum of intensity are always
overestimated with Poisson and clustered NB model. The results in Figure 4(b)
suggest that CBMR framework also provides an accurate estimation of standard
variation (Std) of intensity in each of $x,y,z$ dimension, with relative bias
controlled within $-0.2\%$ to $0.1\%$ for all stochastic models across $20$
meta-analytic datasets, and estimated intensity along $x$ axis are the most
accurate (with smallest Std bias). As shown in Figure 4(c), CBMR with Poisson
model displays the largest negative bias of variance between intensity
estimation and foci count, which suggests that excess variance cannot be
explained by Poisson assumption (only voxels with nonzero foci count are
studied and included in this plot). Clustered NB model displays second largest
negative relative bias in voxel-wise variance estimation (per study), which is
potentially related to the fact that it evaluates a study-specific over-
dispersion parameter over space, but the intensity function is modelled by
Poisson model at voxel level. Small relative bias is found in both NB and
Quasi-Poisson model (with median $-0.28\%$ and $-0.25\%$), with less variation
in relative bias across multiple datasets with Quasi-Poisson model, which
suggests both models are capable of dealing with excess variance in CBMA data.

Overall, we regard these evaluations as an evidence that NB model is
preferred. While it has slight bias for total intensity (Figure 4(c)), it has
much more accurate variance than the Poisson model.

###  4.3 Comparison with ALE

(a) Z-score map generated by ALE

(b) Z-score map generated by CBMR (with Poisson)

(c) Z-score map generated by CBMR (with NB)

(d) Z-score map generated by CBMR (with Clustered NB)

Figure 5: Activation maps (for significant uncorrected p-values under $5\%$
significance level, presented in Z-scores) generated by ALE (with FWHM=$14$)
and CBMR (with a variety of stochastic models) on the Cue Reactivity dataset
with axial slices at $z=-24,-12,0,12,24,36,48$. Under the null hypothesis of
spatial homogeneity, activation regions with z-scores corresponding to
uncorrected p-values below the significance level $0.05$ are highlighted.

We compared our CBMR results to a widely used approach, computing tests for
spatial homogeneity across space with both CBMA and ALE. For simplicity, we
only demonstrate the comparison of detected activation regions on the Cue
Reactivity dataset (total foci count of $6288$) (Hill-Bowen et al., 2021). For
comparison purposes, we show z statistic values at all voxels significant at
$\alpha=0.05$ uncorrected in Figure 5. Here, we choose FWHM=$14$ to obtain
comparative spatial resolution between ALE and CBMR. Consistency in activation
regions is found in left cerebral cortex, frontal orbital cortex, insular
cortex, left and right accumbens, while spatial specificity of activation
regions differ slightly in ALE and CBMR, with ALE detecting slightly more
voxels.

Another criterion of consistency is the dice similarity coefficient (DSC), the
intersection of ALE and CMBR significant voxels divided by the average number
of significant voxels. As shown in Table 3, ALE appears generally more
sensitive than CBMR regardless of foci counts in the datasets, though DSC
varies above $71.89\%$ to $80.33\%$ on the datasets with more than $1200$ foci
counts, which demonstrates good similarity between the methods.

[]

Table 3: Number of voxels in activation regions of ALE (FWHM=$14$) and CBMR
(with Poisson), based on uncorrected p-values with $5\%$ significance level,
as well as Dice similarity coefficient in $20$ meta-analytic datasets
(Datasets are listed in an ascending order according to total number of foci).

Dataset | n_foci | $|AR_{CBMR}|$ | $|AR_{ALE}|$ | $|AR_{CBMR}\cap AR_{ALE}|$ | DSC
---|---|---|---|---|---
14\. Nicotine Use | $77$ | $1312$ | $12431$ | $1154$ | $17.79\%$
2\. PTSD | $154$ | $6306$ | $15866$ | $5067$ | $45.71\%$
13\. Cannabis Use | $314$ | $11841$ | $18390$ | $8235$ | $54.48\%$
17\. Nicotine Administration | $349$ | $11546$ | $18916$ | $8028$ | $52.71\%$
9\. Sleep Deprivation | $454$ | $10250$ | $15461$ | $5732$ | $44.59\%$
20\. n-Back | $640$ | $19404$ | $31512$ | $17627$ | $69.24\%$
3\. Substance Use | $657$ | $19024$ | $26477$ | $13602$ | $59.79\%$
19\. Finger Tapping | $696$ | $19067$ | $33914$ | $17939$ | $67.72\%$
4\. Dementia | $1194$ | $16244$ | $30437$ | $12464$ | $53.41\%$
10\. Naturalistic | $1220$ | $22328$ | $29442$ | $15344$ | $59.28\%$
7\. Decision Making | $1225$ | $28284$ | $36735$ | $23372$ | $71.89\%$
12\. Emotion | $2038$ | $57698$ | $67699$ | $48847$ | $77.91\%$
18\. Executive Function | $2629$ | $33848$ | $46679$ | $31698$ | $78.73\%$
16\. Face Perception | $2920$ | $41682$ | $53109$ | $36710$ | $77.45\%$
11\. Problem Solving | $3043$ | $38466$ | $51315$ | $34757$ | $77.43\%$
5\. Cue Reactivity | $3197$ | $41242$ | $52371$ | $37301$ | $79.69\%$
6\. Emotion Regulation | $3543$ | $36602$ | $48157$ | $31176$ | $73.56\%$
1\. Social Processing | $4934$ | $48376$ | $61136$ | $40740$ | $74.40\%$
15\. Frontal Pole CBP | $9525$ | $53165$ | $65339$ | $47595$ | $80.33\%$
8\. Reward | $6791$ | $43048$ | $51721$ | $37711$ | $79.59\%$

ALE evaluates experimental effect by testing probabilistic maps (generated by
Gaussian kernel) against the null hypothesis, CBMR estimates activation
intensity and conducts hypothesis testing at voxel level, both of them neglect
the effects of testing all voxels simultaneously, and cannot control the rate
of false rejections. Some researchers proposed a conservative threshold
($\alpha=0.0001$) on the uncorrected p-values to reduce type I error
(Turkeltaub et al., 2002), while a more principled approach is to control the
false discovery rate (FDR) via Benjamini-Hochberg (BH) procedure. Figure 6
shows a comparison of results using a $5\%$ FDR threshold, where CMBR
(Poisson) p-values use a $10^{-3}$ truncation, and Table 4 shows a comparison
of number of detected voxels. The sensitivity of ALE and CBMR is comparable,
with sometimes ALE or CBMR being detecting more voxels. The DSC varies between
$70.55\%$ and $79.76\%$ for datasets with more than $1225$ foci counts,
indicating consistency of activation regions between ALE and CBMR approach
after FDR correction.

[]

Table 4: Number of voxels in activation regions of ALE (FWHM=$14$) and CBMR
(with Poisson), based on FDR corrected p-values (using BH procedure) with
$5\%$ significance level , as well as Dice similarity coefficient in $20$
meta-analytic datasets (Datasets are listed in an ascending order according to
total number of foci).

Dataset | n_foci | $|AR_{CBMR}|$ | $|AR_{ALE}|$ | $|AR_{CBMR}\cap AR_{ALE}|$ | DSC
---|---|---|---|---|---
14\. Nicotine Use | $77$ | $209$ | $0$ | $0$ | $0.00\%$
2\. PTSD | $154$ | $0$ | $1201$ | $0$ | $0.00\%$
13\. Cannabis Use | $314$ | $313$ | $152$ | $17$ | $7.31\%$
17\. Nicotine Administration | $349$ | $1338$ | $943$ | $522$ | $45.77\%$
9\. Sleep Deprivation | $454$ | $176$ | $0$ | $0$ | $0.00\%$
20\. n-Back | $640$ | $11456$ | $17725$ | $10212$ | $69.99\%$
3\. Substance Use | $657$ | $3145$ | $2082$ | $1225$ | $46.87\%$
19\. Finger Tapping | $696$ | $12410$ | $23837$ | $11590$ | $63.95\%$
4\. Dementia | $1194$ | $5126$ | $7931$ | $3142$ | $48.13\%$
10\. Naturalistic | $1220$ | $4192$ | $3241$ | $1861$ | $50.07\%$
7\. Decision Making | $1225$ | $15331$ | $20468$ | $12628$ | $70.55\%$
18\. Executive Function | $2629$ | $26039$ | $37797$ | $24690$ | $77.67\%$
16\. Face Perception | $2920$ | $28893$ | $38193$ | $25533$ | $76.12\%$
11\. Problem Solving | $3043$ | $28221$ | $39091$ | $25675$ | $76.29\%$
5\. Cue Reactivity | $3197$ | $30382$ | $38847$ | $27375$ | $78.57\%$
6\. Emotion Regulation | $3543$ | $23388$ | $31620$ | $20056$ | $72.92\%$
1\. Social Processing | $4943$ | $34317$ | $45263$ | $28555$ | $71.76\%$
8\. Reward | $6791$ | $33021$ | $39743$ | $28728$ | $78.96\%$
15\. Frontal Pole CBP | $9525$ | $44030$ | $55251$ | $39594$ | $79.76\%$
12\. Emotion | $22038$ | $50480$ | $57321$ | $41918$ | $77.77\%$

(a) Z-score map generated by ALE

(b) Z-score map generated by CBMR (with Poisson)

(c) Z-score generated by CBMR (with NB)

(d) Z-score map generated by CBMR (with Clustered NB)

Figure 6: Activation maps (for significant FDR corrected p-values under $5\%$
significance level, presented in Z-scores) generated by ALE and CBMR with FDR
correction (by BH procedure) with truncated p-values of Cue Reactivity
dataset. The figure is shown with axial slices at $z=-24,-12,0,12,24,36,48$.
Under the null hypothesis of spatial homogeneity, activation regions with
z-scores corresponding to corrected p-values below the significance level
$0.05$ are highlighted.

###  4.4 Effect of study-level covariates

Unlike ALE, CBMR is based on explicit probabilistic models, and can estimate
the effect of study-level covariates. Here, we integrate two study-level
covariates, study-wise (square root) sample size and year of publication
(after centring and standardisation) into CBMR framework on each of $20$ meta-
analytic datasets. We find, for example, on Cue Reactivity dataset, the year
of publication is not significant ($Z=-0.6880,p=0.4915$), while sample size is
significant ($Z=6.1454,p<10^{-8}$) (See Table S6 for $p-$value and $Z-$score
of study-level covariates on each of $20$ meta-analytic datasets).

##  5 Discussion

In this work we have presented a meta-regression framework with a spatial
model as a general approach for CBMA data, where we have considered multiple
stochastic models and allowed for study-level factors (e.g., sample size and
year of publication). Our approach uses a spline parameterization to model the
smooth spatial distribution of activation foci, and fits a generalised linear
model with different variants of voxelwise (Poisson model, NB model and Quasi-
Poisson model) or study-wise (Clustered NB model) statistical distributions.
Our approach is a computationally efficient alternative to previous Bayesian
spatial regression models, providing the flexibility and interpretability of a
regression model while jointly modelling all of space. For comparison, the
implementation of Bayesian log-Gaussian Cox process regression needed
approximately $30$ hours on an NVIDIA Tesla K20c GPU card Samartsidis et al.
(2019), while our meta-regression runs for roughly $20$ minutes on an NVIDIA
GTX 1080 Graphics Card. Furthermore, as a more intuitive and interpretable
approach derived from generalised linear model, we believe that our meta-
regression framework is more comprehensible to practitioners, relative to the
spatial posterior intensity function. Through simulations on synthetic data
(with simulated foci counts analogous to foci counts in each of $20$ meta-
analytic datasets), we demonstrated valid FDR control for spatial homogeneity
null hypothesis after a truncation of $p-$values below $10^{-3}$. According to
$20$ meta-analytic datasets, we found that NB model is the most accurate
stochastic model in model comparisons via LRT, AIC and BIC, as well as
smallest relative bias in both mean and variance of intensity estimation (per
study), while Poisson and clustered NB model cannot explain over-dispersion
observed in foci count. Meanwhile, we also compare the findings of activation
regions from both ALE and CBMR approach, and justify the validity and
robustness of CBMR, especially on the datasets with relatively high foci
count, e.g., datasets with at least $200$ total foci.

There are a few limitations in our work. Here we have only considered a single
group of studies. In future work, we will extend our method to estimate the
spatial intensity function of multiple groups (e.g., multiple types of stimuli
within a cognitive task), so that we can investigate the consistency and
difference of activation regions by group comparison. Meanwhile, we are
currently not using regularisation term on spatial regression coefficients of
CBMR. We have considered Firth-type penalty which indeed guarantees convergent
estimates (especially on brain regions without any foci) and removes the
first-order asymptotic bias term of maximum likelihood estimates, but it also
causes a significant over-estimation of intensity at the edge of brain mask.
The edge effect induced by Firth-type penalty relates to the structure of the
Jeffreys prior and higher intensity associated with edge and corner basis
elements. However, it’s plausible to consider regularising likelihood
functions with alternative penalty term (e.g., $L_{1}$ or $L_{2}$ norm) in the
future, and figure out the optimal value of hyper-parameter. To estimate the
variance of voxelwise spatial intensity, we approximate the covariance of
spatial regression coefficients by inverting the Fisher Information matrix. It
might give rise to numerical instability because the dimension of Fisher
Information matrix is large (there are hundreds or even thousands elements in
spline bases), and it might even be numerically singular for datasets with low
foci count since most of voxels have near-zero intensity estimation. We have
tried many approaches to improve numerical stability, including adding a
extremely small epsilon ($10^{-6}$) or $1\%$ of the largest diagonal element
on the diagonal of Fisher Information matrix, or computing the Fisher
Information assuming the null hypothesis of homogeneity is true. However, all
of these efforts produced underestimation of the variance of voxelwise spatial
intensity and gave rise to invalid p-values. In future work, we might consider
non-parametric methods to estimate the covariance of spatial regression
coefficient instead of the inverse of Fisher Information, or add a
regularisation term on B-spline roughness to avoid very negative spatial
regression coefficients.

Another potential direction is to conduct meta-analysis with data from
multiple source, specifically, integrate additional information about reported
foci or full statistic map (e.g., p-values or t-scores) if available. Some
researches proposed Markov melding as a fully Bayesian framework for joining
probabilistic sub-models, where evidence from different source is specified in
each sub-model, and sub-models are joined while preserving all information and
uncertainty (Goudie et al., 2019). Such an approach might enrich the inference
obtained from CBMR by integrating the magnitude of CBMA activation or even
image-based meta-analysis data. Finally. it’s worth considering a zero-
inflated stochastic model (e.g., Poisson or NB model) as the current datasets
only consist of studies with at least one focus, there might be inflated zero
foci count than observed. Excess zeros are separated and modelled
independently in zero-inflated models, which might provide a more accurate
approximation for low-rate Binomial data.

## Software

Implementation in the form of Python and Pytorch code can be found in Github
repository. CBMR framework has also been implemented and integrated into
NiMARE python package.

## Acknowledgments

The computational aspects of this research were supported by the Wellcome
Trust Core Award Grant Number 203141/Z/16/Z and the NIHR Oxford BRC. The views
expressed are those of the authors and not necessarily those of the NHS, the
NIHR or the Department of Health.

Conflict of Interest: The authors declare no conflict of interest.

## Funding

This work was supported by the National Institutes of Health (NIH) under Award
Numbers H6R00550_CS00.01.

## References

  * Barndorff-Nielsen and Yeo [1969] Ole Barndorff-Nielsen and GF Yeo.  Negative binomial processes.  _Journal of Applied Probability_ , 6(3):633–647, 1969.
  * Benjamini and Hochberg [1995] Yoav Benjamini and Yosef Hochberg.  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  _Journal of the Royal statistical society: series B (Methodological)_ , 57(1):289–300, 1995.
  * Eickhoff et al. [2012] Simon B Eickhoff, Danilo Bzdok, Angela R Laird, Florian Kurth, and Peter T Fox.  Activation likelihood estimation meta-analysis revisited.  _Neuroimage_ , 59(3):2349–2361, 2012.
  * Eisenberg et al. [1966] Herbert B Eisenberg, Randolph RM Geoghagen, and John E Walsh.  A general use of the poisson approximation for binomial events, with application to bacterial endocarditis data.  _Biometrics_ , pages 74–82, 1966.
  * Gallagher et al. [2019] Anne Gallagher, Christine Bulteau, David Cohen, and Jacques L Michaud.  _Neurocognitive Development: Normative Development_.  Elsevier, 2019.
  * Geoffroy and Weerakkody [2001] Pedro Geoffroy and Govinda Weerakkody.  A poisson-gamma model for two-stage cluster sampling data.  _Journal of Statistical Computation and Simulation_ , 68(2):161–172, 2001.
  * Goudie et al. [2019] Robert JB Goudie, Anne M Presanis, David Lunn, Daniela De Angelis, and Lorenz Wernisch.  Joining and splitting models with markov melding.  _Bayesian analysis_ , 14(1):81, 2019.
  * Hill-Bowen et al. [2021] Lauren D Hill-Bowen, Michael C Riedel, Ranjita Poudel, Taylor Salo, Jessica S Flannery, Julia A Camilleri, Simon B Eickhoff, Angela R Laird, and Matthew T Sutherland.  The cue-reactivity paradigm: An ensemble of networks driving attention and cognition when viewing drug and natural reward-related stimuli.  _Neuroscience & Biobehavioral Reviews_, 130:201–213, 2021\.
  * Kang et al. [2011] Jian Kang, Timothy D Johnson, Thomas E Nichols, and Tor D Wager.  Meta analysis of functional neuroimaging data via bayesian spatial point processes.  _Journal of the American Statistical Association_ , 106(493):124–134, 2011.
  * Kang et al. [2014] Jian Kang, Thomas E Nichols, Tor D Wager, and Timothy D Johnson.  A bayesian hierarchical spatial point process model for multi-type neuroimaging meta-analysis.  _The annals of applied statistics_ , 8(3):1800, 2014.
  * Laird et al. [2005] Angela R Laird, Jack J Lancaster, and Peter T Fox.  Brainmap.  _Neuroinformatics_ , 3(1):65–77, 2005.
  * Lawless [1987] Jerald F Lawless.  Negative binomial and mixed poisson regression.  _The Canadian Journal of Statistics/La Revue Canadienne de Statistique_ , pages 209–225, 1987.
  * Longford [1987] Nicholas T Longford.  A fast scoring algorithm for maximum likelihood estimation in unbalanced mixed models with nested random effects.  _Biometrika_ , 74(4):817–827, 1987.
  * Montagna et al. [2018] Silvia Montagna, Tor Wager, Lisa Feldman Barrett, Timothy D Johnson, and Thomas E Nichols.  Spatial bayesian latent factor regression modeling of coordinate-based meta-analysis data.  _Biometrics_ , 74(1):342–353, 2018.
  * Radua et al. [2012] J Radua, D Mataix-Cols, Mary L Phillips, W El-Hage, DM Kronhaus, N Cardoner, and S Surguladze.  A new meta-analytic method for neuroimaging studies that combines reported peak coordinates and statistical parametric maps.  _European psychiatry_ , 27(8):605–611, 2012.
  * Salimi-Khorshidi et al. [2009] Gholamreza Salimi-Khorshidi, Stephen M Smith, John R Keltner, Tor D Wager, and Thomas E Nichols.  Meta-analysis of neuroimaging data: a comparison of image-based and coordinate-based pooling of studies.  _Neuroimage_ , 45(3):810–823, 2009.
  * Samartsidis et al. [2017a] Pantelis Samartsidis, Silvia Montagna, Angela R Laird, Peter T Fox, Timothy D Johnson, and Thomas E Nichols.  Estimating the number of missing experiments in a neuroimaging meta-analysis.  _BioRxiv_ , page 225425, 2017a.
  * Samartsidis et al. [2017b] Pantelis Samartsidis, Silvia Montagna, Thomas E Nichols, and Timothy D Johnson.  The coordinate-based meta-analysis of neuroimaging data.  _Statistical science: a review journal of the Institute of Mathematical Statistics_ , 32(4):580, 2017b.
  * Samartsidis et al. [2019] Pantelis Samartsidis, Claudia R Eickhoff, Simon B Eickhoff, Tor D Wager, Lisa Feldman Barrett, Shir Atzil, Timothy D Johnson, and Thomas E Nichols.  Bayesian log-gaussian cox process regression: with applications to meta-analysis of neuroimaging working memory studies.  _Journal of the Royal Statistical Society. Series C, Applied statistics_ , 68(1):217, 2019.
  * Shanno [1970] David F Shanno.  Conditioning of quasi-newton methods for function minimization.  _Mathematics of computation_ , 24(111):647–656, 1970.
  * Turkeltaub et al. [2002] Peter E Turkeltaub, Guinevere F Eden, Karen M Jones, and Thomas A Zeffiro.  Meta-analysis of the functional neuroanatomy of single-word reading: method and validation.  _Neuroimage_ , 16(3):765–780, 2002.
  * Ver Hoef and Boveng [2007] Jay M Ver Hoef and Peter L Boveng.  Quasi-poisson vs. negative binomial regression: how should we model overdispersed count data?  _Ecology_ , 88(11):2766–2772, 2007.
  * Wager et al. [2007] Tor D Wager, Martin Lindquist, and Lauren Kaplan.  Meta-analysis of functional neuroimaging data: current and future directions.  _Social cognitive and affective neuroscience_ , 2(2):150–158, 2007.
  * Westfall and Young [1993] Peter H Westfall and S Stanley Young.  _Resampling-based multiple testing: Examples and methods for p-value adjustment_ , volume 279.  John Wiley & Sons, 1993.
  * Yarkoni et al. [2011] Tal Yarkoni, Russell A Poldrack, Thomas E Nichols, David C Van Essen, and Tor D Wager.  Large-scale automated synthesis of human functional neuroimaging data.  _Nature methods_ , 8(8):665–670, 2011.
  * Yue et al. [2012] Yu Ryan Yue, Martin A Lindquist, and Ji Meng Loh.  Meta-analysis of functional neuroimaging data using bayesian nonparametric binary regression.  _The Annals of Applied Statistics_ , pages 697–718, 2012.

## Supplementary Material

##  Appendix S1 Detailed derivation of stochastic models

###  S1.1 Poisson model

We assert that the sum of two independent Poisson random variables is also
Poisson. Let $X\sim\mathrm{Poi}(\lambda_{1})$ and
$Y\sim\mathrm{Poi}(\lambda_{2})$ be two independent random variables, and
$Z=X+Y$, then,

$\begin{split}P(Z=n)=P(X+Y=N)&=\sum\limits_{k=-\infty}^{\infty}P(X=k)P(Y=n-K)\\\
&=\sum\limits_{k=0}^{n}P(X=k)P(Y=n-k)\\\
&=\sum\limits_{k=0}^{n}e^{-\lambda_{1}}\frac{\lambda_{1}^{k}}{k!}e^{-\lambda_{2}}\frac{\lambda_{2}^{n-k}}{(n-k)!}\\\
&=e^{-(\lambda_{1}+\lambda_{2})}\sum\limits_{k=0}^{n}\frac{\lambda_{1}^{k}\lambda_{2}^{n-k}}{k!(n-k)!}\\\
&=\frac{e^{-(\lambda_{1}+\lambda_{2})}}{n!}\sum\limits_{k=0}^{n}\frac{n!}{k!(n-k)!}\lambda_{1}^{k}\lambda_{2}^{n-k}\\\
&=\frac{e^{-(\lambda_{1}+\lambda_{2})}}{n!}(\lambda_{1}+\lambda_{2})^{n}\end{split}$
(28)

Therefore, $Z=X+Y\sim\mathrm{Poi}(\lambda_{1}+\lambda_{2})$ is also a Poisson
variable. The conclusion can be further extended: the sum of multiple Poisson
random variables ($\mathrm{Poi}(\lambda_{i}),i=1,\cdots,n$) also follows
Poisson distribution, with the parameter
$\lambda=\sum\limits_{i=1}^{n}\lambda_{i}$.

Hence, under the assumption of independence of counts across studies, we
believe that the likelihood function is exactly same if we model the voxelwise
total foci count over studies (with length-N) instead of voxelwise foci count
for each study (with expanded length-(NM)). This reformulation can simplify
the computation of the log-likelihood function and reduce the dimensionality
of statistics (never larger than $M$ or $N$ in dimension).

###  S1.2 Negative Binomial (Poisson-Gamma) Model

In this section, we describe the formulation of NB distribution in more
detail. Based on the assumption of NB (Poisson-Gamma) model, there’s a single
parameter $\alpha$ indicates variance in excess of Poisson model, for voxel
$j$ in study $i$, the voxelwise mean of intensity $\lambda_{ij}$ follows a
Gamma distribution with mean $\mathbb{E}(\lambda_{ij})=\mu_{ij}$ and variance
$\mathrm{Var}(\lambda_{ij})=\alpha\mu_{ij}^{2}$

$\lambda_{ij}\sim\mathrm{Gamma}(\alpha^{-1},\frac{\alpha^{-1}}{\mu_{ij}})\Rightarrow\mathbb{E}(\lambda_{ij})=\mu_{ij},\mathrm{Var}(\lambda_{ij})=\alpha\mu_{ij}^{2}$

And $Y_{ij}|\lambda_{ij}$ follows a Poisson distribution with conditional mean
$\mathbb{E}(Y_{ij}|\lambda_{ij})=\lambda_{ij}$

$Y_{ij}|\lambda_{ij}=\mathrm{Poisson}(\lambda_{ij})\Rightarrow
P(Y_{ij}|\lambda_{ij}=k)=\frac{\lambda_{ij}^{k}e^{-\lambda_{ij}}}{k!}$

which gives rise to marginal probability of $Y_{ij}$

$\begin{split}P(Y_{ij}=y_{ij})&=\int\limits_{\lambda_{ij}}P(Y_{ij}|\lambda_{ij})P(\lambda_{ij})d\lambda_{ij}=\int\limits_{\lambda_{ij}=0}^{\infty}\frac{\lambda_{ij}^{y_{ij}}e^{-\lambda_{ij}}}{y_{ij}!}\frac{(\frac{1}{\alpha\mu_{ij}})^{\frac{1}{\alpha}}}{\Gamma(\frac{1}{\alpha})}\lambda_{ij}^{\frac{1}{\alpha}-1}e^{-\frac{\lambda_{ij}}{\alpha\mu_{ij}}}d\lambda_{ij}\\\
&=\frac{\frac{1}{\alpha\mu_{ij}}^{\frac{1}{\alpha}}}{y_{ij}!\Gamma(\frac{1}{\alpha})}\int\limits_{\lambda_{ij}=0}^{\infty}\lambda_{ij}^{y_{ij}}e^{-\lambda_{ij}}\lambda_{ij}^{\frac{1}{\alpha}-1}e^{-\frac{\lambda_{ij}}{\alpha\mu_{ij}}}d\lambda_{ij}=\frac{\frac{1}{\alpha\mu_{ij}}^{\frac{1}{\alpha}}}{y_{ij}!\Gamma(\frac{1}{\alpha})}\frac{\Gamma(y_{ij}+\frac{1}{\alpha})}{(\frac{1}{\alpha\mu_{ij}}+1)^{y_{ij}+\frac{1}{\alpha}}}\\\
&=\frac{\Gamma(y_{ij}+\frac{1}{\alpha})}{\Gamma(y_{ij}+1)\Gamma(\frac{1}{\alpha})}(\frac{\frac{1}{\alpha\mu_{ij}}}{\frac{1}{\alpha\mu_{ij}}+1})^{\frac{1}{\alpha}}(\frac{1}{\frac{1}{\alpha\mu_{ij}}+1})^{y_{ij}}=\frac{\Gamma(y_{ij}+\alpha^{-1})}{\Gamma(y_{ij}+1)\Gamma(\alpha^{-1})}(\frac{1}{1+\alpha\mu_{ij}})^{\alpha^{-1}}(\frac{\alpha\mu_{ij}}{1+\alpha\mu_{ij}})^{y_{ij}}\end{split}$

which satisfies the mathematical form of probability density function of NB
model, $Y_{ij}\sim NB(\alpha^{-1},\frac{\mu_{ij}}{\alpha^{-1}+\mu_{ij}})$,
with mean $\mathbb{E}[Y_{ij}]=\mu_{ij}$ and variance
$\mathbb{V}(Y_{ij})=\mu_{ij}+\alpha\mu_{ij}^{2}$.

###  S1.3 Moment Matching Approach

For the purpose of approximating the sum of multiple independent NB random
variables, we approximate a sum of NB variates with a NB distribution by
moment matching (mean and variance). Suppose the voxelwise count in each
individual study is
$Y_{ij}\sim\mathrm{NB}(\alpha^{-1},\frac{\mu_{ij}}{\mu_{ij}+\alpha^{-1}})$,
and $\alpha$ is global dispersion parameter. Using the independence of studies
at voxel $j$,

$\begin{cases}\mathbb{E}(Y_{\centerdot,j})=\sum\limits_{i=1}^{M}\mathbb{E}(Y_{ij})=\sum\limits_{i=1}^{M}\mu_{ij}\\\
\mathbb{V}(\mathbb{E}(Y_{\centerdot,j}))=\sum\limits_{i=1}^{M}\mathrm{Var}(Y_{ij})=\sum\limits_{i=1}^{M}\mu_{ij}+\sum\limits_{i=1}^{M}\alpha\mu_{ij}^{2}\end{cases}$

To ensure that the proposed NB distribution ($Y_{\centerdot,j}\sim
NB(r^{\prime},p^{\prime})$) matches the mixture of NB distributions, with
regard to both mean and variance, we need

$\begin{cases}\mathbb{E}(Y_{\centerdot,j})=\sum\limits_{i=1}^{M}\mu_{ij}\\\
\mathrm{Var}(Y_{\centerdot,j})=\sum\limits_{i=1}^{M}\mu_{ij}+\sum\limits_{i=1}^{M}\alpha\mu_{ij}^{2}\end{cases}\Rightarrow\begin{cases}p^{\prime}=\frac{\sum\limits_{i=1}^{M}\mu_{ij}^{2}}{\alpha^{-1}\sum\limits_{i=1}^{M}\mu_{ij}+\sum\limits_{i=1}^{M}\mu_{ij}^{2}}\\\
r^{\prime}=\frac{(\sum\limits_{i=1}^{M}\mu_{ij})^{2}}{\alpha\sum\limits_{i=1}^{M}\mu_{ij}^{2}}\end{cases}$

Therefore, approximated NB distribution of sum of foci count at voxel $j$ is,
$Y_{\centerdot,j}\sim\mathrm{NB}\left(\frac{(\sum\limits_{i=1}^{M}\mu_{ij})^{2}}{\alpha\sum\limits_{i=1}^{M}\mu_{ij}^{2}},\frac{\sum\limits_{i=1}^{M}\mu_{ij}^{2}}{v\sum\limits_{i=1}^{M}\mu_{ij}+\sum\limits_{i=1}^{M}\mu_{ij}^{2}}\right)$,
with excess variance in NB approximation $\alpha^{\prime}$,

$\frac{1}{\alpha^{\prime}}=\frac{(\sum\limits_{i=1}^{M}\mu_{ij})^{2}}{\alpha\sum\limits_{i=1}^{M}\mu_{ij}^{2}}\Rightarrow\alpha^{\prime}=\frac{\sum\limits_{i=1}^{M}\mu_{ij}^{2}}{(\sum\limits_{i=1}^{M}\mu_{ij})^{2}}\alpha$

.

###  S1.4 Two-stage hierarchy Poisson-Gamma model

In this section, we propose a two-stage hierarchy Poisson-Gamma model, which
regards random (Gamma) effect as a latent characteristic of each study,
instead of independent voxelwise effects. The name “two-stage hierarchy
Poisson-Gamma model" comes from the modelling procedure: consider the
clustered count data $Y_{ij}$, $i=1,\cdots,M$ (number of studies),
$j=1,\cdots,N$ (number of voxels). Draw $\lambda_{i}$ from a Gamma
distribution with mean $1$ and variance $\alpha$.

$\lambda_{i}\sim\mathrm{Gamma}(\alpha^{-1},\alpha^{-1})\Rightarrow\mathbb{E}(\lambda_{i})=1,\mathbb{V}(\lambda_{i})=\alpha,$
$f(\lambda_{i})=\frac{\frac{1}{\alpha}^{\alpha^{-1}}}{\Gamma(\alpha^{-1})}\lambda_{i}^{\alpha-1}e^{-\alpha^{-1}\lambda_{i}}$

For each study $i$, for each voxel $j$, $Y_{ij}|\lambda_{i}$ are drawn from a
Poisson distribution with mean $\lambda_{i}\mu_{ij}$, where $\mu_{ij}$ is the
spatial mean parameterized by some $\beta$ (B-spline basis coefficients).

$Y_{ij}|\lambda_{i}\sim\mathrm{Poisson}(\lambda_{i}\mu_{ij})\Rightarrow
P(Y_{ij}|\lambda_{i}=k)=\frac{(\lambda_{i}\mu_{ij})^{k}e^{-\lambda_{i}\mu_{ij}}}{k!}$

Therefore, marginal probability of foci count $Y_{ij}$ is,

$\begin{split}P(Y_{ij}=y_{ij})&=\int\limits_{\lambda_{i}}P(Y_{ij}|\lambda_{i})P(\mu_{i})d\lambda_{i}=\int\limits_{\lambda_{i}=0}^{\infty}\frac{(\mu_{ij}\lambda_{i})^{y_{ij}}e^{-\mu_{ij}\lambda_{i}}}{y_{ij}!}\frac{\frac{1}{\alpha}^{\alpha^{-1}}}{\Gamma(\alpha^{-1})}\lambda_{i}^{\alpha^{-1}-1}e^{-\alpha^{-1}\lambda_{i}}d\lambda_{i}\\\
&=\frac{\mu_{ij}^{y_{ij}}\frac{1}{\alpha}^{\alpha^{-1}}}{\Gamma(\alpha^{-1})y_{ij}!}\int\limits_{\lambda_{i}=0}^{\infty}\lambda_{i}^{y_{ij}+\alpha^{-1}-1}e^{-\lambda_{i}(\mu_{ij}+\alpha^{-1})}d\lambda_{i}=\frac{\mu_{ij}^{y_{ij}}\frac{1}{\alpha}^{\alpha^{-1}}}{\Gamma(\alpha^{-1})y_{ij}!}\frac{\Gamma(y_{ij}+\alpha^{-1})}{(\mu_{ij}+\alpha^{-1})^{y_{ij}+\alpha^{-1}}}\\\
&=\frac{\Gamma(y_{ij}+\alpha^{-1})}{\Gamma(y_{ij}+1)\Gamma(\alpha^{-1})}(\frac{\mu_{ij}}{\mu_{ij}+\alpha^{-1}})^{y_{ij}}(\frac{\alpha^{-1}}{\mu_{ij}+\alpha^{-1}})^{\alpha^{-1}}\end{split}$

Therefore, marginal distribution of foci count follows a NB distribution,
$Y_{ij}\sim\mathrm{NB}(\alpha^{-1},\frac{\mu_{ij}}{\mu_{ij}+\alpha^{-1}})$,
with mean $\mathbb{E}(Y_{ij})=\mu_{ij}$ and variance
$\mathrm{Var}(Y_{ij})=\mu_{ij}+\alpha\mu_{ij}^{2}$.

###  S1.5 Covariance structure in Clustered NB model

The two-stage hierarchical clustered NB model also introduces a covariance
structure between foci within a study, specifically, covariance of number of
count in voxel $j$ and voxel $j^{\prime}$ ($Y_{ij}$ and $Y_{ij^{\prime}}$) in
study $i$, modelled by clustered NB model is,

$\begin{split}\mathbb{E}[Y_{ij}Y_{ij^{\prime}}]&=\mathbb{E}_{\lambda_{i}}\left[\mathbb{E}[Y_{ij}Y_{ij^{\prime}}|\lambda_{i}]\right]=\mathbb{E}_{\lambda_{i}}\left[\mathbb{E}[Y_{ij}|\lambda_{i}]\mathbb{E}[Y_{ij^{\prime}}|\lambda_{i}]\right]=\mathbb{E}_{\lambda_{i}}\left[(\lambda_{i}\mu_{ij})(\lambda_{i}\mu_{ij^{\prime}})\right]\\\
&=\mu_{ij}\mu_{ij^{\prime}}\int_{\lambda_{i}}\lambda_{i}^{2}f(\lambda_{i})d\lambda_{i}=\mu_{ij}\mu_{ij^{\prime}}\int_{\lambda_{i}}\frac{\frac{1}{\alpha}^{\alpha^{-1}}}{\Gamma(\alpha^{-1})}\lambda_{i}^{\alpha^{-1}+1}e^{-\alpha^{-1}\lambda_{i}}d\lambda_{i}\\\
&=\mu_{ij}\mu_{ij^{\prime}}\frac{\frac{1}{\alpha}^{\alpha^{-1}}}{\Gamma(\alpha^{-1})}\alpha^{\alpha^{-1}+1}\alpha\int_{\lambda_{i}}(\frac{1}{\alpha}\lambda_{i})^{(\alpha^{-1}+1)}e^{\alpha^{-1}\lambda_{i}}d(\alpha^{-1}\lambda_{i})\\\
&=\mu_{ij}\mu_{ij^{\prime}}\alpha^{2}\frac{1}{\Gamma(\alpha^{-1})}\Gamma(\alpha^{-1}+2)=\mu_{ij}\mu_{ij^{\prime}}\alpha^{2}\alpha^{-1}(\alpha^{-1}+1)\\\
&=(1+\alpha)\mu_{ij}\mu_{ij^{\prime}}\end{split}$

For different voxels $j$ and $j^{\prime}$ within a same study $i$,

$\mathrm{Cov}(Y_{ij}Y_{ij^{\prime}})=\mathbb{E}[Y_{ij}Y_{ij^{\prime}}]-\mathbb{E}[Y_{ij}]\mathbb{E}[Y_{ij^{\prime}}]=(1+\alpha)\mu_{ij}\mu_{ij^{\prime}}-\mu_{ij}\mu_{ij^{\prime}}=\alpha\mu_{ij}\mu_{ij^{\prime}}$

While for different studies $i$ and $i^{\prime}$,

$\mathrm{Cov}(Y_{ij}Y_{i^{\prime}j^{\prime}})=0$

Now, we will look at the total log-likelihood function of clustered NB model,
using dependence between studies. Let
$Y_{i,\centerdot}=\sum\limits_{j=1}^{N}Y_{ij}$ be the sum of foci within study
$i$ regardless of location. The joint probability for the number of count
$Y_{ij}\,(\forall j=1,\cdots,N)$ from the $i^{th}$ study is,

$\begin{split}&f(Y_{i1},Y_{i2},\cdots,Y_{iN})=\int_{\lambda_{i}}f(Y_{i1},Y_{i2},\cdots,Y_{iN}|\lambda_{i})f(\lambda_{i})d\lambda_{i}=\int_{\lambda_{i}}\prod\limits_{j=1}^{N}f(Y_{ij}|\lambda_{i})f(\lambda_{i})d\lambda_{i}\\\
&=\int_{\lambda_{i}}\prod\limits_{j=1}^{N}\frac{(\mu_{ij}\lambda_{i})^{Y_{ij}}e^{-\mu_{ij}\lambda_{i}}}{Y_{ij}!}\frac{\frac{1}{\alpha}^{\alpha^{-1}}}{\Gamma(\alpha^{-1})}\lambda_{i}^{\alpha^{-1}-1}e^{-\alpha^{-1}\lambda_{i}}d\lambda_{i}\\\
&=\frac{\frac{1}{\alpha}^{\alpha^{-1}}\prod\limits_{j=1}^{N}\mu_{ij}^{Y_{ij}}}{\Gamma(\alpha^{-1})\prod\limits_{j=1}^{N}Y_{ij}!}\int_{\lambda_{i}}\exp\\{-\lambda_{i}(\alpha^{-1}+\sum\limits_{j=1}^{N}\mu_{ij})\\}\lambda_{i}^{\sum\limits_{j=1}^{N}Y_{ij}+\alpha^{-1}-1}d\lambda_{i}\\\
&=\frac{\frac{1}{\alpha}^{\alpha^{-1}}\prod\limits_{j=1}^{N}\mu_{ij}^{Y_{ij}}}{\Gamma(\alpha^{-1})\prod\limits_{j=1}^{N}Y_{ij}!}\int_{\lambda_{i}}\exp\\{-\lambda_{i}(\alpha^{-1}+\mu_{i,\centerdot})\\}\lambda_{i}^{Y_{i,\centerdot}+\alpha^{-1}-1}d\lambda_{i}\\\
&=\frac{\frac{1}{\alpha}^{\alpha^{-1}}\prod\limits_{j=1}^{N}\mu_{ij}^{Y_{ij}}}{\Gamma(\alpha^{-1})\prod\limits_{j=1}^{N}Y_{ij}!}\frac{1}{(\alpha^{-1}+\mu_{i,\centerdot})^{Y_{i,\centerdot}+\alpha^{-1}-1}}\frac{1}{\alpha^{-1}+\mu_{i,\centerdot}}\cdot\\\
&\;\;\;\int_{\lambda_{i}}\exp\\{-\lambda_{i}(\alpha^{-1}+\mu_{i,\centerdot})\\}[\lambda_{i}(\alpha^{-1}+\mu_{i,\centerdot})]^{Y_{i,\centerdot}+\alpha^{-1}-1}d[\lambda_{i}(\alpha^{-1}+\mu_{i,\centerdot})]\\\
&=\frac{\prod\limits_{j=1}^{N}\mu_{ij}^{Y_{ij}}}{\Gamma(\alpha^{-1})\prod\limits_{j=1}^{N}Y_{ij}!}\frac{\frac{1}{\alpha}^{\alpha^{-1}}}{(\alpha^{-1}+\mu_{i,\centerdot})^{Y_{i,\centerdot}+\alpha^{-1}}}\Gamma(Y_{i,\centerdot}+\alpha^{-1})\\\
&=\frac{\Gamma(Y_{i,\centerdot}+\alpha^{-1})\frac{1}{\alpha}^{\alpha^{-1}}}{\Gamma(\alpha^{-1})\prod\limits_{j=1}^{N}Y_{ij}!}(\alpha^{-1}+\mu_{i,\centerdot})^{-(Y_{i,\centerdot}+\alpha^{-1})}\exp{[\sum\limits_{j=1}^{N}Y_{ij}\log(\mu_{ij})]}\end{split}$

It gives rise to log-likelihood function for study $i$,

$\begin{split}\log
f(Y_{i1},\cdots,Y_{iN})&=\alpha^{-1}\log(\alpha^{-1})+\log\Gamma(Y_{i,\centerdot}+\alpha^{-1})-\log\Gamma(\alpha^{-1})-\sum\limits_{j=1}^{N}\log
Y_{ij}!\\\
&-(Y_{i,\centerdot}+\alpha^{-1})\log(\alpha^{-1}+\mu_{i,\centerdot})+\sum\limits_{j=1}^{N}Y_{ij}\log(\mu_{ij})\end{split}$

Therefore, using the independence between study $i$ and $i^{\prime}$ ($i\neq
i^{\prime}$), the log-likelihood of $Y_{ij}$($\forall
i=1,\cdots,M,j=1,\cdots,N$) across all studies and voxels is,

$\begin{split}l(\beta,\alpha)&=\sum\limits_{i=1}^{M}\log[f(Y_{i1},Y_{i2},\cdots,Y_{iN})]\\\
&=M\alpha^{-1}\log(\alpha^{-1})-M\log\Gamma(\alpha^{-1})+\sum\limits_{i=1}^{M}\log\Gamma(Y_{i,\centerdot}+\alpha^{-1})-\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}\log
Y_{ij}!\\\
&-\sum\limits_{i=1}^{M}(Y_{i,\centerdot}+\alpha^{-1})\log(\alpha^{-1}+\mu_{i,\centerdot})+\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}Y_{ij}\log(\mu_{ij})\end{split}$

Using the assumption that $Y_{ij}=0\text{ or }1$, so that $\log(Y_{ij}!)=0$
and $\mu_{ij}=\mu_{j}^{X}\mu_{i}^{Z}$,

$\begin{split}l(\beta,\alpha)&=M\alpha^{-1}\log(\alpha^{-1})-M\log\Gamma(\alpha^{-1})+\sum\limits_{i=1}^{M}\log\Gamma(Y_{i,\centerdot}+\alpha^{-1})-\sum\limits_{i=1}^{M}(Y_{i,\centerdot}+\alpha^{-1})\log(\alpha^{-1}+\mu_{i,\centerdot})+\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}Y_{ij}\log(\mu_{ij})\\\
&=M\alpha^{-1}\log(\alpha^{-1})-M\log\Gamma(\alpha^{-1})+\sum\limits_{i=1}^{M}\log\Gamma(Y_{i,\centerdot}+\alpha^{-1})-\sum\limits_{i=1}^{M}(Y_{i,\centerdot}+\alpha^{-1})\log(\alpha^{-1}+\mu_{i,\centerdot})\\\
&+\left(\sum\limits_{j=1}^{N}\sum\limits_{i=1}^{M}Y_{ij}\right)\left(\sum_{k=1}^{P}X_{jk}\beta_{k}\right)\\\
&=M\alpha^{-1}\log(\alpha^{-1})-M\log\Gamma(\alpha^{-1})+\sum\limits_{i=1}^{M}\log\Gamma(Y_{i,\centerdot}+\alpha^{-1})-\sum\limits_{i=1}^{M}(Y_{i,\centerdot}+\alpha^{-1})\log(\alpha^{-1}+\mu_{i,\centerdot})\\\
&+\left[\sum\limits_{j=1}^{N}Y_{\centerdot,j}\sum_{k=1}^{P}X_{jk}\beta_{k}\right]\\\
&=M\alpha^{-1}\log(\alpha^{-1})-M\log\Gamma(\alpha^{-1})+\sum\limits_{i=1}^{M}\log\Gamma(Y_{i,\centerdot}+\alpha^{-1})-\sum\limits_{i=1}^{M}(Y_{i,\centerdot}+\alpha^{-1})\log(\alpha^{-1}+\mu_{i,\centerdot})+Y_{\centerdot,}^{\top}\log(\mu^{X})\end{split}$

##  Appendix S2 Deterministic model

###  S2.1 Model factorisation: Poisson model

We consider model factorisation to replace the full $(MN)-vector$ of foci
count by sufficient statistics (dimension not larger than $M$ or $N$).
Following the total log-likelihood function in Equation 6,

$\begin{split}l&=\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}[Y_{ij}\log(\mu_{ij})-\mu_{ij}-\log(y_{ij}!)]=\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}Y_{ij}\log{\mu_{ij}}-\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}\mu_{ij}-0\\\
&=\left(\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}Y_{ij}\right)\left(\sum\limits_{k=1}^{P}X_{jk}\beta_{k}+\sum\limits_{s=1}^{R}Z_{is}\gamma_{s}\right)-\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}\mu^{X}_{j}\mu^{Z}_{i}\\\
&=\left[\sum\limits_{j=1}^{N}Y_{\centerdot,j}\sum\limits_{k=1}^{P}X_{jk}\beta_{gk}\right]+\sum\limits_{i=1}^{M}Y_{i,\centerdot}\sum\limits_{s=1}^{R}Z_{is}\gamma_{s}-\left[\sum\limits_{j=1}^{N}\mu^{X}_{j}\right]\left[\sum\limits_{i=1}^{M}\mu^{Z}_{i}\right]\\\
&=\left[\sum\limits_{j=1}^{N}Y_{\centerdot,j}\log\mu^{X}_{j}\right]+\sum\limits_{i=1}^{M}Y_{i,\centerdot}\log\mu^{Z}_{i}-\left[\boldsymbol{1}^{\top}\mu^{X}\right]\left[\boldsymbol{1}^{\top}\mu^{Z}\right]\\\
&=Y_{\centerdot,}^{\top}\log(\mu^{X})+Y_{,\centerdot}^{\top}\log(\mu^{Z})-\left[\boldsymbol{1}^{\top}\mu^{X}\right]\left[\boldsymbol{1}^{\top}\mu^{Z}\right]\end{split}$
(29)

###  S2.2 Model factorisation: NB model

Following the log-likelihood function in Equation 8,

$\begin{split}l(\beta,\alpha)&=\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}\left[\log\Gamma(Y_{ij}+\alpha^{-1})-\log\Gamma(Y_{ij}+1)-\log\Gamma(\alpha^{-1})+Y_{ij}\log(\alpha\mu_{ij})-(Y_{ij}+\alpha^{-1})\log(1+\alpha\mu_{ij})\right]\\\
&=\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}\left[\\{\sum\limits_{k=0}^{Y_{ij}-1}\log(k+\alpha^{-1})\\}-\log\Gamma(Y_{ij}+1)+Y_{ij}\log(\alpha\mu_{ij})-(Y_{ij}+\alpha^{-1})\log(1+\alpha\mu_{ij})\right]\\\
&=\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}\left[\\{\sum\limits_{k=0}^{Y_{ij}-1}\log(k+\alpha^{-1})\\}-\log\Gamma(Y_{ij}+1)+Y_{ij}\log(\alpha)+Y_{ij}\log(\mu_{ij})-(Y_{ij}+\alpha^{-1})\log(1+\alpha\mu_{ij})\right]\\\
&=\left(\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}Y_{ij}\log(\alpha^{-1})-\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}\log(1)\right)+\left(\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}Y_{ij}\right)\log(\alpha)\\\
&+\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}Y_{ij}\left(\sum\limits_{k=1}^{P}X_{jk}\beta_{k}+\sum\limits_{s=1}^{R}Z_{is}\gamma_{s}\right)-\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}(Y_{ij}+\alpha^{-1})\log(1+\alpha\mu_{ij})\end{split}$
(30)

Here, the last term
$\sum\limits_{j=1}^{N}(Y_{ij}+\alpha^{-1})\log(1+\alpha\mu_{ij})$ is
impractical to simplify, therefore, we consider moment matching method similar
as Section S1.3,
$Y_{\centerdot,j}\sim\mathrm{NB}(r^{\prime}_{j},p^{\prime}_{j})$ where

$\begin{split}r^{\prime}_{j}&=\alpha^{-1}\frac{(\sum\limits_{i=1}^{M}\mu_{ij})^{2}}{\sum\limits_{i=1}^{M}\mu_{ij}^{2}}=\alpha^{-1}\frac{(\mu^{X}_{j}\sum\limits_{i=1}^{M}\mu^{Z}_{i})^{2}}{\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}}=\alpha^{-1}\frac{(\mu^{X}_{j})^{2}(\sum\limits_{i=1}^{M}\mu^{Z}_{i})^{2}}{\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}}\\\
p^{\prime}_{j}&=\frac{\sum\limits_{i=1}^{M}\mu_{ij}^{2}}{\alpha^{-1}\sum\limits_{i=1}^{M}\mu_{ij}+\sum\limits_{i=1}^{M}\mu_{ij}^{2}}=\frac{\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}}{\alpha^{-1}\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})+\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}}=\frac{\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}}{\alpha^{-1}\mu^{X}_{j}\sum\limits_{i=1}^{M}\mu^{Z}_{i}+\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}}\end{split}$
(31)

And the parameter $\alpha^{\prime}$ of extra variance in NB approximation is

$\frac{1}{\alpha^{\prime}}=\frac{1}{\alpha}\frac{(\sum\limits_{i=1}^{M}\mu_{ij})^{2}}{\sum\limits_{i=1}^{M}\mu_{ij}^{2}}\Rightarrow\alpha^{\prime}=\frac{\sum\limits_{i=1}^{M}\mu_{ij}^{2}}{(\sum\limits_{i=1}^{M}\mu_{ij})^{2}}\alpha=\frac{\sum\limits_{i=1}^{M}(\mu^{X}_{j}\mu^{Z}_{i})^{2}}{(\mu^{X}_{j})^{2}(\sum\limits_{i=1}^{M}\mu^{Z}_{i})^{2}}\alpha$
(32)

###  S2.3 Model factorisation: clustered NB model

Following the total log-likelihood function in Section S1.5, we incorporate
the effect of study-level covariates into the Clustered NB model,

$\begin{split}l(\beta,\alpha)&=Mv\log(v)-M\log\Gamma(v)+\sum\limits_{i=1}^{M}\log\Gamma(Y_{i,\centerdot}+v)-\sum\limits_{i=1}^{M}(Y_{i,\centerdot}+v)\log(v+\mu_{i,\centerdot})+\sum\limits_{i=1}^{M}\sum\limits_{j=1}^{N}Y_{ij}\log(\mu_{ij})\\\
&=Mv\log(v)-M\log\Gamma(v)+\sum\limits_{g=1}^{B}\sum\limits_{i\in
I_{g}}\log\Gamma(Y_{i,\centerdot}+v)-\sum\limits_{g=1}^{B}\sum\limits_{i\in
I_{g}}(Y_{i,\centerdot}+v)\log(v+\mu_{i,\centerdot})\\\
&+\sum_{g=1}^{B}\left(\sum\limits_{j=1}^{N}\sum\limits_{i\in
I_{g}}Y_{ij}\right)\left(\sum_{k=1}^{P}X_{jk}\beta_{g(i)k}+\sum\limits_{s=1}^{R}Z_{is}\gamma_{s}\right)\\\
&=Mv\log(v)-M\log\Gamma(v)+\sum\limits_{g=1}^{B}\sum\limits_{i\in
I_{g}}\log\Gamma(Y_{i,\centerdot}+v)-\sum\limits_{g=1}^{B}\sum\limits_{i\in
I_{g}}(Y_{i,\centerdot}+v)\log(v+\mu_{i,\centerdot})\\\
&+\sum_{g=1}^{B}\left[\sum\limits_{j=1}^{N}Y_{gj}\sum_{k=1}^{P}X_{jk}\beta_{g(i)k}\right]+\sum\limits_{i=1}^{M}Y_{i,\centerdot}\sum\limits_{s=1}^{R}Z_{is}\gamma_{s}\\\
&=Mv\log(v)-M\log\Gamma(v)+\sum\limits_{g=1}^{B}\sum\limits_{i\in
I_{g}}\log\Gamma(Y_{i,\centerdot}+v)-\sum\limits_{g=1}^{B}\sum\limits_{i\in
I_{g}}(Y_{i,\centerdot}+v)\log(v+\mu_{i,\centerdot})\\\
&+\sum\limits_{g=1}^{B}Y_{g}^{\top}\log(\mu_{g}^{X})+Y_{,\centerdot}^{\top}\log(\mu^{Z})\end{split}$
(33)

##  Appendix S3 Statistical inference and generalised linear hypothesis
testing

###  S3.1 Contrasts on regression coefficient of study-level covariates

To investigate the effects of study-level covariates (e.g., sample size, year
of publication) on activation intensity estimation, we conduct generalised
linear hypothesis testing (GLH) on the regression coefficients $\gamma$. For
every study-level covariate $\gamma_{r},\forall r=1,\cdots,s$,

  * •

$H_{0}:C_{\gamma}\gamma=C_{\gamma}[\gamma_{1},\gamma_{2},\cdots,\gamma_{s}]^{T}=\mathbb{0}_{m\times
1}^{T}$ where $C_{\gamma}$ is the contrast matrix of size $m\times s(m\leq s)$

  * •

$H_{1}:C_{\gamma}\gamma=C_{\gamma}[\gamma_{1},\gamma_{2},\cdots,\gamma_{s}]^{T}\neq\mathbb{0}_{m\times
1}^{T}$

The covariance of regression coefficient $\gamma$,
$\mathrm{Cov}_{\gamma}=\mathrm{Cov}([\gamma_{1},\gamma_{2},\cdots,\gamma_{s}]^{T})$
is approximated from the inverse of Fisher Information matrix. According to
asymptotic normality of maximum likelihood estimator,

$\begin{split}\hat{\gamma}-\gamma&\xrightarrow{D}N\left(\mathbb{0}^{T}_{s\times
1},\mathrm{Cov}_{\gamma}\right)\\\
C_{\gamma}(\hat{\gamma}-\gamma)&\xrightarrow{D}N\left(\mathbb{0}^{T}_{s\times
1},C_{\gamma}\mathrm{Cov}_{\gamma}C_{\gamma}^{T}\right)\\\
C_{\gamma}\hat{\gamma}&\xrightarrow{D}N\left(C_{\gamma}\gamma,C_{\gamma}\mathrm{Cov}_{\gamma}C_{\gamma}^{T}\right)\end{split}$

Since a quadratic form of normal distribution has a Chi-square distribution,

$(C_{\gamma}\hat{\gamma})^{T}(C_{\gamma}\mathrm{Cov}_{\gamma}C^{T}_{\gamma})^{-1}(C_{\gamma}\hat{\gamma})\xrightarrow{D}\chi^{2}_{m}$

e.g., contrast matrix $C_{\gamma}=[1,0]$ or $[0,1]$ is for testing if the
regression coefficient of the $1^{st}$ or $2^{nd}$ study-level covariate is
zero.

###  S3.2 PP-plots of spatial homogeneity tests for each $20$ meta-analytic
datasets

Previously, we have only displayed the PP-plots of spatial homogeneity tests
of four representative datasets in Section 4.1, here we will include all PP-
plots on $20$ meta-analytic datasets in Figure S7.

Figure S7: P-P plot of $p$-value (under $-\log_{10}$ scale) with all of $20$
meta-analytic datasets, estimated by CBMR with Poisson model without study-
level covariates, sampled with model-based approach.

###  S3.3 Likelihood-based comparison between Poisson, NB and clustered NB
model

To demonstrate the likelihood-based comparison between Poisson, NB and
clustered NB model, we plot the maximised log-likelihood, AIC and BIC for each
of the $20$ meta-analytic datasets in Figure S8. We also conduct Likelihood
ratio test (LRT) to evaluate the trade-off between model sufficiency and
complexity, here, we only list the p-values of LRT between Poisson and
clustered NB model in Table S5, as $p<10^{-8}$ for LRT between Poisson and NB
model for each of the $20$ meta-analytic datasets.

Figure S8: Likelihood-based comparison (maximised log-likelihood, AIC and BIC)
within CBMR with Poisson, NB and clustered NB model, across all of $20$ meta-
analytic datasets. We found that NB model always gives the most accurate
estimation, and existence of excess variance has been justified among CBMA
data.

[]

Table S5: p-values of Likelihood Ratio test between Poisson and clustered NB
model

Dataset | p-value | Dataset | p-value | Dataset | p-value
---|---|---|---|---|---
Social Processing | $p<10^{-8}$ | PTSD | $p<10^{-8}$ | Substance Use | $p<10^{-8}$
Dementia | $p<10^{-8}$ | Cue Reactivity | $p<10^{-8}$ | Emotion Regulation | $p<10^{-8}$
Decision Making | $p<10^{-8}$ | Reward | $p<10^{-8}$ | Sleep Deprivation | $p<10^{-8}$
Naturalistic | $p<10^{-8}$ | Problem Solving | $p<10^{-8}$ | Emotion | $p<10^{-8}$
Cannabis Use | $1$ | Nicotine Use | $p<10^{-8}$ | Frontal Pole CBP | $p<10^{-8}$
Face Perception | $p<10^{-8}$ | Nicotine Administration | $0.99$ | Executive Function | $p<10^{-8}$
Finger Tapping | $0.99$ | n-Back | $p<10^{-8}$ |  |

###  S3.4 Effect of study-level covariates

Here, we investigate the effect of study-wise (square root) sample size and
year of publication (after centring and standardisation) on each of $20$ meta-
analytic datasets. Under the null hypothesis that regression coefficient of
each covariates is not distinguishable from $0$ ($\gamma_{i}=0$ for $i=1,2$),
we conduct hypothesis testing and summarise $Z-$score and $p-$value in Table
S6.

[]

Table S6: Hypothesis testing on the effect of two study-level covariates on
$20$ meta-analytic datasets

| (Square root) sample size | Year of publication
---|---|---
Dataset |  $Z-$score | p-value |  $Z-$score | p-value
Social Processing | $10.9053$ | $p<10^{-8}$ | $0.4164$ | $0.6771$
PTSD | $2.8789$ | $0.004$ | $0.6029$ | $0.5466$
Substance Use | $4.3887$ | $1.1404\times 10^{-5}$ | $6.8398$ | $p<10^{-8}$
Dementia | $20.7177$ | $p<10^{-8}$ | $-1.3985$ | $0.1620$
Cue Reactivity | $6.1454$ | $p<10^{-8}$ | $-0.6880$ | $0.4915$
Emotion Regulation | $6.8934$ | $p<10^{-8}$ | $-3.9588$ | $7.5329\times 10^{-5}$
Decision Making | $4.1104$ | $3.9499\times 10^{-5}$ | $0.1060$ | $0.9156$
Reward | $-0.1228$ | $0.9022$ | - | -
Sleep Deprivation | $12.8765$ | $p<10^{-8}$ | $0.4201$ | $0.6744$
Naturalistic | $1.7038$ | $0.0884$ | $0.5395$ | $0.5896$
Problem Solving | $4.3079$ | $1.6485\times 10^{-5}$ | $2.2789$ | $0.0227$
Cannabis Use | $3.5915$ | $3.2878\times 10^{-4}$ | $2.2117$ | $0.0270$
Nicotine Use | $5.0024$ | $5.6631\times 10^{-7}$ | $3.1836$ | $0.0015$
Frontal Pole CBP | $5.5190$ | $3.4101\times 10^{-8}$ | $7.4040$ | $p<10^{-8}$
Face Perception | $3.4090$ | $6.5212\times 10^{-4}$ | $5.1018$ | $3.3651\times 10^{-7}$
Nicotine Administration | $1.4594$ | $0.1445$ | $-1.0516$ | $0.2930$
Executive Function | $1.6989$ | $0.0932$ | $0.5047$ | $0.6138$
Finger Tapping | - | - | $0.1764$ | $0.8600$
n-Back | $1.4616$ | $0.1439$ | $0.1239$ | $0.9014$