Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Am Stat Assoc. Author manuscript; available in PMC 2010 December 1.
Published in final edited form as:
J Am Stat Assoc. 2010 September 1; 105(491): 1042–1055.
doi:  10.1198/jasa.2010.tm09129
PMCID: PMC2967047

Correlated z-values and the accuracy of large-scale statistical estimates


We consider large-scale studies in which there are hundreds or thousands of correlated cases to investigate, each represented by its own normal variate, typically a z-value. A familiar example is provided by a microarray experiment comparing healthy with sick subjects' expression levels for thousands of genes. This paper concerns the accuracy of summary statistics for the collection of normal variates, such as their empirical cdf or a false discovery rate statistic. It seems like we must estimate an N by N correlation matrix, N the number of cases, but our main result shows that this is not necessary: good accuracy approximations can be based on the root mean square correlation over all N · (N − 1)/2 pairs, a quantity often easily estimated. A second result shows that z-values closely follow normal distributions even under non-null conditions, supporting application of the main theorem. Practical application of the theory is illustrated for a large leukemia microarray study.

Keywords: rms correlation, non-null z-values, correlation penalty, Mehler's identity, empirical process, acceleration

1 Introduction

Modern scientific studies routinely produce data on thousands of related situations. A familiar example is a microarray experiment in which thousands of genes are being investigated for possible disease involvement. Each gene might produce a z-value, say zi, for the ith gene, by definition a test statistic theoretically having a standard normal distribution


under the null hypothesis H0 of no disease involvement. A great deal of the current literature was developed under the assumption of independence among the zi's. This can be grossly unrealistic in practice, as discussed in Owen (2005) and Efron (2007a), among others. This paper concerns the accuracy of summary statistics of the zi's, for example, their empirical cdf (cumulative distribution function), under conditions of substantial correlation.

Figure 1 concerns a leukemia microarray study by Golub et al. (1999) that we will use for motivation and illustration. Two forms of leukemia are being examined for possible genetic differences: acute lymphoblastic leukemia (ALL) and acute myeloid leukemia (AML). In the version of the data discussed here there are n1 = 47 ALL patients and n2 = 25 AML patients, with expression levels on the same N = 7128 genes measured on each patient.

Figure 1
Histogram of z-values for N = 7128 genes, leukemia study, Golub et al. (1999). Dashed curve f(x), smooth fit to histogram; solid curve “empirical null”, normal density fit from central 50% of histogram, is much wider than ...

A two-sample t-statistic ti comparing AML with ALL expression levels was computed for each gene and converted to a z-value,


where Φ and F70 are the cumulative distribution functions for a standard normal and a Student-t distribution with 70 degrees of freedom. Figure 1 shows a histogram of the zi's, which turns out to be much wider than (1.1) suggests: its central spread is estimated to be [sigma with hat]0 = 1.68 rather than 1, as discussed in Section 3.

Here is an example of the results to be derived in Sections 2 through 4. Let F(x) be the right-sided cdf (“survival curve”) of the z-values,


Then a good approximation for the variance of F(x) is


The first term in (1.4) is the usual binomial variance, while the second term is a correlation penalty accounting for dependence between the zi's. The quantities occurring in the correlation penalty are

  • [sigma with hat]0, the estimate of central spread (1.68 above);
  • [alpha], an estimate of the root-mean-square of the correlations between the N(N − 1)/2 pairs of zi's (equaling about .11 for the leukemia data, as calculated from the simple formula in Section 3);
  • f(1)(x), the first derivative of a smooth fit to the z-value histogram (estimated by a Poisson spline regression in Figure 1).

The row marked sd^ in Table 1 is the square root of formula (1.4) applied to the leukemia data. F(4) = .025 is seen to have sd^=.0040, more than double sd^0=.0018, the binomial standard deviation obtained by ignoring the second term in (1.4). The permutation standard deviation, obtained from repeated permutations of the 72 patients, is only .0001 at x = 4. Permutation methods, which preserve within-microarray correlations, have been advocated for large-scale hypothesis testing (see Westfall and Young, 1993; Dudoit, Shaffer and Boldrick, 2003, Sect. 2.6), but they are inappropriate for the accuracy considerations of this paper.

Table 1
Estimates of standard deviation for right-sided cdf F(x) (1.3); sd^ square root of formula (1.4); sd^0 square root of first term in (1.4); sd^perm permutation standard deviation. Accuracy of False Discovery Rate estimate Fdr^(x) discussed ...

Formula (1.4), and more ambitious versions that include covariances across different values of x, are derived in Section 2 and Section 3: an exact expression is derived first, followed by a series of simplifying approximations and techniques for their estimation. The basic results are extended to provide more general accuracy estimates in Section 4: comparing, for example, the variability of local versus tail-area false discovery rates.

All of our results depend on the assumption that the zi's are normal with possibly different means and variances,


There is no requirement that they be “z-values” in the hypothesis-testing sense of (1.1), and in fact this paper is more concerned with estimation than testing. However, z-values are ubiquitous in large-scale applications, and not only in the two-sample setting of the leukemia study. Section 5 concerns the non-null distribution of z-values. A theorem is derived justifying (1.5) as a good approximation, allowing results like (1.4) to be applied to the leukemia study z-values. Section 6 and Section 7 close with remarks and a brief summary.

The statistics microarray literature has shown considerable interest in the effects of large-scale correlation, some good references being Dudoit, van der Laan and Pollard (2004), Owen (2005), Qiu, Klebanov and Yakovlev (2005b), Qiu, Brooks, Klebanov and Yakovlev (2005a) and Desai, Deller and McCormick (2009). Efron (2007a) used a z-value setting to examine the effects of correlation on false discovery rate analysis; that paper's Section 2 theorem is a null hypothesis version of the general result developed here. A useful extension along different lines appears in Schwartzman and Lin (2009).

