Spatial Confidence Regions for Combinations of Excursion Sets in Image Analysis

In the following sections, we shall need notation to describe the numerous possible low dimensional sub-manifolds which can arise from intersecting the boundaries of $M$ excursion sets. To do so, we first denote the power set of a finite set, $A$, as $\mathcal{P}(A)$, and define $\mathcal{P}^{+}(A)$ as $\mathcal{P}^{+}(A)=\mathcal{P}(A)\setminus\{\emptyset\}$. For each $i\in\mathcal{M}$, we define $\partial\mathcal{A}_{c}^{i}$ as the level set $\{s\in S:\mu^{i}(s)=c\}$. Similarly, we define $\partial\mathcal{F}_{c}$ as the level set $\{s\in S:\min_{i\in\mathcal{M}}\mu^{i}(s)=c\}$. For a spatial set, $B\subseteq S$, its complement in $S$ is denoted $B^{c}=S\setminus B$, its topological closure is denoted $\overline{B}$, its interior is denoted $B^{\circ}$ and its boundary is denoted $\partial B$. In addition, for closed $B$, the notation $B^{1}$ and $B^{-1}$ shall represent $B$ and $\overline{B^{c}}$, respectively. We note here that the notation $\partial\mathcal{A}^{i}_{c}$ and $\partial B$ may conflict with one another. This conflict is resolved, however, by the assumptions of Section 2.2 which ensure that the boundaries of $\{\mathcal{A}_{c}^{i}\}_{i\in\mathcal{M}}$ and $\mathcal{F}_{c}$ are equal to the level sets $\{\partial\mathcal{A}_{c}^{i}\}_{i\in\mathcal{M}}$ and $\partial\mathcal{F}_{c}$ (at least locally, in the vicinity of $\partial\mathcal{F}_{c}$, which is sufficient for our purposes).

We can now partition the level set $\partial\mathcal{F}_{c}$ into sub-manifolds using the set of possible intersections of the $M$ excursion set boundaries. For all $\alpha\in\mathcal{P}^{+}(\mathcal{M})$, we define the boundary segment $\partial^{\alpha}\mathcal{F}_{c}$ as follows;

$$\partial^{\alpha}\mathcal{F}_{c}=\partial\mathcal{F}_{c}\cap\bigg{(}\bigcap_{i\in\alpha}\partial\mathcal{A}_{c}^{i}\bigg{)}\cap\bigg{(}\bigcap_{j\in\mathcal{M}\setminus\alpha}(\partial\mathcal{A}_{c}^{j})^{c}\bigg{)}.$$

We note it is possible that, for some $\alpha\in\mathcal{P}^{+}(\mathcal{M})$, the boundary segment $\partial^{\alpha}\mathcal{F}_{c}$ is empty (in fact, this is often the case in practice). This notation is crucial to the statement of Theorem 2.1 and is illustrated by Fig. 2.

In this section, we outline and discuss the assumptions upon which our theory relies. Additional discussion of why each assumption is required, alongside illustration, is given in Supplementary Theory Section S2.

Each of the examples presented in Fig. 2 satisfy Assumptions 2.2.2 and 2.2.3 and highlight several features worthy of note. Firstly, we do not assume $\partial\mathcal{F}_{c}$ is a $C^{1}$ curve, but rather piece-wise $C^{1}$ (c.f. examples $(d)$ and $(e)$). Secondly, the assumptions allow for the possibility that the excursion sets could be nested within one another (c.f. examples $(c)$ and $(f)$). Thirdly, the assumptions permit the excursion sets to share common boundaries (see examples $(b)$, $(c)$ and $(f)$). The last two observations are of practical relevance as, in many applications, it may be expected that the target functions for different study conditions exhibit a strong similarity. For instance, in the neuroimaging example, it may be expected that some anatomical regions of the brain display evidence of activation for all of the study conditions. We note here that there is no requirement that $\mathcal{F}_{c}$ be (globally) path-connected; $\mathcal{F}_{c}$ may consist of several disconnected components without violating the assumptions of this section.

