NeuroImage 47 (2009) 213–219

Contents lists available at ScienceDirect

NeuroImage j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / y n i m g

Technical Note

Local and covariate-modulated false discovery rates applied in neuroimaging Glenn Lawyer a,⁎, Egil Ferkingstad b, Ragnar Nesvåg c, Katarina Varnäs d, Ingrid Agartz a,c,d a

Department of Psychiatry, University of Oslo, Oslo, Norway Statistics for Innovation, Norwegian Computing Center, Oslo, Norway c Department of Psychiatric Research, Diakonhjemmet Hospital, Oslo, Norway d Department of Clinical Neuroscience, Karolinska Hospital, Stockholm, Sweden b

a r t i c l e

i n f o

Article history: Received 11 September 2008 Revised 11 March 2009 Accepted 18 March 2009 Available online 31 March 2009 Keywords: Inference Empirical Bayes Mixture model Magnetic resonance imaging

a b s t r a c t False discovery rate (FDR) control has become a standard technique in neuroimaging. Recent work has shown that a finer grained estimate of the FDR is obtained by estimating, at a specific value of the test statistic, the scaled ratio of the null density to the observed density of the test statistic. The method can be extended by allowing an external covariate, also measured on the points where the hypothesis was tested, to modulate estimation of this local FDR. The current work, in addition to demonstrating these methods by reanalyzing results from two previously published investigations of cortical thickness, presents a method to test if the covariate modulation differs significantly from chance. The first study compared schizophrenia patients to healthy controls and the second compared genotypes of the −633 T/A polymorphism of the gene coding the brain derived neurotrophic factor (BDNF) protein in a subset of the subjects from the case/control study. Local FDR estimates increased findings over FDR in both studies. Using p-values from the case/control study to modulate local FDR estimation in the BDNF study further increased findings. The relationship between case/control related and BDNF related cortical thickness variation was found to be highly significant, providing support for this gene's involvement in the etiology of the disease. The increased statistical precision from more accurate models of the distribution of the test statistic demonstrates the potential of these methods for neuroimaging and suggests the possibility to test novel hypothesis. © 2009 Elsevier Inc. All rights reserved.

Introduction Vertex based and voxel based morphometry have become important methods in psychiatric and psychological research. Benjamini and Hochberg's (1995) false discovery rate (FDR) is a popular approach to the multiplicity issue inherent in large-scale hypothesis testing. The technique was introduced to neuroimaging by Genovese et al. (2002). Multiple statistical testing is not unique to neuroscience. False discovery rate control has been heavily investigated in a number of fields, primarily genomics. One important development, the local false discovery rate (local FDR) (Efron and Tibshirani, 2002), begins with the premise that, for each test conducted, the null hypothesis (H0) is either true or false, and that the test statistic will have a different distribution if H0 is true rather than false. In other words, the distribution of the test statistic follows a mixture model. Inference proceeds by comparing, for a given value of the test statistic, the density of the null distribution to the mixture distribution. Inference is hence local in terms of the value of the test statistic. The local FDR measure can be interpreted either in ⁎ Corresponding author. Institute of Psychiatry, University of Oslo, P.O. 85, 0319 Oslo, Norway. Fax: +47 2249 5861. E-mail address: [email protected] (G. Lawyer). 1053-8119/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2009.03.047

terms of controlling the rate of false positives or in Bayesian terms as the posterior probability of the null given the value of the test statistic. The local FDR has received widespread attention. A recent report listed four published methods for estimating the parameters of the mixture model, in addition to the authors' own method (Dalmasso et al., 2007). A further development considers the possibility that the mixture model varies over the set of all tests conducted, and that this variation is associated with an observable covariate. In data from one diffusion tensor imaging experiment (Schwartzman et al., 2005) the observed strength of group differences varied with distance from the back of the skull (Efron, 2008). Ignoring this effect resulted in both under- and over-estimation of significance, whereas taking distance into account increased both the accuracy and the power of the inference (Efron, 2008). An approach for incorporating covariate information into local FDR estimates has been proposed by Ferkingstad et al. (2008), who estimate the covariate-modulated posterior probability (CMPP) of the null hypothesis for each test. The possibility that the distribution of a test statistic may vary based on an outside covariate suggests investigations in which the relationship between the covariate and the main test is the question of interest. The current work presents a novel method to test the significance of this relationship. The method is illustrated by testing

