Theoretical framework

The problem of testing statistical hypotheses at each element of a 2 or 3-dimensional field can be conceived as the following general problem:

Given a set L of sites, there is a statistic T(u) defined for each site u∈L for which one wishes to test a null hypothesis (H₀). According to this hypothesis, all elements T(u) are furnished by the distribution P0(T) (the null distribution). H₀ is assumed to be of the form:

where H_0u is a marginal null hypothesis about the probability distribution of the measurements at site u. At H_0u, consequently T(u) is generated by P₀. In the active region, A is defined as the set of sites u where H_0u does not hold, thus affirming the alternative hypothesis H_1u . The problem then is to find an estimate A^ for the set A (note: A and A^ may be empty).

Once a method is selected, its performance must be evaluated with standard tools. Some very common tools that will be used presently are described.

where Ac is the complement of the active region A, l•l denotes the cardinality of a set and E[•] refers to the expected value of a random variable.

It is also possible to define FPR for elements that are not adjacent to the boundary of the active region:

where the r-dilation Dr of A is defined as:

This measure is defined for two reasons. Firstly, the boundary of the active region is usually not well localized, since the activation level often decreases slowly as one moves away from A. Secondly, the methods that consider the neighborhood of each element for estimating A likely show an increased number of false positives close to the boundary of the active region. If most false positives are of this kind, they should have less impact on the measured performance than the false positives disconnected from A.

Thus, FPRr could be a better performance measure in such a case.

denoting the expected proportion of sites correctly estimated as the activation region.

representing the probability Pr of having at least one false positive, given that the activation region is empty.

with W=(|L|-|A|)/|A|,

portraying the proportion of wrongly rejected null hypotheses in A^.

Some of these measures can be combined. For example, a widely accepted way of characterizing the performance of a method is through the receiver operating characteristic (ROC) curve^[13], which indicates, for a fixed A, the maximum attainable TPR for any given maximum allowable FPR. For one-sided tests, the maximum TPR is generally an increasing function of the maximum allowable FPR.

To construct the curve, the true region A and the corresponding activation level a=ETu, u∈A must be known. Then, the study of a method that depends on a parameter θ involves setting a value for θ to obtain a point on the ROC curve. The parameter acts as a cut-off point to distinguish between the sites considered positives or negatives. Take as an example the methods based on the following computation:

where pv(u) depicts the p-value of Tu and therefore pvu=1- P0(Tu). This class of methods, denominated PW, encompasses most of the standard methods. They differ from each other only in the way the threshold θ is computed (note: all PW methods have the same ROC curve, regardless of the way θ is computed). Thus, for the application of the threshold based on a significance level α, without correction for multiple hypotheses (uncorrected method), it holds that θ = α. For the FWER method ^[6,⁸,^14], with significance level α_FWER, it follows that:

where P0M is the distribution of maxu∈L Tu under H₀. Consequently, FWER can be controlled by choosing as the threshold the value located at the (1-α_FWER) portion of the right side of P0M. For an elaborate discussion on the association of FWER with the maximum statistical value, see Pantazis ^[15].

In the case of FDR, with a significance level α_FDR^[16], the procedure for finding the threshold θ begins with ordering the individual p-values of sites u∈L from the largest to the smallest. Accordingly, pvu1≥ pvu2 ≥… pvuN with N=|L|. Let k be the index of the first site on the list, at which the p-value is less than or equal to the desired FDR proportion. Hence,

and θ is set as:

For more details on the method, consult Benjamini ^[17]. Alternatively, given a desired value ϵ for FPR_max, θ can simply be set as ϵ (e.g., ϵ= 10^-5). If the local hypotheses H_0u are independent, this is equivalent to the application of the standard PW method without correction for multiple hypotheses, but with low level of significance α = ϵ. Since (as mentioned) the maximum TPR is an increasing function of FPR, the value of θ* = ϵ will be the one that maximizes the TPR while keeping FPR under control:

