A novel characterization of the generalized family wise error rate using empirical null distributions

Monotonicity of the kFWER estimators

The proposed kFWER estimator is defined as a function of a set Z which is usually an interval. For the left hand side (LHS) kFWER, for a point z, we define Z=(−∞, z), while for the right hand side (RHS) kFWER, we define Z=(z, ∞). In short, the parametric and empirical (LHS or RHS) kFWER estimates are monotone. To see this we can re-examine either equations (21) and (22). Without loss of generality, consider (21) restated here:

(35)kFWER¯(Z)=∑η=kN(1−FB(k−1|η, F0¯(Z))) fB(η|N, π0¯). (35)

The only term in (35) involving z or LHS version Z=(−∞, z) is (1−FB(k−1|η, F0¯(Z))), which is monotone (non-decreasing) as a function of z, hence kFWER¯(Z) is monotone (non-decreasing) as a function of z. A similar argument will show that the (RHS) kFWER is monotone (non-increasing) as a function of z.

The bias of kFWER¯(Z)

Under certain constraints (zero assumption and large π₀) we can show E(kFWER¯(Z))=kFWER(Z). To see this, we examine the convexity of kFWER¯(Z) as a function of π₀. We assume that F₀ is given and so π₀ is the only parameter to estimate in kFWER (Z). To emphasize this dependence on π₀ we will denote our estimate kFWER¯(Z) as kFWER¯π0(Z). According to (21) in our manuscript we have

(36)kFWER¯π0(Z)=∑η=kN(1−FB(k−1|η, F¯0(Z))) fB(η|N, π0). (36)

Figure 8 shows (36) as a function of π₀ for certain fixed values of k, N, and F¯0(Z). As shown kFWER¯ appears nearly linear (and constant) for large values of N and π₀. This can be further confirmed by examining the second derivative as a function of π₀. Computing the second derivative of kFWER π0(Z) as a function of π₀ we have

Figure 8

kFWER¯π0(Z) for various values of N, k, and F0(Z). For large N, the function appears linear and nearly constant for large values of π₀.

(37)kFWER¯π0′​′(Z)=∑η=kN[(Nη)(1−FB(k−1|η, F0(Z)))]((η−1)ηπ0η−2(1−π0)N−η−2(N−η)(1−π0)N−η−1ηπ0η−1+π0η(N−η)(N−η−1)(1−π0)N−η−2) (37)

The second derivative in (37) also depends on k, N, and F0(Z). In Figure 9 we explore the second derivative as a function of π₀ and these other variables. We see that for large values of N and π₀ the second derivative is near 0 confirming that kFWER is nearly linear. If we assume linearity for large π₀, by Jensen’s inequality we have that,

Figure 9

Second Derivative of kFWER¯π0(Z): for various values of N, k, and F0(Z). The second derivative is nearly zero for large values of N and k suggesting linearity of the kFWER¯ estimator.

(38)E(kFWER¯π0¯(Z))=k FWER E(π0¯)(Z) (38)

(39)=k FWER π0(Z) (39)

(40)=k FWER (Z) (40)

Note, we assume E(π0¯)=π0. Since F₀ is known, we have by definition π0¯=N+(A0)/(N⋅F0(A0)) and so

(41)E(π0¯)=E(N+(A0)N⋅F0(A0)) (41)

(42)=1N⋅F0(A0)E(N+(A0)) (42)

(43)=1N⋅F0(A0)π0N⋅F0(A0)  (by zero assumption)  (43)

(44)=π0 (44)

We can also examine kFWER^(Z) with the same assumptions. That is, we assume the following:

• each z_i follows the two group model independently,

• the zero assumption where f₁(z) is near zero for a subset of the sample space near zero, say A0,

• and that F₀ is the standard normal CDF.

Since we are using the nonparametric form kFWER^(Z), we assume that we need to estimate F₀ (non parametrically). We use 25,000 Monte Carlo simulations to study the performance of the MLE method to estimate F₀ and subsequently π₀. The bias of kFWER^π^0,F0^ is shown in Figure 10(A)–(C). Note the bias depends on the direction of significance for the non-null genes. This is reasonable since the direction of the non-null genes will skew F0^ in that direction. That is, if the non-null genes have positive z-values, F0^ will have a positive mean and thus the RHS version of kFWER^ will be an overestimate while the LHS version of kFWER^ will be an underestimate (Figure 10(A)). The opposite will occur with non-null genes with negative z-values (Figure 10(B)). With an equal balance of negative and positive z-scores, F0^ will be similar to a N(0, 1) CDF and thus kFWER^ will be similar kFWER¯, that is, a somewhat conservative estimate of kFWER (Figure 10(C)).