214

G. Lawyer et al. / NeuroImage 47 (2009) 213–219

for a relationship between thinner brain cortex associated with variation in the gene coding for the brain derived neurotrophic factor protein (BDNF) and the thinner cortex frequently observed in patients with schizophrenia. In addition to presenting a method to test for relationships between a covariate and a main effect, this work demonstrates the use of local FDR and CMPP in neuroimaging studies by re-analyzing data from two previous publications. One compared cortical thickness in patients with schizophrenia with healthy controls (Nesvåg et al., 2008), and the other investigated BDNF polymorphisms and cortical thickness in schizophrenia (Varnäs et al., 2008). The subjects in the BDNF study were a subset of the patients from the case/control study. Overview of the methods General overview The local FDR estimates the probability that the null is true for a given value of the test statistic, considering the overall distribution of the test statistic in the experiment. It is most easily conceived of graphically. Consider a curve fit to a histogram of the test statistics from a large-scale testing experiment, say p-values from an investigation of cortical thickness differences (see Fig. 2). The height of the curve at a specific value of p can be modeled as a mixture of the null and some alternative hypothesis: f ðpÞ = π0 f0 ð pÞ + ð1 − π0 Þf1 ðpÞ

ð1Þ

where the mixture parameter π0 is the probability that the null is true, f0(p) is the density under the null hypothesis, and f1(p) the density under the alternative hypothesis. For any given value of p, the local FDR is the scaled ratio of the height predicted under the null to the observed height of the curve (Efron and Tibshirani, 2002): π f ð pÞ : local FDR ðpÞu 0 0 f ð pÞ

ð2Þ

A straightforward application of Bayes' law shows that the local FDR is the posterior probability of the null hypothesis, given the test statistic. π0 f0 ð pÞ prðnullÞprð pjnullÞ u = prðnull jpÞ: f ð pÞ pr ðpÞ

ð3Þ

In this general overview, p can refer to any type of test statistic (i.e. p, z, t, F, χ2, …). Local FDR was developed for large-scale testing situations, in which the goal is to identify cases (vertices or voxels) where the alternative hypothesis is likely to be true, and where a controlled proportion of false positive results is tolerable. Efron (2005) proposed declaring all hypothesis tests with tolerable local FDR scores “interesting,” reserving the word “significant” for true significance tests, i.e. investigations seeking to establish one (or a few) strongly motivated hypothesis by showing that the null is unlikely. It can sometimes happen that a factor which is observed at, but not included in, each hypothesis test is believed to influence the distribution of the test statistic. Genetic effects, for example, would be more likely at highly heritable regions than at regions with low heritability. The CMPP approach addresses such covariate influence by allowing π0 and f1(p) to vary across the different hypothesis tests, with this variation dependent on a covariate x also observed at each hypothesis test. Inference is based on the pair (pi, xi), where i indexes a specific hypothesis test. The mixture model describing the distribution of the test statistic (Eq. (1)) becomes: f ðp jxÞ = π0 ðxÞf0 ðpÞ + ð1 − π0 ðxÞÞf1 ðp jxÞ:

ð4Þ

The CMPP is defined as CMPP ð pjxÞuπ0 ðxÞf0 ðpÞ = f ðpj xÞ:

ð5Þ

One simple and transparent method to resolve the dependence on x when estimating the CMPP is to bin the paired data (pi, xi) into M bins, B1,B2,…BM, increasing in x. Bins should be chosen small enough that the influence of x is nearly constant in each bin; in practice between 10 and 20 bins generally suffices. The local FDR is estimated in each bin, possibly with smoothing across the bins. The estimated CMPP of (pi, xi) is the estimated local FDR of pi in the bin containing xi. It is possible that the covariate has no effect on the distribution, despite the investigator's belief. If the local FDR has in fact no dependence on x, then the true local FDR would be the same in each bin. The method reverts to a slightly less efficient estimate of the local FDR. It may be that the relationship between the main test and the proposed covariate is itself a question of interest. This allows the framing of novel hypotheses which can be tested against the null hypothesis of no dependence. The above observation of CMPP's behavior when there is no dependence suggests the following significance test. In computing the CMPP, each of the M bins has its own estimate for the scalar quantity π0. Under the null hypothesis this value would be nearly the same in each bin. Collect the estimates in a vector π = [π10,π20,…πM 0 ] (the superscript indexes the bin). Compare the observed range of π range ðπ Þ = maxðπ Þ − minðπ Þ