where ϵ is a user-specified small positive number. Although the actual ROC curve for a given problem is unknown (because it depends on the values of Tu for u∈A), it is possible to specify the value of ϵ. Thus, the optimal PW method, according to Eq. (11) with θ = θ*, will correspond to the UC method with α = 1- ϵ. It may be applied if the field of p-values, or equivalently the field P0Tu, u ∈L, is known. Such fields can be estimated either theoretically or by non-parametric empirical means (e.g., permutation or re-sampling procedures)^[18,^19]. The optimal PW method will be denoted by PW* (ϵ).

The signal-to-noise ratio (SNR) ^[20] is defined as:

where σ₀ is the variance of T under H₀ and the minimum activation level in T has been arbitrarily assigned as the information signal. This represents the worst-case scenario in which the entire region A has a minimum value. On the other hand, the mean activation level or the amplitude could also be employed, as discussed by Welvaert and Rosseel ^[17] and Acosta-Franco et al.^[5] For low values of SNR, the value of TPR obtained with PW* (ϵ) in the corresponding PW ROC curve is usually low for reasonable values of ϵ.

Consequently, PW methods have low sensitivity and the estimated A^ is too conservative.

Hence, more sensitive methods must be developed that are able to accurately reflect the active regions, which in most cases consist of clusters of several contiguous elements of L and not of isolated elements. Customarily, these methods use the field T as a starting point for the generation of a new field T^. The new statistic takes the spatial correlation of A into account, and then processes T^ with a PW method as before. For example, T^ may refer to the size or mass of a cluster of arc-connected elements (u) with values of T above a certain threshold ^[21].

However, this approach has some drawbacks, one being that the results depend strongly on the value of the selected threshold, and in general no principled way exists to make such a selection. Some variants of the method alleviate the problem to a certain extent by computing T^ as a weighted combination of cluster sizes obtained with different thresholds. Another problem is the loss of interpretability, since in many cases the significance of the statistic sought is directly related to the activation level (i.e., the value of T), and to the extension of supra-threshold clusters.

A distinct approach was developed to address these problems, being the morphology-based hypothesis testing (MBHT) method^[12]. It involves the computation of T^ as a combination of the results found when applying a set of K morphological erosions with different sizes of structuring elements to the field T. The statistics are calculated as:

where W_k(u) denotes the structuring element k (e.g., a circle of radius r_k) centered at u. Then, these values are integrated into the statistic T^:

where P_0K is the null distribution of the statistic Tk. Afterwards, the optimal PW procedure can be applied to the field T^ as detailed above.

The results of this approach are competitive with those based on supra-threshold cluster statistics^[12] and allow a clearer interpretation. Once again, the disadvantage is that the results may depend on the shape of the structuring elements, and for their selection no principled solution exists.

In the current contribution, the approach introduced is capable of overcoming these difficulties. It formulates the hypothesis testing task as a Bayesian estimation problem, with an MRF previously applied to the active regions, thus implementing a prior constraint on the spatial contiguity of A. For the application of a Gaussian Markov random field (GMRF) for fMRI data analysis, see Mejia et al. and the references therein^[3].

The new scheme proved to have better performance than PW methods, while maintaining interpretability. It also has better performance than MBHT, and does not require the selection of any particular shape for the structuring elements.

The regularized hypothesis testing method

In RHT, the hypothesis testing problem is written in terms of an image segmentation problem, which is solved by using a Bayesian estimation framework. First, the hypothesis testing formulation is explained. The procedure is laid out for calculating the prior distributions of sites and the likelihood that they belong to an active region. Subsequently, an approximation algorithm for finding the maximum a posterior probability (MAP) estimate is established. Finally, the parameter selection problem is addressed.

Hypothesis testing formulation

Hypothesis testing may be formulated as a binary segmentation problem. Accordingly, given the set of sites L, the problem is reduced to partitioning L into non-overlapping cohesive regions A0, A1, such that the activation level in region A0 is zero, and in region A1 is greater than or equal to a known constant α₁. Hence, A0 ∩ A1= ∅ and L= A0 ∪ A1, where A=defA1 is the active region, and therefore A0= Ac.

Let c be an unknown discrete label field that identifies the partitions of L, defining c(u)∈K=def{0,1} and cu=k if u∈Ak, for each u∈L. Thus, ac(u) depicts the activation level at site u (note: α₀=0).

