Spatial confidence regions for combinations of excursion sets in image analysis.

The analysis of excursion sets in imaging data is essential to a wide range of scientific disciplines such as neuroimaging, climatology, and cosmology. Despite growing literature, there is little published concerning the comparison of processes that have been sampled across the same spatial region but which reflect different study conditions. Given a set of asymptotically Gaussian random fields, each corresponding to a sample acquired for a different study condition, this work aims to provide confidence statements about the intersection, or union, of the excursion sets across all fields. Such spatial regions are of natural interest as they directly correspond to the questions 'Where do all random fields exceed a predetermined threshold?', or 'Where does at least one random field exceed a predetermined threshold?'. To assess the degree of spatial variability present, our method provides, with a desired confidence, subsets and supersets of spatial regions defined by logical conjunctions (i.e. set intersections) or disjunctions (i.e. set unions), without any assumption on the dependence between the different fields. The method is verified by extensive simulations and demonstrated using task-fMRI data to identify brain regions with activation common to four variants of a working memory task.

Using the atmospheric CO2 growth rate to constrain the CO2 flux from land use and land cover change since 1900.

We explore the ability of the atmospheric CO2 record since 1900 to constrain the source of CO2 from land use and land cover change (hereafter "land use"), taking account of uncertainties in other terms in the global carbon budget. We find that the atmospheric constraint favors land use CO2 flux estimates with lower decadal variability and can identify potentially erroneous features, such as emission peaks around 1960 and after 2000, in some published estimates. Furthermore, we resolve an offset in the global carbon budget that is most plausibly attributed to the land use flux. This correction shifts the mean land use flux since 1900 across 20 published estimates down by 0.35 PgC year-1 to 1.04 ± 0.57 PgC year-1 , which is within the range but at the low end of these estimates. We show that the atmospheric CO2 record can provide insights into the time history of the land use flux that may reduce uncertainty in this term and improve current understanding and projections of the global carbon cycle.

Simultaneous confidence bands for functional data using the Gaussian Kinematic formula.

We propose a construction of simultaneous confidence bands (SCBs) for functional parameters over arbitrary dimensional compact domains using the Gaussian Kinematic formula of t-processes (tGKF). Although the tGKF relies on Gaussianity, we show that a central limit theorem (CLT) for the parameter of interest is enough to obtain asymptotically precise covering even if the observations are non-Gaussian processes. As a proof of concept we study the functional signal-plus-noise model and derive a CLT for an estimator of the Lipshitz-Killing curvatures, the only data-dependent quantities in the tGKF. We further discuss extensions to discrete sampling with additive observation noise using scale space ideas from regression analysis. Our theoretical work is accompanied by a simulation study comparing different methods to construct SCBs for the population mean. We show that the tGKF outperforms state-of-the-art methods with precise covering for small sample sizes, and only a Rademacher multiplier-t bootstrap performs similarly well. A further benefit is that our SCBs are computational fast even for domains of dimension greater than one. Applications of SCBs to diffusion tensor imaging (DTI) fibers (1D) and spatio-temporal temperature data (2D) are discussed.

Confidence Sets for Cohen's d effect size images.

Current statistical inference methods for task-fMRI suffer from two fundamental limitations. First, the focus is solely on detection of non-zero signal or signal change, a problem that is exacerbated for large scale studies (e.g. UK Biobank, N=40,000+) where the 'null hypothesis fallacy' causes even trivial effects to be determined as significant. Second, for any sample size, widely used cluster inference methods only indicate regions where a null hypothesis can be rejected, without providing any notion of spatial uncertainty about the activation. In this work, we address these issues by developing spatial Confidence Sets (CSs) on clusters found in thresholded Cohen's d effect size images. We produce an upper and lower CS to make confidence statements about brain regions where Cohen's d effect sizes have exceeded and fallen short of a non-zero threshold, respectively. The CSs convey information about the magnitude and reliability of effect sizes that is usually given separately in a t-statistic and effect estimate map. We expand the theory developed in our previous work on CSs for %BOLD change effect maps (Bowring et al., 2019) using recent results from the bootstrapping literature. By assessing the empirical coverage with 2D and 3D Monte Carlo simulations resembling fMRI data, we find our method is accurate in sample sizes as low as N=60. We compute Cohen's d CSs for the Human Connectome Project working memory task-fMRI data, illustrating the brain regions with a reliable Cohen's d response for a given threshold. By comparing the CSs with results obtained from a traditional statistical voxelwise inference, we highlight the improvement in activation localization that can be gained with the Confidence Sets.