ð6Þ

to the null distribution of π's range. A null distribution can be computed by permutation testing, i.e. measuring the range of π under repeated random assignment of the xi to the pi. Applications The overview presented three inferential concepts: the local FDR; covariate-modulated local FDR; and testing if such modulation is statistically significant. We demonstrate these using results from two previous studies of cortical thickness variation published by our group. The first study compared patients with schizophrenia with healthy controls. Widespread areas of strongly significant difference between the two groups were observed (Nesvåg et al., 2008). The second study, using the same subject group, compared subjects carrying different variants of single-nucleotide polymorphisms (SNPs) on the BDNF gene. Analysis of the −633 T/A SNP yielded weak findings of difference between the AA and TT carriers in a small region of the frontal lobe in patients with schizophrenia but not in healthy controls (Varnäs et al., 2008); we re-analyze only this section of the study here. The analysis plan is as follows. We begin by estimating the local FDR for each of the studies. Comparison is made to Benjamini and Hochberg's FDR by counting the number of cortical locations considered interesting under each technique. Our example studies show how the two approaches compare both in situations when findings are widespread (case/control) and when findings are weak (BDNF). CMPP is demonstrated by using the p-values from the case/ control study to modulate inference in the BDNF study. We count the number of results found interesting using CMPP, local FDR, and FDR. The covariate choice is motivated by the twin observations that the BDNF findings occurred only in the patients (Varnäs et al., 2008) and that BDNF gene variants may be associated with risk of schizophrenia (Jönsson et al., 2006). Note that the covariate here is used as a covariate of the inferential model. The dependent variables are parameters of the local FDR equation. This is quite different from the more familiar use of covariates (age, sex, etc.) to

G. Lawyer et al. / NeuroImage 47 (2009) 213–219

215

Fig. 1. Uncorrected p-values from the case/control study (left) and BDNF study (right). Colors are coded by −log10(p), i.e. red (1.3) is p b 0.05. Blue indicates p b 0.05 regions with opposite direction of effect; most such regions fell outside of the regions considered in the current analysis.

explain a dependent variable (cortical thickness, haemodynamic response, etc.) in the model on which inference is performed. Finally, the use of CMPP to frame and test novel hypothesis is demonstrated by measuring the statistical significance of the influence of schizophrenia related cortical thinning on inference of BDNF related cortical thinning. This putative relationship implied that BDNF related effects were more likely at locations with strong patient/control differences, or, in terms of the CMPP approach, the estimated value of π0 would be low when patient/control difference were strong and high when patient/control differences were minimal. Subject material and preprocessing Subject demographics, scan acquisition, and scan processing are fully described in the relevant publications. Subjects in the case/ control study were 203 unrelated Caucasian individuals living in Stockholm county in Sweden. Of these, 96 were patients with chronic schizophrenia recruited from outpatient clinics and 107 were healthy control subjects. The section of the BDNF study reanalyzed here relied on a subset of the patients from the case/ control study consisting of carriers of the TT (n = 32) and AA (n = 24) variant of the BDNF −633 T/A polymorphism. Patients carrying the AT variant were excluded, as were all healthy subjects. Subject recruitment and scan acquisition was conducted as part of

the HUBIN project (Hall et al., 2000; Arnborg et al., 2000) at Karolinska Institutet, Stockholm, Sweden. Both studies were conducted in accordance with the Declaration of Helsinki and approved by the Ethics Committee of the Karolinska Hospital and the Swedish Data Inspection Board (“Datainspektionen”). All subjects participated after giving informed written consent, including consent for the data to be used for further analysis. T1-weighted MR images were acquired from each subject using a 3D spoiled gradient pulse recall sequence on a 1.5 Tesla GE Signa (GE, Milwaukee, Wis, USA) system at the Magnetic Resonance Research Center, Karolinska Hospital, Stockholm, Sweden, between 1999 and 2003. Cortical thickness of each subject was measured as the shortest distance between the pial and white matter surfaces at each vertex, as determined by FreeSurfer stable release 1.2 (Dale et al., 1999; Fischl et al., 1999; Fischl and Dale, 2000). Thickness maps were registered to a common coordinate system aligned across subjects according to cortical folding patterns using a nonrigid high dimensional spherical method (Fischl et al., 1999) and based on a template brain made by averaging the 203 subjects in the case/control study. The same template was used for both reanalyses. Images were smoothed with a 10 mm Gaussian kernel. Group differences were measured as p-values generated via contrast analysis of a general linear model fitted independently at each vertex in the aligned and co-registered cortical thickness

