# Random field theory-based p-values: a review of the SPM implementation

Dirk Ostwald, Sebastian Schneider, Rasmus Bruckner, Lilla Horvath
Arbeitsbereich Computational Cognitive Neuroscience, Department of Education
and Psychology,
Freie Universität Berlin, Berlin, Germany11footnotemark: 1

###### Abstract

P-values and null-hypothesis significance testing are popular data-analytical
tools in functional neuroimaging. Sparked by the analysis of resting-state
fMRI data, there has recently been a resurgence of interest in the validity of
some of the p-values evaluated with the widely used software SPM. The default
parametric p-values reported in SPM are based on the application of results
from random field theory to statistical parametric maps, a framework we refer
to as RFP. While RFP was established almost two decades ago and has since been
applied in a plethora of fMRI studies, there does not exist a unified
documentation of the mathematical and computational underpinnings of RFP as
implemented in current versions of SPM. Here, we provide such a documentation
with the aim of contributing to contemporary efforts towards higher levels of
computational transparency in functional neuroimaging.

###### keywords:

fMRI , GLM , SPM , p-values , random field theory

##  1 Introduction

Despite their debatable value in a scientific context, p-values and
statistical hypothesis testing remain popular data-analytical techniques
(Fisher,, 1925; Neyman and Pearson,, 1933; Benjamin et al.,, 2018). Given
recent discussions about the reproducibility of quantitative empirical
findings (e.g. Ioannidis,, 2005; Button et al.,, 2013), there has been a
resurgence of interest in the validity of p-values and statistical hypothesis
tests involved in the analysis of fMRI data with the popular software packages
SPM and FSL (Eklund et al.,, 2016; Mumford et al.,, 2016; Brown and Behrmann,,
2017; Eklund et al.,, 2017; Cox et al.,, 2017; Kessler et al.,, 2017; Flandin
and Friston,, 2017; Bowring et al.,, 2018; Geuter et al.,, 2018; Turner et
al.,, 2018; Slotnick, 2017b, ; Slotnick, 2017a, ; Mueller et al.,, 2017;
Eklund et al.,, 2018; Gopinath et al., 2018b, ; Gopinath et al., 2018a, ). The
default parametric p-values reported in SPM and FSL are based on the
application of results from random field theory to statistical parametric maps
and derive from the fundamental aim to control the multiple testing problem
entailed by the mass-univariate analysis of fMRI data (Friston et al., 1994a,
; Poline and Brett,, 2012; Monti,, 2011; Ashburner,, 2012; Jenkinson et al.,,
2012; Nichols,, 2012). In brief, the p-values reported by SPM and FSL are
exceedance probabilities of topological features of a data-adapted null
hypothesis random field model. Historically, this framework was established in
a series of landmark papers in the 1990’s, with some additional contributions
in the early 2000’s (Friston et al.,, 1991; Worsley et al.,, 1992; Friston et
al., 1994b, ; Worsley et al.,, 1996; Friston et al.,, 1996; Kiebel et al.,,
1999; Jenkinson,, 2000; Hayasaka and Nichols,, 2003; Hayasaka et al.,, 2004).
In the following, we will collectively refer to this framework as “random
field theory-based p-value evaluation” (RFP).

While RFP was established almost two decades ago and has since been applied in
a plethora of functional neuroimaging studies, there does not exist a unified
documentation of the mathematical and computational underpinnings of RFP. More
specifically, the early landmark papers on the approach are characterized by
an evolution of ideas which were often superseded by later developments. In
these later developments, however, the RFP framework was not constructed from
first principles again, such that many important implementational details
became rendered opaque. As concrete examples, the conceptual papers by Friston
et al., 1994b  and Friston et al., (1996) mention the expected Euler
characteristic, which is central to the current SPM implementation of RFP,
only in passing. They are also concerned with Gaussian random fields only,
rather than the T- and F-statistical random fields that are routinely assessed
using RFP. Similarly, the technical papers by Kiebel et al., (1999), Hayasaka
and Nichols, (2003), and Hayasaka et al., (2004) present variations on the
estimation of RFP’s null model parameters without a definite commitment to the
approach that is currently implemented in SPM. Finally, reviews of the
approach, such as provided by Nichols and Hayasaka, (2003), Friston, 2007b ,
Worsley, (2007), and Flandin and Friston, (2015) only cover certain aspects of
RFP and usually omit many computational details that would be required for its
de-novo implementation. Given this unsatisfying state of the RFP literature
and the recent commotion with regards to the approach, we reasoned that a
comprehensive documentation of the mathematical and computational
underpinnings of RFP may help to alleviate uncertainties about the conceptual
and implementational details of the approach for neuroimaging analysis
practitioners and theoreticians alike. Here, we thus aim to provide such a
documentation.

We chose to constrain our treatment by what probably constitutes the most
important practical aspect of the RFP approach: the SPM results table (Figure
1). We focus on the RFP implementation in SPM rather than FSL, because SPM
remains the most commonly used software tool in the neuroimaging community and
because the origins of RFP are closely linked to the development of SPM
(Carp,, 2012; Borghi and Van Gulick,, 2018; Ashburner,, 2012). More
specifically, we focus on the mathematical theory and statistical evaluation
of the p-values reported in the set-level, cluster-level, and peak-level
columns, as well as the footnote of the results table of SPM12 Version 7219.
Note that in the following, all references to SPM imply a reference to SPM12
Version 7219. We particularly address the application of RFP in the context of
GLM-based fMRI data analysis, and do not discuss other applications, such as
in the analysis of M/EEG data (e.g., Kilner and Friston,, 2010). Our guiding
example is a second-level one-sample $T$-test design as shown in Figure 1.
With respect to multiple testing control, we focus on the family-wise error
rate (FWER)-corrected p-values and do not cover the later addition of false
discovery rate-corrected q-values (Chumbley and Friston,, 2009; Chumbley et
al.,, 2010). With respect to FWER control, we focus on the evaluation of
p-values rather than statistical hypothesis tests. This is warranted by the
fact that SPM de-facto only evaluates p-values, while statistical hypothesis
tests are performed by the SPM user, who typically decides to report only
activations at a given level of significance. Finally, from an
implementational perspective, we primarily focus on the functionality of the
spm_list.m, spm_P_RF.m, spm_est_smoothness.m, and spm_resel_vol.c functions
and cover computational aspects of more low-level SPM routines only in
passing. With these constraints in place, we next review the outline of our
treatment and preview the central points of each section.

Figure 1: The SPM results table for a second-level one-sample $T$-test design
relating to the positive main effect of visual stimulus coherence in a
perceptual decision making fMRI study (Georgie et al.,, 2018). The central aim
of this work is to document the mathematical and computational underpinnings
of the p-values listed in the SPM results table. As indicated by the whitened
parts of the table, we only consider corrected p-values related to FWER
control. For implementational details of the evaluation of the SPM results
table using the SPM12 Version 7219 distribution, please see rfp_1.m.

#### Overview

Overall, we proceed in three main parts, Section 2: Background, Section 3:
Theory, and Section 4: Application. In Section 2: Background, we establish a
few key definitions from the two main pillars of RFP: random field theory and
geometry (Adler,, 1981; Adler and Taylor,, 2007). Specifically, in Section
2.1: Real-valued random fields, we define the concept of a real-valued random
field and discuss some subtle but important aspects of this definition. Our
definition rests on mathematical probability theory, in which random variables
are understood as mappings between probability spaces. For readers unfamiliar
with these measure-theoretic underpinnings of probability theory, we provide a
short synopsis of the required concepts in Supplement S1. After discussing a
few analytical features of random fields (expectation and covariance
functions, stationarity), we introduce the central classes of random fields
for RFP in Section 2.2: Gaussian and statistical random fields. As a toy
example that we will return to repeatedly, we introduce two-dimensional
Gaussian random fields with Gaussian covariance functions. Such Gaussian
random fields have the feature that their smoothness in the RFP-sense can be
captured by a single parameter of their covariance function. We also introduce
the central notions of $Z\text{-},T\text{-}$ and $F\text{-}$ fields, which we
will often refer to as statistical random fields and consider as the
theoretical analogues to statistical parametric maps that derive from the
practical analysis of fMRI data. In Section 2.3: Volumetrics we then introduce
two geometric foundations of RFP, Lebesgue measure and intrinsic volumes.
These geometric concepts are of central importance in RFP, because the
probability distributions of topological features of statistical random fields
(for example, the probability of the global maximum to exceed a given value)
depend not only on the characteristics of the statistical random field, but
also on the volume of the space it extends over. In RFP, this dual dependency
is condensed in the notion of resel volumes. In brief, resel volumes are
smoothness-adjusted intrinsic volumes, which in turn are coefficients in a
formula for the Lebesgue measure of tubes. The treatment in Section 2.3 is
necessarily shallow and emphasizes intuition over formal rigour.

Section 3: Theory is devoted to the theoretical development of RFP. The
central part of this section is the delineation of the parametric probability
distributions used by RFP to capture the stochastic behaviour of topological
features of a statistical random field’s excursion set under the null
hypothesis. The exceedance probabilities of these parametric distributions for
observed data are the p-values reported in the SPM results table. We discuss
the theory of RFP against the background of a continuous space, discrete data
point model. This model is conceived as the theoretical analogue to the
familiar mass-univariate GLM and is formally introduced at the outset of
Section 3. Section 3.1: Excursion sets, smoothness, and resel volumes then
introduces three foundational concepts of RFP: excursion sets of statistical
random fields are subsets of the field’s domain on which the field exceeds a
pre-specified value $u$. In the applied neuroimaging literature, this value
$u$ (or, more specifically, its equivalent p-value) is known as the cluster-
forming threshold (Nichols et al.,, 2016) and has been at the center stage of
the recent debate on the validity of parametric p-values in fMRI data analysis
(e.g., Flandin and Friston,, 2017). The topological constitution of excursion
sets, and hence the probability distributions of their features, depends in a
predictable manner on the value of $u$ and on the smoothness of the underlying
statistical random field. In the applied neuroimaging literature, the
smoothness measure employed by RFP is known as “FWHM” (e.g., Eklund et al.,,
2016). In fact, the parameterization of a statistical random field’s
smoothness in terms of the full widths at half maximum (FWHMs) of isotropic
Gaussian convolution kernels applied to hypothetical white-noise Gaussian
random fields conforms to a reparameterization of a more fundamental
smoothness measure and was introduced early on in the RFP literature (Friston
et al.,, 1991; Worsley et al.,, 1992). We introduce this more fundamental
smoothness measure in the Smoothness subsection of Section 3.1. To obtain an
intuition about this measure, we show analytically in Supplement S3 that for
Gaussian random fields with Gaussian covariance functions it evaluates to a
scalar multiple of the covariance function’s length parameter. We close
Section 3.1 by introducing the FWHM reparameterization of smoothness and
finally conjoining the probabilistic and geometric threads of our exposition
in the notion of resel volumes. Equipped with these foundations, we then
proceed to the core of RFP theory in Section 3.2: Probabilistic properties of
excursion set features. More specifically, upon establishing a set of expected
values, we review the probability distributions of the following excursion set
features as evaluated by the SPM implementation of RFP: (1) the global
maximum, (2) the number of clusters, (3) the cluster volume, (4) the number of
clusters of a given volume, and (5) the maximal cluster volume. We supplement
this discussion with a review of the relevant parametric distributions in
Supplement S1 and number of proofs in Supplement S3. The distributions of the
global maximum and the maximal cluster volume are the distributions that endow
the SPM implementation of RFP with FWER control at the peak- and cluster-
level, respectively. To this end, we review the principal idea of controlling
the FWER by means of maximum statistics in Supplement S2. With all theoretical
aspects in place, we are then in the position to consider the application of
RFP theory in a GLM-based fMRI data-analytical setting.

We thus begin Section 4: Application by reconsidering the continuous space GLM
introduced in Section 3 from the discrete spatial sampling perspective of
fMRI. In Section 4.1: Parameter estimation we then discuss how RFP furnishes a
data-adaptive null hypothesis model by estimating the smoothness and
approximating the intrinsic volumes of the statistical random field that
underlies the observed data. The estimated FWHM smoothness parameters and
intrinsic volumes are combined in estimated resel volumes, which form the
pivot point between the data-driven and theory-based aspects of RFP. We then
close our treatment by considering the p-values reported in the SPM results
table. As will become evident, the p-values reported for set-, cluster-, and
peak-level inferences can actually be evaluated using a single (Poisson)
cumulative distribution function that is adapted for different topological
statistics and data-based null model characteristics (Friston et al.,, 1996).
Finally, in Section 5: Discussion, we briefly explore potential future avenues
for the further refinement of RFP.

In summary, we make the following novel contributions to the literature. From
a scientific perspective, we provide a unified review of RFP, which is both
mathematically comprehensive and computationally explicit. As such, our
treatment allows for the detailed statistical interpretation of the
experimental effects reported using RFP in the last decades. Further, given
the recent discussions on the validity of SPM’s cluster-level p-values, our
treatment has the potential to serve as a readily accessible starting point
for the further refinements of RFP. Finally, from an educational perspective,
we provide a novel resource to introduce newcomers to the field of
computational cognitive neuroscience to one of functional neuroimaging’s
equally most basic (univariate cognitive process mapping using the GLM) and
sophisticated (random field-based modelling) analysis tools. All custom
written Matlab code (The MathWorks, NA) implementing simulations and
visualizations reported on herein is available from https://osf.io/3dx9w/. The
code repository also contains the group-level fMRI data evaluated in Figure 1,
and documented and revised versions of SPM’s parameter estimation and p-value
evaluation routines (rfp_spm_est_smoothness.m and rfp_spm_table.m,
respectively). For details, please refer to the OSF project documentation.

#### Prerequisites and notation

We assume throughout that the reader is familiar with the theoretical and
practical aspects of GLM-based fMRI analysis as presented for example in
Kiebel and Holmes, (2007); Monti, (2011); Poline and Brett, (2012) and
Poldrack et al., (2011). We presume that a familiarity with basic concepts of
the multiple testing problem entailed by the mass-univariate GLM-based
analysis of fMRI data, as covered for example in Brett et al., (2007);
Nichols, (2012) and Ashby, (2011) is beneficial. As mentioned above, we
provide a synopsis of essential aspects from mathematical probability theory
and the theory of multiple testing in Supplement S1 and Supplement S2,
respectively. In Table 1, we list mathematical symbols that will be used
throughout.

Symbol | Meaning
---|---
$\mathbb{N}^{0}$, $\mathbb{N}_{m}$, $\mathbb{N}_{m}^{0}$ |  The sets $\mathbb{N}\cup\\{0\\}$, $\\{1,2,...,m\\}$, and $\\{0,1,...,m\\}$, respectively
$|S|$ |  Cardinality of a set $S$
$e_{i}$ |  $i$th standard basis vector for $\mathbb{R}^{d}$
$||x||_{2}$ |  Euclidean norm of $x\in\mathbb{R}^{d}$
$|A|$ |  Determinant of a matrix $A\in\mathbb{R}^{d\times d}$
p.d. | Positive-definite
$\mbox{inf},\mbox{sup}$ | Infimum, supremum
$\nabla f(x)$ |  Gradient of a function $f:\mathbb{R}^{d}\to\mathbb{R}$ at $x$
$\frac{\partial}{\partial x_{i}}f(x)$ |  $i$th partial derivative of a function $f:\mathbb{R}^{d}\to\mathbb{R}$ at $x$, $i=1,...,d$
$\Gamma(x)$ |  Gamma function evaluated at $x\in\mathbb{R}$
$\mathbb{P},\mathbb{E},\mathbb{C}$ | Probability, expectation, covariance
$N(\mu,\Sigma)$ | Gaussian distribution
$\bullet$ | End of definition
$\Box$ | End of proof
Table 1: List of mathematical symbols. Note that the meaning of the symbol
$|\cdot|$ is usually clear from the context.

##  2 Background

###  2.1 Real-valued random fields

We commence by discussing some basic aspects of real-valued random fields. For
comprehensive introductions to the theory of random fields, see for example
Christakos, (1992), Abrahamsen, (1997), Adler, (1981) and Adler and Taylor,
(2007) and for the required background material from mathematical probability
theory, see Supplement S1. We use the following definition:

######  Definition 1 (Real-valued random field).

A real-valued random field $\\{X(x)|x\in S\\}$ on a domain
$S\subset\mathbb{R}^{D}$ is a set of random variables $X(x)$ on a probability
space $(\Omega,\mathcal{A},\mathbb{P})$, i.e., for each $x\in
S,X(x):\Omega\to\mathbb{R}$ is a random variable.