In this section, we present the central result of this work, Theorem 2.1, alongside two corollaries, Corollary 2.1 and Corollary 2.2. To do so, we now define the nested sets $\hat{\mathcal{F}}^{-}_{c}$ and $\hat{\mathcal{F}}^{+}_{c}$ by thresholding the statistic $\tau_{n}^{-1}\min_{i\in\mathcal{M}}\hat{g}^{i}_{n}(s)$:

$$\hat{\mathcal{F}}^{-}_{c}:=\hat{\mathcal{F}}^{-}_{c}(a)=\{s\in S:\tau_{n}^{-1}\min_{i\in\mathcal{M}}\hat{g}^{i}_{n}(s)\geq-a\},$$

$$\hat{\mathcal{F}}^{+}_{c}:=\hat{\mathcal{F}}^{+}_{c}(a)=\{s\in S:\tau_{n}^{-1}\min_{i\in\mathcal{M}}\hat{g}^{i}_{n}(s)\geq+a\},$$

using some constant $a\in\mathbb{R}^{+}$. To find an appropriate value of $a$, such that Equation (1) holds for a desired tolerance level (e.g. $\alpha=0.05$), Theorem 2.1 is required.

Conceptually, the event $\{H>a\}$ may be thought of as the situation in which, at some point, $s$, inside some boundary segment, $\partial^{\alpha}\mathcal{F}_{c}$, a large value was observed for the statistic $|\min_{i\in\alpha}(G^{i}(s))|$. Therefore, Theorem 2.1 tells us the probability that the statement $\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}$ is violated is directly determined by such points. If we can find a value of $a$ such that $\mathbb{P}[H\leq a]=1-\alpha$ then asymptotically, the inclusion probability, $\mathbb{P}[\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}]$, will also equal $1-\alpha$. In other words, if the $(1-\alpha)^{th}$ quantile of the distribution of $H$ is known, then confidence regions $\hat{\mathcal{F}}^{-}_{c}$ and $\hat{\mathcal{F}}^{+}_{c}$ may be constructed which satisfy Equation (1). We shall take up the task of evaluating the quantiles of $H$ in Section 3.

Theorem 2.1 is key to generating confidence regions for conjunction inference (intersections of excursion sets). However, Theorem 2.1 may also be extended to allow investigation of other logical statements involving negations or disjunctions. Below we provide two corollaries, each of which is illustrated via a worked example. Corollary 2.1 extends Theorem 2.1 to account for logical negations, whilst Corollary 2.2 provides an equivalent statement for inference performed upon logical disjunctions (unions of excursion sets).

For each $i\in\mathcal{M}$ (i.e. for each study condition), we assume that the data follow a spatially-varying Linear Model (LM) defined at location $s\in S$ as:

$$Y^{i}(s)=X^{i}\beta^{i}(s)+\epsilon^{i}(s),\quad\epsilon^{i}(s)\sim N(0,\Sigma^{i}(s)).$$

(3)

The known quantities in the model are; the $(n\times 1)$ vector of responses, $Y^{i}(s)$, and the $(n\times p)$ design matrix, $X^{i}$. The unknown model parameters are; the $(p\times 1)$ vector of regression coefficients, $\beta^{i}(s)$, and the $(n\times n)$ random error covariance matrix, $\Sigma^{i}(s)$. For each spatial location, $s\in S$, the Generalized Least Squares (GLS) estimator of the parameter vector $\beta^{i}(s)$, $\hat{\beta}^{i}(s)$, is given by:

$$\hat{\beta}^{i}(s)=(X^{i^{\prime}}\Sigma^{i}(s)^{-1}X^{i})^{-1}X^{i^{\prime}}\Sigma^{i}(s)^{-1}Y^{i}(s).$$