Peak p-values and false discovery rate inference in neuroimaging.

Peaks are a mainstay of neuroimage analysis for reporting localization results. The current peak detection procedure in SPM12 requires a pre-threshold for approximating p-values and a false discovery rate (FDR) nominal level for inference. However, the pre-threshold is an undesirable feature, while the FDR level is meaningless if the null hypothesis is not properly defined. This article provides: 1) a peak height distribution for smooth Gaussian error fields, which does not require a screening pre-threshold; 2) a signal-plus-noise model where FDR of peaks can be controlled and properly interpreted. Matlab code for calculation of p-values using the exact peak height distribution is available as an SPM extension.

Spatial confidence sets for raw effect size images.

The mass-univariate approach for functional magnetic resonance imaging (fMRI) analysis remains a widely used statistical tool within neuroimaging. However, this method suffers from at least two fundamental limitations: First, with sufficient sample sizes there is high enough statistical power to reject the null hypothesis everywhere, making it difficult if not impossible to localize effects of interest. Second, with any sample size, when cluster-size inference is used a significant p-value only indicates that a cluster is larger than chance. Therefore, no notion of confidence is available to express the size or location of a cluster that could be expected with repeated sampling from the population. In this work, we address these issues by extending on a method proposed by Sommerfeld et al. (2018) (SSS) to develop spatial Confidence Sets (CSs) on clusters found in thresholded raw effect size maps. While hypothesis testing indicates where the null, i.e. a raw effect size of zero, can be rejected, the CSs give statements on the locations where raw effect sizes exceed, and fall short of, a non-zero threshold, providing both an upper and lower CS. While the method can be applied to any mass-univariate general linear model, we motivate the method in the context of blood-oxygen-level-dependent (BOLD) fMRI contrast maps for inference on percentage BOLD change raw effects. We propose several theoretical and practical implementation advancements to the original method formulated in SSS, delivering a procedure with superior performance in sample sizes as low as N=60. We validate the method with 3D Monte Carlo simulations that resemble fMRI data. Finally, we compute CSs for the Human Connectome Project working memory task contrast images, illustrating the brain regions that show a reliable %BOLD change for a given %BOLD threshold.

Expected Number and Height Distribution of Critical Points of Smooth Isotropic Gaussian Random Fields.

We obtain formulae for the expected number and height distribution of critical points of smooth isotropic Gaussian random fields parameterized on Euclidean space or spheres of arbitrary dimension. The results hold in general in the sense that there are no restrictions on the covariance function of the field except for smoothness and isotropy. The results are based on a characterization of the distribution of the Hessian of the Gaussian field by means of the family of Gaussian orthogonally invariant (GOI) matrices, of which the Gaussian orthogonal ensemble (GOE) is a special case. The obtained formulae depend on the covariance function only through a single parameter (Euclidean space) or two parameters (spheres), and include the special boundary case of random Laplacian eigenfunctions.

MULTIPLE TESTING OF LOCAL MAXIMA FOR DETECTION OF PEAKS IN RANDOM FIELDS.

A topological multiple testing scheme is presented for detecting peaks in images under stationary ergodic Gaussian noise, where tests are performed at local maxima of the smoothed observed signals. The procedure generalizes the one-dimensional scheme of [31] to Euclidean domains of arbitrary dimension. Two methods are developed according to two different ways of computing p-values: (i) using the exact distribution of the height of local maxima, available explicitly when the noise field is isotropic [9, 10]; (ii) using an approximation to the overshoot distribution of local maxima above a pre-threshold, applicable when the exact distribution is unknown, such as when the stationary noise field is non-isotropic [9]. The algorithms, combined with the Benjamini-Hochberg procedure for thresholding p-values, provide asymptotic strong control of the False Discovery Rate (FDR) and power consistency, with specific rates, as the search space and signal strength get large. The optimal smoothing bandwidth and optimal pre-threshold are obtained to achieve maximum power. Simulations show that FDR levels are maintained in non-asymptotic conditions. The methods are illustrated in the analysis of functional magnetic resonance images of the brain.

Characterizing brain tissue by assessment of the distribution of anisotropic microstructural environments in diffusion-compartment imaging (DIAMOND).