$\hfill\bullet$

A variety of notations for real-valued random fields exist. In the RFP
literature, random fields are most commonly denoted by “$X(x),x\in S$” (e.g.,
Worsley,, 1994; Taylor and Worsley,, 2007; Flandin and Friston,, 2015), a
convention which we shall also follow herein. With regards to this notation,
it is important to realize that the symbol $X(x)$ denotes a random variable,
while it may look suspiciously like the value of a function $X$ evaluated for
an input argument $x$. This is of course intentional: because each $X(x)$ of a
real-valued random field is a real-valued random variable, it can be written
as

$X(x):\Omega\to\mathbb{R},\omega\mapsto X(x)(\omega).$ (1)

In expression (1), $X(x)(\omega)$ denotes the value in $\mathbb{R}$ that
$X(x)$ takes on for the input argument $\omega\in\Omega$. Because the double
bracket notation $X(x)(\omega)$ is rather unconventional, an alternative
notation mainly used in the mathematical literature is to write a real-valued
random field as

$X:S\times\Omega\to\mathbb{R},(x,\omega)\mapsto X(x,\omega),$ (2)

or to denote the constituent random variables of a real-valued random field by

$X_{x}:\Omega\to\mathbb{R},\omega\mapsto X_{x}(\omega).$ (3)

The three different notations of expressions (1), (2) and (3) make the same
point: for a fixed $\omega$, $X(x)$, $X$, or $X_{x}$ is a deterministic, real-
valued function with domain $S$. The notation $X(x)$ suppresses the dependency
of this function on the random elementary outcome $\omega$, while the other
two notations do not. The clearest notation in this regard is perhaps offered
by expression (2), but as mentioned above, the notation $X(x),x\in S$ seems to
be preferred in the RFP literature. Crucially, and these notational subtleties
aside, the domain $S$ of random fields is commonly an uncountable infinite set
(such as a $D$-dimensional interval), which implies that a random field
commonly comprises an uncountable infinite number of random variables.

Two further aspects of Definition 1 are worth discussing. First, we chose the
same letter for both the random field (capital $X$) and the domain values
(lowercase $x$), because on the one hand, we will need to define specific
random fields later on, and $X$ is for the current purposes perhaps the most
generic choice, while on the other hand, $x$ invokes a clearer notion of a
spatial domain than, say, $t$. In fact, $t$ is a popular letter for the
elements of the domain of a random field (e.g., Worsley, (1994)). However,
because random fields will be used in the context of RFP for GLM-based fMRI
data analysis as spatial models, we prefer $x$. Second, we defined the domain
$S$ to be a subset of $\mathbb{R}^{D}$. We use the letter $S$, because in the
context of RFP, the domains of random fields of interest are typically
referred to as search spaces. Intuitively, the search spaces of interest in
GLM-based fMRI data analysis correspond to the brain (in whole-brain analyses)
or brain regions of interest (in small volume correction analyses). The $D$ we
have primarily in mind is $D=3$, i.e., the case of random fields over three-
dimensional space, such as a volume of observed test statistics. For
visualization purposes we will also consider the case $D=2$. In the case of
$D=1$, random fields are commonly referred to as random processes, or perhaps
even more often as stochastic processes. Finally, a linguistic remark:
because, as the pinnacle of mass-univariate GLM-based fMRI data analysis, RFP
concerns only univariate statistical values, we are in fact only concerned
with real-valued (as opposed to, say, vector-valued) random fields. We will
therefore use the terms “real-valued random field” and “random field”
interchangeably henceforth.

#### Expectation, covariance, and variance functions

Like random variables, random fields have certain analytical features that
express their expected (or “average”) behaviour over many realizations. These
features, which generalize the concepts of a random variable’s expectation and
variance, are known as expectation, covariance, and variance functions. We use
the following definitions for these entities (Abrahamsen,, 1997):

######  Definition 2 (Expectation, covariance, and variance functions).

For a real-valued random field $X(x),x\in S$, the expectation function is
defined as

$m:S\to\mathbb{R},x\mapsto m(x):=\mathbb{E}(X(x)),$ (4)

the covariance function is defined as

$c:S\times S\to\mathbb{R},(x,y)\mapsto c(x,y):=\mathbb{C}(X(x),X(y))$ (5)

and the variance function is defined as

$s^{2}:S\to\mathbb{R}_{\geq 0},x\mapsto s^{2}(x):=c(x,x).$ (6)

$\hfill\bullet$

Note that the expectation and covariance expressions (4) and (5) are evaluated
for the random variables $X(x)$ and $X(x),X(y)$, respectively, with regards to
the probability measure of the underlying probability space. Further note that
the variance function is defined directly in terms of the expectation and
covariance functions.

#### Stationarity

An important property of real-valued random fields, and a fundamental
assumption for the random fields considered in RFP, is stationarity. We define
stationarity in the wide sense as follows (Abrahamsen,, 1997):

######  Definition 3 (Wide-sense stationarity).

A real-valued random field $X(x),x\in S$ is called stationary in the wide
sense, if the expectation function (4) of $X(x)$ is a constant function,

$m:S\to\mathbb{R},x\mapsto m(x):=\bar{m}\mbox{ for }\bar{m}\in\mathbb{R}$ (7)

and the covariance function of $X(x)$ is a function of the separation of its
arguments only, i.e., for all $x,y\in S$ and $d:=x-y$, the covariance function
(5) can be written as

$c:S\to\mathbb{R},d\mapsto c(d).$ (8)

$\hfill\bullet$

In addition to stationarity in the wide sense, there exists the notion of
stationarity in the strict sense (Abrahamsen,, 1997). Stationarity in the
strict sense means that all finite-dimensional distributions of a real-valued
random field are invariant under arbitrary translations of their support
arguments in the domain of the real-valued field. Because stationarity in the
wide and the strict sense are equivalent for Gaussian and Gaussian-related
random fields, which are our primary concern herein, and because the concept
of wide-sense stationarity is more readily applicable in the computational
simulation of random fields, we only formalize stationarity in the wide sense.
Note, however, that in general, stationarity in the strict sense implies
stationarity in the wide sense, but not vice versa. Henceforth, we shall only
use the term “stationarity” to mean stationarity in the wide sense. Finally,
note that stationary random fields are also sometimes referred to as
“homogeneous” random fields (e.g., Adler, (1981), Worsley, (1995)).

The role of stationarity for RFP is fundamental, because, intuitively, it
corresponds to the assumption that the null hypothesis holds true at every
location of the search space $S$, and the null hypotheses over the entire
search space are identical. More explicitly, defining $\bar{m}:=0$ for a
stationary real-valued random field corresponds to the assumption that the
expected value of the random variable at each location is $0$, which, from the
perspective of statistical testing is identical with the expectation of no
effect, i.e., the null hypothesis.

###  2.2 Gaussian and statistical random fields

#### Gaussian random fields

We next introduce a central class of random fields for RFP, Gaussian random
fields, and use these to make the abstract definitions of random fields, their
associated expectation, covariance, and variance functions, as well as the
concept of stationarity more concrete.

######  Definition 4 (Gaussian random field).

A Gaussian random field (GRF) on a domain $S\subset\mathbb{R}^{D}$ is a random
field $G(x),x\in S$, such that every finite-dimensional vector of random
variables

$g:=(G(x_{1}),G(x_{2}),...,G(x_{n}))^{T}$ (9)

is distributed according to a multivariate Gaussian distribution.

$\hfill\bullet$

This definition implies that the random vector $g$ is distributed according to
a multivariate Gaussian distribution for every choice of the $x_{i},i=1,...,n$
and $n\in\mathbb{N}$, such that we can write

$g\sim N(\mu,\Sigma)\mbox{ for
}\mu\in\mathbb{R}^{n},\Sigma\in\mathbb{R}^{n\times n}\mbox{ p.d. }.$ (10)

Moreover, if we denote the expectation function and the covariance function of
a GRF by $m$ and $c$, respectively, the expectation parameter
$\mu\in\mathbb{R}^{n}$ and covariance matrix parameter
$\Sigma=(\Sigma_{ij})\in\mathbb{R}^{n\times n}$, p.d. in expression (10) have
the entries

$\mu_{i}=m(x_{i})\mbox{ and }\Sigma_{ij}=c(x_{i},x_{j})\mbox{ for
}i,j=1,...,n.$ (11)

Notably, expressions (10) and (11) achieve two things: first, they allow for
reducing the fairly abstract concept of a random field to the multivariate
Gaussian distribution as a more familiar entity. Second, they immediately
imply a computational approach to obtain realizations from a GRF: if one
defines a set of discrete support points $x_{1},...,x_{n}$ in the domain of a
GRF of interest and evaluates the expectation and covariance functions of the
GRF at these points, one obtains an expectation parameter and a covariance
matrix that can be used as parameters for a multivariate Gaussian distribution
random vector generator. Of course, this requires that the expectation and
covariance functions of the GRF are defined as well.

Defining appropriate covariance functions and deciding whether a given
function is actually a suitable covariance function, such that the
corresponding covariance matrix is positive-definite, is a mathematical
problem on its own (see e.g., Rasmussen and Williams, (2006, Chapter 4),
Abrahamsen, (1997)), and we will not delve further into this issue here. We
do, however, provide an important example of a valid covariance function, the
so-called Gaussian covariance function (GCF) given by (e.g., Christakos,
(1992, p. 71) and Powell et al., (2014, p. 5))

$\gamma:S\times
S\to\mathbb{R}_{>0},(x,y)\mapsto\gamma(x,y):=v\exp\left(-\frac{||x-y||_{2}^{2}}{\ell^{2}}\right)\mbox{
with }v,\ell>0.$ (12)

Note that the variance function of a GRF with a GCF evaluates to

$s^{2}(x)=\gamma(x,x)=v,$ (13)

and that the GCF is de-facto a function of the scalar Euclidean distance

$\delta:=||x-y||_{2}\in\mathbb{R}_{\geq 0}$ (14)

between $x$ and $y$ only. It may thus equivalently be expressed as

$\tilde{\gamma}:\mathbb{R}_{\geq
0}\to\mathbb{R}_{>0},\delta\mapsto\tilde{\gamma}(\delta):=v\exp\left(-\frac{\delta^{2}}{\ell^{2}}\right)\mbox{
with }v,\ell>0,$ (15)

for which $v=\tilde{\gamma}(0)$. The parameter $\ell$ in eqs. (12) and (15) is
a length constant and determines the strength of the covariation of the random
variable $X(x)$ and $X(y)$ at a distance $\delta=||x-y||_{2}$ (Figure 2A and
B). Intuitively, the larger the value of $\ell$, the stronger the covariation
at a given distance, and hence the less variable the profile of GRF
realizations over space. In Figure 2C, we visualize the dependency of
realizations of GRFs with GCF on the parameter $\ell$, for GRFs with domain
$S=[0,1]\times[0,1]\subset\mathbb{R}^{2}$ and constant zero mean function.
Note that the GCF is called “Gaussian” because of its functional form, and not
because we consider it in the context of Gaussian random fields. Further note
that for a constant mean function, GRFs with GCF are stationary. Finally, note
that GRFs with $D=1$ are known as Gaussian processes and enjoy some popularity
in the machine learning community (Rasmussen,, 2004; Rasmussen and Williams,,
2006).

Figure 2: Realizations of Gaussian random fields with Gaussian covariance
functions (GCFs). (A) Panel A visualizes a GCF for a random field with two-
dimensional domain $D=2$ of the form of expression (12) with $y=(0,0)^{T}$
with parameters $v=1$ and $\ell=1$. (B) Panel B visualizes GCFs for a random
field with two-dimensional domain $D=2$ of the form of expression (15) with
parameter $v=1$ and varying length constants $\ell$. Note that the height of
$\tilde{\gamma}$ at a distance $\delta$ indicates the amount of covariation
between two random variables at distance $\delta$ in the random field. (C)
Panel C depicts nine realizations of a Gaussian random field with domain
$S=[0,1]\times[0,1]$ with $2^{6}$ support points in each dimension,
corresponding to spatially arranged samples from a $2^{6}$ dimensional
multivariate Gaussian distribution of the form given by expressions (10) and
(11). Each realization is based on a Gaussian random field with Gaussian
covariance function of varying length parameter $\ell$, as indicated in the
subpanel titles. Note that for larger values of $\ell$, and hence larger
covariation of two arbitrary random variables of the field at distance
$\delta$, the realizations show less variability over space, and appear
“smoother” than for smaller values of $\ell$. For the full implementational
details of these simulations, please see rfp_2.m.

#### Statistical random fields

Three additional random fields that derive from Gaussian random fields are of
central importance in RFP. These fields, often referred herein under the
umbrella term statistical random fields, can be considered the random field
analogues of the $Z\text{-},T\text{-},$ and $F\text{-}$distributions (Shao,,
2003). Following Worsley, (1994) and Adler, (2017), we define these fields as
follows:

######  Definition 5 (Statistical random fields: $Z\text{-}$,$T\text{-}$, and
$F\text{-}$fields).

A $Z\text{-}$field on a domain $S\subset\mathbb{R}^{D}$ is a stationary
Gaussian random field with constant expectation function $m(x):=0$ and
variance function $s^{2}(x)=1$. Let

$Z_{1}(x),...,Z_{n}(x),x\in S$ (16)

be $n$ independent $Z\text{-}$fields, i.e.,

$\mathbb{C}(Z_{i}(x),Z_{j}(x))=0\mbox{ for all }1\leq i,j\leq n,i\neq j,x\in
S$ (17)

and let

$\mathbb{C}(\nabla Z_{i}(x))=\mathbb{C}(\nabla Z_{j}(x))\mbox{ for all }1\leq
i,j\leq n,i\neq j,x\in S.$ (18)

Then

$T(x):=\frac{Z_{1}(x)\sqrt{n-1}}{\sqrt{\sum_{i=2}^{n}Z_{i}^{2}(x)}},x\in S$
(19)

is called a $T\text{-}$field with $n-1$ degrees of freedom. Finally, let

$Z_{1}(x),...,Z_{n}(x),Z_{n+1}(x),...,Z_{n+m}(x),x\in S$ (20)

be $n+m$ independent $Z$-fields, i.e.,

$\mathbb{C}(Z_{i}(x),Z_{j}(x))=0\mbox{ for all }1\leq i,j\leq n+m,i\neq j,x\in
S$ (21)

and let

$\mathbb{C}(\nabla Z_{i}(x))=\mathbb{C}(\nabla Z_{j}(x))\mbox{ for all }1\leq
i,j\leq n+m,i\neq j,x\in S.$ (22)

Then

$F(x):=\frac{m\sum_{i=1}^{n}Z_{i}^{2}(x)}{n\sum_{i=n+1}^{n+m}Z_{i}^{2}(x)},x\in
S$ (23)

is called an $F\text{-}$field with $n,m$ degrees of freedom.

$\hfill\bullet$

In the context of RFP, we conceive of $Z\text{-},T\text{-}$ and
$F\text{-}$fields as the theoretical, continuous-space analogues of
statistical parametric maps, i.e., discrete-space spatial maps of realized
$Z\text{-}$, $T\text{-}$ and $F\text{-}$statistics (cf. Friston, 2007a, ).

###  2.3 Volumetrics

A fundamental aspect of RFP is the fact that the probability distributions of
topological features of random fields, such as the maximum of a random field
to exceed a given threshold value, depend on (1) the stochastic
characteristics of the random field and (2) the size of the search space over
which the random field extends. Intuitively, the larger the space over which
the random field extends and the more variable the random field per unit
space, the higher the probability that, for example, the maximum of the field
exceeds a given threshold value. Because random fields are mathematically
developed over continuous space, the question of the size of a subset of a
random field’s domain is not trivial and rests on a large body of mathematical
work that falls into the realms of differential and integral geometry (Adler
and Taylor,, 2007, Part II). We here provide a bare minimum of terminology for
describing volumes in continuous space by focussing on the notions of Lebesgue
measure and intrinsic volumes. Lebesgue measure is a fundamental volume
measure in mathematical measure theory and forms the basis for the concept of
intrinsic volumes. Intrinsic volumes in turn are the foundation for the
concept of resel volumes as discussed in Section 3.1: Excursion sets,
smoothness, and resel volumes.

#### Lebesgue measure