In this context, under the $i^{th}$ study condition, at each spatial location $s$, our interest lies in assessing whether linear relationships hold between the elements of the parameter vector $\beta^{i}(s)$. Such relationships may be expressed using a contrast vector $L^{i}$ of dimension $(p\times 1)$. In the notation of the previous sections, the target and estimator functions for the spatially-varying LM are given as $\mu^{i}(s)=L^{i^{\prime}}\beta^{i}(s)$ and $\hat{\mu}^{i}_{n}(s)=L^{i^{\prime}}\hat{\beta}^{i}(s)$ respectively, and $\tau_{n}^{-1}\sigma^{i}(s)$ is the contrast standard error, given as:

$$\tau_{n}^{-1}\sigma^{i}(s)=\big{[}L^{i^{\prime}}(X^{i^{\prime}}\Sigma^{i}(s)^{-1}X^{i})^{-1}L^{i}\big{]}^{\frac{1}{2}}.$$

Whilst in the above notation we have presented $\{\hat{\beta}^{i}(s)\}_{i\in\mathcal{M}}$ as being estimated separately using $M$ different models, we note that nothing in our theory prevents the same model being used across all $M$ study conditions. For example, it is possible for the model matrices, $\{Y^{i}(s)\}_{i\in\mathcal{M}}$ and $\{X^{i}\}_{i\in\mathcal{M}}$, and parameter vector, $\{\beta^{i}(s)\}_{i\in\mathcal{M}}$, to be equal for all $i\in\mathcal{M}$, with the only distinction between study conditions being represented by the contrast vector $L^{i}$. This observation is noteworthy as, in many applications, it is common for researchers to compile all study conditions into a single LM to obtain a heightened statistical power.

To apply Theorem 2.1 to the above LM, we now assume that Assumptions 2.2.1-2.2.3 hold. Such assumptions are commonly satisfied by standard conditions placed on the continuity and boundedness of the increments and moments of $\{\beta^{i}(s)\}_{s\in S}$ and $\{\epsilon^{i}(s)\}_{s\in S}$ (see Sommerfeld et al. (2018) for further detail). We note here that the covariance matrix, $\Sigma^{i}(s)$, is usually unknown in practice, meaning the function $\tau_{n}^{-1}\sigma^{i}(s)$ cannot be calculated. As demonstrated in Sommerfeld et al. (2018), however, $\Sigma^{i}(s)$ can be replaced by any consistent estimator $\hat{\Sigma}^{i}(s)$ and the assumptions given in Section 2.2 are still satisfied. Such estimation is common in the statistics literature and is typically achieved by assuming that some fixed correlation structure applies to $\Sigma^{i}(s)$ (e.g. diagonal independence, auto-regressive, etc.) so that the number of independent variance parameters which must be estimated is small.

In practice, the distribution of $H$ is unknown and must be estimated. In this section, for the model specification described in Section 3.1, we employ a wild $t-$bootstrap resampling scheme to obtain the quantile, $a$, of $H$ which satisfies $\mathbb{P}[H\leq a]=1-\alpha$.

To do so, we first define the $(n\times 1)$-dimensional residual vector, $R^{i}$, for the LM of the $i^{th}$ study condition ($i\in\mathcal{M}$) as:

$$[R^{i}_{1}(s),...,R^{i}_{n}(s)]^{\prime}=R^{i}(s)=[\hat{\Sigma}^{i}(s)]^{-\frac{1}{2}}(Y^{i}(s)-X^{i}\hat{\beta}^{i}(s))$$

The wild $t-$bootstrap proceeds by, in each bootstrap instance, generating $n$ i.i.d. Rademacher random variables, $\{r_{1},...,r_{n}\}$, (random variables which take the values $-1$ and $+1$ with equal probability) independently of the data. For the $i^{th}$ study condition, at spatial location $s$, a new bootstrap sample is then obtained by multiplying the elements of the residual vector by the Rademacher variables (i.e. the new sample is given by $\{r_{1}R^{i}_{1}(s),...,r_{n}R^{i}_{n}(s)\}$). Once the boostrap sample has been generated, the standard deviation of the bootstrap sample, $\hat{\sigma}^{*,i}(s)$, is calculated. A bootstrap instance of the variable $G^{i}$, $\tilde{G}^{i}$, is now generated as follows;