Fig. 2. Histograms of the p-values from the case/control study (left) and BDNF study (right). The red vertical line marks p = 0.05. The blue curve is the estimated density f(p). Note that the case/control histogram only shows the region p ∈ (0,0.1), beyond which the histogram was level.

216

G. Lawyer et al. / NeuroImage 47 (2009) 213–219

was made of the number of vertices found interesting at each rate for each method. A further count was made of the number of vertices passing the arbitrary threshold of p ≤ 0.001. To illustrate the behavior of the algorithms under differing conditions, comparison was made considering vertices from brain regions of decreasing size (see Fig. 3): 1. the left hemisphere rostral to the superior portion of the central sulcus (RLH) 2. a large region of the left prefrontal cortex (PFC) 3. a small subregion of the left prefrontal cortex (subPFC). Regions were defined by manually tracing on the study's template brain and were not intended to perfectly represent standard anatomical distinctions.

Fig. 3. The regions examined were the rostral left hemisphere (RLH, yellow + green + red), the prefrontal cortex (PFC, green + red π), and a subsection of the prefrontal cortex (subPFC, red).

maps which tested for a difference between the two groups after controlling for age. Fig. 1 shows results from this re-analysis of the two studies, presented as −log10(p) thresholded to show values less than p = 0.05. Fig. 2 shows histograms of the p-values. Estimating local FDR and CMPP Local FDR was estimated independently in each study, using the method suggested by Ferkingstad et al. (2008). The method applies to p-values, thus the density f0(p) was uniform. The density f1(p) was modeled as a convex decreasing beta distribution, which allowed the joint posterior distribution of π0 and the parameters ξ and θ of the beta distribution be written as f ðπ0 ; n; θ jp1 ; N ; pn Þ~

n Y

½π0 + ð1 − π0 Þf1 ðpi Þ × f ðπ 0 Þf ðnÞf ðθÞ

ð7Þ

i=1

where p1,…,pn were the observed p-values from the experiment, f1 ð p Þ =

Cðn + θÞ n − 1 θ−1 p ð1− pÞ CðnÞCðθÞ

ð8Þ

Assessing significance of relationship We hypothesized that BDNF related cortical thickness differences were more likely at locations with strong patient/control differences. The significance of this hypothesis was measured by assessing the probability of the observed range of π under the (null) hypothesis of no relationship between these two factors. The null distribution was simulated by computing, 10,000 times, the range of π after random reassignment of the covariate values to the p-values. A twenty bin model was used. Tests were made for both the PFC and the subPFC. Results In the case/control study, local FDR showed a 10% increase in findings over FDR at the 5% threshold, an approximate 2% increase at the 10% threshold, and a 1% decrease at the 20% threshold, over both the RLH and the PFC. Very little difference in findings was seen for the subPFC, with both methods concluding that this region consisted almost exclusively of interesting cortical thickness differences. Both methods provided substantial gains over arbitrarily thresholding the p-values at p ≤ 0.001, for all regions. Table 1 shows the percentage of each region found interesting at each threshold. Fig. 4 shows the estimated local FDR over the RLH in the case/control study. The original BDNF paper reported no findings as surviving FDR correction at a rate of 5% applied to the entire cortical hemisphere. Similar results were observed here, with FDR reporting no interesting results for the RLH and only finding interesting results in the PFC at the 20% rate. The local FDR maintained this skeptical view of the data at the 5% rate, but found 4% of the PFC interesting at the 10% rate and twice as many interesting results as FDR at the 20% rate. Including covariate information further increased findings. CMPP thresholded at a 5% rate found that −633 T/A BDNF gene variation may influence cortical thickness over 12% of the PFC. Results for local FDR are shown in Table 2 and for CMPP in Table 3. Fig. 4 shows the estimated CMPP over the PFC in the BDNF study. Schizophrenia related cortical thinning proved to have a strongly significantly effect on the local FDR estimates in the BDNF study. In the PFC, the simulated null distribution of π's range varied between 0.02 and 0.17. The observed range was 0.96. In the subPFC, the null