Lebesgue measure is a fundamental building block of mathematical measure
theory. Intuitively, Lebesgue measure can be understood as deriving from the
desire to allocate a meaningful notion of volume to geometric objects that are
modelled as subsets of $\mathbb{R}^{D}$. Specifically, the measure should (a)
be compatible with the usual notion of length, area, and volume for lines,
rectangles, and cuboids, respectively, (b) be smaller for an object A than for
an object B, if A can be fit into object B, (c) be invariant under
translations of the object, and (d) be additive in the sense that if two
objects do not overlap, their joint volume is the sum of their individual
volumes (Meisters,, 1997). Lebesgue measure is a measure that fulfils these
properties for a large class of subsets of $\mathbb{R}^{D}$. A full
development of Lebesgue measure from first principles is beyond our current
scope and can be found in many books (e.g., Billingsley,, 1978; Cohn,, 1980;
Stein and Shakarchi,, 2009). Instead, we here provide a general definition of
Lebesgue measure based on Hunter, (2011), discuss some of the intuition
associated with this definition, and finally list some values for the Lebesgue
measure of some familiar geometric objects.

######  Definition 6 (Lebesgue outer measure, Lebesgue measure).

Let

$R=[a_{1},b_{1}]\times\dots\times[a_{D},b_{D}]\subset\mathbb{R}^{D},\mbox{
with }-\infty<a_{d}\leq b_{d}<\infty,\,d=1,...,D$ (24)

denote a $D$-dimensional closed rectangle with sides oriented parallel to the
coordinate axes, let $\mathcal{R}\big{(}\mathbb{R}^{D}\big{)}$ denote the set
of all such rectangles in $\mathbb{R}^{D}$, and let

$\varphi:\mathcal{R}\big{(}\mathbb{R}^{D}\big{)}\to[0,\infty[,R\mapsto\varphi(R):=\prod_{d=1}^{D}(b_{d}-a_{d})$
(25)

denote the volume of a rectangle in $\mathbb{R}^{D}$. Then the Lebesgue outer
measure is defined as

$\lambda^{*}:\mathcal{P}\big{(}\mathbb{R}^{D}\big{)}\to[0,\infty],S\mapsto\lambda^{*}(S):=\inf\left\\{\sum_{i=1}^{\infty}\varphi(R_{i})|S\subset\cup_{i=1}^{\infty}R_{i},R_{i}\in\mathcal{R}\big{(}\mathbb{R}^{D}\big{)}\
\right\\},$ (26)

where the infimum is taken over all countable collections of rectangles whose
union contains $S$, a set $S\subset\mathbb{R}^{D}$ is called Lebesgue
measurable, if for every $A\subset\mathbb{R}^{D}$

$\lambda^{*}(A)=\lambda^{*}(A\cap S)+\lambda^{*}(A\textbackslash S),$ (27)

and the Lebesgue measure is defined as the function

$\lambda:\mathcal{P}\big{(}\mathbb{R}^{D}\big{)}\to[0,\infty],S\mapsto\lambda(S):=\lambda^{*}(S).$
(28)

$\hfill\bullet$

The definition of Lebesgue measure thus rests on a familiar concept: the
measure of a $D$-dimensional rectangle corresponds to the product of its side
lengths. For $D=1$, this corresponds to the length of a line, for $D=2$ to the
area of a rectangle (Figure 3A), and for $D=3$ to the volume of a cuboid. This
familiar measure of the content of $D$-dimensional rectangles is used in
Definition 6 to define an approximation of the measure of arbitrary geometric
objects modelled by a subset $S\subset\mathbb{R}^{D}$ in expression (26).
Specifically, the Lebesgue outer measure of $S$ is given by finding the
smallest possible value of the countable sum of rectangle volumes, the
rectangles of which cover the set $S$ of interest (Figure 3B). Lebesgue
measure itself is then defined as the Lebesgue outer measure restricted to a
subset of subsets of $\mathbb{R}^{D}$ which fulfil the condition of eq. (27).

The definition of Lebesgue measure by Definition 6 is fairly abstract. For
concrete subsets of $\mathbb{R}^{D}$ of geometrical interest, for example a
circle in $\mathbb{R}^{2}$ or a cuboid in $\mathbb{R}^{3}$, it is not
immediately clear how to compute their Lebesgue measure based on this
definition. Lebesgue measure has, however, many mathematically desirable
properties. For example, it can be shown that Lebesgue measure indeed fulfils
the desiderata (a) - (d) discussed above, and, for $D$-dimensional rectangles,
simplifies to the volume of rectangles (25). Some other noteworthy values of
Lebesgue measure are the Lebesgue measures of intervals in $\mathbb{R}$,

$\lambda(]a,b[)=\lambda([a,b[)=\lambda(]a,b])=\lambda([a,b])=b-a,$ (29)

and the Lebesgue measure of a closed ball $B$ with radius $r>0$ in
$\mathbb{R}^{D}$,

$\lambda(B)=\lambda\left(\\{x\in\mathbb{R}^{D}|\,||x||_{2}\leq
r\\}\right)=\frac{\pi^{D/2}}{\Gamma(D/2+1)}r^{D}.$ (30)

Finally, Lebesgue measure forms the foundation for the notion of the intrinsic
volumes of a subset of $\mathbb{R}^{D}$ as discussed next.

#### Intrinsic volumes

In general, intrinsic volumes can be thought of as measures of the $0$\- to
$D$-dimensional size of a set $S\subset\mathbb{R}^{D}$. Historically, Worsley
et al., (1992) introduced the RFP framework for three-dimensional excursion
sets only, for which it was assumed that the probability of them to touch the
search space boundary was negligible. For these excursion sets, Lebesgue
measure as a measure of volume was sufficient. Worsley et al., (1996) then
introduced corrections for these “boundary effects” which require more refined
measures of volume. These more refined measures of volume are furnished by
intrinsic volumes, which appear in the mathematical literature under a variety
of names, such as Quermassintegrale, Minkowski functionals, Dehn-Steiner
functionals, integral curvatures, and Euclidean Lipschitz-Killing curvatures
(Adler and Taylor,, 2007; Adler et al.,, 2015; Adler,, 2017). As for the
Lebesgue measure, a development of the concept of intrinsic volumes from first
principles is beyond our scope. We thus content with a general definition of
intrinsic volumes, for which we attempt to provide some intuition, and with
listing some values for the intrinsic volumes of familiar geometric entities.

As a general definition of the intrinsic volumes of a $D$-dimensional search
volume $S\in\mathbb{R}^{D}$, we use the following implicit definition provided
by Taylor and Worsley, (2007, Appendix A.1.1) and Adler et al., (2015, pp. 99
- 101)):

######  Definition 7 (Tube, intrinsic volumes).

Let $S\subset\mathbb{R}^{D}$ and let a tube $T_{S,r}$ of radius $r$ around $S$
be defined as

$T_{S,r}:=\\{x\subset\mathbb{R}^{D}|\delta(x,S)\leq r\\}\mbox{ with
}\delta(x,S):=\inf_{y\in S}||x-y||_{2}.$ (31)

Then, for a convex set $S\in\mathbb{R}$ the Lebesgue measure of $T_{S,r}$ is
given by

$\lambda(T_{S,r})=\sum_{i=0}^{D}b_{D-i}r^{D-i}\mu_{d}(S)\mbox{ with
}b_{j}:=\frac{\pi^{j/2}}{\Gamma(j/2+1)}\mbox{ for }j=D-i\mbox{ and
}i=0,...,D.$ (32)

and $\mu_{0}(S),\mu_{1}(S),...,\mu_{D}(S)$ are the intrinsic volumes of $S$.

$\hfill\bullet$

Note that the definition of the intrinsic volume $\mu_{0}(S),...,\mu_{D}(S)$
of a set $S\subset\mathbb{R}^{D}$ is implicit here in the sense that the
quantities $\mu_{0}(S),...,\mu_{D}(S)$ appear as coefficients in the sum
formula for the Lebesgue measure of a tube, and no direct way of how to
evaluate these coefficients is provided. Also note that the $b_{j}$ correspond
to the Lebesgue measures of unit balls in $\mathbb{R}$ (cf. (30) with $r:=1$).
To gain some insight into the meaning of this implicit definition of intrinsic
volumes, we follow (Adler et al.,, 2015, pp. 100-101) and discuss its
application to a two-dimensional tube around a triangle (Figure 3C). Here, the
triangle corresponds to the set $S\subset\mathbb{R}^{2}$ in eq. (31). A two-
dimensional tube $T_{S,r}$ of radius $r$ around $S$ is then formed by
considering all points $x\in\mathbb{R}^{2}$ for which the distance
$\delta(x,S)$ between the point $x$ and the subset $S$ is smaller than $r\geq
0$. The characteristic shape of a tube is afforded by the measure of distance
$\delta(x,S)$: for its evaluation, all points $y\in S$ are considered, their
Euclidean distance to $x$ evaluated, and the distance between $x$ and $S$ then
corresponds to the smallest distance between any $y\in S$ and $x$. This can be
made more concrete by considering the tube around the triangle in Figure 3C.
First, consider the three sides of the triangle. Here, the smallest distance
between a point $x\in\mathbb{R}^{2}$ and a point $y\in S$ is given by moving
away from the triangle in a perpendicular manner. However, at the three
corners of the triangle, the smallest distance from the cornerpoint is given
by the corresponding arc of the circle of radius $r$ centred on the triangle’s
corner. By inspection of Figure 3C it can be inferred that the Lebesgue
measure of the tube $T_{S,r}$, which in this two-dimensional scenario
corresponds to the area of the tube, comprises three principal contributions:
the area of the original triangle, the areas of the three rectangles at the
sides of the triangle, and the contributions of the disc sections at the three
corners of the triangle. With a little geometric intuition, it can also be
inferred that these three disc sections in fact together form a full circle of
radius $r$. Steiner’s formula of expression (32) in the current context then
states

$\displaystyle\begin{split}\lambda(T_{S,r})&=\pi
r^{2}\mu_{0}(S)+2r\mu_{1}(S)+\mu_{2}(S).\\\ \end{split}$ (33)

The first term in eq. (33) corresponds to the area of a circle with radius
$r$, if $\mu_{0}(S)=1$. Thus, the $0$th-order intrinsic volume of a triangle
appears to be 1. More generally, the $0$th-order intrinsic volume corresponds
to the Euler characteristic of a subset $S\subset\mathbb{R}^{D}$, which will
be discussed in Section 3.2: Probabilistic properties of excursion set
features. The second term in eq. (33) corresponds to the area of the three
rectangles of the tube, if $\mu_{1}(S)$ measures the perimeter of the triangle
(i.e., the sum of its side lengths) and is multiplied by $1/2$. Finally, it
follows that $\mu_{2}(S)$ must measure the area of the triangle. Note that it
is also possible to form a three- or higher dimensional tube around a
triangle. More generally, it can be shown that the intrinsic volumes of two-
dimensional subsets $S\subset\mathbb{R}^{2}$ are given by

$\displaystyle\begin{split}\mu_{0}(S)&:\mbox{Euler characteristic of }S\\\
\mu_{1}(S)&:0.5\cdot\mbox{Circumference of }S\\\ \mu_{2}(S)&:\mbox{Area of
}S,\end{split}$ (34)

and for three-dimensional subsets $S\subset\mathbb{R}^{3}$ by

$\displaystyle\begin{split}\mu_{0}(S)&:\mbox{Euler characteristic of }S\\\
\mu_{1}(S)&:2\cdot\mbox{Caliper diameter of }S\\\
\mu_{2}(S)&:0.5\cdot\mbox{Surface area of }S\\\ \mu_{3}(S)&:\mbox{Volume of
}S.\end{split}$ (35)

Here, the caliper diameter of a convex object is evaluated by placing the
solid between two parallel planes (or calipers), measuring the distance
between the planes, and averaging over all rotations of object. This kind of
measurement can be formed using a caliper, a tool that is shown in Figure 3D.
As mentioned above, the concept of the Euler characteristic will be elaborated
on below. As for Lebesgue measure, it should be evident by now that evaluating
the intrinsic volumes for a given geometric object that is modelled as a
subset of $\mathbb{R}^{D}$ is not trivial. In Figure 3E, we collect the
intrinsic volumes of a number of elementary geometrical objects.

Finally, we note that in the context of RFP, the intrinsic volumes of the
search space of interest are evaluated numerically based on an algorithm
proposed by (Worsley et al.,, 1996). This algorithm implies that the first and
second order intrinsic volumes in the case $S\subset\mathbb{R}^{3}$ can be
decomposed into contributions from one- and two-dimensional subspaces along
the cardinal axes, i.e., that

$\displaystyle\mu_{1}(S)$ $\displaystyle=\mu_{1}^{x_{1}}(S)$
$\displaystyle\,+\mu_{1}^{x_{2}}(S)$ $\displaystyle\,+\mu_{1}^{x_{3}}(S)$ (36)
$\displaystyle\mu_{2}(S)$ $\displaystyle=\mu_{2}^{x_{1}x_{2}}(S)$
$\displaystyle\,+\mu_{2}^{x_{1}x_{3}}(S)$
$\displaystyle\,+\mu_{2}^{x_{2}x_{3}}(S),$ (37)

where $\mu_{1}^{x_{d}}(S)\in\mathbb{R},\,d=1,2,3$ denotes a contribution from
the respective one-dimensional subspace of $\mathbb{R}^{3}$, and
$\mu_{2}^{x_{i}x_{j}}(S),1\leq i,j\leq 3,i\neq j$ denotes a contribution from
the respective two-dimensional subspace of $\mathbb{R}^{3}$. Inspection of the
analytical intrinsic volumes in Figure 3E indicates that this decomposition
clearly holds for box-shaped geometric objects.

Figure 3: Lebesgue measure and intrinsic volumes. (A) A rectangle $R$ with
sides oriented parallel to the coordinate axes in $\mathbb{R}^{2}$. (B) A
visualization of Lebesgue outer measure. Consider the central object with the
strong outline as a subset $S\subset\mathbb{R}^{2}$. Then $S$ can be covered
by a (possibly infinite) union of rectangles, such that
$S\subset\cup_{i=1}^{\infty}R_{i}$. The panel shows a coverage of $S$ afforded
by a finite number of rectangles. Clearly, in general the volume of the
coverage overestimates the volume of the object of interest. (C) A two-
dimensional tube around a triangle. The central grey area shows the triangle.
The union of the triangle, the white rectangles, and the dark grey disc
sections forms the two-dimensional tube around the triangle. (D) A caliper.
(E) Intrinsic volumes for a selected set of basic geometric shapes as listed
in Worsley et al., (1996, Table I)). Note that the intrinsic volumes of orders
higher than the dimensionality of the geometric object are zero.

##  3 Theory

We develop the theory of RFP against the background of the continuous space,
discrete data point model

$Y_{i}(x)=m_{i}\beta(x)+\sigma Z_{i}(x)\quad\mbox{ for }i=1,...,n\mbox{ and
}x\in S\subseteq\mathbb{R}^{3}.$ (38)

Here, $Y_{i}(x)$ denotes a random variable that models the $i$th of $n$ data
observations at spatial location $x\in\mathbb{R}^{3}$,
$m_{i}\in\mathbb{R}^{p}$ is a known space-independent vector of $p$ model
constraints (commonly the $i$th row of a design matrix), $\beta(x)$ is an
unknown data point-independent value of a space-dependent effect size
parameter function $\beta:\mathbb{R}^{3}\to\mathbb{R}^{p}$, $\sigma>0$ is an
unknown standard deviation parameter, and $Z_{i}(x)$ is a $Z$-field modelling
observation error. In line with Definition 5 we assume that the
$Z_{i}(x),i=1,...,n$ are independent and of identical gradient covariances.
The model equation (38) should thus be read as the generalization of the
familiar mass-univariate GLM-based fMRI data analysis equation (e.g., Kiebel
and Holmes,, 2007, eq. 8.6) to continuous space, which entails that the data
observations $Y_{i}(x)$ are (usually non-stationary) random fields. In the
context of GLM-based fMRI data analysis, eq. (38) may represent a first-level
model, in which case the $m_{i}$ are typically derived from the convolution of
condition-specific trial onset functions with a haemodynamic response function
and it is assumed that temporal error correlations have been accounted for by
pre-whitening (e.g., Glaser and Friston,, 2007). Equivalently, eq. (38) may
represent a second-level model in the summary statistics approach (e.g.
Mumford and Nichols,, 2006), in which case the $m_{i}$ typically represent
categorical statistical designs and comprise primarily ones and zeros, and the
assumption of independent error contributions is natural.