Figure 10

Bias: The bias in kFWER^(Z) when assuming k=1, N=500 and π₀=0.95. The dotted line is kFWER^(Z) with the LHS version for x<0 and the RHS version for x≥0 with the non null genes possessing large positive z-scores (a), large negative z-scores (b) and a (roughly) equal mixture of large positive and negative z-scores (c). The black line is the true kFWER. The bias of kFWER^(Z) depends on the direction of the non-null z-scores with kFWER^(Z) being conservative when there is roughly an equal number of non-null genes in each direction.

Sensitivity of estimators to A0 and π₀

For our simulations we used the program locfdr default choice of A0 which centers A0 at the median of z₁, z_s, …, z_N, with half-width about 2 times a preliminary estimate of σ₀ based on the interquartile range. For our simulations this results in A0≈(−2, 2). In Figure 11 we provide a simulation showing the sensitivity of our methods to the choice of A0 and π₀. In Table 2 we provide a simulation showing the sensitivity of our parameter estimates to the choice of π₀ and N. In short, the empirical kFWER estimator is generally robust to the choice of A0 but does not estimate kFWER well when π₀ and N are small.

Figure 11

A₀ and π₀ sensitivity: The dotted line is the mean kFWER^(Z) with Monte Carlo estimated 95 percent confidence interval as the shaded region. The mean kFWER¯(Z) is shown as the dashed line with the true kFWER (Z) in solid black. In this simulation we examine the LHS version of our estimators with N=500 and the non-null genes with large negative z-scores. In general, the estimators are robust to the choice of A0 but the accuracy of kFWER^(Z) depends greatly on the choice of π₀ with large values of π₀ required for accurate F₀ estimation. The kFWER¯(Z) is an accurate estimator of the true kFWER (Z) under all settings for π₀ and A₀ considered in this simulation.

In the following subsections we present results related to the adjusted Bonferroni, Šidàk, and Holm procedures. Additionally we present several methods to empirically estimate null distributions.

Adjusted Bonferroni method

Theorem 1Using the two sample mixture model defined in (1), zbon=F0−1(kα/n)is such that kFWER(z_bon)≤α when using the (LHS) kFWER definition wherek FWER (z)=kFWER(Z)whereZ=(−∞, z).

Proof. Let Z=(−∞, zbon). Then,

 (LHS)  k FWER (zbon)=Pr(N0(Z)≥k)≤E(N0(Z))k  (by Markov Inequality) =E(E(N0(Z)|N0))/k  (Law of total expectation) =E(N0F0(zbon))/k  (since N0(Z) |N0~Bin(N0, F0(zbon))) =Nπ0F0(F0−1(kα/N))/k  (since N0~Bin(N, π0)) =π0α≤α.

□

Corollary 9.1Usingzbon=F0−1(1−kα/n)and the (RHS) kFWER definition, we show that (RHS) kFWER(z_bon)≤α where (RHS)k FWER (z)=k FWER (Z) whereZ=(z, ∞).

Proof. A similar argument to Theorem 1 establishes the (RHS) kFWER result.□

Šidàk method

Theorem 2Consider the two sample mixture model defined in (1), zsid=F0−1(psid)where p_sid is such that

(45)FB(k−1|N, psid)=1−α. (45)

Then (LHS)kFWER(z_sid)≤α when using the (LHS) kFWER definition where k FWER (z)=k FWER (Z) where Z=(−∞, z).

Proof. We can assume that p_sid=F₀(z_sid). Then we have

(46)(LHS) kFWER(zsid)=∑i=kN(1−FB(k−1|i, psid)) fB(i|N, π0). (46)

Note from (45) we can assume that F_B(k–1|N, p_sid)=1–α and so 1–F_B(k–1|N, p_sid)=α. Hence the last term of the sum in (46) is απ0N. Importantly, we further note that F_B(k–1|n, p)>F_B(k–1|N, p) for any n<N. Thus,