$$\tilde{G}^{i}(s)=n^{-\frac{1}{2}}\sum_{l=1}^{n}\frac{r_{l}R^{i}_{l}(s)}{\hat{\sigma}^{*,i}(s)}.$$

A bootstrap instance of the variable $H$, $\tilde{H}$, may then be computed using the bootstrap variables $\{\tilde{G}^{i}(s)\}_{s\in S,i\in\mathcal{M}}$ as follows:

$$\tilde{H}=\max_{\alpha\in\mathcal{P}^{+}(\mathcal{M})}\bigg{(}\sup_{s\in\partial^{\alpha}\hat{\mathcal{F}}_{c}}\big{|}\min_{i\in\alpha}(\tilde{G}^{i}(s))\big{|}\bigg{)}$$

(4)

Quantiles of the distribution of $H$ may now be estimated empirically from the observed sampling distribution of $\tilde{H}$. Discussion of how the above supremum, taken across continuous space, is evaluated in practice is deferred to Supplementary Results Section S$2.1$. However, we do briefly note here that, as the location of the true boundary segment $\partial^{\alpha}\mathcal{F}_{c}$ is in practice unknown, the boundary segment $\partial^{\alpha}\mathcal{F}_{c}$ has been replaced by an estimate $\partial^{\alpha}\hat{\mathcal{F}}_{c}$. It should be noted that, as a result of this substitution, our proposed method may not be applied when $\hat{\mathcal{F}}_{c}$ is empty.

Several strong references exist that argue and verify the correctness of using the wild $t-$bootstrap for estimating quantiles of maxima distributions in the above manner. Of particular note are the works of Chang et al. (2017) and Bowring et al. (2019), in which the wild $t-$bootstrap is proposed for the estimation of Gaussian maxima distributions and the generation of confidence regions, respectively. Discussion of the wild $t-$bootstrap in an imaging context may also be found in Telschow and Schwartzman (2021). Other significant references which serve as precursors to this work include Chernozhukov et al. (2013) and Sommerfeld et al. (2018), in which the Gaussian multiplier bootstrap was investigated for estimating Gaussian maxima distributions and generating confidence regions, respectively. More general introductions to such bootstrapping procedures may be found in Wu (1986), Efron and Tibshirani (1994) and Hesterberg (2014).

In the form that it has been presented here, the wild $t-$bootstrap requires that $\{G^{i}\}_{i\in\mathcal{M}}$ be symmetric. Further, as noted above, the wild $t-$bootstrap has been extensively verified for estimating maxima distributions only in settings in which $\{G^{i}\}_{i\in\mathcal{M}}$ are Gaussian. Here we stress that, whilst the theory presented in Section 2 makes no assumptions on the distribution of $\{G^{i}\}_{i\in\mathcal{M}}$, the bootstrap theory presented throughout Section 3 requires $\{G^{i}\}_{i\in\mathcal{M}}$ to be Gaussian. In order to utilise the results of Section 2 when the random variables $\{G^{i}\}_{i\in\mathcal{M}}$ are not Gaussian, alternative methods must be employed for estimating the quantiles of $H$.

To assess the correctness and performance of the method, a range of simulations were conducted using synthetic data. Each simulation was designed to investigate how the empirical coverage (i.e. the proportion of simulation instances in which the inclusion $\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}$ held) was affected by various secondary factors of practical interest. In this section, we describe three simulations designed to investigate how the methods’ performance was influenced by; $(1)$ the degree of overlap between excursion sets, $(2)$ the presence of correlation between the noise fields for different study conditions, and $(3)$ the magnitude of the spatial rate of change of the signal.