The fundamental aim of RFP is to evaluate the probabilities of topological
features of statistical parametric maps under the null hypothesis of no
activation. With regards to the model of eq. (38), such a null hypothesis
corresponds to $\beta(x)=0$ for all $x\in\mathbb{R}^{3}$. The null hypothesis
thus entails that the standardized data observations
$\sigma^{-1}Y_{i}(x),i=1,...,n$ are $Z$-fields. Consequently, the evaluation
of location-specific $T\text{-}$ or $F\text{-}$ratios (cf. the right-hand
sides of eqs. (19) and (23), respectively) results in $T$\- or $F$-fields. In
the current section, we are concerned with the probabilistic theory of
topological features of these statistical random fields, which we will denote
generically by $X(x),x\in S$. As outlined in the Introduction, we first
consider the fundamental RFP concepts of excursion sets, smoothness, and resel
volumes of statistical random fields. Building on these concepts, we then
consider the theory of the probabilistic properties of excursion set features.
It should be noted that throughout this section we assume a continuous model
of space and the validity of the null hypothesis.

###  3.1 Excursion sets, smoothness, and resel volumes

#### Excursion sets

Excursion sets of statistical random fields above a level $u\in\mathbb{R}$ are
an elementary building block in the theory of RFP, because the topological
features of statistical random fields that RFP is concerned with are
topological features of its excursion sets. We use the following definition
(Worsley,, 1994, Section 2, p. 14):

######  Definition 8 (Excursion set).

Let $X(x),x\in\mathbb{R}^{D}$ denote a statistical random field, let
$S\subset\mathbb{R}^{D}$ denote a subset of $\mathbb{R}^{D}$ referred to as
the search space, and let $u\in\mathbb{R}$ denote a level. Then the excursion
set of $X(x)$ inside the search space $S$ above the level $u$ is defined as

$E_{u}:=\\{x\in S|X(x)\geq u\\}.$ (39)

$\hfill\bullet$

In words, the excursion set $E_{u}$ comprises all points $x$ in the search
space $S$ for which the statistical random field $X$ takes on values larger or
equal to $u$. In the neuroimaging literature, the level $u$ is referred to as
the cluster-forming (or cluster-defining) threshold (Nichols et al.,, 2017;
Flandin and Friston,, 2017; Eklund et al.,, 2016). Clearly, because $X$ is a
random entity, the excursion set is also a random entity. From a probabilistic
perspective, this entails that the characterization of excursion sets is
concerned with expectations and distributions of excursion set features,
because the precise properties of an excursion set vary from one realization
of a statistical random field to another. Figure 4A visualizes the excursion
set above a level of $u:=1$ for four realizations of a $Z\text{-}$field with
Gaussian covariance function.

The probabilistic properties of topological features of excursion sets depend
on many things. One observation, which is exploited repeatedly in the theory
of RFP, is the following insight: as the level $u$ increases, the constitution
of an excursion set follows a predictable path (Figure 4B): at low values of
$u$, excursion sets typically display complex topological shapes, comprising
multiple clusters, which in turn may have parts that are not part of the
excursion set and appear as cluster holes. At higher values of $u$, holes tend
to disappear and the number of clusters decreases. Finally, at very high
values of $u$ around the level of the statistical random field realization’s
global maximum, the excursion set comprises only a single cluster (for $u$
slightly smaller than the realization’s global maximum), a single point (for
$u$ being equal to the global maximum), or no points at all (for all values of
$u$ that are larger than the global maximum). As will be discussed below, this
observation is used frequently in the approximation of probabilistic
properties of excursion sets. A second observation is that, naturally, the
topological properties of the excursion set at a fixed level $u$ depend on the
characteristics of the underlying statistical random field. Clearly, for
general nonstationary random fields, regions characterized by higher values of
the expectation function than others have a higher probability of being
contained in the excursion set. For the stationary statistical random fields
that are of interest in RFP, the topological properties of excursion sets
depend strongly on the field’s smoothness: smoother random fields tend to have
a more slowly and shallowly varying profile, while less smooth (rougher)
random fields typically display a more “peaky” profile. This results in the
observation that excursion sets of smooth statistical random fields either
contain few, but larger clusters, while excursion sets of rough statistical
random fields contain more, but smaller clusters (Figure 4C).

Figure 4: Excursion sets. (A) This panel visualizes the random nature of
excursion sets. The upper row of four panels depicts four realizations of a
$Z\text{-}$field with a Gaussian covariance function with parameters
$v=1,\ell=0.15$. The lower four panels depict the realization-specific
excursion set for a fixed level of $u=1$. Points included in the excursion
sets are marked white, points not included in the excursion sets are black.
(B) This panel visualizes the topological properties of the excursion set of a
single $Z\text{-}$field realization (leftmost subpanel) as a function of the
level value $u$. As $u$ increases from $-0.5$ to $2$, the topological
constitution of the excursion set becomes less complex, and, in the limit of
high levels $u$ comprises a single cluster comprising the global maximum of
the statistical random field’s realization or no points at all. (C) This panel
visualizes the dependency of the topological properties of the excursion set
on the smoothness of the statistical random field. The left two subpanels show
a $Z\text{-}$field realization and its corresponding excursion set for $u=1$
for a rough (non-smooth) field. The excursion set comprises many small
clusters. The right two subpanels show the same entities for a smoother field.
Here, the excursion set comprises a smaller number of clusters, which are
larger than in the case of the rough $Z\text{-}$ field. For implementational
details of these simulations, please see rfp_4.m.

#### Smoothness

As visualized in Figure 2 and Figure 4, realizations of random fields can be
smooth, i.e., varying little over a given spatial distance, or rough, i.e.,
varying strongly over a given spatial distance. This intuitive notion of
smoothness is an important determinant of the probabilistic behaviour of
topological features of random fields and of central importance in RFP
(Friston et al.,, 1991). Thus, in order to establish quantitative
relationships between a statistical random field’s smoothness and the
probabilistic behaviour of its topological features, a quantitative measure of
smoothness is required. The principle measure used to quantify smoothness in
RFP is the reciprocal unit-square Lipschitz-Killing curvature (Taylor and
Worsley,, 2007, eq. 3), defined as

$\varsigma:=\frac{1}{\int_{[0,1]^{D}}|\mathbb{V}\left(\nabla
X(x)\right)|^{\frac{1}{2}}\,dx},$ (40)

where

$\mathbb{V}(\nabla X(x)):=\left(\mathbb{C}\left(\frac{\partial}{\partial
x_{i}}X(x),\frac{\partial}{\partial x_{j}}X(x)\right)\right)_{1\leq i,j\leq
D}$ (41)

denotes the $D\times D$ variance matrix of the gradient components (i.e.,
partial derivatives) of the random field at $x\in[0,1]^{D}$. Note that we
refer to $\mathbb{V}(\nabla X(x))$ as a variance matrix (as opposed to a
covariance matrix) despite the fact that entries of $\mathbb{V}(\nabla X(x))$
are given by covariances. This is motivated by the fact that these covariances
refer to the partial derivatives of the identical random variable $X(x)$ and
not to partial derivatives of covariances of two different random variables
(Friston et al.,, 1991; Worsley et al.,, 1992; Taylor and Worsley,, 2007). In
the following, we first attempt to elucidate in which sense eqs. (40) and (41)
capture the intuitive notion of a random field’s smoothness. We then discuss,
how eq. (40) is reformulated in RFP to give rise to the “full width at half
maximum” parameterization of smoothness. Throughout, we make the assumption
that the random field’s gradient exists at all locations of the random field’s
domain, i.e., that the random field is differentiable with respect to space.

To obtain a first intuition in which sense eqs. (40) and (41) provide a
measure of smoothness, we consider the case of a stationary GRF on a one-
dimensional domain ($D=1$) with a GCF $\gamma$ and parameters $v:=1$ and
$\ell>0$. In this case, the gradient of $X$ at the location $x\in[0,1]$
simplifies to the derivative of $X$ with respect to $x$, which, in line with
the conventions of spatial statistics, we denote by $\dot{X}(x)$. Further, the
variance matrix of the gradient components simplifies to the variance of this
spatial derivative, and the determinant to the absolute value, which in turn
is redundant for a non-negative variance. We are thus led to consider

$\varsigma=\frac{1}{\int_{[0,1]}\mathbb{V}(\dot{X}(x))^{\frac{1}{2}}\,dx}.$
(42)

We may first note that, regardless of the type of random field, the reciprocal
of the square of the variance of the spatial derivative of the random field at
a location $x$ is a sensible measure of the field’s smoothness, if we assume
that this value is constant over space: naturally, the spatial derivative
quantifies the rate of change of the values of the random field, and if this
is high, the random field changes a lot, if it is low, the random field does
not change much. In addition, the variance describes how variable this spatial
variability is over realizations of the random field, and, if it is small,
most realizations of the random field with small spatial derivatives will
appear smooth. Furthermore, for the current example of a one-dimensional GRF
with a GCF, one may analytically evaluate eq. (42) and, as shown in Supplement
S3.1, obtains

$\varsigma=\frac{\ell}{\sqrt{2}}.$ (43)

Thus, for the current scenario, the smoothness measure defined by eqs. (40)
and (41) is directly related to the spatial covariation parameter of the
covariance function of the GRF. This is, of course, a special case. An
important aspect of eq. (40) is that it is defined as long as the partial
derivatives of the random field of interest are defined, but for many random
fields it may not directly map onto a single parameter of the field’s
covariance function. This is analogous to the difference between a variance of
a random variable and the variance parameter of a univariate Gaussian
distribution: while the former concept applies to any random variable, the
case that it directly maps onto a parameter of the probability density
function of a random variable, like the variance parameter of the univariate
Gaussian, is rather an exception.

More generally, the smoothness measure defined by eqs. (40) and (41) involves
the determinants of the spatial gradients at locations $x\in[0,1]^{D}$.
Assuming again that the gradients are constant over space, we next consider
the two-dimensional scenario ($D=2$). In this case, the determinant is given
by

$|\mathbb{V}(\nabla X(x))|=\mathbb{V}\left(\frac{\partial}{\partial
x_{1}}X(x)\right)\mathbb{V}\left(\frac{\partial}{\partial
x_{2}}X(x)\right)-\mathbb{C}\left(\frac{\partial}{\partial
x_{1}}X(x),\frac{\partial}{\partial x_{2}}X(x)\right)^{2}.$ (44)

As is evident from eq. (44), the variances and covariances of the field’s
partial derivatives make opposing contributions to the field’s smoothness: the
variances of the partial derivatives contribute to roughness, while the
covariances of the partial derivatives contribute to the field’s smoothness.
The former was already observed for the one-dimensional scenario. The latter
can be interpreted as follows: a systematic simultaneous change of values in
the direction of both the $x_{1}$ and $x_{2}$ ordinates indicates smoothness.
Finally, an intuition about the three-dimensional case is afforded by the
intuitive view of the determinant of a $3\times 3$ matrix as the magnification
factor of the volume of a unit cuboid under the linear transform furnished by
the matrix (Hannah,, 1996). In this case, large values for the off-diagonal
elements with respect to the diagonal elements imply a smaller magnification
and vice versa, thus implicating large variances of the partial derivatives in
low smoothness and large correlations of the partial derivatives in high
smoothness. In summary, if a random field is differentiable with respect to
space, then eqs. (40) and (41) provide a sensible scalar measure of the
intuitive notion of a random field’s smoothness.

#### Smoothness reparameterization

While $\varsigma$ thus constitutes a perfectly fine scalar measure of a
statistical random field’s smoothness, Friston et al., (1991) and Worsley et
al., (1992) proposed to reparameterize $\varsigma$ in the three-dimensional
case ($D=3$) as

$\varsigma=(4\ln 2)^{-\frac{3}{2}}f_{x_{1}}f_{x_{2}}f_{x_{3}}\mbox{ with
}f_{x_{1}},f_{x_{2}},f_{x_{3}}>0.$ (45)

In the context of RFP, the values $f_{x_{1}},f_{x_{2}},f_{x_{3}}$ are referred
to as the full widths at half maximum (FWHMs) in direction $x_{1},x_{2}$ and
$x_{3}$, respectively. The reparameterization of smoothness in terms of FWHMs
oriented along the cardinal axes of three-dimensional space was motivated by
assuming that the realization of a statistical random field of interest had
been created by convolving a GRF with a white-noise covariance function with a
Gaussian convolution kernel (see Supplement S3.2 for formal definitions of
these entities). In this case, the smoothness of the resulting statistical
random field depends on the spatial width of the Gaussian convolution kernel.
The spatial width of a Gaussian convolution kernel, in turn, can be described
in terms of the covariance matrix of the kernel. If it is assumed that the
covariance matrix of the Gaussian convolution kernel is diagonal, then the
width of the convolution kernel is determined by the three diagonal parameters
$\sigma_{x_{1}}^{2},\sigma_{x_{2}}^{2}$ and $\sigma_{x_{3}}^{2}$, which govern
the widths of the respective one-dimensional marginal Gaussian functions of
the convolution kernel. An alternative representation of the width of a
Gaussian function in turn is afforded by its FWHM: in general, the FWHM of a
function $f$ is the value $x$ for which $f(-x/2)=f(x/2)=f(0)/2$. As shown in
Supplement S3.2, for a univariate Gaussian function extending over the $x_{1}$
domain with variance parameter $\sigma^{2}_{x}$, this value is given by

$f_{x}=\sqrt{8\ln 2}\sigma_{x}.$ (46)

In other words, the FWHM of a univariate Gaussian function is a scalar
multiple of the square root of its variance parameter. Eq. (45) was then
motivated by approximating the variance matrix of partial derivatives
$\mathbb{V}(\nabla X(x))$ of a given statistical random field $X(x),x\in S$ by
the variance matrix of partial derivatives of a hypothetical random field
$Y(x),x\in S$ that was imagined to have been created by the convolution of a
white-noise GRF with an isotropic Gaussian convolution kernel parameterized in
terms of its FWHMs $f_{x_{1}},f_{x_{2}}$, and $f_{x_{3}}$. For such a field,
the variance matrix of partial derivatives was proposed to be independent of
space and of the form (Worsley et al.,, 1992, eq. 2)

$\mathbb{V}(\nabla Y(x))=4\ln 2\left(\begin{matrix}f_{x_{1}}^{-2}&0&0\\\
0&f_{x_{2}}^{-2}&0\\\ 0&0&f_{x_{3}}^{-2}\end{matrix}\right)=:\Lambda.$ (47)

For this hypothetical field, it then follows directly that

$\varsigma=\left(\int_{[0,1]^{D}}|\Lambda|^{\frac{1}{2}}\,dx\right)^{-1}=\left((4\ln
2)^{3}(f_{x_{1}}f_{x_{2}}f_{x_{3}})^{-2}\right)^{-\frac{1}{2}}=(4\ln
2)^{-\frac{3}{2}}f_{x_{1}}f_{x_{2}}f_{x_{3}}.$ (48)

In Supplement S3.2 we show that the diagonal elements of the variance matrix
of gradient components of a three-dimensional white-noise GRF convolved with
an isotropic Gaussian convolution kernel with FWHMs $f_{x_{1}},f_{x_{2}}$, and
$f_{x_{3}}$ are indeed of the form implied by the right-hand side of eq. (47)
(see Holmes, (1994) and Jenkinson, (2000) for similar work). Note, however,
that this only (at least partially) validates the construction of the
approximation, but not the approximation

$\mathbb{V}(\nabla X(x))\approx\mathbb{V}(\nabla Y(x)),$ (49)

which is implicit in the reparameterization of the smoothness parameter of a
given statistical random field $X(x),x\in S$ in terms of FWHMs, perse.

#### Resel volumes

Building on the parameterization of smoothness in terms of FWHMs, the concept
of resel volumes was introduced by Worsley et al., (1992). Intuitively, resel
volumes are the smoothness-adjusted intrinsic volumes of subsets of domains of
random fields. Stated differently, the smoother a random field, the larger the
resel volumes of subsets of its domain at constant intrinsic volumes. Based on
the general definition of intrinsic volumes for the three-dimensional case
(35), the decomposition of the first- and second-order intrinsic volumes (36),
and the FWHM parameterization of a random field’s smoothness (48), the resel
volumes of a set $S\subset\mathbb{R}^{D}$ in the domain of a random field are
then defined as follows (Worsley et al.,, 1996, pp. 63, eq. (3.2)):

######  Definition 9 (Resel volumes).

Let $X(x),x\in\mathbb{R}^{3}$ denote a statistical random field with
smoothness

$\varsigma=(4\ln 2)^{-\frac{3}{2}}f_{x_{1}}f_{x_{2}}f_{x_{3}}$ (50)