(47)1−FB(k−1, N, p)>1−FB(k−1, n, p). (47)

Thus, we can rewrite (46) with the understanding that (LHS) k FWER (zsid)=k FWER (Z) with Z=(−∞, zsid) as follows,

(48)(LHS) k FWER (zsid)=∑i=kN(1−FB(k−1|i,psid))⋅fB(i|N, π0)≤α[fB(k|N, π0))+fB(k+1|N, π0)+…+   fB(N|N, π0)]≤α(1−FB(k−1|N, π0))≤α   (since 1−FB(k−1|N, π0) ≤ 1) . (48)

□

Corollary 9.2Similarly, if we definezsid=F0−1(1−psid) where p_sid is such that

(49)FB(k−1|N, psid)=1−α. (49)

Then (RHS)kFWER(z_sid)≤α when using the (RHS) kFWER definition where (RHS) k FWER (z)=k FWER (Z) where Z=(z, ∞).

Proof. Similar to the proof for Theorem 2. □

Holm method

Theorem 3Consider the two sample mixture model defined in (1) and zholm=F0−1(p(r)) where r is the largest index satisfying (28). Then (LHS) kFWER(z_holm)≤α when using the (LHS) kFWER definition with one sided p-values defined by p_i=F₀(z_i).

Proof. If r≤k, then our α adjustment is the same as the Bonferroni α adjustment and we can use the Markov inequality as employed in the proof for the adjusted Bonferroni argument.

Assume r>k, then we use the technique employed in Lehmann and Romano (2005). Let y1, y2, …, yN0 denote the ordered z statistics for the true null hypotheses where y1≤y2≤y3…≤yN0. Then let z_j=y_k where z₁≤z₂≤…≤z_N denote the ordered z statistics. Thus, the following probability statements hold,