and where π0, ξ, and θ were transformed to allow unconstrained parametrization. The parametric model was justified by Allison et al's (2002) demonstration that the density underlying a set of observed p-values can be well approximated by a mixture of beta distributions, and that generally one uniform and one beta distribution suffice. In the context of local FDR, using only one beta is a conservative choice in that it avoids overestimating the proportion of false null hypothesis (Ferkingstad et al., 2008). The joint posterior was approximated as a Gaussian Markov random field (Rue and Martino, 2007), providing significant times savings compared with Markov Chain Monte Carlo approaches. CMPP was estimated in the BDNF study using the p-values from the case/control study as the covariate. Since a common template brain was used to generate each set of p-values, pairing the main pvalue with the covariate was straightforward. The BDNF p-values were binned based on the covariate and the local FDR was estimated in each bin. This estimation used a modified form of Eq. (7) which included smoothing parameters for π0, ξ, and θ on the right hand side. This encouraged parameter values in neighboring bins to be similar. CMPP was calculated with both ten and twenty bins. Results were similar for both models. Only the twenty bin results are reported.

RLH PFC subPFC

Comparison of local FDR and CMPP to FDR The efficacy of local FDR and of CMPP was assessed by comparison to FDR. Comparison was made at three rates: 5%, 10%, and 20%. A count

Percentage of each region found interesting at FDR and local FDR of 5%, 10%, and 20%. For comparison, the percentage of the region with p ≤ 0.001 is also given. Abbreviations: FDR = false discovery rate; RLH = rostral left hemisphere; PFC = prefrontal cortex; subPFC = subsection of the PFC.

Table 1 Case/control data. Region

p ≤ 0.001 20% 32% 42%

5%

10%

20%

FDR

Local FDR

FDR

Local FDR

FDR

Local FDR

50% 81% 98%

54% 89% 98%

62% 90% 99%

64% 91% 98%

73% 95% 100%

72% 93% 99%

G. Lawyer et al. / NeuroImage 47 (2009) 213–219

217

Fig. 4. Local FDR rates for case/control differences in schizophrenia, calculated over the RLH (left), and CMPP rates for the BDNF study, as calculated over the PFC (right). The yellow line on the right figure marks the boundary of the PFC. Colors are coded by −log10(X), where X is either the local FDR (left) or the CMPP (right). At this scale 1.3 (red) is X b 0.05.

distribution varied between 0.05 and 0.35. The observed value was 0.52. Fig. 5 shows the simulated nulls for the PFC and the subPFC along with the observed ranges of π. Fig. 6 shows the value of π0 in each covariate bin as estimated over the PFC. Discussion Use of the local FDR increased findings compared to FDR in both studies when the false discovery proportion was tightly controlled. Applied to the numerous and strong findings of the case/control study, the difference between the two methods decreased as the allowed proportion of false discoveries increased. Applied to the few and weak findings of the BDNF study, local FDR showed substantial increase in findings over FDR at all thresholds examined. Use of CMPP further increased findings. Studies on simulated data suggest these increases reflect greater statistical power (Ferkingstad et al., 2008). The test of the strength of covariate modulation found the relationship between the BDNF results and the thinner cortex observed in schizophrenia strongly significant, a previously untested relationship. The local FDR estimate used here was based on p-values. Current neuroimaging software typically generates p-values by assuming the theoretical null for the underlying test statistic. The large number of tests conducted, however, allows empirical estimation of this underlying null. Some evidence suggests such estimation is desirable. Correlations in the underlying data could contract or dilate the theoretical null (Efron, 2004; Schwartzman et al., 2009). As permutation and/or randomization methods would not necessarily reveal such features, Schwartzman et al. (2009) provided parametric forms of four distributions commonly used in neuroimaging (z, t, χ2, and F) which can be fit to histograms of observed test statistics with Poisson regression. These parametric forms could alternately be combined with an estimate f(p) to compute the local FDR directly. The