and let $S\subset\mathbb{R}^{3}$ denote a subset of the domain of the
statistical random field. Then the $0$th- to $3$rd-order resel volumes of $S$
are defined in terms of the intrinsic volumes of $S$ and the smoothness
parameters $f_{x_{1}},f_{x_{2}}$ and $f_{x_{3}}$ by

$\displaystyle R_{0}(S)$ $\displaystyle:=\mu_{0}(S)$ (51) $\displaystyle
R_{1}(S)$ $\displaystyle:=\frac{1}{f_{x_{1}}}\mu_{1}^{x_{1}}(S)$
$\displaystyle\,+\frac{1}{f_{x_{2}}}\mu_{1}^{x_{2}}(S)$
$\displaystyle\,+\frac{1}{f_{x_{3}}}\mu_{1}^{x_{3}}(S)$ (52) $\displaystyle
R_{2}(S)$ $\displaystyle:=\frac{1}{f_{x_{1}}f_{x_{2}}}\mu_{2}^{x_{1}x_{2}}(S)$
$\displaystyle\,+\frac{1}{f_{x_{1}}f_{x_{3}}}\mu_{2}^{x_{1}x_{3}}(S)$
$\displaystyle\,+\frac{1}{f_{x_{2}}f_{x_{3}}}\mu_{2}^{x_{2}x_{3}}(S)$ (53)
$\displaystyle R_{3}(S)$
$\displaystyle:=\frac{1}{f_{x_{1}}f_{x_{2}}f_{x_{3}}}\mu_{3}(S)$ (54)

$\hfill\bullet$

The resel volumes of a given set $S\subset\mathbb{R}^{3}$ can thus be readily
evaluated if the intrinsic volumes of $S$ and their subspace contributions, as
well as the FWHM smoothness parameters $f_{x_{1}}$,$f_{x_{2}}$, $f_{x_{3}}$ of
the statistical random field extending over $S$, are known.

###  3.2 Probabilistic properties of excursion set features

The probabilistic properties of the following five topological features of
excursion sets of statistical random fields form the core of RFP theory: (P1)
the global maximum of the statistical random field, (P2) the number of
clusters within an excursion set, (P3) the volume of clusters within an
excursion set, (P4) the number of clusters of a given volume within an
excursion set, and (P5) the maximal volume of a cluster within an excursion
set. In the theory of RFP, these topological features are described by random
variables, the probability distributions of which are governed by the
stochastic and geometric characteristics of the underlying statistical random
field and the level $u$ of the respective excursion sets. The dependencies of
the respective distributions on the properties of the underlying statistical
random field and the excursion set level are mitigated by four expected
values: (E1) the expected volume of an excursion set, (E2) the expected number
of local maxima within an excursion set, (E3) the expected number of clusters
within an excursion set, and (E4) the expected volume of a cluster within an
excursion set. In the following, we first discuss the parametric dependencies
of these expectations on the type of statistical random field and its resel
volumes, i.e., its smoothness and intrinsic volumes. We then discuss the
probability distributions of the five topological features that form the core
of RFP theory. As will become evident throughout this section, a defining
characteristic of RFP theory is that the parametric distributional forms of
virtually all topological features of interest derive from approximations,
rather than analytical transforms of the probability density functions
describing the statistical random fields under functions that instantiate the
topological features. For the presentation of RFP theory herein, this entails
that we will most often only report the parametric distributional form of a
given topological feature under RFP and discuss some background about its
motivation, rather than analytically deriving the respective parametric
distributional form from first principles. The approximative-parametric nature
is inherent to RFP and is rooted in the seminal works by Adler and Hasofer,
(1976) and Adler, (1976) and has more recently seen major advances by Taylor
et al., (2005, 2006) and Adler and Taylor, (2007) (see also Adler, (2000) for
a historical perspective).

#### Expected values

(E1) The expected volume of an excursion set

Recall that an excursion set $E_{u}$ is a subset of the search space $S$ for
which a statistical random field takes on values larger or equal to
$u\in\mathbb{R}$. As discussed above, the standard analytical measure that
assigns a notion of volume to subsets of $\mathbb{R}^{D}$ is the Lebesgue
measure $\lambda$, which returns the Lebesgue volume $\lambda(E_{u})$ of the
excursion set. As discussed by Adler, (1981, Section 1.7, pp. 18 - 19), the
expected value of the Lebesgue volume of an excursion set is given by

$\mathbb{E}\left(\lambda(E_{u})\right)=\lambda(S)(1-F_{X}(u)),$ (55)

where $\lambda(S)$ denotes the Lebesgue volume of the search space and $F_{X}$
denotes the distribution function of the statistical random field at the
origin (see also Hayasaka and Nichols, (2003, eq. 4) and, for
$Z\text{-}$fields, Friston et al., 1994b (, eq. 3)). Because for stationary
random fields, this distribution function equals the distribution functions at
all other locations of the random field, we denote it by $F_{X}$ without
reference to its spatial position. A brief derivation of eq. (55) is provided
in Supplement S3.3. As it stands, eq. (55) has strong intuitive appeal: the
expected value of the excursion set $E_{u}$ is a proportion of the search
space $S$, and the proportion factor is given by the probability of the
statistical random field to take on values larger than $u$ at each point in
space. For example, for small values of $u$, the probability that the
statistical random field takes on values larger than $u$ is always higher than
for larger values of $u$. Accordingly, the expected value of the volume of the
excursion set is larger in the former than in the latter case. From eq. (55),
it follows directly that if the three-dimensional volume of the search space
is expressed in terms of the third-order resel volume $R_{3}(S)$, then the
expected third-order resel volume of the excursion set is given by

$\mathbb{E}(R_{3}(E_{u}))=R_{3}(S)(1-F_{X}(u)).$ (56)

In the following, we will denote the random variable modelling the expected
volume of an excursion set independent of its unit of measurement by $V_{u}$,
where $V_{u}:=\lambda(E_{u})$ or $V_{u}:=R_{3}(E_{u})$ depending on the
specified context. We visualize the expected volume of the excursion set as a
function of $u$ and $R_{3}(S)$ in Figure 5B. As is evident from this
visualization, the expected volume of an excursion set $\mathbb{E}(V_{u})$
decreases with higher cluster-forming thresholds $u$ and higher statistical
random field smoothness, i.e., decreasing third-order resel volume $R_{3}(S)$.

(E2) The expected number of local maxima within an excursion set

Let $M_{u}$ denote the number of local maxima of a statistical random field
$X(x),x\in\mathbb{R}^{D}$ above $u$ inside a set $S\subset\mathbb{R}^{D}$.
Then $M_{u}$ is also the number of local maxima of the excursion set $E_{u}$.
In RFP, the expected number of local maxima of a statistical random field is
approximated by the expected Euler characteristic $\chi(\cdot)$ of the
excursion set (Worsley,, 1994, Section 2, p. 15), i.e.

$\mathbb{E}(M_{u})\approx\mathbb{E}(\chi(E_{u})).$ (57)

To unpack expression (57), we first discuss the notion of the Euler
characteristic as a topological measure. We then consider the motivation for
the approximation of expression (57), and finally discuss the parametric
closed form for the expected Euler characteristic used in RFP.

The Euler characteristic is a function and is referred to as a topological
invariant. Topological invariants measure the geometrical characteristics of
geometrical objects irrespective of the way they are bent. A formal
description of the Euler characteristic is beyond our current scope, because,
as Adler, (2000, p.16) puts it, the Euler characteristic “(…) is one of those
unfortunate cases in which what is easy for the human visual system to do
quickly and effectively requires a lot more care when mathematicised”. This
care is not warranted here, because, as will be seen below, the expected Euler
characteristic of the excursion set $E_{u}$ will itself be approximated by a
function of $u$ without recourse to an actual evaluation of $\chi(E_{u})$. We
thus only provide a brief intuition for the Euler characteristic itself: for a
three-dimensional volume of interest, the Euler characteristic counts, in an
alternating sum, the number of three types of topological features: (1)
connected components, (2) visible open holes, often referred to as handles,
and (3) invisible voids. Semi-formally, the Euler characteristic of an
excursion set $E_{u}$ is hence given by

$\displaystyle\begin{split}\chi(E_{u})=\quad&\,\mbox{Number of connected
components of }E_{u}\\\ -&\,\mbox{Number of handles of }E_{u}\\\
+&\,\mbox{Number of voids in }E_{u}.\end{split}$ (58)

The following examples for familiar geometric objects provided by Worsley,
1996a  are instructive:

  * $\chi(E_{u})=1$ for a single solid excursion set (1 connected component - 0 handles + 0 voids),

  * $\chi(E_{u})=0$ for a doughnut-shaped excursion set (1 connected component - 1 handles + 0 voids),

  * $\chi(E_{u})=-2$ for a pretzel-shaped excursion set (1 connected component - 3 handles + 0 voids),

  * $\chi(E_{u})=2$ for a tennisball-shaped excursion set (1 connected component - 0 handles + 1 voids).

More qualitatively, if an excursion set comprises many disconnected
components, each with very few holes (in astrophysics referred to as a
“meatball” topology), the Euler characteristic is positive; if the clusters of
an excursion set are connected by many bridges, thus creating many holes (in
astrophysics referred to as a “sponge” topology), the Euler characteristic is
negative, and if the excursion set comprises many surfaces that enclose
hollows (in astrophysics referred to as a “bubble” topology), the Euler
characteristic is again positive.

In the theory of RFP, it is, however, not the value of the Euler
characteristic for a specific realization of an excursion set that is of
primary interest, but rather its expected value - as evident from the right-
hand side of expression (57). This approximation (or “heuristic”) is based on
the following intuition: as the threshold value $u$ that defines the excursion
set $E_{u}$ increases, voids and handles of the excursion set tend to
disappear, and what remains are isolated connected components of the excursion
set (in the neuroimaging literature referred to as clusters). In this
scenario, $\chi(E_{u})$ thus counts the number of connected components (cf.
eq. (58)). Each of the clusters in the excursion set naturally contains a
local maximum of the excursion set, which thus motivates to approximate the
expected number of local maxima by the expected Euler characteristic over
realizations of the statistical random field $X(x)$. As for the accuracy of
this approximation, it has been demonstrated that it becomes exact in the
limit of very high values of $u$, in which $\chi(E_{u})$ evaluates to either 1
(one remaining cluster) or 0 (no remaining clusters) (Taylor et al.,, 2005).
For lower values of $u$, an analytical validation does not seem to exist so
far. Crucially though, the approximation in expression (57) is helpful from a
mathematical perspective, because there exists a parametric closed form for
the expected Euler characteristic. This parametric closed form was introduced
by (Worsley et al.,, 1996, eq. 3.1) and is given by

$\mathbb{E}(\chi(E_{u}))=\sum_{d=0}^{D}R_{d}(S)\rho_{d}(u),$ (59)

where the $R_{d}(S),d=0,...,D$ are the $d$-dimensional resel volumes of the
search space and the $\rho_{d}(u)$ are referred to as Euler characteristic
(EC) densities. The formula on the right-hand side of eq. (59) is powerful
because it decomposes the expected value of the Euler characteristic into (1)
smoothness-adjusted volumetric contributions from each of the $d=0$ to $d=D$
dimensions of the search space, and (2) unit volume-smoothness densities
$\rho_{d}$ that only depend on the value of $u$ and the type of statistical
random field on the other hand. We have already discussed the meaning of the
$R_{d}$ in Section 3.1 and hence focus on the EC densities in the following.
Intuitively, the functions

$\rho_{d}:\mathbb{R}\to\mathbb{R},u\mapsto\rho_{d}(u)\mbox{ for }d=0,...,D$
(60)

describe the contribution of unit smoothness resel volumes to the expected
value of the Euler characteristic. The values taken by these functions depend
on $u$, and in general, as the value of $u$ increases (i.e., the excursion set
is formed based on a higher threshold), these values decrease. For constant
resel volumes, this implies that the value of the expected Euler
characteristic, and thus the approximated expected number of local maxima,
decreases. The exact parametric forms of the function $\rho_{d}$ depend on the
type of statistical random field and have been analytically evaluated by
Worsley, (1994). We list their functional forms in Table 2 and visualize their
form based on their SPM implementation in Figure 5A. Note that the $0$th order
EC density of each field corresponds to $1-F_{X}(u)$, where $F_{X}$ denotes
the distribution function of the respective field. This allows for evaluating
the expected volume of the excursion set as expressed in eq. (56) in terms of
the third order resel volume and the $0$th order EC density. In summary, the
expected number of local maxima of an excursion set can thus be approximated
for a given statistical random field in a straight-forward manner, if the
field’s resel volumes are known (or have been estimated).

$Z$-field Euler characteristic densities
---
$\rho_{0}(u)=\int_{u}^{\infty}\frac{1}{(2\pi)^{1/2}}\exp\left(-t^{2}/2\right)dt$
$\rho_{1}(u)=\frac{(4\ln 2)^{1/2}}{2\pi}\exp\left(-u^{2}/2\right)$
$\rho_{2}(u)=\frac{4\ln 2}{(2\pi)^{3/2}}\exp\left(-u^{2}/2\right)u$
$\rho_{3}(u)=\frac{(4\ln
2)^{3/2}}{(2\pi)^{2}}\exp\left(-u^{2}/2\right)(u^{2}-1)$
$T$-field Euler characteristic densities
$\rho_{0}(u;\nu)=\int_{u}^{\infty}\frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\sqrt{\nu\pi}\Gamma\left(\frac{\nu}{2}\right)\left(1+\frac{t^{2}}{\nu}\right)^{-1/2(\nu+1)}}dt$
$\rho_{1}(u;\nu)=\frac{(4\ln
2)^{1/2}}{2\pi}\left(1+\frac{u^{2}}{\nu}\right)^{-1/2(\nu-1)}$
$\rho_{2}(u;\nu)=\frac{(4\ln
2)\Gamma\left(\frac{\nu+1}{2}\right)}{(2\pi)^{3/2}\left(\frac{\nu}{2}\right)^{1/2}\Gamma\left(\frac{\nu}{2}\right)}\left(1+\frac{u^{2}}{\nu}\right)^{-1/2(\nu-1)}u$
$\rho_{3}(u;\nu)=\frac{(4\ln
2)^{3/2}}{(2\pi)^{2}}\left(1+\frac{u^{2}}{\nu}\right)^{-1/2(\nu-1)}\left(\frac{\nu-1}{\nu}u^{2}-1\right)$
$F$-field Euler characteristic densities
$\rho_{0}(u;\nu_{1},\nu_{2})=\int_{u}^{\infty}\frac{\Gamma((\nu_{1}+\nu_{2})/2)}{\Gamma(\nu_{1}/2)\Gamma(\nu_{2}/2)}\frac{\nu_{1}}{\nu_{2}}\left(\frac{\nu_{1}}{\nu_{2}}t\right)^{1/2(\nu_{1}-2)}\left(1+\frac{\nu_{1}t}{\nu_{2}}\right)^{-1/2(\nu_{1}+\nu_{2})}\,dt$
$\rho_{1}(u;\nu_{1},\nu_{2})=\frac{(4\ln
2)^{1/2}}{(2\pi)^{1/2}}\frac{\Gamma((\nu_{1}+\nu_{2}-1)/2)2^{1/2}}{\Gamma(\nu_{1}/2)\Gamma(\nu_{2}/2)}\left(\frac{\nu_{1}u}{\nu_{2}}\right)^{1/2(\nu_{1}-1)}\left(1+\frac{\nu_{1}u}{\nu_{2}}\right)^{-1/2(\nu_{1}+\nu_{2}-2)}$
$\rho_{2}(u;\nu_{1},\nu_{2})=\frac{4\ln
2}{2\pi}\frac{\Gamma((\nu_{1}+\nu_{2}-2)/2)}{\Gamma(\nu_{1}/2)\Gamma(\nu_{2}/2)}\left(\frac{\nu_{1}u}{\nu_{2}}\right)^{1/2(\nu_{1}-2)}\left(1+\frac{\nu_{1}u}{\nu_{2}}\right)^{-1/2(\nu_{1}+\nu_{2}-2)}\left((\nu_{2}-1)\frac{\nu_{1}u}{\nu_{2}}-(\nu_{1}-1)\right)$
$\rho_{3}(u;\nu_{1},\nu_{2})=\frac{(4\ln
2)^{3/2}}{(2\pi)^{3/2}}\frac{\Gamma((\nu_{1}+\nu_{2}-3)/2)2^{-1/2}}{\Gamma(\nu_{1}/2)\Gamma(\nu_{2}/2)}\left(\frac{\nu_{1}u}{\nu_{2}}\right)^{1/2(\nu_{1}-3)}\left(1+\frac{\nu_{1}u}{\nu_{2}}\right)^{-1/2(\nu_{1}+\nu_{2}-2)}$
$\quad\quad\quad\quad\quad\quad\quad\quad\quad\times\left((\nu_{2}-1)(\nu_{2}-2)\left(\frac{\nu_{1}u}{\nu_{2}}\right)^{2}-(2\nu_{1}\nu_{2}-\nu_{1}-\nu_{2}-1)\frac{\nu_{1}u}{\nu_{2}}+(\nu_{1}-1)(\nu_{1}-2)\right)$
Table 2: Euler characteristic densities for $Z\text{-}$fields and for
$T\text{-}$ and $F\text{-}$fields with $\nu$ and $\nu_{1},\nu_{2}$ degrees of
freedom, respectively. Note that the $0$th-order Euler characteristic
densities correspond to complementary distribution functions of $Z,T$ and $F$
random variables. Reproduced from Table II in Worsley et al., (1996).