In all simulations, for $i\in\mathcal{M}$, the model used to generate the synthetic data took the form;

$$Y^{i}(s)=\mu^{i}(s)+\epsilon^{i}(s),\quad\epsilon^{i}(s)\sim N(0,\sigma^{i}(s)I_{n}).$$

(5)

where $n$ was allowed to vary from $40$ to $500$ (in increments of $20$). The aim, across all simulations, was to assess the performance of the method when employed for conjunction inference on the true signal $\{\mu^{i}\}_{i\in\mathcal{M}}$ (i.e. when asked the question “Where do all of the $\{\mu^{i}\}_{i\in\mathcal{M}}$ exceed a predefined threshold?”).

The mechanism used to generate the true signal, $\{\mu^{i}\}_{i\in\mathcal{M}}$, varied across simulations. In Simulations 1 and 2, $\{\mu^{i}\}_{i\in\mathcal{M}}$ were generated using two binary images of circles and squares, respectively, positioned in a ‘Venn diagram’ arrangement. Each binary image was scaled by a predefined amount and then smoothed using an isotropic Gaussian filter with a Full-Width Half Maximum (FWHM) of $5$ pixels. In Simulation $3$, $\mu^{1}$ and $\mu^{2}$ were simulated as a horizontal and vertical linear ramp, respectively, with predefined gradients.

Each simulation was performed twice, once using noisier ‘low-Signal-to-Noise Ratio (SNR)’ synthetic data and again using less noisy ‘high-SNR’ synthetic data. To generate the high-SNR synthetic data for Simulations $1$ and $2$, all simulated $\{\mu^{i}\}_{i\in\mathcal{M}}$ were scaled by a factor of $3$ prior to smoothing. This data generation process is identical to that employed in Bowring et al. (2019), in which it is noted that synthetic data generated using these parameters strongly resembles that of an fMRI analysis when voxels are of dimension $2$mm³. In Simulation $3$, for the high-SNR data, the linear ramps were initially simulated with a gradient of $8$ per $50$ pixels. Across all simulations, the low-SNR data was generated by reducing the signal magnitude of the high-SNR data by a factor of $4$. In all simulations, thresholds of $c=1/2$ and $c=2$ were employed for the low-SNR data and high-SNR data, respectively.

As Simulation $1$ aimed to investigate the method’s performance as the overlap between excursion sets increased, in this simulation, the distance between the centre of the circles was varied, ranging from $0$ pixels (full overlap) to $50$ pixels (barely overlapping) in increments of $2$ pixels. Similarly, as Simulation $2$ aimed to assess the effect of correlation between $\epsilon^{1}$ and $\epsilon^{2}$, in this simulation, this correlation was varied between $-1$ and $1$ in increments of $0.1$ (in all other simulations, the noise fields were generated independently). As Simulation $3$ served to assess how sensitive the empirical coverage was to the rate of change of $\mu^{1}$ and $\mu^{2}$ over space, in this simulation, the initial ramp gradients were multiplied by a factor of $k$, where $k$ was varied from $0.25$ to $1.75$ in $0.05$ increments.

All simulation results were obtained as averages taken across $2500$ simulation instances and compared to the nominal coverage using $95\%$ binomial confidence intervals. In all simulation instances, the image dimensions were $(100\times 100)$ pixels, confidence regions were computed for a tolerance level of $\alpha=0.05$ using $5000$ bootstrap realizations, $\sigma^{i}(s)$ was computed using the ordinary least squares estimator and $\tau_{n}$ was computed as $\tau_{n}=n^{-0.5}$. In each simulation instance, the assessment of whether the inclusion statement, $\hat{\mathcal{F}}^{+}_{c}\subseteq\mathcal{F}_{c}\subseteq\hat{\mathcal{F}}^{-}_{c}$, held was verified using the interpolation-based methods of Bowring et al. (2019). In this approach, the inclusion statement was deemed to hold if, and only if, the sets of pixels which were identified as $\hat{\mathcal{F}}_{c}^{-}$, $\mathcal{F}_{c}$ and $\hat{\mathcal{F}}_{c}^{+}$ were appropriately nested and, in addition, interpolation along the boundary $\partial\mathcal{F}_{c}$ confirmed that no violations of the inclusion statement had occurred.