Without loss of generality, the activation level in A1 is assumed to be equal to α₁, since this represents a worst case in terms of the desired false positive control in the estimation of A. Considering T as the field formed by the statistics T(u) for all u∈L, the following observation model is proposed:

where n is a noise field showing distribution P_n(n), with:

where Z_n is a normalization constant and U_n(n) is a so-called energy function. Then, the likelihood P(T|c) is furnished as:

where ZT|c is a normalization constant. It is important to find an estimate that considers the spatial correlation in both the noise field n and the label field c .

As the real statistic T is obtained through a regression analysis of data provided by a model, it could take on a positive or negative value (depending on the choice of the statistic), resulting in positive or negative differences− activation regions− of the data with respect to the model. Zero-value regions would still portray sites where there is no activation, and the distribution of values of field n would be determined by the selection of the statistic T. Hereafter, a noise model is presented as a GMRF. In later sections, we remove this assumption by introducing a procedure for standardizing the noise field in order to map noise fields ranging from arbitrary distributions to standard normal distributions.

Noise field model

For the noise field n, a GMRF model was employed, taking the form of Eq. (16) with:

where γ > 0 is a scale parameter, n(u) is the value of the field n in site u, and n(v) is the value of the field n in a neighboring site v. The parameter τ₁ ≥ 0 is related to the spatial correlation of the field n and the second term added to the first is taken from all pairs of neighboring sites ⟨u,v⟩ in the image.

When τ₁ = 0 and γ = 1, the formulation is reduced to a Gaussian white noise model with zero mean and unit variance (the “Parameter Selection” section explains other types of noise fields). With such a model, the likelihood takes the form of Eq. (17), being:

In the above equation, c(u) is used as an indicator function. Accordingly, if its value is 1 (i.e., u ∈ A1), the term (1-c(u))T(u)2 becomes zero, simplifying Eq. (19) to an expression that is equivalent to substituting n (u) for T(u)-a1c(u). In Eq. (18), the activation level in T(u) is eliminated, keeping only the noise component, as in the observation model proposed in (15).

Moreover, if c(u)=0 (i.e., u∈A0), the term c(u)(T(u)-a1c(u))2 becomes zero and the term T(u)-a1c(u) in the second summation is simplified to T(u). This is equivalent to substituting n(u) for T(u) in Eq. (18), since T(u) already corresponds to purely noise.

Eq. (19) may be rewritten as a quadratic function of a vector-valued field b defined as b₁(u) = c(u) and b₀(u) = 1 - c(u),

with a0=0 and:

Bayesian formulation of the hypothesis testing problem

For field b, we propose the previous application of an MRF model ^[22,^23]. In order to introduce the prior constraint that the active regions are spatially cohesive, b is modeled with Ising potentials, furnishing the following distribution:

where τ₂ ≥ 0 is a parameter that controls the granularity of field b.

With this equation and the Gauss-Markov model (16) and (18) for noise, it is possible to estimate field b by employing a Bayesian formulation, affording the posterior probability distribution:

where the likelihood is:

with UT|b given by (20) and field b satisfying (21)-(22). Finally, the posterior distribution can be written as:

where UT|bb and Ubb are furnished by (20) and (24) respectively.

The MAP estimate for field b is obtained by minimizing (28), subject to the constraints (21)-(22). This is a combinatorial optimization problem that is, in general, very difficult to solve. Several methods have been proposed, such as the iterated conditional modes (ICM) algorithm^[24], stochastic relaxation^[25], graph cuts procedure^[26], etc.

Here, because of the form of the cost function that includes the noise correlation term, we prefer option^[22], which is based on the relaxation of constraint (22). Hence, a new formulation is made:

Dividing (28) by γ and utilizing ν=defτ1/γ and λ=defτ2/γ, the functional is minimalized as:

If the noise has a mean of zero, then α₀ = 0. The resulting optimization problem is:

This is a quadratic minimization problem, subject to linear constraints, which can be solved efficiently by employing a projected gradient descent method. Consequently, each b(u) (in order to simplify the notation we remove the dependency of ν, α₁, λ) may now be interpreted as an approximation for the posterior marginal distribution of the indicator functions  1Ak(u).