(50)(LHS) kFWER(zholm)=Pr({N0(−∞, zholm)≥k})=Pr(#  of null z-values  ∈(−∞, zholm)≥k). (50)

The event {# of null z-values ∈(–∞, z_holm)≥k} is equal to the event {y_k=z_j≤z_holm}. In order to reject at least k true nulls, the largest possible value of j is N–N₀+k, namely, the situation where the N–N₀ true alternatives are the smallest z statistics. Hence, we have that yk=zj≤zN−N0+k. Hence, we have Pr[#  of null z-values ∈(−∞, zholm)≥k)=Pr({yk=zj≤zN−N0+k}]. Now apply F₀ to the event {yk=zj≤zN−N0+k} which is a non decreasing function in order to obtain, F0(yk)=F0(zj)=pj≤F0(zN−N0+k)=pN−N0+k. However, since z_j≤z_holm we must have that p_j≤α_j. Also, αj≤αN−N0+k=kα/N0 since α is an increasing function, see (29). Hence,

(51)(LHS) k FWER (zholm)=Pr(#  of null z-values ∈(−∞,zholm)≥k)=Pr({yk=zj≤zN−N0+k})=Pr({F0(yk)=F0(zj)≤pN−N0+k})≤Pr(F0(yk)=F0(zj)=pj≤αN−N0+k=kα/N0.=Pr(Uk≤kαN0)  where Uk ≡ F0(yk); the k-th null p-value =Pr(W>k)  where W ~ Bin(N0, kα/N0)≤E(W)/K   By Markov inequality =α (51)

Corollary 9.3Consider the two sample mixture model defined in (1) and zholm=1−F0−1(p(r)) where r is the largest index satisfying (28). Then kFWER(z_holm)≤α when using the (RHS) kFWER definition where k FWER (z)=k FWER (Z) where Z=(z, ∞).

Proof. Similar to the proof for Theorem 3.□

Theorem 4We assume that each z_i follows the two group model independently and that F₀is given as the standard normal CDF. We also assume the zero assumption where f₁(z) is near zero for a subset of the sample space near zero, say, A0. Then E(kFWER¯(Z))≥k FWER (Z).

Proof. Since F₀ is given, π₀ is the only parameter to estimate in kFWER. To emphasize the dependence on π₀, we will denote our estimator kFWER¯(Z) as kFWER¯π0¯(Z). According to (21) we have

(52)kFWER¯π0¯(Z)=∑η=kN(1−FB(k−1|η, F0¯(Z))) fB(η|N, π0¯) (52)

(53)=∑η=→kNC⋅fB(η|N, π0¯) (53)

(54)=∑η=kNC⋅(Nη)ηπ0¯(N−η)1−π0¯ (54)

(55)=∑η=kNC′ηπ0¯(N−η)1−π0¯, (55)

where C and C′ do not depend on π₀. By straightforward calculus the function in (55) is convex in terms of π0¯. Thus by Jensen’s inequality we have

(56)E(kFWER¯π0¯(Z))≥kFWERE(π0¯)(Z) (56)

(57)=kFWERπ0(Z) (57)

(58)=kFWER(Z) (58)

Note, in going from (6) to (7) we assume E(π0¯)=π0. Since F₀ is known, we have by definition π0¯=N+(A0)/(N⋅F0(A0)) and so

(59)E(π0¯)=E(N+(A0)N⋅F0(A0)) (59)

(60)=1N⋅F0(A0)E(N+(A0)) (60)

(61)=1N⋅F0(A0)π0N⋅F0(A0)  (by zero assumption)  (61)

(62)=π0 (62)

Methods to estimate the null distribution

In this section, we paraphrase two methods described in Efron (2010) for estimating the null distribution. We assume that f₀(z) is normal but not necessarily N(0, 1) say,

(63)f0(z)~N(δ0, σ02), (63)

and we define f_π0(z)₌π₀f₀(z). This implies that

(64)log(fπ0(z))=[log(π0)−12{δ02σ02+log(2πσ02)}]+δ0σ02z−12σ02z2, (64)

is a quadratic function of z.

The MLE method for empirical null distribution

This method was first introduced in Efron (2007). The maximum likelihood estimator (MLE) method starts with the zero assumption, where we assume that f₁(z) is zero for a certain subset A0 of the sample space. In other words,

(65)f1(z)=0  for  z∈A0. (65)

We assume N₀ is the number of z_i in A0 and ℐ0 their indices, ℐ0={i:zi∈A0} and N0=#ℐ0. We define z₀ as the corresponding collection of z-values,

(66)z0={zi, i∈ℐ0}. (66)

Also, let φδ0,σ0(z) be the N(δ0, σ02) density function,

(67)φδ0,σ0(z)=12πσ02 exp  {−12(z−δ0σ0)2} (67)

and

(68)H0(δ0, σ0)≡∫A0φδ0,σ0(z)dz, (68)

this being the probability that a N(δ0, σ02) variate falls in A0.

We suppose that the N z_i values follow the two-group model (1) with f0~N(δ0, σ02) and f₁(z)=0 for z∈A0. Then z_o has density and likelihood function

(69)fδ0,σ0,π0(z0)=[(NN0)θN0(1−θ)N−N0][∏ℐ0φδ0,σ0(zi)H0(δ0, σo)] (69)

when θ=π0H0(δ0, σ0)=Pr({zi∈A0}).

Computations can produce maximum likelihood estimators (δ^0, σ^0, π^0);fδ0,σ0,π0(z0) is the product of two exponential families which can be solved separately (the two bracketed terms). The binomial term gives θ^=N0/N while δ^0 and σ^0 are the MLEs from a truncated normal family, obtained by familiar iterative calculations, finally yielding

(70)π^0=θ^/Ho(δ^0, σ^0). (70)

The log of (69) is concave in (δ₀, σ₀, π₀) guaranteeing that the MLE solutions are unique. This is described more fully in Section 6.3 of Efron (2010).

The central matching method for empirical null distribution

This method was first introduced in Efron (2004). In this method, we define y_k as the number of observations z_i in the kth bin,

(71)yk=#{zi∈Zk}, (71)

where we partition the range Z of z_i values into K bins of equal width d with

(72)Z=∪k=1K Zk. (72)

Then, with the central matching method, we estimate f₀(z) and π₀ by assuming that log(f(z)) is quadratic near 0 and equal to (64) with,

(73)log(f(z))≈β0+β1z+β2z2. (73)

Estimating (β₀, β₁, β₂) can be done using least squares with the histogram counts y_k around z=0 and matching coefficients between (64) and (73). In other words, via matching, we obtain,

(74)σ02=−1/(2β2), (74)

(75)δ0=−β1/(2β2), (75)

(76)logπ0=β0−β124β2+log(−π/β2). (76)