PURPOSE: To develop a statistical model for the tridimensional diffusion MRI signal at each voxel that describes the signal arising from each tissue compartment in each voxel. THEORY AND METHODS: In prior work, a statistical model of the apparent diffusion coefficient was shown to well-characterize the diffusivity and heterogeneity of the mono-directional diffusion MRI signal. However, this model was unable to characterize the three-dimensional anisotropic diffusion observed in the brain. We introduce a new model that extends the statistical distribution representation to be fully tridimensional, in which apparent diffusion coefficients are extended to be diffusion tensors. The set of compartments present at a voxel is modeled by a finite sum of unimodal continuous distributions of diffusion tensors. Each distribution provides measures of each compartment microstructural diffusivity and heterogeneity. RESULTS: The ability to estimate the tridimensional diffusivity and heterogeneity of multiple fascicles and of free diffusion is demonstrated. CONCLUSION: Our novel tissue model allows for the characterization of the intra-voxel orientational heterogeneity, a prerequisite for accurate tractography while also characterizing the overall tridimensional diffusivity and heterogeneity of each tissue compartment. The model parameters can be estimated from short duration acquisitions. The diffusivity and heterogeneity microstructural parameters may provide novel indicator of the presence of disease or injury. Magn Reson Med 76:963-977, 2016. © 2015 Wiley Periodicals, Inc.

Lognormal Distributions and Geometric Averages of Symmetric Positive Definite Matrices.

This article gives a formal definition of a lognormal family of probability distributions on the set of symmetric positive definite (SPD) matrices, seen as a matrix-variate extension of the univariate lognormal family of distributions. Two forms of this distribution are obtained as the large sample limiting distribution via the central limit theorem of two types of geometric averages of i.i.d. SPD matrices: the log-Euclidean average and the canonical geometric average. These averages correspond to two different geometries imposed on the set of SPD matrices. The limiting distributions of these averages are used to provide large-sample confidence regions and two-sample tests for the corresponding population means. The methods are illustrated on a voxelwise analysis of diffusion tensor imaging data, permitting a comparison between the various average types from the point of view of their sampling variability.

Clinical trial of a farnesyltransferase inhibitor in children with Hutchinson-Gilford progeria syndrome.

Hutchinson-Gilford progeria syndrome (HGPS) is an extremely rare, fatal, segmental premature aging syndrome caused by a mutation in LMNA that produces the farnesylated aberrant lamin A protein, progerin. This multisystem disorder causes failure to thrive and accelerated atherosclerosis leading to early death. Farnesyltransferase inhibitors have ameliorated disease phenotypes in preclinical studies. Twenty-five patients with HGPS received the farnesyltransferase inhibitor lonafarnib for a minimum of 2 y. Primary outcome success was predefined as a 50% increase over pretherapy in estimated annual rate of weight gain, or change from pretherapy weight loss to statistically significant on-study weight gain. Nine patients experienced a ≥50% increase, six experienced a ≥50% decrease, and 10 remained stable with respect to rate of weight gain. Secondary outcomes included decreases in arterial pulse wave velocity and carotid artery echodensity and increases in skeletal rigidity and sensorineural hearing within patient subgroups. All patients improved in one or more of these outcomes. Results from this clinical treatment trial for children with HGPS provide preliminary evidence that lonafarnib may improve vascular stiffness, bone structure, and audiological status.

Distribution of the Height of Local Maxima of Gaussian Random Fields.

Let {f(t) : t ∈ T} be a smooth Gaussian random field over a parameter space T, where T may be a subset of Euclidean space or, more generally, a Riemannian manifold. We provide a general formula for the distribution of the height of a local maximum [Formula: see text] is a local maximum of f(t)} when f is non-stationary. Moreover, we establish asymptotic approximations for the overshoot distribution of a local maximum [Formula: see text] is a local maximum of f(t) and f(t0) > v} as v → ∞. Assuming further that f is isotropic, we apply techniques from random matrix theory related to the Gaussian orthogonal ensemble to compute such conditional probabilities explicitly when T is Euclidean or a sphere of arbitrary dimension. Such calculations are motivated by the statistical problem of detecting peaks in the presence of smooth Gaussian noise.

Voxelwise multivariate analysis of multimodality magnetic resonance imaging.

Most brain magnetic resonance imaging (MRI) studies concentrate on a single MRI contrast or modality, frequently structural MRI. By performing an integrated analysis of several modalities, such as structural, perfusion-weighted, and diffusion-weighted MRI, new insights may be attained to better understand the underlying processes of brain diseases. We compare two voxelwise approaches: (1) fitting multiple univariate models, one for each outcome and then adjusting for multiple comparisons among the outcomes and (2) fitting a multivariate model. In both cases, adjustment for multiple comparisons is performed over all voxels jointly to account for the search over the brain. The multivariate model is able to account for the multiple comparisons over outcomes without assuming independence because the covariance structure between modalities is estimated. Simulations show that the multivariate approach is more powerful when the outcomes are correlated and, even when the outcomes are independent, the multivariate approach is just as powerful or more powerful when at least two outcomes are dependent on predictors in the model. However, multiple univariate regressions with Bonferroni correction remain a desirable alternative in some circumstances. To illustrate the power of each approach, we analyze a case control study of Alzheimer's disease, in which data from three MRI modalities are available.