To obtain the latter it is possible to set 1Ak(u)=1 if b_k(u) > b_j(u) for all j ≠ k. Thus, the estimation of the active region is:

Parameter selection

A general method can now be introduced to allow the model (31)-(33) to be used for any type of noise. The parameters α₁, ν, λ will be able to be automatically selected by controlling the false positive proportion while maximizing the performance of the method with respect to the TPR.

Standardization of the noise field

The derivation of functional (31) is based on the observation model (15) with the Gauss-Markov model (16)-(18) for noise. However, the model is not always valid in practice, since the T field may have a different distribution. The first step, therefore, is to generate a transformation capable of mapping the arbitrary distribution to a standard one. The proposed method is based on two properties:

Combining both properties, it follows that for any random variable x, the random variable G ^-1 (Q(x)) is normal. Consequently, the only requirement for transforming a random variable x into a normal variable is to know the corresponding CDF Q(x) or to estimate it by obtaining the empirical CDF (ECDF) Q^(x) from samples of the particular distribution.

Assuming that random field samples T0 can be generated under H₀, it is possible to estimate the corresponding null distribution P₀ (⋅). Given a field T derived from the observations (e.g., the F-test for the generalized linear model (GLM) computed at each voxel), the sample can be mapped to a standard normal distribution Φ:

where erf(y) denotes the error function^[30], the probability that a random variable Z normally distributed with mean μ = 0 and variance σ = 1/2 falls in the range [-y,y]. The field  T^, with a standard normal distribution, is afforded by the following transformation:

Hyper-parameter calibration

The hyper-parameters that control the behavior of the method are the spatial correlation parameter ν, the minimum activation level α₁, and the regularization parameter λ. The general idea is to estimate ν from the observations, and subsequently compute α₁, λ from the estimated parameter ν and a parameter ϵ representing the level of FPR control.

Estimation of parameter ν: The parameter ν = τ ₁/γ in (31) controls the noise spatial correlation and can be estimated by minimizing the negative logarithm of the pseudo-likelihood of samples of the field T0 generated under H₀ (i.e., samples of the noise field n).

The pseudo-likelihood is defined as the product of the conditional distributions for n(u) at each pixel u, given the values at its neighboring fields ^[31]. Using Eqs. (17)-(18), the result is:

where t is a constant, N is the neighborhood of pixel u and:

The minimization of Eq. (37) yields a closed formula for the estimation of the parameters:

Calibration of parameters α₁, λ: For a given value ϵ for FPR control and a fixed ν= ν^ calculated with Eq. (44), the present proposal is to find the optimal parameter settings α ₁*, λ* as the maximizers of the TPR, subject to FPR ≤ ϵ. However, the solution of such an optimization problem is very difficult to obtain because both TPR and FPR depend on the unknown region A. To find an approximate solution, the following observations are made:

where P(α) is taken as the uniform distribution in the set S. With these simplifications, the parameter estimation algorithm simply consists of sweeping the plausible λ values in a given set Λ (a discretization of an interval [0, λ_max]), followed by finding, for each λ, the value a^1(λ) (such that FPR(a^1(λ),λ) =ϵ) and the corresponding TPR-(λ), as explained above. The final λ* is then found as the maximizer of TPR-(λ), affording a1*= a^1(λ*). This is summarized in Algorithm 1.

The optimal parameters α₁*, λ* depend on the data only through the estimated parameter ν and the desired FPR⁰ control ϵ. Thus, the optimal values may be pre-computed and stored in a table. When a new data set arrives, one only needs to standardize the noise, estimate ν, read α₁*, λ* from the table for the desired FPR control, and minimize U (b; ν, α₁*, λ*) furnished by Eq. (30). The tables are represented as false color images in Figure 1.

The final algorithm for estimating the active region A^ for a given data set is summarized in Algorithm 2.

Data description

To study the behavior of the present method, we include various experiments based on both synthetic and real data.