A further six simulations which investigated the influence of other secondary factors of interest were also conducted. The secondary factors considered by these simulations include; the presence of spatial structure in the noise, the effect of $\partial^{\{1\}}\mathcal{F}_{c}$ and $\partial^{\{2\}}\mathcal{F}_{c}$ having differing lengths, the effect of $\epsilon^{1}$ having a much larger variance than $\epsilon^{2}$ and the impact of varying the value of $M$. A full discussion of these simulations is deferred to Supplementary Results Section S$4$. In addition, in keeping with the previous literature, tolerance levels of $\alpha=0.1$ and $\alpha=0.2$, as well as the substitution of $\partial\mathcal{F}_{c}$ for $\partial\hat{\mathcal{F}}_{c}$ in Equation (4), were also considered for simulation (c.f. Sommerfeld et al. (2018), Bowring et al. (2019), Bowring et al. (2021)). Again, the results of these simulations may be found in the Supplementary Results document in Sections S$5$ and S$6$.

For a majority of the simulations conducted, the empirical coverage estimates were tightly distributed around the expected nominal coverage. In particular, for the high-SNR data, across all simulations, the $95\%$ binomial confidence intervals consistently captured the nominal coverage at the predicted rate. However, when the simulations employed the low-SNR data, it can be seen that the confidence regions produced by the method were conservative when the data were generated under certain unfavourable conditions.

For example, in Fig. 3, it can be seen that some over-coverage was observed for Simulation $1$ when the data was noisy, and the number of subjects was low. This observation is corroborated by the results of the remaining low-SNR synthetic data simulations (c.f. Supplementary Results Section S$5$). However, in all cases, the degree of agreement between the empirical and nominal coverage improved as the sample size increased. Furthermore, no such over-coverage was observed for the equivalent high-SNR simulations. This observation matches expectation, as Theorem 2.1 holds only asymptotically, and is supported by the previous literature on confidence regions, in which similar results were observed for the single-‘study condition’ setting (c.f. Sommerfeld et al. (2018), Bowring et al. (2019)).

The presence of strong negative correlation between the study conditions and the rate at which the signal varied across space both also appeared to influence the methods’ performance (c.f. Fig 3, Simulations 2 and 3). In both instances, it is likely that the observed over-coverage is due to difficulties in estimating the true location of the boundary $\partial\mathcal{F}_{c}$, as each of these factors make it harder to identify the locations at which small changes in $\min_{i\in\mathcal{M}}\mu^{i}$ occur. When high-SNR data was used in the place of low-SNR data, both of these factors had a substantially reduced impact on the empirical coverage.

Despite the above observations, the simulation results overwhelmingly supported the claim that the proposed method is robust to a range of secondary factors. Included in the list of such factors are; variation in the shape and size of $\mathcal{F}_{c}$, spatial correlation in the noise, variation in the lengths of individual boundary segments, between-‘study condition’ correlation, realistic levels of variation in the magnitude of the noise, large numbers of study conditions and situations in which boundary segments are shared by many excursion sets. For further detail, see Supplementary Results Sections S$5$ and S$6$.

In terms of time efficiency, the wild $t-$bootstrap provided extremely fast computational performance. For instance, using an Intel(R) Xeon(R) Gold 6126 2.60GHz processor with 16GB RAM, and averaging over the $2500$ simulation instances conducted for $n=500$ high-SNR observations, in Simulation $1$, the time taken to generate confidence regions for a circle separation of $20$ pixels was $5.98$ seconds. For a comprehensive summary of computation times, see Supplementary Results Section S$7$.