(E3) The expected number of clusters within an excursion set

In the SPM implementation of RFP, the expected number of clusters is
approximated by the expected Euler characteristic (Friston et al.,, 2007,
Chpt. 18, pp. 233-234, Chpt. 19, pp. 240-241). This means that the
approximations of the expected number of local maxima and the approximation of
the expected number of clusters are identical. This identity entails the
assumption that each cluster contains only a single local maximum of the
excursion set. Cast more formally, let $M_{u}$ denote the number of local
maxima of an excursion set $E_{u}$, and let

$\mathcal{M}:=\\{m_{1},...,m_{M_{u}}\\}\subset E_{u}$ (61)

denote the set of local maxima of $E_{u}$. Clusters are unconnected components
of the excursion set and thus form a partition $\mathcal{P}_{u}$ of $E_{u}$,
i.e.,

$\mathcal{P}_{u}=\\{\mathcal{C}_{1},...,\mathcal{C}_{C_{u}}\\}\mbox{ with
}\mathcal{C}_{i}\cup\mathcal{C}_{j}=\emptyset\mbox{ for }i\neq j\mbox{ and
}\cup_{i=1}^{C_{u}}\mathcal{C}_{i}=E_{u},$ (62)

where

$C_{u}:=|\mathcal{P}_{u}|$ (63)

denotes the number of clusters. While each cluster contains at least one local
maximum, in general it is very well possible that a given cluster comprises
multiple local maxima. This is also a common finding in the actual analysis of
fMRI data using SPM, in which case the local maxima for large clusters are
indeed listed in the SPM results table based on a distance selection criterion
of typically larger than 8 mm apart. The validity of the approximation chain

$\mathbb{E}(C_{u})\approx\mathbb{E}(M_{u})\approx\mathbb{E}(\chi(E_{u}))=\sum_{d=0}^{D}R_{d}(S)\rho_{d}(u),$
(64)

thus strongly rests on the assumption of high values for the excursion set-
defining threshold $u$: in the case that only the global maximum of the random
field is contained in the single connected component of the field’s excursion
set, the number of local maxima of the excursion set and the number of
clusters are identical with probability one. In Figure 5B, we visualize the
expected number of clusters $\mathbb{E}(C_{u})$ within an excursion set over a
range of values for $u$ and interpolations of the resel volumes
$R_{0}(S),...,R_{3}(S)$ of the exemplary data set underlying Figure 1. As is
evident from this visualization, the expected number of clusters within an
excursion set decreases with increasing values for the cluster-forming
threshold and increasing smoothness of the statistical random field,
corresponding to decreasing resel volumes.

(E4) The expected volume of clusters within an excursion set

The expected volume of a single cluster is evaluated in RFP based on the
expressions for (1) the expected volume of the excursion set, and (2) the
expected number of clusters in such an excursion set. For a value
$u\in\mathbb{R}$, let the cluster volume be defined by the random variable

$K_{u}:=\frac{V_{u}}{C_{u}}.$ (65)

Then, if $V_{u}$ and $C_{u}$ are independent random variables with finite
expected values, the expected cluster volume is given by

$\mathbb{E}(K_{u})=\frac{\mathbb{E}(V_{u})}{\mathbb{E}(C_{u})}$ (66)

(e.g., Friston et al., (1996, eq. 2), Hayasaka and Nichols, (2003, eq. 3)).
Naturally, if the expected volume of the excursion set is expressed in terms
of Lebesgue measure, i.e., $V_{u}=\lambda(E_{u})$, eq. (66) denotes the
expected cluster Lebesgue volume, and if the expected volume of the excursion
set is expressed in terms of third-order resel volume, i.e.,
$V_{u}=R_{3}(E_{u})$, eq. (66) denotes the expected cluster resel volume. We
visualize the expected cluster resel volume as a function of the cluster-
forming threshold and the smoothness of the statistical random field in Figure
5B. As is evident from this visualization, the expected cluster resel volume
decreases with increasing values of the cluster-forming threshold and
increasing statistical random field smoothness.

Figure 5: Probabilistic properties of excursion set features. (A) Euler
characteristic densities. The panels depict, from left to right, the $Z$-field
Euler characteristic (EC) densities, the $T$-field EC densities for a
$T$-field with 15 degrees of freedom and the $F$-field EC densities for an
$F$-field with (1,20) degrees of freedom. (B) Excursion set feature
expectations for a $T$-field with $15$ degrees of freedom. From left to right,
the panels show the expected excursion set volume $\mathbb{E}(V_{u})$ in units
of resels, the expected number of clusters $\mathbb{E}(C_{u})$, and the
expected cluster volume $\mathbb{E}(K_{u})$ in resels as a function of the
excursion set level $u$ and a parameterization of smoothness. For the expected
excursion set volume, this smoothness parameter is given by the third-order
resel volume. For the remaining panels, including the panels in (C), we
selected the resel volumes of the exemplary data set (cf. Figure 1) given by
$R_{0}(S)=6.0,R_{1}(S)=32.8,R_{2}(S)=353.6,$ and $R_{3}(S)=704.6$ and
evaluated (59) over a range of scalar multiples $wR_{d},d=0,1,2,3$. The value
of $w\in[0,5]$ is displayed on the $x$-axis. Note that smoothness increases
with decreasing resel volumes. (C) Excursion set feature probability
distributions for a $T$-field with $15$ degrees of freedom as functions of
smoothness. From the left to right, the panels show the probability for the
maximum of the statistical random field to exceed a value $u$,
$\mathbb{P}(X_{m}\geq u)$, the probability for the number of clusters within
an excursion set to assume a value $c$, $\mathbb{P}(C_{u}=c)$, the log
probability density for a cluster to assume a volume $k$, $f_{K_{u}}(k)$. (D)
The distribution of the number of clusters of a given volume within an
excursion set. The panels depict $\mathbb{P}(C_{\geq k,u}=c)$ as a function of
$u$ and $k$.

#### Probability distribution of excursion set features

(P1) The distribution of the global maximum of a statistical random field

Historically, the distribution of the global maximum of a statistical random
field has been one of the seeds of the RFP framework (Worsley et al.,, 1992).
This is because maximum statistics afford FWER control in multiple testing
scenarios, as we show in Supplement S2. To introduce the distribution of the
global maximum used in RFP, let

$X_{m}:=\max_{x\in S}{X(x)}$ (67)

denote the global maximum of the statistical random field $X(x)$ on a set $S$.
In RFP, the probability distribution of the global maximum is specified in
terms of its distribution function, which is assumed to be of the following
form (e.g., Worsley et al., (1992, eqs. 1,6), Worsley et al., (1996, eq. 3.1))

$F_{X_{m}}:\mathbb{R}\to[0,1],u\mapsto
F_{X_{m}}(u):=1-\mathbb{E}(\chi(E_{u})).$ (68)

The reasoning behind this form is similar to that justifying the approximation
of the expected number of local maxima by the expected Euler characteristic:
first, if the excursion set defining value $u$ is high enough, $E_{u}$ will
not comprise any handles and voids any more. Second, in the limit of even
higher values $u$, there will remain only one connected component comprising
the global maximum or no connected component at all. In this case, the Euler
characteristic of the excursion set can assume only the values $1$ or $0$: if
the global maximum is equal to or larger than the excursion set defining
value, i.e., if $X_{m}\geq u$, the Euler characteristic will assume the value
$1$, i.e., $\chi(E_{u})=1$. If the global maximum is smaller than the
excursion set defining value, i.e., if $X_{m}<u$, the Euler characteristic
will assume the value $0$, i.e., $\chi(E_{u})=0$. That is, in the limit of
high values for $u$, the Euler characteristic $\chi(E_{u})$ of the excursion
set can be conceived as a random variable taking on values in the set
$\\{0,1\\}$ with a probability mass function defined in terms of the mutually
exclusive event probabilities

$\mathbb{P}(\chi(E_{u})=0):=\mathbb{P}(X_{m}<u)\mbox{ and
}\mathbb{P}(\chi(E_{u})=1):=\mathbb{P}(X_{m}\geq u).$ (69)

Evaluating the expected value $\mathbb{E}(\chi(E_{u}))$ in this scenario then
yields

$\mathbb{E}(\chi(E_{u}))=\mathbb{P}(\chi(E_{u})=0)\cdot
0+\mathbb{P}(\chi(E_{u})=1)\cdot 1=\mathbb{P}(X_{m}\geq u).$ (70)

With the definition of the cumulative distribution function of a continuous
random variable (cf. Supplement S1), expression (68) then follows directly.

It is instructive to consider the dependency of the probability for the global
maximum to exceed a value $u$ as a function of the statistical random field’s
smoothness: an increase in the degree of smoothness at constant intrinsic
volumes results in a decrease of the statistical random field’s resel volumes
(cf. Definition 9). A decrease in resel volumes at constant value $u$ in turn
leads to a decrease in the value of the expected Euler characteristic (cf. eq.
(59)). A decrease in the expected Euler characteristic in turn results in an
increase of the probability $\mathbb{P}(X_{m}<u)$ (cf. eq. (68)), which in
turn is equivalent to a decrease of the probability $\mathbb{P}(X_{m}\geq u)$.
In summary, an increase in a statistical random field’s smoothness thus yields
a lower probability for the field’s global maximum to exceed a given value
$u$. We visualize this dependency in Figure 5C. Note that $\mathbb{E}(C_{u})$
and $\mathbb{P}(X_{m}\geq u)$ are identically equal to
$\mathbb{E}(\chi(E_{u}))$.

(P2) The distribution of the number of clusters within an excursion set

The RFP approximation for the distribution of the number of clusters within an
excursion set was introduced by Friston et al., 1994b (, eq. 8) and is of the
form

$C_{u}\sim\mbox{Poiss}(\lambda_{C_{u}}),$ (71)

where

$\lambda_{C_{u}}:=\mathbb{E}(C_{u}).$ (72)

That is, the number of clusters is assumed to be distributed according to a
Poisson distribution with parameter $\mathbb{E}(C_{u})$, which in turn is
available from eq. (64). Friston et al., 1994b  base the approximation for the
distribution of the number of clusters in an excursion set on Adler, (1981,
Theorem 6.9.3, p. 161) and the fact that a Poisson distribution is fully
specified in terms of its expectation parameter. Theorem 6.9.3 in Adler,
(1981) is concerned with the limit probability of an Euler characteristic
number for high values of $u$ and remarks that “Since this number is (…)
essentially the same as the number of components of the excursion set (…) and
the number of local maxima above $u$, these should have the same limiting
distribution (…). This, however, has never been rigorously proven”.
Intuitively, eq. (71) can be understood as reflecting the fact that the
occurrence of clusters in an excursion set under the null hypothesis should
not depend on spatial location, but only depend on an “average spatial
occurrence rate” given by $\mathbb{E}(C_{u})$. By means of the approximation
(64), an increase in the smoothness of a statistical random field results in a
decrease of the expected number of clusters in an excursion set at a fixed
level $u$. In Figure 5D, we visualize the dependency of $\mathbb{P}(C_{u}=c)$
as a function of the statistical random field’s smoothness for a fixed level
$u=6.5$. As evident from the visualization, increases in smoothness result in
a shift of probability mass towards smaller values of $c$.

(P3) The distribution of the volume of clusters within an excursion set

As in definition (65), let $K_{u}$ denote a real-valued random variable that
models the volume of clusters within an excursion set measured either in terms
of Lebesgue measure or third-order resel volume. In the SPM implementation of
RFP, the probability distribution of this random variable is assumed to be
specified in terms of the probability density function (cf. Friston et al.,
1994b (, p. 213, eq. (12)))

$f_{K_{u}}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0},k\mapsto
f_{K_{u}}(k):=\frac{2\kappa}{D}k^{\frac{2}{D}-1}\exp\left(-\kappa
k^{\frac{2}{D}}\right),$ (73)

where

$\kappa:=\left(\Gamma\left(\frac{D}{2}+1\right)\frac{1}{\mathbb{E}(K_{u})}\right)^{\frac{2}{D}}.$
(74)

These expressions are motivated as follows. First, eq. (73) has been derived
by Friston et al., 1994b  based on a result by Nosko, 1969b ; Nosko, 1969a ;
Nosko, (1970). Nosko’s result, which is reported in Adler, (1981, p. 158),
states that the size of the distribution of connected components of a GRF’s
excursion set to the power of $2/D$ has an exponentional distribution with
expectation parameter $\kappa$. As shown in Supplement S3.3, eq. (73) then
follows directly from Nosko’s result by the change of variables theorem (e.g.,
Fristedt and Gray, (2013, p. 136)). Note that Adler, (1981, p. 158) remarks
that “no detailed proofs are given” for Nosko’s result. Second, eq. (74)
follows from the expression for the expected cluster volume eq. (66) and the
fact that an exponential distribution is fully specified in terms of the
reciprocal value of its expectation. It should be noted that the probability
density function eq. (73) applies to GRFs and, in the limit of high degrees of
freedom, to $T$-fields (cf. comments in spm_P_RF.m). Refined results for the
distribution of cluster volumes in excursion sets of $T$\- and $F$-fields were
derived by Cao, (1999), but are not implemented in SPM.

For later reference, we note two formal consequences of the parametric form
eq. (73) for the probability distribution of the cluster volume $K_{u}$.
First, as shown in Supplement S3.3, eq. (73) implies that the cumulative
density function of the cluster volume is given by

$F_{K_{u}}:\mathbb{R}\to[0,1],k\mapsto F_{K_{u}}(k)=1-\exp\left(-\kappa
k^{\frac{2}{D}}\right).$ (75)

Second, from eq. (75) it follows directly, that the probability of the cluster
volume to exceed a constant $k$ under the RFP framework evaluates to (cf.
Friston et al., 1994b (, eq. 11))

$\mathbb{P}(K_{u}\geq k)=\exp\left(-\kappa k^{\frac{2}{D}}\right).$ (76)

In Figure 5C, we visualize the log probability density function $f_{K_{u}}$
for a fixed value $u=6.5$ as a function of the statistical random field’s
smoothness. The allocation of probability density mirrors the effect of
smoothness on the expected volume of clusters within the excursion set
$\mathbb{E}(K_{u})$ in Figure 5B.

(P4) The distribution of the number of clusters of a given volume within an
excursion set

Let $C_{\geq k,u}$ denote a random variable that models the number of clusters
within an excursion set of level $u$ that have a volume equal to or larger
than some constant $k$. Under the RFP framework, it is assumed that (cf.
Friston et al., (1996, eq. 4))

$C_{\geq k,u}\sim\mbox{Poiss}(\lambda_{C_{\geq k,u}}),$ (77)

where

$\lambda_{C_{\geq k,u}}:=\mathbb{E}(C_{u})\mathbb{P}(K_{u}\geq k).$ (78)

That is, the random variable $C_{\geq k,u}$ is assumed to be distributed
according to a Poisson distribution, the expectation parameter of which is
given by the product of the expected number of clusters and the probability of
a cluster volume to assume a value equal to or larger than the constant $k$.
Note again that parametric forms for $\mathbb{E}(C_{u})$ and
$\mathbb{P}(K_{u}\geq k)$ are available from eqs. (64) and (76), respectively.
Further note that eq. (77) specifies probabilities of the form
$\mathbb{P}(C_{\geq k,u}=i)$, which should be read as representing the
probability that the number of clusters within an excursion set defined by the
level value $u$ with a volume larger than or equal to $k$ is $i$. For later
reference, we note that the sum $\sum_{i=0}^{c-1}\mathbb{P}(C_{\geq k,u}=i)$
thus represents the probability that the number of clusters within an
excursion set defined by level $u$ with volume larger than or equal to $k$ is
either smaller than or equal to a constant $c\in\mathbb{N}$ minus $1$, and
hence the expression $1-\sum_{i=0}^{c-1}\mathbb{P}(C_{\geq k,u}=i)$ represents
the probability that the number of clusters within an excursion set defined by
level $u$ having volume equal to or larger than $k$ is equal to or larger than
$c$.