Synthetic data: The synthetic data was generated dynamically and consisted of two components. Firstly, S₀ = {n₁, n₂, …, n₄₀} is a set of noise fields of a regular lattice L of 50×50 pixels, without any activation region. That is, the images were produced under H₀, using Eqs. (16) and (18) with n(u) as random independent variables showing standard Gaussian distribution. Parameter γ = 1 and ν were to be specified in each experiment. Secondly, S₁ = {c₁, c₂, …, c₄₀} consists of 40 binary labeled fields of a regular lattice of 50×50 pixels, representing indicator functions for activation regions of different shapes (see Figure 2). Then, the simulated observed fields are generated by the following model:

where α denotes the activation level and its value is specified in each experiment.

Real data: Experiments with real data were based on the auditory dataset^[32]. The data corresponded to 96 volumes of a single subject (each volume composed of 64×64×64 voxels of 1×1×3 mm), which were acquired in blocks of 6 volumes. Since the repetition time between scanning was set to 7 seconds, there were a total of 16 blocks of 42 seconds (although due to the effects related to T1, the first two blocks were discarded). The sequence of volumes alternated between blocks of rest and stimulation, starting with rest. Auditory stimulation consisted of two-syllable words presented binaurally at a rate of 60 per minute.

Localization of false positives for the RHT algorithm

This experiment was designed to test the assumption that in the proposed method the false positive errors are concentrated close to the boundary of the active region. Consequently, FPR₂ can be well approximated by FPR⁰ (i.e., FPR computed over images produced under H₀).

A total of 1000 independent runs were conducted for the procedure. Briefly, set S₀ was generated with parameter ν = 0.75, a value estimated from fMRI images by utilizing Eq. (44). Fields T₁, T₂, …, T₄₀ were obtained by using model (47), sets S₀ and S₁, and randomly selected activation level α in the interval ^[2,^4].

Algorithm 2 was applied for levels of FPR control ϵ ∈ {0.001, 0.0001, 0.00001, 0.000001} with the optimally calibrated parameters, computing the average value of FPR₂ for the set S₁and FPR for the set S₀. The results are shown in Table 1. As can be appreciated, by calibrating the parameters of the method to control FPR⁰, the appropriate control over FPR₂ is also achieved.

Performance of the method

The first set of the following experiments was designed to make a comparison of RHT to FWER ^[8] and MBHT ^[12], the latter two being the state-of-the-art methods that take the spatial expanse of the active region into account. Accordingly, the set S ₀ was produced with parameter ν ∈ {0, 0.5, 1, 1.5}. A total of 1000 runs were performed, randomly selecting the active regions A from S₁ with an activation level a ∈ {0.0, 0.25, 0.5, ⋯,5}. The average TPR was computed. The RHT method was carried out by employing Algorithm 2 and a level of false positive control of ϵ^=0.000001.

Based on the results from each method, RHT proved capable of providing better TPR values than the other two algorithms for most activation levels (Figure 3). The improvement is likely due to the consideration in the new model of the spatial cohesion of the activation regions (controlled by the regularization parameter λ) as well as the spatial correlation in the noise field (controlled by the parameter ν). Contrarily, the rest of the algorithms assume that these two components of statistic T are formed by independent variables.

Experiments with fMRI data

In the second set of experiments, Algorithm 2 of the proposed method was applied to the real fMRI data described at the beginning of this section. RHT was not applied directly to the original data, but rather to a field T (an F-test field computed from the data). The first step was to standardize the field in order to obtain T^ with Eq. (36), which is the input of the segmentation algorithm. Hence, it was necessary to calculate the ECDF from the data. Hereafter, some details about the aforementioned steps are explained.