approach critically requires that the observed statistics are mainly null cases and is less reliable when π0 b 90% (Efron, 2004; Schwartzman et al., 2009). This condition was violated for some of the regions investigated here, making an empirical null problematic for the current investigation. One very natural question regarding the CMPP is why covariates which are believed to have an effect are not simply included in the statistical model? Neuroimaging software toolkits which allow investigations with spatially varying design matrices are available (Casanova et al., 2007; Oakes et al., 2007). The answer is that a covariate in the inferential model, as in CMPP, serves a fundamentally different purpose than a covariate in the statistical model. Intuitively, the CMPP approach can be thought of as using a covariate to segment the analysis into regions of (increasing) interest. For example, the likelihood that a strong observed genetic effect represented a true positive would be larger in a region with high heritability than in a region with limited heritability. CMPP accommodates this by varying π0 and f1(p) in the local FDR estimation (see Eqs. (4) and (5)). A covariate in the statistical model, by contrast, modulates the data. One may, to continue the example, wish to remove age and gender effects from the data before testing for a genetic effect. The biological observation that true effects in the brain have spatial extent suggests using a localized measure of the spatial correlation structure as a covariate. Such a measure would make the CMPP approach applicable to almost any neuroimaging study. It is a problematic hope. If the brain images were smoothed, one would be “double-dipping” as the smoothing would influence both the outputs of the statistical tests and the inference based on those outputs. Even ignoring this, some evidence suggests that the benefit of incorporating voxel-wise correlation information into FDR-styled inference procedures is marginal (Logan and Rowe, 2004). An appropriate FDR threshold for declaring results “interesting” has yet to be agreed on. The standard in neuroimaging seems to be settling on 5%, though this appears to be primarily based on familiarity with the 5% level from traditional significance testing

Table 2 BDNF data. Region

p ≤ 0.001

5% FDR

Local FDR

FDR

Local FDR

FDR

Local FDR

RLH PFC subPFC

0% 0% 1%

0% 0% 0%

0% 0% 22%

0% 0% 22%

0% 4% 47%

0% 13% 53%

0% 27% 61%

10%

20%

Percentage of each region found interesting at FDR and local FDR of 5%, 10%, and 20%. For comparison, the percentage of the region with p ≤ 0.001 is also given. Abbreviations: FDR = False Discovery Rate; RLH = rostral left hemisphere; PFC = prefrontal cortex; subPFC = subsection of the PFC.

Table 3 Percentage of the region found interesting in the BDNF data at a CMPP of 5%, 10%, and 20%. Region

CMPP ≤ 5%

CMPP ≤ 10%

CMPP ≤ 20%

RLH PFC subPFC

0% 12% 37%

0% 19% 56%

1% 35% 69%

Abbreviations: CMPP = Covariate-modulated posterior probability; RLH = rostral left hemisphere; PFC = prefrontal cortex; subPFC = subsection of the PFC.

218

G. Lawyer et al. / NeuroImage 47 (2009) 213–219

rather than on mathematical argument. Efron suggests 20% for local FDR in cases where π0 ≈ 0.90 (Efron, 2005). This level is motivated by comparison to the weight of evidence in favor of the null required in traditional statistical testing (Efron and Gous, 2001). FDR can be compared mathematically to local FDR by recalling that FDR represents the expected proportion of false positives given p ≤ thresh. FDR, then, is based on tail areas. The local FDR, by contrast, is an estimate of the proportion of false positives for a specific p-value. FDR can be expressed in terms of local FDR as the expectation of all local FDR scores within the tail area (Efron, 2005): FDR ðpÞ = E½local FDRðP Þ jP V p:

ð9Þ