Expressions (77) and (78) can be motivated as follows. Assume the existence of
a random variable $C_{u}$ that models the number of clusters within an
excursion set at level $u$ and the existence of a parametric expression for
its distribution (provided by eq. (71) and eq. (72) in the theory of RFP).
Assume further that one would like to specify the distribution
$\mathbb{P}(C_{\geq k,u})$ of a second random variable $C_{\geq k,u}$ that
models the number of clusters in an excursion set $E_{u}$ at level $u$ that
have a volume equal to or larger than some constant $k$. By definition, this
distribution is a marginal distribution of the joint distribution

$\mathbb{P}(C_{u},C_{\geq k,u})=\mathbb{P}(C_{u})\mathbb{P}(C_{\geq
k,u}|C_{u}).$ (79)

Define the conditional distribution of $C_{\geq k,u}$ given $C_{u}$ by the
binomial distribution

$C_{\geq k,u}|C_{u}\sim\mbox{BN}_{j}(\mu),\mbox{ with
}\mu:=\mathbb{P}(K_{u}\geq k)$ (80)

such that by definition

$\mathbb{P}(C_{\geq k,u}=i|C_{u}=j)=\left(\begin{matrix}j\\\
i\end{matrix}\right)\mu^{i}(1-\mu)^{j-i},\mbox{ for
}i\in\mathbb{N}_{j}^{0}\mbox{ and }j\in\mathbb{N}.$ (81)

Intuitively, the event that $i$ of $j$ realized clusters (i.e., $i=0,1,...,j$)
have a volume equal or larger to $k$ is thus conceived as one realization of
$j$ independent Bernoulli events with “success probability”
$\mathbb{P}(K_{u}\geq k)$. Based on eq. (81), it can then be shown that (cf.
Friston et al., (1996, Appendix))

$\mathbb{P}(C_{\geq k,u}=i)=\sum_{j=1}^{\infty}\mathbb{P}(C_{u}=i,C_{\geq
k,u}=j)=\frac{\lambda_{C_{\geq k,u}}^{i}\exp\left(-\lambda_{C_{\geq
k,u}}\right)}{i!},$ (82)

i.e., that such a random variable $C_{\geq k,u}$ has a Poisson distribution
with expectation parameter $\lambda_{C_{\geq k,u}}$. We document the
derivation of eq. (82) in Supplement S3.3. The importance of the random
variable $C_{\geq k,u}$ in the SPM implementation of RFP will be elaborated on
in Section 4.2. Finally, in Figure 5D, we visualize the distribution of
$C_{\geq k,u}$ as a function of smoothness and the values of $u$ and $k$. At
constant smoothness, increasing $u$ and $k$ leads to a shift of probability
mass to lower values of $c$.

(5) The distribution of the maximal cluster volume within an excursion set

Finally, to afford FWER control in multiple testing scenarios relating to the
volume of clusters within excursion sets, a parametric form of the respective
maximum statistic is required (cf. Supplement S2). To this end, let, for a set
of $c\in\mathbb{N}$ of clusters within an excursion set $E_{u}$ at level $u$,
$k_{i},i=1,...,c$ denote the volumes of clusters $i=1,...,c$. Further, let

$K_{u}^{m}:=\max\\{k_{1},...,k_{c}\\}$ (83)

denote the maximal cluster volume within the excursion set. Then, with the
definition of the random variable $C_{\geq k,u}$, we have for $k\in\mathbb{R}$

$\displaystyle\begin{split}\mathbb{P}(K_{u}^{m}<k)=\mathbb{P}(C_{\geq
k,u}=0)\mbox{ and }\mathbb{P}(K_{u}^{m}\geq k)=\mathbb{P}(C_{\geq k,u}\geq
1).\end{split}$ (84)

In words, first, the probability that the maximal cluster volume within an
excursion set $E_{u}$ at level $u$ is smaller than a constant $k$ is the same
as the probability that the number of clusters in the excursion set $E_{u}$
which have a volume equal to or larger than $k$ is zero. Second, the
probability that the maximal cluster volume in an excursion set $E_{u}$ is
equal to or larger than a constant $k$ is the same as the probability that the
number of clusters in the excursion set $E_{u}$ which have a volume equal to
or larger than $k$ is equal to or larger than $1$. The distribution of the
maximal cluster volume $K_{u}^{m}$ is thus fully determined by the
distribution of $C_{\geq k,u}$. For later reference, we note that with eqs.
(77), (78), and (82) (cf. (Friston et al., 1994b, , eq. 14))

$\displaystyle\begin{split}\mathbb{P}(K_{u}^{m}\geq k)&=\mathbb{P}(C_{\geq
k,u}\geq 1)\\\ &=1-\sum_{i=0}^{0}\mathbb{P}(C_{\geq k,u}=i)\\\
&=1-\exp\left(-\mathbb{E}(C_{u})\mathbb{P}(K_{u}\geq k)\right).\end{split}$
(85)

##  4 Application

The general theory of random fields - and the theory of RFP in particular -
model three-dimensional space by $\mathbb{R}^{3}$. In reality, however, fMRI
data points are only acquired at a finite number of discrete spatial
locations, commonly known as voxels. This discretization of space has to be
accounted for when considering the application of RFP theory to GLM-based fMRI
data analysis. In the current section, we thus consider the discrete space,
discrete data-point linear model

$Y_{iv}=m_{i}\beta_{v}+\sigma Z_{iv}\mbox{ for }i=1,...,n\mbox{ and
}v\in\mathcal{V},$ (86)

where

$\mathcal{V}:=\left\\{v|v=(v_{1},v_{2},v_{3})\mbox{ with
}v_{1}\in\mathbb{N}_{m_{1}},v_{2}\in\mathbb{N}_{m_{2}},v_{3}\in\mathbb{N}_{m_{3}}\right\\}$
(87)

denotes a voxel index set. Eq. (86) should be read as the discrete-space
analogue of eq. (38) for which the continuous space coordinates
$x\in\mathbb{R}^{3}$ have been replaced by voxel indices $v\in\mathcal{V}$. We
assume that these voxel indices correspond to the weighting parameters of
three orthogonal, but not necessarily orthonormal, basis vectors
$b_{d}\in\mathbb{R}^{3},d=1,2,3$ of $\mathbb{R}^{3}$, such that the set

$L:=\left\\{\sum_{d=1}^{3}v_{d}b_{d}|v_{1}\in\mathbb{N}_{m_{1}},v_{2}\in\mathbb{N}_{m_{2}},v_{3}\in\mathbb{N}_{m_{3}}\right\\}$
(88)

is a finite lattice covering the subset of $\mathbb{R}^{3}$ of interest.
Assuming an MNI space bounding box of -78 to 78 mm in the coronal, -112 mm to
76 mm in the sagittal, and -70 mm to 95 mm in the axial direction and a 2
$\times$ 2 $\times$ 2 mm voxel size, typical values for the number of indices
are $m_{1}=79$, $m_{2}=95$, and $m_{3}=69$. With these definitions, we thus
have the correspondences

$Y_{iv}=Y_{i}(\tilde{x}),\beta_{v}=\beta_{i}(\tilde{x}),\mbox{ and
}Z_{iv}=Z_{i}(\tilde{x})\mbox{ for }\tilde{x}:=\sum_{d=1}^{3}v_{d}b_{d}\mbox{
and }v\in\mathcal{V}$ (89)

between the discrete space and the continuous space linear model
representations of eqs. (86) and (38), respectively.

Building on these correspondences between the theoretical scenario of a
continuous space and the applied scenario of discrete spatial sampling, we are
now in a position to discuss the application of the theoretical framework of
the previous section to GLM-based fMRI data analysis. As outlined in the
Introduction, we first consider the parameter estimation of the RFP null
model, which comprises the estimation of the statistical random field’s
smoothness and the approximation of the search space’s intrinsic and resel
volumes. We then consider the evaluation of the parametric p-values reported
in the SPM results table.

###  4.1 Parameter estimation

As discussed in Section 3.2, the probability distributions of all excursion
set features of a statistical random field depend on the resel volumes of the
search space via the expected Euler characteristic (cf. eq. (59)). The resel
volumes in turn depend on the FWHM parameterization of the statistical random
field’s smoothness and the search space’s intrinsic volumes (cf. Definition
(9)). Application of the parametric forms for the probability distributions of
excursion set features discussed in Section 3.2 in a data-analytical context
thus necessitates the estimation of the statistical random field’s smoothness,
the approximation of its search space’s intrinsic volumes, and the subsequent
evaluation of its resel volumes. In effect, the estimated resel volumes are
intended to furnish a data-adaptive null hypothesis model against which the
actually observed data can be evaluated. Because the data typically do not
fully conform to the null hypothesis $\beta_{v}=0$ for all $v\in\mathcal{V}$,
the SPM implementation of RFP bases its parameter estimation scheme on the so-
called standardized residuals

$r_{iv}:=\frac{Y_{iv}-m_{i}\hat{\beta_{v}}}{\sqrt{\frac{1}{n-p}\sum_{i=1}^{n}\left(Y_{iv}-m_{i}\hat{\beta_{v}}\right)^{2}}},i=1,...,n,v\in\mathcal{V},$
(90)

where $\hat{\beta}_{v}$ denotes the Gauss-Markov estimator for the effect size
parameter at voxel $v$. Practically, these values are stored in ResI_000*.nii
files, which are created and deleted by SPM’s spm_spm.m function during SPM’s
GLM estimation. Note that the denominator of the standardized residuals in eq.
(90) corresponds to the estimate $\hat{\sigma}$ of the standard deviation
parameter $\sigma$. Under the assumption that $m_{i}\hat{\beta_{v}}$
veridically reflects the expectation of $Y_{iv}$, the standardized residuals
can thus be conceived as realizations of the $Z$-fields $Z_{i}(x)$ at discrete
spatial locations (cf. eq. (38) and the introduction of Section 3). Stated
differently, the RFP framework assumes that the standardized residuals can be
considered observations of the null hypothesis model used to evaluate the
probabilities of observed excursion set features. RFP parameter estimation as
implemented in SPM’s spm_est_smoothness.m function then comprises three steps:
first, the estimation of the FWHM smoothness parameters based on the
standardized residuals, second, the approximation of the intrinsic volumes of
the search space, and finally, the evaluation of the resel volumes.

#### FWHM smoothness parameter estimation

SPM’s estimation scheme for the FWHM smoothness parameters can be conceived as
a two-step procedure. In the first step, an estimate of the smoothness
parameter $\varsigma$ (cf. eq. (40)) of the $Z$-fields $Z_{i}(x)$ is
constructed based on the standardized residuals. In a second step, this
estimate is used together with an estimate of the approximating matrix
$\Lambda$ to afford the estimation of FWHM smoothness parameters that conform
to the equality of eq. (48) and, in their relative magnitude, reflect the
diagonal entries of the variance matrix of partial derivatives (cf. eq. (41)).
More specifically, the smoothness parameter estimate of the first step is
evaluated as

$\hat{\varsigma}:=\frac{1}{\frac{1}{m}\sum_{v=1}^{m}|\frac{n}{\tilde{n}(n-p)}V_{v}|^{1/2}},$
(91)

where $V_{v}\in\mathbb{R}^{3\times 3}$ denotes a scaled voxel-wise estimate of
the variance matrix of gradient components with elements

$(V_{v})_{jk}:=\sum_{i=1}^{\tilde{n}}{(r_{i(v+e_{j})}-r_{iv})(r_{i(v+e_{k})}-r_{iv})}\,\,\mbox{
for }j,k=1,2,3\mbox{ and }v\in\mathcal{V},$ (92)

$m=m_{1}+m_{2}+m_{3}$ denotes the cardinality of the voxel index set, and
$\tilde{n}\leq n$ denotes the size of a data point subset potentially used by
SPM for computational efficiency. The forms of eqs. (91) and (92) are
motivated as follows (Hayasaka and Nichols,, 2003; Hayasaka et al.,, 2004):
the $(V_{v})_{jk}$ are products of the first-order finite differences between
the standardized residuals of adjacent voxel indices. These first-order finite
differences may be interpreted as scaled versions of numerical first-order
derivatives. Notably, SPM thus foregoes the direct approximation of the
gradient components of the statistical random field. The first order
differences in the factors in eq. (92) are evaluated by spm_sample_vol.c. The
formation of the sum of products of these first-order finite differences of
standardized residuals in eq. (92) and its subsequent division by $\tilde{n}$
can be interpreted as an estimation of the (co)variances of the gradient
components, if it is assumed that the sample average over available data
points $\tilde{n}$ is zero. The subsequent multiplication by $n$ and division
by the degrees of freedom $n-p$ is meant to rescale the respective estimates
to the entire set of data points and render the estimate unbiased by
effectively converting standardized to normalized residuals (cf.
spm_est_smoothness.m, comments ll. 144 - 169 and Worsley, 1996b ). In effect,
the sum terms in the denominator of eq. (91) can be interpreted as scaled
voxel-specific numerical estimates of the square root of the determinants of
$\mathbb{V}(\nabla X(x)),x\in S$. Finally, the summation of these terms and
division by the number of voxel indices $m$ can be interpreted as the
numerical integration over the unit cuboid for first-order finite differences,
rather than numerical derivative approximations to the statistical random
fields’ gradient components. In Figure 6 we demonstrate that a numerically
simplified version of SPM’s smoothness parameter estimation approach indeed
results in a valid recovery of the smoothness parameter of a two-dimensional
Gaussian random field with Gaussian covariance function.

Figure 6: Smoothness parameter estimation. To exemplify the smoothness
estimation scheme of SPM’s RFP implementation, we capitalize on the fact that
for one-dimensional GRFs with Gaussian covariance function, the reciprocal
unit-square Lipschitz-Kiling curvature $\varsigma$ is identical to the
Gaussian covariance function’s length parameter $\ell$ up to a factor of
$2^{-1/2}$ (cf. eq. (43)). A straightforward estimator for $\ell$ is thus
$\hat{l}:=\sqrt{2}\varsigma$. For the simulations depicted in the figure, we
evaluated $\varsigma$ using numerical differentiation and sample variances to
estimate $\ell$ for one-dimensional GRFs on the interval $[0,1]$. We varied
$\ell$ in the interval $[0.04,0.2]$ and examined sample sizes in the range
from $n=100$ to $n=2000$. (A) Panel A visualizes a sample of $n=100$ rough GRF
realizations with $\ell=0.04$ . (B) Panel B visualizes a sample of $n=100$
smooth GRF realizations with $\ell=0.2$. (C) Panel C visualizes the absolute
percent relative estimation error $|100\cdot(\hat{\ell}-\ell)/\ell|$ for the
ensuing estimates. In general, this error is small ($\lesssim 3$ % on average)
and decreases for larger sample sizes uniformly over the space of $\ell$. For
the full implementational details of these simulations, please see rfp_6.m.

The reparameterization of smoothness in terms of FWHM smoothness parameters in
the second step then is achieved by setting

$\tilde{f}_{x_{j}}:=\frac{1}{4\ln
2}\left(\frac{1}{m}\sum_{v=1}^{m}\frac{n}{\tilde{n}(n-p)}(V_{v})_{jj}\right)^{-\frac{1}{2}}$
(93)

and

$\hat{f}_{x_{j}}=\frac{\hat{\varsigma}^{\frac{1}{3}}}{(4\ln
2)^{-\frac{1}{2}}}\frac{\tilde{f}_{x_{j}}}{\left(\prod_{k=1}^{3}\tilde{f}_{x_{j}}\right)^{\frac{1}{3}}}$
(94)

for $j=1,2,3$. Eqs. (93) and (94) implement SPM’s assumption that the variance
matrix of gradient components is of diagonal form and thus that the smoothness
parameter $\varsigma$ can be reparameterized in terms of the FWHM smoothness
parameters $f_{x_{j}}$ (cf. eq. (47)). More specifically, these equations are
motivated as follows (cf. spm_est_smoothness.m, comments ll. 246 - 248): the
$\tilde{f}_{x_{j}}$ can be interpreted as estimates of the diagonal elements
of $\Lambda$ based on the $\mathbb{V}(\nabla X(x))$ matrix estimate discussed
above. By normalizing these values using their geometric average in the second
factor of eq. (94) and scaling the resulting normalized $\tilde{f}_{x_{j}}$’s
by the first factor, it is then ensured that the $\hat{f}_{x_{j}}$’s reproduce
the smoothness parameter estimate $\hat{\varsigma}$ if they are substituted as
estimates of the $f_{x_{j}}$’s in the matrix $\Lambda$ (cf. eq. (48)). The
thus evaluated estimates are then used for the evaluation of the resel volumes
and are reported in the footnote of the SPM results table as FWHM (voxels).