Gene expression profiling of human breast tissue samples using SAGE-Seq.

We present a powerful application of ultra high-throughput sequencing, SAGE-Seq, for the accurate quantification of normal and neoplastic mammary epithelial cell transcriptomes. We develop data analysis pipelines that allow the mapping of sense and antisense strands of mitochondrial and RefSeq genes, the normalization between libraries, and the identification of differentially expressed genes. We find that the diversity of cancer transcriptomes is significantly higher than that of normal cells. Our analysis indicates that transcript discovery plateaus at 10 million reads/sample, and suggests a minimum desired sequencing depth around five million reads. Comparison of SAGE-Seq and traditional SAGE on normal and cancerous breast tissues reveals higher sensitivity of SAGE-Seq to detect less-abundant genes, including those encoding for known breast cancer-related transcription factors and G protein-coupled receptors (GPCRs). SAGE-Seq is able to identify genes and pathways abnormally activated in breast cancer that traditional SAGE failed to call. SAGE-Seq is a powerful method for the identification of biomarkers and therapeutic targets in human disease.

Evaluation of resampling-based inference for topological features of neuroimages.

Many recent studies have demonstrated the inflated type 1 error rate of the original Gaussian random field (GRF) methods for inference of neuroimages and identified resampling (permutation and bootstrapping) methods that have better performance. There has been no evaluation of resampling procedures when using robust (sandwich) statistical images with different topological features (TF) used for neuroimaging inference. Here, we consider estimation of distributions TFs of a statistical image and evaluate resampling procedures that can be used when exchangeability is violated. We compare the methods using realistic simulations and study sex differences in life-span age-related changes in gray matter volume in the Nathan Kline Institute Rockland sample. We find that our proposed wild bootstrap and the commonly used permutation procedure perform well in sample sizes above 50 under realistic simulations with heteroskedasticity. The Rademacher wild bootstrap has fewer assumptions than the permutation and performs similarly in samples of 100 or more, so is valid in a broader range of conditions. We also evaluate the GRF-based pTFCE method and show that it has inflated error rates in samples less than 200. Our R package, pbj , is available on Github and allows the user to reproducibly implement various resampling-based group level neuroimage analyses.

A Deep Learning Approach to Improve Retinal Structural Predictions and Aid Glaucoma Neuroprotective Clinical Trial Design.

PURPOSE: To investigate the efficacy of a deep learning regression method to predict macula ganglion cell-inner plexiform layer (GCIPL) and optic nerve head (ONH) retinal nerve fiber layer (RNFL) thickness for use in glaucoma neuroprotection clinical trials. DESIGN: Cross-sectional study. PARTICIPANTS: Glaucoma patients with good quality macula and ONH scans enrolled in 2 longitudinal studies, the African Descent and Glaucoma Evaluation Study and the Diagnostic Innovations in Glaucoma Study. METHODS: Spectralis macula posterior pole scans and ONH circle scans on 3327 pairs of GCIPL/RNFL scans from 1096 eyes (550 patients) were included. Participants were randomly distributed into a training and validation dataset (90%) and a test dataset (10%) by participant. Networks had access to GCIPL and RNFL data from one hemiretina of the probe eye and all data of the fellow eye. The models were then trained to predict the GCIPL or RNFL thickness of the remaining probe eye hemiretina. MAIN OUTCOME MEASURES: Mean absolute error (MAE) and squared Pearson correlation coefficient (r2) were used to evaluate model performance. RESULTS: The deep learning model was able to predict superior and inferior GCIPL thicknesses with a global r2 value of 0.90 and 0.86, r2 of mean of 0.90 and 0.86, and mean MAE of 3.72 µm and 4.2 µm, respectively. For superior and inferior RNFL thickness predictions, model performance was slightly lower, with a global r2 of 0.75 and 0.84, r2 of mean of 0.81 and 0.82, and MAE of 9.31 µm and 8.57 µm, respectively. There was only a modest decrease in model performance when predicting GCIPL and RNFL in more severe disease. Using individualized hemiretinal predictions to account for variability across patients, we estimate that a clinical trial can detect a difference equivalent to a 25% treatment effect over 24 months with an 11-fold reduction in the number of patients compared to a conventional trial. CONCLUSIONS: Our deep learning models were able to accurately estimate both macula GCIPL and ONH RNFL hemiretinal thickness. Using an internal control based on these model predictions may help reduce clinical trial sample size requirements and facilitate investigation of new glaucoma neuroprotection therapies. FINANCIAL DISCLOSURE(S): Proprietary or commercial disclosure may be found after the references.