This implies that FDR serves as a lower bound on local FDR, that is, local FDR at worst provides the same level of multiple comparison control as FDR, and potentially better control. The conception of local FDR is generally credited to Efron and Tibshirani (2002), with much development of the idea taking place in the context of DNA microarray experiments. An earlier neuroimaging paper, however, estimated the posterior probability of a voxel being active, defined as in Eq. (2), using a parametric form of Eq. (1) (Everitt and Bullmore, 1999). The approach appeared to identify activated regions far more distinctly than inference based on pvalues alone. Posterior Probability Mapping (Friston and Penny, 2003) also provides an empirical Bayesian estimation of the posterior probability of activation. Estimation is based on the general linear model applied independently at a given voxel, and cannot be considered a form of local FDR. The article mentions that thresholding posterior probabilities is similar in spirit to FDR control. It should be noted that FDR, unlike local FDR, only has a Bayesian interpretation under special circumstances. Schwartzman et al. (2009) argues, however, that when p-values are generated with an empirical null, FDR can be interpreted as a posterior probability even when the test statistics show dependency. Inference based on local FDR cannot be directly compared to cluster-based inference as the two methods are fundamentally different in their approach. The two methods can, however, be combined. Inference in clustering techniques is independent of the actual measure used to determine the cluster. If empirical methods such as those proposed here do offer increased sensitivity, nothing prevents their use as inputs to a clustering scheme. The current study found strong support for association between the BDNF −633 T/A gene variant and the thinner cortex observed in

Fig. 6. Estimated value of π0, with 95% credible interval, in each bin of the CMPP in the BDNF study, as calculated over the PFC.

schizophrenia. The relationship, however, is unlikely to be direct. A direct relationship would have given stronger results in Varnäs et al. (2008). It is more likely that BDNF is involved in a regulatory network which has become disturbed in schizophrenia, making patients with the TT variant of the −633 T/A allele more vulnerable to thinning. The decision not to use synthetic data in the current study was motivated by two factors. The algorithms presented, while new to neuroimaging, have been extensively tested in other domains. The critically different feature of neuroimaging data is the spatial characteristics of the signal, which the algorithms do not incorporate. Given this, it was not clear that the presentation and discussion of synthetic data would justify the additional space required. Conclusions False discovery rates estimated at each specific value of the test statistic are more precise than estimates based on tail areas of the statistic's distribution. Here the local FDR increased findings in two brain cortical thickness investigations. In situations where an additional relevant measure is available, precision is further increased by allowing this covariate to modulate local FDR estimation. Here using CMPP markedly strengthened the evidence that the BDNF gene influenced cortical thickness in prefrontal regions. Novel hypotheses can be tested by measuring if the covariate's modulating effect differs significantly from chance. Here this provided evidence that the BDNF gene may be involved in the etiology of schizophrenia. Acknowledgments

Fig. 5. Simulated null distributions of the range of π for both the PFC and the subPFC. The observed range is indicated by the diamond.

This study was financially supported by the Wallenberg Foundation, the Swedish Research Council (2007–3687), the Norwegian Research Council (160181/V50), and South-Eastern Norway Regional Health Authority (2005;A135). None of the funding organizations took part in the collection, management, analysis or interpretation of the data. Monica Hellberg, Emma Bonnet and Lilian Frygnell have assisted in recruitment and handling of patients at the Karolinska Hospital throughout the study period. The MRI processing and data analysis was performed at the departments of Psychiatry and Psychology at the University of Oslo. Our sincere gratitude is given to Arnoldo Frigessi for many interesting statistical

G. Lawyer et al. / NeuroImage 47 (2009) 213–219

discussions and general support during the preparation of this manuscript. The manuscript profited greatly from the dedication of the reviewers. References Allison, D., Gadbury, G., Heo, M., Fernández, J., Lee, C., Prolla, T., Weindruch, R., 2002. A mixture model approach for the analysis of microarray gene expression data. Comput. Stat. Data Anal. 1, 1–20. Arnborg, S., Agartz, I., Nordström, M., Hall, H., Sedvall, G., 2000. Human Brain Informatics — understanding causes of mental illness. ERCIM News 3, 2s4–25. Benjamini, Y., Hochberg, Y., 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. Royal Stat. Soc. B 57 (1), 289–300. Casanova, R., Srikanth, R., Baer, A., Laurienti, P.J., Burdette, J.H., Hayasaka, S., Flowers, L., Wood, F., Maldjian, A., 2007. Biological parametric mapping: a statistical toolbox for multimodality brain image analysis. NeuroImage 34 (1), 137–143. Dale, A., Fischl, B., Sereno, M., February 1999. Cortical surface-based analysis — I: segmentation and surface reconstruction. NeuroImage 9 (2), 179–194. Dalmasso, C., Bar-Hen, A., Broët, P., 2007. A constrained polynomial regression procedure for estimating the local false discovery rate. BMC Bioinformatics 8, 229. Efron, B., 2004. Large-scale simultaneous hypothesis testing: the choice of a null hypothesis. J. Am. Stat. Assoc. 99, 96–104. Efron, B., 2005. Local false discovery rates. http://www-stat.stanford.edu/ckirby/brad/ papers/2005LocalFDR.pdf. Efron, B., 2008. Simultaneous inference: when should hypothesis testing problems be combined? Ann. Appl. Stat. 2 (1), 197–223. Efron, B., Gous, A., 2001. Scales of evidence for model selection: Fisher versus Jeffreys. Model Selection IMS Monograph 38, 208–246. Efron, B., Tibshirani, R., 2002. Empirical Bayes methods and false discovery rates for microarrays. Genet. Epidemiol. 23, 70–86. Everitt, B., Bullmore, E., 1999. Mixture model mapping of brain activation in functional magnetic resonance images. Hum. Brain Mapp. 7 (1), 1–14. Ferkingstad, E., Frigessi, A., Rue, H., Thorleifsson, G., Kong, A., 2008. Unsupervised empirical Bayesian multiple testing with external covariates. Ann. Appl. Stat. 2 (2), 714–735.