Data pre-processing. The aim of the pre-processing was to remove artifacts in the data as well as to prepare the data to maximize the statistical analysis. Here we use spatial pre-processing provided in the script auditory_spm12_batch.m, which implies realignment, co-registration, segmentation and normalization. Although the original script (available online at http://www.fil.ion.ucl.ac.uk/spm/data/auditory/) includes a smoothing step to ensure that some assumptions about noise distribution are fulfilled^[11,^18], RHT omits this step because it is capable of handling different noise distributions. Actually, the inclusion of the step would not be beneficial. The data for these calculations was processed on SPM software version 12 (available at http://www.fil.ion.ucl.ac.uk/spm/) and MATLAB 2018a. A slice of the pre-processed data was selected to perform the following experiments.

F-test field. After the pre-processing stage, the statistical analysis was carried out to determine the active voxels that correspond to a given stimulus. To obtain the active regions, a voxel-wise analysis is typically conducted by fitting models to a single voxel time course. The data at each voxel is modeled in a univariate way with a linear model:

where y_t, the dependent variable, is a vector formed by the intensity values at each time point. β₀ and β₁ are the intercept and the slope of the linear model, respectively. Meanwhile, ϵ_t is the error term in the model fitting, which captures factors other than x_t, such as signal noise, capable of affecting y_t. The explanatory variable x_t corresponds to the model of neural activity. It is also a vector comprised of entries depicting a real value at each time point:

where h_t stands for the hemodynamic response function and d_t ∈ {0,1} is an impulse train that indicates whether the stimulus was present at time t. Then, the neural activity is modeled as a convolution,* in which the hemodynamic response function acts as a filter.

The F-test represents the ratio between the variance described in a reduced model y_t = β₀ + ϵ_t (without any additional effect) and the full model (48)-(49) that includes the effect of interest. Finally, the F-test was computed for every voxel in the volume. The resulting field was the input for the RHT algorithm, calculated by using the script auditory_spm12_batch.m with the default parameters.

Null distribution. Before the segmentation step of the RHT algorithm, it was necessary to normalize the F-test field and thus to estimate the null distribution. More specifically, a sample of the F-test was obtained under H₀, allowing for the computation of the corresponding ECDF. This was carried out by permuting the order of the stimulus labels d(t) for the volumes (see ^[33] for more details), and by calculating the F-test field for each permutation, as explained in the previous section. After completion of these procedures, it was possible to apply Eq. (36) to transform the original F-test field (i.e., the original order of the labels d(t), without permutation) to one that follows a standard normal distribution. Finally, T^ fields generated under H₀, the value of ν was estimated by utilizing Eq. (44), finding ν^= 0.75.

RHT algorithm robustness with respect to SNR

In order to investigate the stability of the present method, the algorithm was tested by modifying the SNR in the data from which the activation region would be established. Accordingly, blocks of observations were eliminated from the full experiment. As aforementioned, the first two of the 16 original blocks were eliminated due to effects related to T1. Hence, 14 blocks (84 volumes) were taken initially to perform the procedure described earlier in this section to identify the activation regions. Subsequently, the number of blocks was reduced by two, taking 12 blocks (72 volumes) and performing the procedure again to determine the activation regions, and so on until reaching 4 blocks.

In each case, a tolerance for false positives to ϵ = 0.0001 was established and the results of RHT were compared to those obtained with MBHT and FWER using the same number of blocks. It can be verified by visual inspection (Figure 4) that the three methods− RHT, MBHT and FWER− detected activation regions corresponding to the primary auditory cortex, located at the upper sides of the temporal lobes, specifically on the transverse temporal gyri ^[34,^35].

The MBHT and RHT methods afforded better performance and robustness for the reduction of information in most cases. MBHT and FWER displayed less accuracy when reducing the signal information for ϵ = 0.00001. While FWER exhibited difficulties in the identification of activation regions with 6 and 4 blocks (barely finding any region), RHT and MBHT gave better performance. However, when the number of blocks was decreased to only 4, MBHT was not able to detect any activation region, while RHT still had a proficient outcome. RHT outperformed the other models two by successfully recovering both regions.

Degree of false positive control. To avoid having to specify a particular value for ϵ, it is possible to estimate A for a set S of ϵ values and produce a map in which the degree of false positive control (log ϵ for each A^(ϵ) = RHT (T,ϵ^,a1*(ν,ϵ),λ*(ν,ϵ)) is color coded. The result is very similar to the p-value maps that are created for the classic PW method. Specifically, we define the degree of false positive control at site u as:

where 1A^ϵ(u) is the indicator function for the set A^(ϵ). An example of this type of map is shown in Figure 5, with S = {0.01, 0.001, 0.0001, 0.00001}. The map portrays important differences between the first and second levels, but in the last levels the detected regions are more similar to each other. The slight variations exist only at the border of the active region.