Paradoxical results of adaptive false discovery rate procedures in neuroimaging studies.

Adaptive false discovery rate (FDR) procedures, which offer greater power than the original FDR procedure of Benjamini and Hochberg, are often applied to statistical maps of the brain. When a large proportion of the null hypotheses are false, as in the case of widespread effects such as cortical thinning throughout much of the brain, adaptive FDR methods can surprisingly reject more null hypotheses than not accounting for multiple testing at all-i.e., using uncorrected p-values. A straightforward mathematical argument is presented to explain why this can occur with the q-value method of Storey and colleagues, and a simulation study shows that it can also occur, to a lesser extent, with a two-stage FDR procedure due to Benjamini and colleagues. We demonstrate the phenomenon with reference to a published data set documenting cortical thinning in attention deficit/hyperactivity disorder. The paper concludes with recommendations for how to proceed when adaptive FDR results of this kind are encountered in practice.

Functional delta residuals and applications to simultaneous confidence bands of moment based statistics.

Given a functional central limit (fCLT) for an estimator and a parameter transformation, we construct random processes, called functional delta residuals, which asymptotically have the same covariance structure as the limit process of the functional delta method. An explicit construction of these residuals for transformations of moment-based estimators and a multiplier bootstrap fCLT for the resulting functional delta residuals are proven. The latter is used to consistently estimate the quantiles of the maximum of the limit process of the functional delta method in order to construct asymptotically valid simultaneous confidence bands for the transformed functional parameters. Performance of the coverage rate of the developed construction, applied to functional versions of Cohen's d, skewness and kurtosis, is illustrated in simulations and their application to test Gaussianity is discussed.

Impact of Changing Winds on the Mauna Loa CO2 Seasonal Cycle in Relation to the Pacific Decadal Oscillation.

Long-term measurements at the Mauna Loa Observatory (MLO) show that the CO2 seasonal cycle amplitude (SCA) increased from 1959 to 2019 at an overall rate of 0.22  ±  0.034 ppm decade-1 while also varying on interannual to decadal time scales. These SCA changes are a signature of changes in land ecological CO2 fluxes as well as shifting winds. Simulations with the TM3 tracer transport model and CO2 fluxes from the Jena CarboScope CO2 Inversion suggest that shifting winds alone have contributed to a decrease in SCA of -0.10  ±  0.022 ppm decade-1 from 1959 to 2019, partly offsetting the observed long-term SCA increase associated with enhanced ecosystem net primary production. According to these simulations and MIROC-ACTM simulations, the shorter-term variability of MLO SCA is nearly equally driven by varying ecological CO2 fluxes (49%) and varying winds (51%). We also show that the MLO SCA is strongly correlated with the Pacific Decadal Oscillation (PDO) due to varying winds, as well as with a closely related wind index (U-PDO). Since 1980, 44% of the wind-driven SCA decrease has been tied to a secular trend in the U-PDO, which is associated with a progressive weakening of westerly winds at 700 mbar over the central Pacific from 20°N to 40°N. Similar impacts of varying winds on the SCA are seen in simulations at other low-latitude Pacific stations, illustrating the difficulty of constraining trend and variability of land CO2 fluxes using observations from low latitudes due to the complexity of circulation changes.

MULTIPLE TESTING OF LOCAL MAXIMA FOR DETECTION OF PEAKS IN 1D.

A topological multiple testing scheme for one-dimensional domains is proposed where, rather than testing every spatial or temporal location for the presence of a signal, tests are performed only at the local maxima of the smoothed observed sequence. Assuming unimodal true peaks with finite support and Gaussian stationary ergodic noise, it is shown that the algorithm with Bonferroni or Benjamini-Hochberg correction provides asymptotic strong control of the family wise error rate and false discovery rate, and is power consistent, as the search space and the signal strength get large, where the search space may grow exponentially faster than the signal strength. Simulations show that error levels are maintained for nonasymptotic conditions, and that power is maximized when the smoothing kernel is close in shape and bandwidth to the signal peaks, akin to the matched filter theorem in signal processing. The methods are illustrated in an analysis of electrical recordings of neuronal cell activity.