219

Fischl, B., Dale, A., 2000. Measuring the thickness of the human cerebral cortex from magnetic resonance images. Proc. Natl. Acad. Sci. U. S. A. 97, 11044–11049. Fischl, B., Sereno, M., Dale, A., February 1999. Cortical surface-based analysis — II: inflation, flattening, and a surface-based coordinate system. NeuroImage 9 (2), 195–207. Friston, K., Penny, W., 2003. Posterior probability maps and SPMs. NeuroImage 19 (3), 1240–1249. Genovese, C., Lazar, N., Nichols, T., 2002. Thresholding of statistical maps in functional neuroimaging using the false discovery rate. NeuroImage 15 (4), 870–878. Hall, H., McNeil, T., Arnborg, S., Agartz, I., Ösby, U., Linder, J., Sedvall, G., 2000. HUBIN — Human Brain Informatics: a clinical database project for multidisciplinary research in schizophrenia. Eur. Psychiatry 15 (suppl 2), 299s. Jönsson, E., Edman-Ahlbom, B., Sillén, A., Gunnar, A., Kulle, B., Frigessi, A., Vares, M., Ekholm, B., Wode-Helgodt, B., Schumacher, J., Cichon, S., Agartz, I., Sedvall, G., Hall, H., Terenius, L., 2006. Brain-derived neurotrophic factor gene (BDNF) variants and schizophrenia: an association study. Prog. Neuropsychopharmacol. Biol. Psychiatry 30 (5), 924–933. Logan, B.R., Rowe, D.B., May 2004. An evaluation of thresholding techniques in fMRI analysis. NeuroImage 22 (1), 95–108. Nesvåg, R., Lawyer, G., Varnäs, K., Fjell, A., Walhovd, K., Frigessi, A., Jönsson, E., Agartz, I., 2008. Regional thinning of the cerebral cortex in schizophrenia: effects of diagnosis, age and antipsychotic medication. Schizophr Res. 98 (1–3), 16–28. Oakes, T.R., Fox, A.S., Johnstone, T., Chung, M.K., Kalin, N., Davidson, R.J., Jan 2007. Integrating VBM into the general linear model with voxelwise anatomical covariates. NeuroImage 34 (2), 500–508. Rue, H., Martino, S., 2007. Approximate Bayesian inference for hierarchical Gaussian Markov random field models. J. Stat. Plan Inference 10 (1), 3177–3192. Schwartzman, A., Dougherty, R., Taylor, J., 2005. Cross-subject comparison of principal diffusion direction maps. Magn. Reson. Med. 53 (6), 1423–1431. Schwartzman, A., Dougherty, R., Lee, J., Ghahremani, D., Taylor, J., 2009. Empirical null and false discovery rate analysis in neuroimaging. NeuroImage 44 (1), 71–82. Varnäs, K., Lawyer, G., Jönsson, E., Kulle, B., Nesvåg, R., Hall, H., Terenius, L., Agartz, I., 2008. BDNF polymorphisms and frontal cortex morphology in schizophrenia. Psychiatr. Genet. 18 (4), 177–183.