Clarke and Hall's (2009) asymptotic calculations support the use of the independence standard deviation sd^0 in Table 1, even in the face of correlation. The situations they consider are low-correlation by the standard here, with the root-mean-square value [alpha] of (1.4) approaching zero (from their assumption (3.2)). Since [alpha] is often easy to estimate, formulas such as (1.4) provide a quantitative check on the use of sd^0.

2 The distribution of correlated normal variates

Given N correlated normal variates z1, z2, … , zN, with possibly different means and standard deviations, let F(x) denote their right-sided empirical cdf1


This section presents tractable formulas for the mean and covariance of the process {F(x), −∞ < x < ∞}, and a simpler approximation that we will see is nicely suited for applications.

Rather than work directly with cdfs, it will be easier, and in a sense more basic, to first derive results for a discretized version of the empirical density of the zi values. We partition the range Ƶ into K bins Ƶk,


each bin being of width Δ. Let xk indicate the midpoint of Ƶk, and yk the number of zi's in Ƶk,


We will derive expressions for the mean and covariance of the vector y = (y1, y2, … , yK)′. In effect, y is the order statistic of z = (z1, z2, … , zN)′, becoming exactly that as the bin width Δ → 0. (In which case the yk values go to 1 or 0, with the non-zero bin xk values indicating the locations of the ordered zi's, assuming no ties.) Familiar statistical applications, of the type described in Section 4, depend on z only through y.

Suppose that the zi's are divided into a finite number of classes, with members of the cth class Cc having mean μc and standard deviation σc,


Let Nc be the number of members of Cc, with pc the proportion


so Σc Nc = N and Σc pc = 1. The use of model (2.4) for z-values is supported by the results of Section 5.

If x is the K-vector of bin midpoints, let xkc = (xkμc)/σc and


Likewise, for any real-valued function h(x) we define hc to be the K-vector of function values


also denoted by h(xc) in what follows.

It is easy to calculate the expectation of the count vector y under the multi-class model (2.4)-(2.5). Let πkc equal the probability that zi from class Cc falls into the kth bin,


Here [var phi](x) = exp(−x2/2)/√2π, the standard normal density. The approximation πkc [equals, single dot above] Δ[var phi](xkc)/σc from (2.4) becomes arbitrarily accurate for Δ sufficiently small, and we will take it as exact in what follows. Then


The K × K covariance matrix of the count vector y depends on the N × N correlation matrix of z, but in a reasonably simple way discussed next. Two important definitions are needed to state the first result: there are M = N(N − 1)/2 correlations ρii between pairs (zi, zi) of members of z, and we denote by “g(ρ)” the distribution putting weight 1/M on each ρii. Also, for [var phi]ρ(u, v) the bivariate normal density having zero means, unit standard deviations, and correlation ρ, we define




(the integral notation being shorthand for summing over M discrete points).

Lemma 1

Under the multi-class model (2.4)-(2.5), the covariance of the count vector y (2.3) has two components,






Here diag(πc) is the K × K diagonal matrix having diagonal elements πkc, similarly diag(πd), while λcd is the K × K matrix with klth element λ(xkc, xld); the summations are over all classes.

Note. Equation (2.14) assumes that the correlation distribution g(ρ) is the same across all classes Cc. The proof of Lemma 1, which is similar to that for the simpler situation of Efron (2007a), appears in Remark C of Section 6.

The cov0 term in (2.12)-(2.13) is the sum of the multinomial covariance matrices that would apply if the zi's were mutually independent with fixed numbers drawn from each class; cov1 is a penalty for correlation, almost always increasing cov(y). The N2 factor in (2.14) makes the correlation penalty more severe as N increases, assuming g(ρ) stays the same.

Expression (2.14) for the correlation penalty can be considerably simplified. Mehler's identity for λρ(u, v) (2.10) is


where hj is the jth Hermite polynomial. (See Lancaster, 1958 for an enlightening discussion of (2.15), also known as the “tetrachoric series”, and its connections to the singular value decomposition, canonical correlation, Pearson's coefficient of contingency, and correspondence analysis.) Denoting the jth moment of the correlation distribution g(ρ) by αj,


(2.11) becomes


so λcd in (2.14) can be written in outer product notation as


Making use of (2.8), taken as exact,


where φc(j) indicates the jth derivative of [var phi](u) evaluated at each component of xc (using [var phi](j)(u) = (−1)j[var phi](u)hj(u)).

Rearranging (2.14) then gives a simplified formula.

Lemma 2.



(2.14) for the correlation penalty becomes


A convenient approximation to cov1 is based on three reductions of (2.21):

  • The second term in (2.21) is neglible for large N.
  • Common standardization methods for large-scale data sets often make α1, the expectation of g(ρ), exactly or nearly zero, as illustrated in Section 3 for the leukemia data; see Section 3 of Efron (2009).
  • This leaves α2 of (2.16) as the leading term in (2.21). With ρ confined to [−1, 1], the higher-order moments αj = Eg{ρj} often decrease quickly to zero.

The root mean square (rms) correlation


featured in Efron (2007a) (where it is called the total correlation), is a single-number summary of zi's entire correlation structure. Carrying out the three reductions above produces a greatly simplified form of (2.21),


with [phi with macron](2) in (2.20) depending on the second derivative of the normal density, [var phi](2)(u) = [var phi](u) · (u2 − 1).

Figure 2 compares the exact formulas (2.12)-(2.14) for cov(y) with the simplified formula based on the rms approximation (2.23); for a numerical example having N = 6000, α = .10, and two classes (2.4)-(2.5), initially with


but then recentered as in the leukemia example; see Remark D of Section 6 for more detail. The plotted curves show the standard deviations sd{yk}=covkk(y)12 from (2.12), the corresponding rms approximation (2.23), and also sd0{yk}=(cov0kk)12 from (2.13). We can see there is a substantial correlation penalty over most of the range of z, and also that the rms approximation is quite satisfactory here.

Figure 2
Comparison of exact formula for standard deviation of yk from (2.12) (heavy curve) with rms approxmation from (2.23) (dotted curve); N = 6000, α = .10 in (2.22), two classes as in (2.24). Dashed curve is standard deviation from (2.13) ignoring ...

Returning to right-sided cdfs (2.1), let B be the K × K matrix




is a K-vector with kth component the proportion of zi's in bins indexed ≥ k,


(B would be transposed if we were dealing with left-sided cdfs.) The expectation of F is both obvious and easy to obtain from (2.9),


where Φ+(u) = 1 − Φ(u). Now that we are working with tail areas rather than densities we can let Δ → 0, making (2.28) exact.

F has covariance matrix Bcov(y)B/N2. The same kind of calculations as in (2.28) applied to Lemma 1 gives the following theorem.

Theorem 1

Under the multiclass model (2.4)-(2.5),


where Cov0 has klth entry




Here pc is from (2.5), xkc and xlc from (2.6), αj is as in (2.16) and


(Notice the distinction between [phi] and [phi with macron] (2.20), and between Cov and cov etc., Lemma 1.)

The three-step reduction leading to (2.31) also can be applied to Cov1: for α as in (2.22),


with [phi](1) depending on the first derivative of the normal density, [var phi](1)(u) = −[var phi](u)u. Section 3 shows that (2.33) is especially convenient for applications.

Figure 3 is the version of Figure 2 applying to F: the heavy curve tracks sd(Fk) from (2.29), the dotted curve is from Rms approximation (2.33), and the dashed curve shows the standard deviations from Cov0 (2.30), ignoring the correlation penalty. Once again the simple approximation formula performs well, particularly for extreme values of z, which are likely to be the ones of interest in applications. The correlation penalty is more severe here than in Figure 2, especially in the tails.

Figure 3
Comparison of exact formula for sd{Fk} from Theorem 1 (heavy curve) with Rms approximation using (2.33) (dotted curve); same example as in Figure 2. Dashed curve shows standard deviation estimates ignoring the correlation penalty.

The Cov0 formula (2.30) is essentially the covariance function for a Brownian bridge. Results related to Theorem 1 can be found in the empirical process literature; see equation (2.2) of Csörgő and Mielniczuk (1996) for example, which applies to the “one-class” case when all the zi's are [mathematical script N](0, 1). Desai et al. (2009) extend the covariance calculations in Efron (2007a) to include skewness corrections.

3 Estimation of the correlation parameters

Application of Section 2's theory requires us to estimate several parameters: the rms correlation α (2.22), and the class components (pc, μc, σc) in (2.4)-(2.5) (though we will see that the latter task can be avoided under some assumptions). This section illustrates the estimation process in terms of the leukemia study of Section 1. X, the data matrix for the study, has N = 7128 rows, one for each gene, and n = 72 columns, one for each patient; the n1 = 47 ALL patients precede the n2 = 25 AML patients. Entry xij of X is the expression level for gene i on patient j. The columns of X were individually standardized to have mean 0 and variance 1; see Remark E.

The ith row of X gives ti, the two-sample t-statistic comparing expression levels on gene i for AML versus ALL patients. These are converted to z-values zi = Φ−1(F70(ti)) (1.2), whose histogram appears in Figure 1. As noted before, the histogram is much wider near its center than a theoretical [mathematical script N](0, 1) null distribution: analysis using the locfdr program described in Efron (2007b, 2008) estimated that proportion p0 = .93 of the genes were “null” (i.e., identically distributed for ALL and AML), and that z-values for the null genes followed a [mathematical script N](.09, 1.682) distribution.

We wish to estimate the rms correlation α (2.22). Let X0 indicate an N × n0 subset of X pertaining to a single population of subjects, for example the 47 ALL patients. There are N · (N − 1)/2 sample correlations [rho with circumflex]ii between rows i and i′ of X0. Computing all of these, or a sufficiently large random sample, yields the empirical mean and variance (m, v) of the [rho with circumflex] distribution,


(m, v) = (.002,.1902) for the ALL patients. As discussed in Section 3 of Efron (2009), standardizing the columns of X0 to have mean 0 forces m [equals, single dot above] 0, and we will assume m = 0 in what follows. (This is equivalent to taking α1 = 0 as we did following (2.21).)

The obvious choice α=v12 tends to greatly overestimate α : each [rho with circumflex]ii is nearly unbiased for its true correlation ρii, a normal-theory approximation for mean and variance being


(Johnson and Kotz, 1970), but the considerable variances in (3.2) can greatly broaden the empirical distribution of the [rho with circumflex]'s. Two corrected estimates of α are developed in Efron (2009). The simpler correction formula is


based on an identity between the row and column correlations of X0. The second approach uses an empirical Bayes analysis of the variance term in (3.3) to justify a more elaborate formula,


The first three columns of Table 2 compare [alpha] with [alpha] for X0 based on the ALL patients, the AML patients, and both. The final column reports mean ± standard deviation for [alpha] and [alpha] in 100 simulations of model (2.24): N = 6000, n1 = n2 = 40 patients in each class, true α = .10; see Remark D. The two estimates are effectively linear functions of each other for typical values of v; [alpha], the simpler choice, is preferred by the author.

Table 2
Estimates [alpha] and [alpha], (3.3) and (3.4), for rms correlation α (2.22) of leukemia data; also 100 simulations of model (2.24), N = 6000, n1 = n2 = 40, true α = .10, showing mean ± standard deviation. ...

It seems that we need to estimate the class components (pc, μc, σc) in (2.4)-(2.5) in order to apply the theory of Section 2, but under certain assumptions this can be finessed, as discussed next.

The marginal density f(z) under model (2.4)-(2.5) is


so, letting f = f(x) (the density evaluated at the K-vector of bin midpoints), we have Δ · f = Σcpcπc as in (2.8). Formula (2.13) can be expressed as


Here we are assuming, as in (2.5), that the class sample sizes Nc are fixed. A more realistic assumption might be that the numbers N1, N2, … ,NC are a multinomial sample of size N, sampled with probabilities p1, p2, … ,pC, in which case (3.6) becomes the usual multinomial covariance matrix


A smooth curve f fit to the histogram heights2 as in Figure 1 then yields cov^0 by substitution into (3.7), without requiring knowledge of the class structure (2.4)-(2.5). In the same way, we can estimate the Cov0 for F in (2.30) by the standard multinomial formula


Under some circumstances, a similar tactic can be applied to estimate the correlation penalties cov1 and Cov1, (2.23) and (2.33). The first and second derivatives f(1)(z) and f(2)(z) of (3.5) are


Suppose we make the homogeneity assumption that all σc values are the same, say σc = σ0. Comparison with definitions (2.32) and (2.20) then gives


with f(j) = (f(j)(xk))′. This leads to the convenient covariance penalty formulas,


from (2.33) and (2.23).

A smooth estimate f(z) of f(z) can be differentiated to give estimated values of Cov1 and cov1, for example,


for the correlation penalty standard deviation of F(xk) (2.27). (This provides the second term in formula (1.4).) The heavy curve in Figure 4 shows (3.12) for the leukemia data, using [sigma with hat]0 = 1.68, [alpha] = .114, and f(z) from Figure 1.

Figure 4
Leukemia data; two estimates of correlation penalty standard deviation sd1 {Fk} for Fk (2.27). Solid curve formula (3.12); dashed curve Rms approximation (2.33) using class estimates from Table 3. Dotted curve is independence ...

Suppose we are unwilling to make the homogeneity assumption. A straightforward approach to estimating Cov1 or cov1 requires assessments of the parameters (pc, μc, σc) in (2.4)-(2.5). These can be based on the “non-null counts” (Efron, 2007b), the small bars plotted negatively in Figure 1; see Remark B. The figure suggests three classes, left, center and right, with parameter values as estimated in Table 3.

Table 3
Three-class model (2.4)–(2.5) for leukemia data. Parameter estimates based on non-null counts, Remark B.

The dashed curve in Figure 4 shows sd1(Fk) estimated directly from (2.32)-(2.33) using the values in Table 3. It is similar to the homogeneity estimate (3.12) except in the extreme tails.

Formula (1.4) for the standard deviation of F(x) was tested in a simulation experiment. The specifications were the same as in Figure 3, with N = 6000, α = .10, and two classes of z-values (2.24). One hundred X matrices were generated as in the simulation for Table 2, each yielding a vector of 6000 correlated z-values, followed by [sigma with hat]0, [alpha] and f(1)(x) for use in (1.4); see Remark D for further details. Finally, sd^, the square root of (1.4), was calculated along with sd^0, the square root of just the first term.

The solid curve in Figure 5 shows the average of the sd^ values for x between −4 and 4.5, with solid bars indicating standard deviations of the 100 sd^s. There is a good match of the average with the exact sd curve from Figure 3. The error bars indicate moderate variability across the replications. The average for sd^0, dashed curve, agrees with the corresponding curve in Figure 3 and shows that correlation cannot be ignored in this situation.

Figure 5
Simulation experiment for formula (1.4). Solid curve average of sd^, square root of (1.4), 100 replications, with bars indicating standard deviation of sd^ at x = −4, −3, …, 4; dotted curve exact sd from Figure 3; dashed curve ...

4 Applications

Correlation usually degrades statistical accuracy, an important question for the data analyst being the severity of its effects on the estimates and tests at hand. The purpose of Section 2 and Section 3 was to develop practical methods for honestly assessing the accuracy of inferences made in the presence of large-scale correlation. This section presents a few examples of the methodology in action.

We have already seen one example: in Table 1 the accuracy of F(x), the right-sided empirical cdf for the leukemia data, computed from the usual binomial formula that assumes independence among the z-values,


was compared with sd^ from formula (1.4) in which the correlation penalty term was included: sd^ more than doubled sd^0 over most of the range.

Suppose we assume, as in Efron (2008), that each of the N cases (the N genes in the leukemia study) is either null or non-null with prior probability p0 or p1 = 1 − p0, and with the corresponding z-values having density either f0(z) or f1(z),

p0=Pr{null}f0(z)density if null,p1=Pr{non-null}f1(z)density if non-null.

Let F0 and F1 be the right-sided cdfs of f0 and f1, and F the mixture cdf


The probability of a case being null given that z exceeds x is


according to Bayes theorem, “Fdr” standing for false discovery rate.

If p0 and F0 are known then Fdr has the obvious estimate


(2.1). Benjamini and Hochberg's celebrated 1995 algorithm uses Fdr^(x) for simultaneous hypothesis testing, but it can also be thought of as an empirical Bayes estimator of the Bayesian probability Fdr(x). The bottom row of Table 1 shows Fdr^(x) for the leukemia data, taking p0 = .93 and F0 ~ N(.09, 1.682) as in Figure 1. (Later we will do a more ambitious calculation taking into account the estimation of p0 and F0.)

The coefficient of variation for Fdr^(x) approximately equals that for F(x) (when p0F0(x) is known in (4.5)). At x = 5 we have Fdr^(5)=.15, with coefficient of variation about .19. An Fdr^ of .15 might be considered small enough to trigger significance in the Benjamini–Hochberg algorithm, but in any case it seems clear that the probabiltiy of being null is quite low for the 71 genes having zi above 5. Even taking account of correlation effects, we have a rough upper confidence limit of .21 (i.e., .15 · (1 + 2 · .19)) for Fdr(5).

Next we consider accuracy estimates for a general class of statistics Q(y), where Q is a q-dimensional function of the count vector y (2.3). As in Section 5 of Efron (2007b), we assume that a small change dy in the count vector (considered as varying continuously) produces change dQ in Q according to


If cov^(y) is a covariance estimate for y, obtained perhaps as in (2.12), (3.8), or (3.11), then the usual delta-method estimate for cov(Q) is


In a theoretical context, where cov(y) is known, we might instead use


now with the derivative matrix D evaluated at the expectation of y.

Model (4.2) yields the local false discovery rate


f(x) being the mixture density


fdr(x) is inferentially more appropriate than Fdr(x) from a Bayesian point of view, but it is not as immediately available since it involves estimating the density f(x). However, because z-value densities are mixtures of near-normals as shown in Section 5, it is usually straightforward to carry out the estimation.

Locfdr, the algorithm discussed in Efron (2007a, 2008), estimates f(x) by means of Poisson regression of the counts yk as a spline function of the xk, the bin midpoints in (2.2)-(2.3). The structure matrix M for the Poisson regression is K × d, where K is the number of bins and d is degrees of freedom (e.g., the number of free parameters of the spline fit; see Remark G for details). Let f be the vector of fitted values f(xk), and [l with hat] the vector with components [l with hat]k = log(f(xk)). Then, as discussed in (Efron, 2007b, Sect. 5), (4.6) takes the form


and we can use (4.7) or (4.8) to approximate cov([l with hat]).

For any function v(x) define the vector


as with f and [l with hat] above. If




implying, if p0 and f0 are known, that


with D = M(M′ diag(NΔf)M)−1M(4.8).

The solid curves in Figure 6 plot standard deviations for log(fdr^(x)), obtained as square root of the diagonal elements of cov(lfdr^) (4.15), for model (2.24) with N = 6000 and rms correlation α equal 0, .1, or .2; see Remark D. The horizontal axes are plotted in terms of the upper percentiles of F(x), the right end of each plot corresponding to the far right tail of the z-value distribution. For α = 0, sd(logfdr^(x)) increases from .03 to .08 as we move from the fifth to the first percentile of F. The coefficient of variation (CV) of fdr^(x) nearly equals sd(logfdr^(x)), so fdr^(x) is quite accurately estimated for α = 0, but substantially less so for α = .2. Reducing N to 1500 doubles the standard deviation estimates for α = 0, but has less effect in the correlated situations: for α = .1 for example, the increase is only 20% at percentile .025. Simulations confirmed the correctness of these results.

Figure 6
Solid curves show standard deviation of log(fdr^(x)) as a function of x at the upper percentiles of the z-value distribution for model (2.24), N = 6000 and α = 0, .1, .2. Dotted curves (green) same for log(Fdr^(x)) (4.5), nonparametric Fdr estimator. ...

Intuitively it seems that fdr should be harder to estimate than Fdr, but that is not what Figure 6 shows. Let Lk = log(F(xk)), with corresponding vector L. Then D in (4.6) has


with B as in (2.25), giving an estimate of cov(L) from (4.7) or (4.8). The same argument as (4.13)-(4.15) shows that this also estimates cov(logFdr^), the log of vector (4.5), assuming p0F0(x) is known. The dotted curves in Figure 6 show standard deviations for log(Fdr^(x)). If anything, Figure 6 suggests that fdr^ is less variable than Fdr^, particularly at the smaller percentiles.

Here we are comparing the nonparametric estimator Fdr^(x) (4.5) with the parametric estimator fdr^(x). The Poisson spline estimate f(x) that gave fdr^(x) can be summed to give parametric estimates of F(x) and Fdr(x), say Fdr~(x). Straightforward calculations show that the derivative matrix D for Fdr~(x) is


with B from (2.25) and Df equaling D in (4.11). Standard deviations for log(Fdr~), shown by the dashed curves in Figure 6, indicate about the same accuracy for Fdr~(x) as for Fdr^.

All of these calculations assumed that p0 and f0(z) (or F0(z)) in (4.2) were known. This is unrealistic in situations like the leukemia study, where there is clear evidence that a textbook N(0, 1) theoretical null distribution is too narrow. Estimating an “empirical null” distribution, such as N(.09, 1.682) in Figure 1, is both necessary and feasible (see Efron, 2008) but can greatly increase variability, as discussed next.

Formula (4.14) becomes


when p0 and f0 are themselves estimated. The corresponding derivative matrix D^=dlfdr^dy in (4.6) appears as equation (5.8) in Efron (2007b), this formula applying to the central matching method for estimating p0f0(z). The second row of Table 4 shows {logfdr^(x)} obtained from Dcov(y)D′ for the same situation as in the middle panel of Figure 6. Comparison with the theoretical null standard deviations (from the solid curve in the middle panel) shows that estimating the null distribution greatly increases variability.

Table 4
Comparison of {log(fdr^(x))} using empirical null versus theoretical null for the situation in the middle panel of Figure 6. The empirical null standard deviations are much larger, as seen also in Efron (2007b).

Here are some points to note:

  • Accuracy is worse for log(Fdr^) then for log(fdr^) in the top line of Table 4.
  • Accuracy is somewhat better when p0f0(z) is estimated by the MLE option in locfdr (Lemma 2 of Efron, 2007b).
  • The big empirical null standard deviations in Table 4 are at least partially misleading: some of the variability in 1+O(n12) is “signal” rather than “noise”, tracking conditional changes in the appropriate value of fdr(x). See Figure 2 of Efron (2007a) and the discussion in that paper.

Remark H of Section 6 describes a parametric bootstrap resampling scheme that avoids the Taylor series computations of (4.7), but which has not yet been carefully investigated.

5 The non-null distribution of z-values

The results of the previous sections depend on the variates zi having normal distributions (1.5). By definition, a z-value is a statistic having a N(0, 1) distribution under a null hypothesis H0 of interest (1.1): but will it still be normal for non-null conditions? This section shows that under repeated sampling the non-null distribution of z will typically have mean O(1), standard deviation 1+O(n12), and non-normality Op(n−1) (as measured by the magnitude of skewness and kurtosis). In other words, normality degrades more slowly than unit standard deviation as we move away from the null hypothesis.

Figure 7 illustrates the phenomenon for the case of non-central t distributions,


the notation indicating a non-central t variable with ν degrees of freedom and non-centrality parameter δ (not δ2), as described in Chapter 31 of Johnson and Kotz (1970). Here, as in (1.2), Fν is the cdf of a central tν distribution. The standard deviation of z decreases as |δ| increases; for δ = 5, ν = 20, z has (mean, sd) equal (4.01, 0.71). The useful and perhaps surprising observation is that normality holds up quite well even far from the null case δ = 0. We tacitly used this fact to justify application of our theoretical results to the leukemia study.

Figure 7
Density of the z-value statistic (5.1) when t has a noncentral t distribution with ν = 20 degrees of freedom; for non-centrality parameter δ = 0, 1, 2, 3, 4, 5. The densities are seen to be nearly normal; dashed curves are exact normal ...

To begin the theoretical development, suppose that y1, y2, … , yn are independent and identically distributed (iid) observations sampled from Fθ, a member of a one-parameter family of distributions,


having its moment parameters {mean, standard deviation, skewness, kurtosis}, denoted


defined differentiably in θ. The results that follow are heuristic in the sense that they only demonstrate second-order Cornish–Fisher expansion properties, with no attempt to provide strict error bounds.

Under the null hypothesis H0 : θ = 0, which we can write as


the standardized variate




Normality can be improved to second order by means of a Cornish–Fisher transformation,


which reduces the skewness in (5.6) from O(n12) to O(n−1),


See Chapter 1 of Johnson and Kotz (1970) or, for much greater detail, Section 2.2 of Hall (1992). We can interpret (5.8) as saying that Z0 is a second-order z-value,


e.g., a test statistic giving standard normal p-values accurate to O(n−1).

Suppose now that H0 is false, and instead H1 is true, with y1, y2, … , yn iid according to


rather than (5.4). Setting


makes Z1 second-order normal under H1,


We wish to calculate the distribution of Z0 (5.7) under H1. Define


Some simple algebra yields the following relationship between Z0 and Z1.

Lemma 3

Under definitions (5.7), (5.11) and (5.13),




The asymptotic relationships claimed at the start of this section are easily derived from Lemma 3. We consider a sequence of alternatives θn approaching the null hypothesis value θ0 at rate n12,


The parameter d = √n(μθnμ0)/σ0 defined in (5.13) is then of order O(1), as is


while standard Taylor series calculations give


the dot indicating differentiation with respect to θ.

Theorem 2

Under model (5.2), (5.16), and the assumptions of Lemma 3,


with M and S as given in (5.17)-(5.18). Moreover,


Proof. The proof of Theorem 2 uses Lemma 3, with θn playing the role of H1 in (5.14). Both 1 − c2 and (γ1/γ0)Sc2 are of order O(n12); the former from (5.18) and the latter using γ1/γ0 = 1 + ([gamma with dot above]0/γ0)(θnθ0) + O(n−1). Since Y121 is Op(1), this makes the bracketed term in (5.14) Op(n12); multiplying by g0 = γ0/(6 √n) reduces it to Op(n−1), and (5.19) follows from (5.12). Differentiating M and S in (5.17)-(5.18) with respect to d verifies (5.20). ■

Theorem 2 supports our claim that, under non-null alternatives, the null hypothesis normality of Z0 degrades more slowly than its unit standard deviation, the comparison being Op(n−1) versus O(n12).

One-parameter exponential families are an important special case of (5.2). With θ the natural parameter of F and y its sufficient statistic, i.e., with densities proportional to exp{θy}g0(y), (5.20) reduces to


The parameter γ0=(6√n) is called the acceleration in Efron (1987), interpreted as “the rate of change of standard deviation with respect to expectation on the normalized scale,” which agrees with its role in (5.21).

As an example, suppose y1,y2,,yniidθΓ1, Γn indicating a standard gamma distribution with n degrees of freedom, so (y) = (1/θ)exp(y/θ) for y ≤ 0. (Equivalently, y ~ θΓn/n.) This is an exponential family having skewness γ0 = 2 for any choice of θ0. An exact z-value for testing H0 : θ = θo is


where Gn is the cdf of Γn. Table 5 shows the mean, standard deviation, skewness and kurtosis of Z0 for n = 10, θ0 = 1, evaluated for several choices of the alternative θ1. The standard deviation of Z0 increases steadily with θ1; here γ0/(6√n) = .1054, matching to better than three decimal places the observed numerical derivative dS/dM. Skewness and kurtosis are both very small; in the equivalent of Figure 7, there is no visible discrepancy at all between the density curves for Z0 and their matching normal equivalents.

Table 5
Gamma example, n = 10, θ0 = 1, indicating the distribution of z-value (5.22) for various non-null choices of θ. Standard deviation increases with θ1 in accordance with (5.21), while maintaining near-perfect normality for Z0.

So far we have considered z-values obtained from an average y of iid observations, but the results of Theorem 2 hold in greater generality. Section 5 of Efron (1987) considers one-parameter families where [theta w/ hat], an estimator of θ , has MLE-like asymptotic properties in terms of its bias, standard deviation, skewness and kurtosis,


Letting [theta w/ hat] play the role of y and μθ = θ + βθ/n in definitions (5.5)-(5.12), Lemma 3 and Theorem 2 remain true, assuming only the validity of the Cornish–Fisher transformations (5.9)-(5.12). Ignoring the bias βθ, i.e., taking Y0 = √n([theta w/ hat]θ0)/σ0 at (5.5), adds an O(n12) term to M in (5.17).

Moving beyond one-parameter families, suppose F is a p-parameter exponential family, having densities proportional to exp{η1x1+η2x2}g0(x1,x2), where η1 and x1 are real-valued while η2 and x2 are (p − 1)-dimensional vectors, but where we are only interested in η1, not the nuisance vector η2. The conditional distribution of x1 given x2 is then a one-parameter exponential family with natural parameter η1, which puts us back in the context of Theorem 2. Remark H of Section 6 suggests a further extension where the parameter of interest “ θ” can be a general real-valued function of η, not just a coordinate such as η1.

The non-central t family does not meet the conditions of Lemma 3 or Theorem 2: (5.1) is symmetric in δ around zero, causing γ0 in (5.14) to equal zero and likewise the derivative in (5.20). Nevertheless, as Figure 7 shows, it does exhibit impressive non-null normality. Table 6 displays the moment parameters of z = Φ−1(Fν(t)) (1.3), for t ~ tν(δ), ν = 20 and δ = 0, 1, 2, 3, 4, 5. The non-null normality isn't quite as good as in the gamma example of Table 5, but is still quite satisfactory for its application in Section 4.

Table 6
Non-central t example t ~ tν(δ) for ν = 20, δ = 0,1,2,3,4,5; moment parameters of z = Φ−1(Fν(t)) (1.3) indicate near-normality even for δ far from 0. (Moments calculated using (6.10).)

Microarray studies can be more elaborate than two-sample comparisons. Suppose that in addition to the N × n expression matrix X we have measured a primary response variable yj and covariates wj1, wj2, … , wjp on each of the n subjects. Given the observed expression levels xi1, xi2, … , xin for gene i, we could calculate ti, the usual t-value for yj as a function of xij, in a linear model that includes the p covariates. Then


is a z-value (1.1) under the usual Gaussian assumption, showing behavior like that in Table 6 for non-null genes.

6 Remarks

Some remarks, proofs, and details relating to the previous sections are presented here.

A. Poisson regression

The curve f(z) in Figure 1 is a Poisson regression fit to the counts yk, as a natural spline function of the bin centers xk. Here the xk ranged from −7.8 to 7.8 in steps of Δ = .2, while the spline had five degrees of freedom, so M in (4.11) was 79 × 6 (including the intercept column). See Section 5 of Efron (2007b).

B. Table 3

Section 3 of Efron (2008) defines the non-null counts yk(1)=(1fdr^(xk))yk. Since, under model (4.2), 1fdr^(xk) estimates the proportion of non-null z-values in bin k, yk(1) estimates the number of non-nulls. The yk(1) values are plotted below the x axis in Figure 1, determining the “left” and “right” distribution parameters in Table 3. “Center” was determined by the empirical null fit from locfdr, using the MLE method described in Section 4 of Efron (2007b). This method tends to underestimate the non-null counts near z = 0, and also the σc values for the left and right classes, but increasing then to 1.68 had little effect on the dashed curve in Figure 4.

C. Proof of Lemma 1

Let Ik(i) denote the indicator function of the event zi [set membership] Ƶk (2.2) so that the number of zi's from class Cc Ƶk is


the boldface subscript indicating summation over the members of Cc. We first compute E{ykcyld} for bins k and l, kl, and classes c and d,


following notation (2.5)-(2.10), with χij the indicator function of event i = j (which can only occur if c = d). This reduces to


under the assumption that the same correlation distribution g(ρ) applies across all class combinations. Since yk = Σc ykc (2.3), we obtain


the non-bold subscripts indicating summation over classes.



from (6.4) results, after some rearrangement, in


Using πkc = Δ · [var phi](xkc)/σc as in (2.8), expression (6.6) is seen to equal the klth element of cov(y) in Lemma 1, when kl.

The case k = l proceeds in the same way, the only difference being that NΔpcχcd[var phi](xkc)/σc must be added to formula (6.3). This adds NΔΣ cpc[var phi](xkc)/σc to (6.6), again in agreement with cov(y) in Lemma 1.

The assumption that g(ρ) is the same across all classes can be weakened for the rms approximations (2.23) and (2.33), where we only need the second moments α2 to be the same. In fact, the class structure can disappear entirely for rms formulas, as seen in (3.12).

D. Model (2.24)

Specifications (2.24) were recentered to give overall expectation 0 in (2.4), (2.5):


these being the parameter values used in Figures Figures2,2, ,3,3, ,55 and and6.6. Recentering overall expectations to zero is common in practice, a consequence of the data matrix X having its column-wise means subtracted off.

The 6000 × 80 data matrices X used in the simulations for Table 2 and Figure 5 had entries xij ~ [mathematical script N](δij, 1) independent across columns j: δij = 0 for j ≤ 40, while for columns j > 40,


z-values based on the difference of means between the last and first 40 “patients” then satisfy (2.4), (6.7). The correlation distribution g(ρ) was supported on two points, 20% ρ = .20 and 80% ρ = −.05, giving α = .10.

E. Leukemia data standardization

The original entries xij of the leukemia data matrix were genetic expression levels obtained using Affymetrix oligonucleotide microarrays. For the analyses here, each column of X was replaced by its normal score values xij = Φ−1((rij−.5)/7128), where rij was the rank of xij in its column. Transformations such as this reduce the disturbing effects of sensitivity differences between microarrays; see Bolstad, Irizarry, Astrand and Speed (2003).

F. A parametric bootstrap method

Section 3 of Efron (2007a) discusses a hierarchical Poisson simulation scheme that can be adapted to the more general context of this paper. Following the notation in (3.11), we first simulate a vector u,


and then take y ~ Poisson(u), that is ykindPoisson(uk). The simulated y vectors can then be used to assess the variability of any function Q(y), obviating the need for the derivative matrix D. This amounts to a parametric bootstrap approach to accuracy estimation. It produced similar answers to (4.11) when applied to the leukemia data but seemed prone to biases in other applications. (Nonparametric bootstrapping, resampling columns of the data matrix X, can produce erratic results for the kind of large-scale accuracy problems considered in Section 4.)

G. z-value densities

Suppose test statistic t has possible densities {fθ(t), θ [set membership] Θ}, with corresponding cdfs Fθ(t), and we wish to test H0 : θ = θ0. The z-value statistic z = Φ−1{Fθ0(t)} then has densities


The density curves in Figure 7 were obtained from (6.10), with fθ(t) the noncentral tν(θ) density, ν = 20 and θ = 0, 1, 2, 3, 4, 5.

H. Extensions of Theorem 2

In some circumstances, Theorem 2 can be extended to multi-parameter families F = {Fη} where we wish to test θ = θ0 for θ a real-valued function of η. This is straightforward to verify in the context of Efron (1985), which includes for example Fieller's problem, and is conjectured to be true in general exponential families.

7 Summary

The paper considers studies where a large number N of cases are under investigation, N perhaps in the hundreds or thousands, each represented by its own z-value zi, and where there is the possibility of substantial correlation among the zi's. Our main result is a simple approximation formula for the accuracy of summary statistics such as the empirical cdf of the z-values or an estimated false discovery rate. The argument proceeds in five steps:

  • Exact formulas for the accuracy of correlated z-value cdfs are derived under normal distribution assumptions (Section 2).
  • Simple approximations to the exact formulas are developed in terms of the root mean square correlation of all N · (N − 1)/2 cases (Section 2, (2.23) and (2.33)).
  • Practical estimates for the approximation formulas are derived and demonstrated through simulations and application to a microarray study (Section 3 and Section 4).
  • Delta-method arguments are used to extend the cdf results to more general summary statistics (Section 3 and Section 4).
  • Under reasonable assumptions, it is shown that z scores tend to have nearly normal distributions, even in non-null situations (Section 5), justifying application of the theory to studies in which the individual variates are z-values.

Our main conclusion is that by dealing with normal variates, a practical assessment of large-scale correlation effects on statistical estimates is possible.


This work was supported in part by NIH grant 8R01 EB002784 and NSF grant DMS0505673.


1It is convenient for the applications of Section 4 to deal with right-sided cdfs or survival curves instead of the usual left-sided ones in (1.2), and we will use this definition in what follows.

2The estimate used here is a Poisson spline regression as described following (4.10).


  • Bolstad B, Irizarry R, Astrand M, Speed T. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003;19:185–193. [PubMed]
  • Clarke S, Hall P. Robustness of multiple testing procedures against dependence. Ann. Statist. 2009;37:332–358.
  • Csörgő S, Mielniczuk J. The empirical process of a short-range dependent stationary sequence under Gaussian subordination. Probab. Theory Related Fields. 1996;104:15–25.
  • Desai K, Deller J, McCormick J. The distribution of number of false discoveries for highly correlated null hypotheses. Ann. Appl. Statist. 2009 Submitted, under review.
  • Dudoit S, Laan M. J. van der, Pollard KS. Multiple testing. I. Single-step procedures for control of general type I error rates. Stat. Appl. Genet. Mol. Biol. 2004;3:71. Art. 13. electronic. [PubMed]
  • Dudoit S, Shaffer JP, Boldrick JC. Multiple hypothesis testing in microarray experiments. Statist. Sci. 2003;18:71–103.
  • Efron B. Bootstrap confidence intervals for a class of parametric problems. Biometrika. 1985;72:45–58.
  • Efron B. Better bootstrap confidence intervals. J. Amer. Statist. Assoc. 1987;82:171–200. with comments and a rejoinder by the author.
  • Efron B. Correlation and large-scale simultaneous significance testing. J. Amer. Statist. Assoc. 2007a;102:93–103.
  • Efron B. Size, power and false discovery rates. Ann. Statist. 2007b;35:1351–1377.
  • Efron B. Microarrays, empirical Bayes and the two-groups model. Statist. Sci. 2008;23:1–22.
  • Efron B. Are a set of microarrays independent of each other? Ann. Appl. Statist. 2009 To appear. [PMC free article] [PubMed]
  • Golub TR, Slonim DK, Tamayo P, Huard C, Gaasenbeek M, Mesirov JP, Coller H, Loh ML, Downing JR, Caligiuri MA, Bloomfield CD, Lander ES. Molecular Classification of Cancer: Class Discovery and Class Prediction by Gene Expression Monitoring. Science. 1999;286:531–537. [PubMed]
  • Hall P. The bootstrap and Edgeworth expansion. Springer-Verlag; New York: 1992. (Springer Series in Statistics).
  • Johnson NL, Kotz S. Distributions in statistics. Continuous univariate distributions. 1. Houghton Mifflin Co; Boston, Mass: 1970.
  • Lancaster HO. The structure of bivariate distributions. Ann. Math. Statist. 1958;29:719–736.
  • Owen AB. Variance of the number of false discoveries. J. R. Stat. Soc. Ser. B Stat. Methodol. 2005;67:411–426.
  • Qiu X, Brooks A, Klebanov L, Yakovlev A. The effects of normalization on the correlation structure of microarray data. BMC Bioinformatics. 2005a;6:120. [PMC free article] [PubMed]
  • Qiu X, Klebanov L, Yakovlev A. Correlation between gene expression levels and limitations of the empirical Bayes methodology for finding differentially expressed genes. Stat. Appl. Genet. Mol. Biol. 2005b;4:32. Art. 34. electronic. [PubMed]
  • Schwartzman A, Lin X. The effect of correlation in false discovery rate estimation. Harvard University; 2009. (Biostatistics Working Paper Series number 106).
  • Westfall P, Young S. Resampling-based Multiple Testing: Examples and Methods for p-Value Adjustment. Wiley-Interscience; New York, NY: 1993.