#### Intrinsic volume approximation

To approximate the intrinsic volumes of the search space, SPM uses an
algorithm proposed by Worsley et al., (1996, pp. 62 -63) and implemented in
spm_resels_vol.c. This algorithm corresponds to a summation of the
contributions of zero- to three-dimensional subspaces along the cardinal axes
to the respective intrinsic volumes of the search space. More specifically,
the algorithm evaluates

  * •

the total number of voxels

$m:=m_{x_{1}}m_{x_{2}}m_{x_{3}},$ (95)

  * •

the number of edges in the $x_{1},x_{2}$ and $x_{3}$ directions, i.e.,

$\displaystyle\begin{split}m_{1}^{e}&:=|\\{(v,v+e_{1})|v,v+e_{1}\in\mathcal{V}\\}|,\\\
m_{2}^{e}&:=|\\{(v,v+e_{2})|v,v+e_{2}\in\mathcal{V}\\}|,\\\
m_{3}^{e}&:=|\\{(v,v+e_{3})|v,v+e_{3}\in\mathcal{V}\\}|,\\\ \end{split}$ (96)

respectively,

  * •

the number of faces in the $x_{1}\text{-}x_{2}$-, $x_{1}\text{-}x_{3}$-, and
$x_{2}\text{-}x_{3}$\- directions, i.e.,

$\displaystyle\begin{split}m_{12}^{f}&:=|\\{(v,v+e_{1},v+e_{2},v+e_{1}+e_{2})|v,v+e_{1},v+e_{2},v+e_{1}+e_{2}\in\mathcal{V}\\}|,\\\
m_{13}^{f}&:=|\\{(v,v+e_{1},v+e_{3},v+e_{1}+e_{3})|v,v+e_{1},v+e_{3},v+e_{1}+e_{3}\in\mathcal{V}\\}|,\\\
m_{23}^{f}&:=|\\{(v,v+e_{2},v+e_{3},v+e_{2}+e_{3})|v,v+e_{2},v+e_{3},v+e_{2}+e_{3}\in\mathcal{V}\\}|,\\\
\end{split}$ (97)

respectively, and

  * •

the number of cubes, i.e.,

$\displaystyle\begin{split}m^{c}:=|\\{(&v,v+e_{1},v+e_{2},v+e_{3},v+e_{1}+e_{2},v+e_{2}+e_{3},v+e_{1}+e_{3},v+e_{1}+e_{2}+e_{3})\\\
|&v,v+e_{1},v+e_{2},v+e_{3},v+e_{1}+e_{2},v+e_{2}+e_{3},v+e_{1}+e_{3},v+e_{1}+e_{2}+e_{3}\in\mathcal{V}\\}|.\\\
\end{split}$ (98)

Based on these counts and the physical voxel sizes $\delta_{1},\delta_{2}$ and
$\delta_{3}$ in the $x_{1}$-, $x_{2}$-, and $x_{3}$-directions, respectively,
the intrinsic volumes of the search space are then approximated by

$\displaystyle\begin{split}\tilde{\mu}_{0}(S)&=m-(m_{1}^{e}+m_{2}^{e}+m_{3}^{e})+(m_{12}^{f}+m_{13}^{f}+m_{23}^{f})-m^{c},\\\
\tilde{\mu}_{1}(S)&=\delta_{1}(m_{1}^{e}-m_{12}^{f}-m_{13}^{f}+m^{c})+\delta_{2}(m_{2}^{e}-m_{12}^{f}-m_{23}^{f}+m^{c})+\delta_{3}(m_{3}^{e}-m_{13}^{f}-m_{23}^{f}+m^{c}),\\\
\tilde{\mu}_{2}(S)&=\delta_{1}\delta_{2}(m_{12}^{f}-m^{c})+\delta_{1}\delta_{3}(m_{13}^{f}-m^{c})+\delta_{2}\delta_{3}(m_{23}^{f}-m^{c}),\\\
\tilde{\mu}_{3}(S)&=\delta_{1}\delta_{2}\delta_{3}m^{c}.\end{split}$ (99)

Note, for example, that $\tilde{\mu}_{3}(S)$ corresponds to a summation of the
volume $\delta_{1}\delta_{2}\delta_{3}$ of each of the $m^{c}$ voxels
constituting the search space. For a whole-brain GLM-based fMRI data analysis,
this quantity thus typically corresponds to an approximation of the physical
brain volume by many small cuboids.

#### Resel volume estimation

De-facto, SPM does not evaluate an approximation of the search space’s
intrinsic volumes perse. Instead, it capitalizes on the subspace decomposition
assumption of eq. (36) and directly estimates the resel volumes in
spm_resel_vol.m by

$\displaystyle\begin{split}\hat{R}_{0}(S)&=\tilde{\mu}_{0}(S),\\\
\hat{R}_{1}(S)&=\left(\begin{matrix}\frac{\delta_{1}}{\hat{f}_{x_{1}}}&\frac{\delta_{2}}{\hat{f}_{x_{2}}}&\frac{\delta_{3}}{\hat{f}_{x_{3}}}\end{matrix}\right)\left(\begin{matrix}m_{1}^{e}-m_{12}^{f}-m_{13}^{f}+m^{c}\\\
m_{2}^{e}-m_{12}^{f}-m_{23}^{f}+m^{c}\\\
m_{3}^{e}-m_{13}^{f}-m_{23}^{f}+m^{c}\end{matrix}\right),\\\ \\\
\hat{R}_{2}(S)&=\left(\begin{matrix}\frac{\delta_{1}\delta_{2}}{\hat{f}_{x_{1}}\hat{f}_{x_{2}}}&\frac{\delta_{1}\delta_{3}}{\hat{f}_{x_{1}}\hat{f}_{x_{3}}}&\frac{\delta_{2}\delta_{3}}{\hat{f}_{x_{2}}\hat{f}_{x_{3}}}\end{matrix}\right)\left(\begin{matrix}m_{12}^{f}-m^{c}\\\
m_{13}^{f}-m^{c}\\\ m_{23}^{f}-m^{c}\\\ \end{matrix}\right),\mbox{and }\\\
\hat{R}_{3}(S)&=\frac{1}{\hat{f}_{x_{1}}\hat{f}_{x_{2}}\hat{f}_{x_{3}}}\tilde{\mu}_{3}(S).\end{split}$
(100)

The estimated resel volumes $\hat{R}_{d}(S),d=0,1,2,3$ are the pivot points
between the practical perspective of the current section and the theoretical
perspective of Section 3: Theory. As previously noted, virtually all
probability distributions of the topological features of interest in RFP
depend on the resel volumes via the expected Euler characteristic (cf. eq.
(59)). Substitution of $\hat{R}_{d}(S),d=0,1,2,3$ for the resel volumes in the
expected Euler characteristic formula thus results in a set of probability
distributions that reflect a data-adapted statistical random field null model.
The exceedance probabilities of the observed topological features of interest
under this data-adapted null model are then reported in the SPM results table.

This concludes our discussion of the RFP parameter estimation scheme
implemented in SPM. We provide a reproduction of this estimation scheme in a
simplified and annotated version of spm_est_ smoothness.m as
rfp_spm_est_smoothness.m in the accompanying OSF project.

###  4.2 The SPM results table

The majority of p-values in the SPM results table relate to the random
variable $C_{\geq k,u}$ introduced in Section 3.2 (Friston et al.,, 1996).
Recall that $C_{\geq k,u}$ models the number of clusters within a level
$u$-dependent excursion set $E_{u}$ that have a volume equal to or larger than
a constant $k\in\mathbb{R}$. As discussed in Section 3.2, the random variable
$C_{\geq k,u}$ takes on values in $\mathbb{N}^{0}$ and is distributed
according to a Poisson distribution with parameter $\lambda_{C_{\geq
k,u}}:=\mathbb{E}(C_{u})\mathbb{P}(K_{u}\geq k)$, where $\mathbb{E}(C_{u})$
represents the expected number of clusters within an excursion set $E_{u}$ and
is given by eq. (64), while $\mathbb{P}(K_{u}\geq k)$ denotes the probability
for a cluster within an excursion set $E_{u}$ to assume a volume equal to or
larger than the constant $k$ and is given by eq. (76). Formally, we have

$C_{\geq k,u}\sim\mbox{Poiss}(\lambda_{C_{\geq k,u}}),$ (101)

and thus

$\mathbb{P}(C_{\geq k,u}=c)=\exp\left(-\mathbb{E}(C_{u})\mathbb{P}(K_{u}\geq
k)\right)\frac{\left(\mathbb{E}(C_{u})\mathbb{P}(K_{u}\geq
k)\right)^{c}}{c!}\mbox{ for }c=0,1,...\mbox{.}$ (102)

For ease of notation, we denote the distribution function of $C_{\geq k,u}$ in
the current section by

$F_{\lambda_{C_{\geq k,u}}}:\mathbb{N}^{0}\to[0,1],c\mapsto
F_{\lambda_{C_{\geq k,u}}}(c):=\sum_{j=0}^{c}\mathbb{P}(C_{\geq k,u}=j)$ (103)

(cf. Supplement S1). The majority of p-values in the SPM results table are
then special cases of probabilities related to $C_{\geq k,u}$ for the values
$c\in\mathbb{N}^{0},k\in\mathbb{R}$, and $u\in\mathbb{R}$ and are evaluated
using the implementation of the Poisson distribution function (103) in
spm_Pcdf.m. To distinguish these special cases, we make the following
assumptions and use the following notation. We assume that

  1. (1)

a statistical parametric map, for example a three-dimensional voxel map
comprising $T$ value realizations, has been obtained,

  2. (2)

for this map, a cluster-forming threshold value $\tilde{u}$ and a cluster-
extent threshold value $\tilde{k}$ in resel volume have been defined,

  3. (3)

the number of clusters $\tilde{c}$ in the resulting excursion set
$E_{\tilde{u}}$ has been evaluated,

  4. (4)

the volumes of these clusters $\tilde{k}_{i},i=1,...,\tilde{c}$ have been
assessed, and

  5. (5)

for each cluster, the maximum test statistic within this cluster has been
evaluated and is denoted by $\tilde{t}_{i},i=1,...,\tilde{c}$.

A few remarks from a practical perspective on these assumptions are in order.
In SPM, statistical parametric maps for $T$ tests are evaluated based on GLM
beta and variance parameter estimates for user-defined contrast weight vectors
and are saved under the filename spmT_000*.nii or the like. When interacting
with SPM’s Results graphical user interface, or specifying options in SPM’s
batch editor for the Results Report module, the user is prompted to specify
the cluster-forming threshold value $\tilde{u}$ by selecting a Threshold Type
with options “None” and “FWE” for p-value adjustment, and a Threshold Value in
the form of a p-value. The combination of these two options is evaluated by
SPM to define the value $\tilde{u}$ as follows. In the case of “p-value
adjustment: none” and the subsequent specification of a p-value (e.g. 0.001),
SPM uses the function spm_u.m to evaluate the inverse cumulative density
function for the probabilistic complement of the specified value and the
respective statistic. For example, in the case of a statistical parametric map
comprising $T$-values, $\tilde{u}$ thus corresponds to the “critical
$T$-value”, i.e., the value $t_{c}$ for which values of the $T$-statistic
equal to or larger than $t_{c}$ have a probability of equal to or smaller than
the specified p-value under the standard $T$\- distribution with the
appropriate degrees of freedom (e.g. for a $T$ distribution with 30 degrees of
freedom and p-values of 0.05, 0.025, and 0.001, the respective $\tilde{u}$
values evaluate to 1.697,2.042, and 3.385, respectively). In the case of
“p-value adjustment: FWE” and the subsequent specification of a p-value (e.g.
0.05), SPM uses the function spm_uc.m to evaluate a critical threshold using
the RFP framework. This option leads, in general, to higher values for
$\tilde{u}$. The cluster-extent threshold $\tilde{k}$ is specified by the user
in units of voxels and is transformed into units of resels by SPM by
multiplying the user-specified value by $1/\prod_{d=1}^{3}\hat{f}_{x_{d}}$.
Finally, the excursion set $E_{\tilde{u}}$, the number of clusters
$\tilde{c}$, and the cluster-specific volume and maxima statistics
$\tilde{k}_{i}$ and $\tilde{t}_{i},i=1,...,\tilde{c}$ are evaluated using low-
level SPM functionality, such as spm_bwlabel.c, spm_clusters.m, and spm_max.m.
For the computational details of these routines, the reader is referred
directly to the documentation in the functions’ headers. Based on an
interaction between the spm_list.m and spm_P_RF.m functions, SPM then
documents p-values in the SPM results table for set-level, cluster-level and
peak-level inferences (Friston et al.,, 1996). In the following, we specify
for each level of inference which probabilities these p-values correspond to.

#### Set-level inference

Assume we are interested in the probability

$\mathbb{P}(C_{\geq\tilde{k},\tilde{u}}\geq\tilde{c}),$ (104)

i.e., the probability that for the specified cluster-forming threshold
$\tilde{u}$, the number of clusters within the excursion set $E_{\tilde{u}}$
that have a volume equal to or larger than the specified cluster-extent
threshold value $\tilde{k}$ is equal to or larger than the observed number of
clusters $\tilde{c}$. Clearly, we have

$\displaystyle\begin{split}\mathbb{P}(C_{\geq\tilde{k},\tilde{u}}\geq\tilde{c})&=1-\mathbb{P}(C_{\geq\tilde{k},\tilde{u}}<\tilde{c})\\\
&=1-F_{\lambda_{C_{\geq\tilde{k},\tilde{u}}}}(\tilde{c})\end{split}$ (105)

and eq. (104) can thus be evaluated using the Poisson distribution function
(103). In the SPM results table, the probability
$\mathbb{P}(C_{\geq\tilde{k},\tilde{u}}\geq\tilde{c})$ is reported as the
“set-level p value” and the value $\tilde{c}$ is documented as the number of
clusters c (Figure 1 and Table 3). The idea behind using (104) is to enable
inferences about “network-activations” (Friston et al.,, 1996, pp. 229-230):
if the number of actually observed clusters $\tilde{c}$ exceeds the number
$\mathbb{E}(C_{u})$ of clusters expected under the RFP null hypothesis model
by a large margin, eq. (104) evaluates to a small value. In this case, it may
thus be concluded that the observed pattern of distributed clusters reflects
the engagement of a distributed network of brain regions in the GLM contrast
of interest. Note that the set-level inference probability (104) is the only
p-value of the SPM results table that directly depends on the user specified
cluster-extent threshold $\tilde{k}$.

#### Cluster-level inference

Assume that for a cluster of the excursion set $E_{\tilde{u}}$ indexed by $i$,
where $i=1,...,\tilde{c}$, we are interested in the probability

$\mathbb{P}(C_{\geq\tilde{k_{i}},\tilde{u}}\geq 1),$ (106)

i.e., the probability that for the specified cluster-forming threshold
$\tilde{u}$ the number of clusters within the excursion set that have a volume
equal to or larger than the observed cluster volume $\tilde{k}_{i}$ is equal
to or larger than $1$. For this probability, we have

$\displaystyle\begin{split}\mathbb{P}(C_{\geq\tilde{k}_{i},\tilde{u}}\geq
1)&=1-\mathbb{P}(C_{\geq\tilde{k_{i}},\tilde{u}}<1)\\\
&=1-F_{\lambda_{C_{\geq\tilde{k}_{i},\tilde{u}}}}(0)\\\
&=1-\mathbb{P}(C_{\geq\tilde{k}_{i},\tilde{u}}=0)\\\
&=1-\exp\left(-\mathbb{E}(C_{\tilde{u}})\mathbb{P}(K_{\tilde{u}}\geq\tilde{k}_{i})\right)\\\
&=\mathbb{P}\left(K_{\tilde{u}}^{m}\geq\tilde{k}_{i}\right).\end{split}$ (107)

From the second equality of eq. (107), we thus see that (106) can be evaluated
using the Poisson distribution function (103). Moreover, from the last
equality, we see that $\mathbb{P}(C_{\geq\tilde{k_{i}},\tilde{u}}\geq 1)$