Code
set.seed(1)
m <- 1000
p <- c(runif(950), 2 * pnorm(-abs(rnorm(50, mean = 3.5)))) # 950 nulls, 50 non-nulls (z-tests, effect 3.5)
sum(p.adjust(p, "bonferroni") <= 0.05)[1] 9Code
sum(p.adjust(p, "holm") <= 0.05)[1] 9Multiple testing when the tests are not independent
FWER, FDR, what Benjamini-Hochberg assumes, and permutation calibration for networks.
Before reading: Individual p-values and Type I error from Uncertainty, confidence intervals, and tests, and the null hypotheses for the analysis. This page moves from one prespecified test to a family of claims.
1 Define the family of claims
An omics paper often makes a list of claims: these genes differ, these edges change, or these pathways are associated with an outcome. Even when each individual test is well calibrated, inspecting many tests creates many opportunities for false claims. The aim is to control an explicitly chosen error rate for the whole reporting procedure.
List the hypotheses before choosing an adjustment: all measured genes, a prespecified pathway, all candidate network edges, or several contrasts per feature. These are different families. For example, H_{0g}:\psi_g=0 might concern a standardized expression contrast for gene g. Testing two contrasts for each of 1,000 genes creates 2,000 hypotheses if all belong to the same reported family. Separating them into families needs a scientific reporting rationale decided before looking at the results.
The starting point is a valid marginal null p-value. For every true null i, validity means
\begin{aligned} \Pr(p_i\le t)&\le t,\\ 0&\le t\\ &\le1. \end{aligned}
In words, a null test rejects at threshold t no more than a fraction t of repeated datasets. Equality gives the familiar uniform null p-value; discrete tests can be conservative and still satisfy the inequality. “Marginal” refers to one test at a time. It does not mean that tests are independent, and it does not mean that the scientific target must be a marginal association.
Multiplicity adjustment starts after this calibration argument. It does not repair pseudoreplication, an invalid permutation scheme, or a test of the wrong estimand. Filtering or removing failed tests can also change the family. State which hypotheses were eligible and how uncomputable tests were handled; do not quietly redefine the family as the tests that produced attractive p-values.
2 The problem, sized for our work
A differential expression analysis tests 15,000 genes. A differential network analysis on 1,000 genes tests 1000(999)/2=499{,}500 unordered edges. Under a global null with exactly uniform marginal p-values, unadjusted 0.05 testing gives expected false-positive counts of 750 and 24,975, respectively.
To see why dependence does not change those expectations, let \mathcal H_0 index the m_0 true nulls. At a common unadjusted threshold \alpha,
\begin{aligned} V&=\sum_{i\in\mathcal H_0}1\{p_i\le\alpha\},\\ E[V]&=\sum_{i\in\mathcal H_0}\Pr(p_i\le\alpha)\\ &\le m_0\alpha. \end{aligned}
Linearity of expectation does not require independence. Equality holds for exact uniform nulls. Dependence changes how the false positives occur together, and therefore can greatly change the variability of a discovery list. Shared genes, pathways, people, and preprocessing can all induce it.
3 Counts, realized fractions, and error rates
Let D be the total number of discoveries and V the number of those discoveries whose null hypotheses are actually true. Both vary with the sampled dataset. Define
\begin{aligned} \mathrm{FDP}&=\frac{V}{\max(D,1)},\\ \mathrm{FDR}&=E[\mathrm{FDP}],\\ \mathrm{FWER}&=\Pr(V\ge1). \end{aligned}
The false discovery proportion (FDP) describes one list. It is zero by convention when nothing is rejected. In real data, we know D but generally do not know V, so we cannot inspect the list and calculate its true FDP. Simulation gives us that opportunity because we supply the truth.
The family-wise error rate (FWER) is the probability of making even one false rejection. It is useful when a single false claim is costly. The false discovery rate (FDR) averages the false fraction over repetitions of the entire procedure. It can suit an exploratory discovery list, but it allows some individual lists to have much larger false fractions than the target q. Neither criterion is automatically right for every confirmatory or exploratory analysis; choose from the consequences of the claims.
An expectation of fractions is different from a fraction of expected counts. Imagine two equally likely outcomes: a list with (D,V)=(1,1) and a list with (D,V)=(99,0). Their FDPs are 1 and 0, so FDR is 1/2. Pooling the lists gives one false claim in 100 discoveries, or 0.01. That pooled fraction weights large lists more heavily and is not FDR. In a simulation, average one FDP per dataset.
Under the global null every discovery is false. FDP is then 1 if there is any rejection and 0 otherwise, so FDR equals FWER. Away from that special case, the two criteria ask different questions.
4 FWER: Bonferroni and Holm
With m tests, Bonferroni rejects when p_i\le\alpha/m. The union bound explains its guarantee:
\begin{aligned} \Pr(V\ge1)&\le\sum_{i\in\mathcal H_0}\Pr(p_i\le\alpha/m)\\ &\le m_0\alpha/m\\ &\le\alpha. \end{aligned}
The probability of a union is at most the sum of its component probabilities, regardless of how the events overlap. That is why valid individual p-values suffice under any dependence. Holm also controls FWER under arbitrary dependence and can reject more hypotheses by progressively reducing the number still under consideration.
These FWER procedures remain valid with arbitrary dependence if the marginal null p-values are valid. They can be conservative when rejection events overlap, but positive correlations do not supply a universal effective-test count. With m=5\times10^5, the Bonferroni threshold is 10^{-7}; whether it has useful power depends on the effects and information available.
5 FDR: Benjamini-Hochberg
BH uses a rank-dependent threshold: a larger discovery set is allowed a larger raw-p-value cutoff. Its error criterion is FDR, so it generally rejects more than a procedure designed to protect against even one false claim.
The difference from Holm is essential: BH does not stop at the first failed comparison. It looks for the last passing rank and then includes the entire prefix of the sorted list.
5.1 Walk through one list
Take m=6 and q=0.05. The p-values below are already sorted.
| Rank i | p_{(i)} | iq/6 | Own comparison |
|---|---|---|---|
| 1 | 0.004 | 0.00833 | Pass |
| 2 | 0.018 | 0.01667 | Fail |
| 3 | 0.022 | 0.02500 | Pass |
| 4 | 0.070 | 0.03333 | Fail |
| 5 | 0.300 | 0.04167 | Fail |
| 6 | 0.800 | 0.05000 | Fail |
Rank 3 is the last passing rank, so BH rejects the first three hypotheses. Rank 2 is included because its p-value is below the selected cutoff 3q/6=0.025. Stopping at rank 2 would implement a different procedure.
Code
BH <- function(p, q = 0.05) {
m <- length(p); o <- order(p)
k <- max(c(0, which(p[o] <= q * seq_len(m) / m)))
rej <- logical(m); if (k > 0) rej[o[seq_len(k)]] <- TRUE
rej
}
sum(BH(p, 0.05)); sum(p.adjust(p, "BH") <= 0.05) # the same decisions[1] 28[1] 28The function implements the five steps for a complete vector of finite p-values. Package routines also offer conventions for missing p-values and a specified family size; choose those deliberately rather than letting dropped rows determine the scientific family.
5.2 What an adjusted p-value means
For a fixed family, the BH-adjusted p-value is the smallest target level q at which that hypothesis enters the BH list. In sorted order it is
\widetilde p_{(i)} =\min\left\{1,\ \min_{j\ge i}\frac{m\,p_{(j)}}{j}\right\}.
Start from the largest rank, take a running minimum of the scaled values mp_{(j)}/j, and cap at one. The running minimum enforces the step-up rule: an earlier rank can enter because a later rank passes. In the six-test example the adjusted values are 0.024, 0.044, 0.044, 0.105, 0.360, and 0.800. Rank 2’s own scaled value is 0.054, but rank 3 brings it down to 0.044. Comparing these adjusted values with 0.05 reproduces the three rejections.
An adjusted p-value of 0.044 is not a 4.4% probability that this gene’s null is true. It indexes a decision rule for the family. It also does not quantify effect size or supply a confidence interval with automatic simultaneous or post-selection coverage.
5.3 The dependence assumption and the BY alternative
A convenient sufficient condition for BH is independence of the entire p-value vector. Benjamini and Yekutieli extended control to positive regression dependence on the subset of true nulls (PRDS). This is a joint monotonicity condition: conditioning on a larger value of a null coordinate must not reduce the probability of any coordinatewise increasing event in the full vector. It constrains more than pairwise correlations. Certain one-sided Gaussian tests satisfy it; replacing each signed statistic by its absolute value for a two-sided test does not automatically preserve it. Positive correlations among measured genes therefore do not, by themselves, justify BH for arbitrary two-sided gene or edge tests. The BY paper gives the precise condition and sufficient examples.
When that justification is unavailable, Benjamini–Yekutieli (BY) supplies a general dependence correction.
For 1,000 tests, H_m\approx7.49: a BY target of 0.10 uses BH thresholds at an effective level of about 0.0134. This can cost substantial power. It buys a dependence guarantee; it does not improve invalid individual tests.
Under BH’s stated conditions the bound is q m_0/m, although m_0 is unknown. Adaptive approaches such as Storey’s q-values estimate \pi_0=m_0/m using the p-value distribution. Their estimation and dependence assumptions need their own justification. They are not interchangeable with ordinary BH-adjusted p-values, and an estimated null fraction is not permission to rescale any dependent family arbitrarily.
6 Dependence: what it does to BH
We will hold each test’s marginal distribution fixed and change only how the tests vary together. Each simulated dataset contains 1,000 test statistics. Tests 1–100 are nonnull with mean 3; tests 101–1,000 are null with mean zero. Every statistic has variance one. These are already-standardized test statistics, not expression measurements from 1,000 people.
In the independent world, every statistic gets its own independent standard-normal noise. In the dependent world, divide the tests into 20 blocks of 50: indices 1–50 form block 1, 51–100 block 2, and so on. The 100 signals occupy the first two complete blocks. A block is a set of related tests, such as tests sharing a source of variation; it is not a set of participants.
Write b(j) for the block containing test j. Generate
\begin{aligned} Z_j^{\rm ind}&=\mu_j+\varepsilon_j,\\ Z_j^{\rm block}&=\mu_j+\sqrt{\rho}\,F_{b(j)} +\sqrt{1-\rho}\,\varepsilon_j, \end{aligned}
where all block factors F_b and individual errors \varepsilon_j are independent standard normals. Each new dataset gets fresh factors and errors. Two tests in the same block share the same F_b, giving covariance \rho. Tests in different blocks share nothing, giving covariance zero. The total variance is \rho+(1-\rho)=1, so these covariances are also correlations. We use \rho=0.8.
Figure 1 shows the exact correlation structures used to generate the statistics. The signal locations and the dependence pattern are separate design choices.
In both worlds, a null statistic is exactly N(0,1), so 2\Phi(-|Z_j|) is a valid, uniform, two-sided null p-value. The nonnull marginal distributions also agree between worlds. In the block world a large positive or negative factor can make many null p-values small together. False discoveries can therefore arrive in clusters. This mechanism can change the spread of FDP even when its mean stays below q; a wide spread alone is not a failure of FDR control.
Code
set.seed(2)
m <- 1000; m1 <- 100; blocks <- 20; rho <- 0.8; n_sim <- 1000
mu <- c(rep(3, m1), rep(0, m - m1))
FDP <- function(z) { r <- BH(2 * pnorm(-abs(z)), q = 0.10); sum(r[-(1:m1)]) / max(sum(r), 1) }
fdp <- replicate(n_sim, {
z_ind <- rnorm(m) + mu
f <- rnorm(blocks)
z_dep <- sqrt(rho) * rep(f, each = m / blocks) + sqrt(1 - rho) * rnorm(m) + mu
c(independent = FDP(z_ind), block_dependent = FDP(z_dep))
})
rbind(mean_FDP = rowMeans(fdp), MCSE = apply(fdp, 1, sd) / sqrt(n_sim),
SD_FDP = apply(fdp, 1, sd), q95_FDP = apply(fdp, 1, quantile, 0.95)) independent block_dependent
mean_FDP 0.089981420 0.072636385
MCSE 0.001113368 0.004050366
SD_FDP 0.035207774 0.128083834
q95_FDP 0.149253731 0.341967495For each dataset, the code applies BH at q=0.10, counts rejections among the 900 known nulls, and divides by the total discoveries (or one if there are none). It repeats that complete analysis 1,000 times under each world. These are independently generated datasets, not 1,000 permutations of one dataset.
The mean FDP estimates FDR. Its Monte Carlo SE is the sample SD of the dataset-level FDPs divided by \sqrt{n_{\rm sim}}. A general FDP is a fraction, not a Bernoulli outcome, so a binomial MCSE is inappropriate. The SD and 95th percentile describe how much individual lists vary; they answer a different question from the mean. In the independent world the theoretical FDR bound is 0.10(900/1000)=0.09.
This finite experiment concerns the two specified two-sided Gaussian constructions. The one-sided Gaussian PRDS result does not itself establish the needed condition for the two-sided block construction, and one simulation cannot prove a general dependence guarantee. Compare the estimated mean with its MCSE before interpreting small differences from a theoretical level.
In the seeded run, estimated FDR is about 0.090 under independence (MCSE 0.0011) and 0.073 under block dependence (MCSE 0.0041). Yet the 95th percentile of FDP rises from about 0.15 to 0.34. Thus this example does not demonstrate BH exceeding its target FDR: it demonstrates that a smaller mean false fraction can coexist with much more erratic individual lists. A paper should state which of those properties it evaluated.
Code
plot(ecdf(fdp[1, ]), xlim = c(0, 1), verticals = TRUE, do.points = FALSE,
main = "FDP varies across discovery lists", xlab = "Realized FDP", ylab = "Fraction of datasets",
col = "#2F80C8", lwd = 2)
lines(ecdf(fdp[2, ]), verticals = TRUE, do.points = FALSE, col = "#D55E00", lwd = 2)
abline(v = 0.10, lty = 2, col = "gray50")
legend("bottomright", c("Independent", "Block dependent"),
col = c("#2F80C8", "#D55E00"), lty = 1, lwd = 2, bty = "n")
To read the empirical CDF, choose an FDP value on the horizontal axis. The curve’s height is the fraction of datasets whose realized FDP was at most that value. At the dashed line, one minus the height estimates how often a list exceeded 10% false discoveries. FDR is the mean of all the FDP values; it is neither that height nor the 95th percentile. This is why controlling FDR does not promise that 95% of lists have FDP below q.
Shared-edge tests also need a distributional model before a covariance formula can be used. For jointly Gaussian standardized features, if \rho_{jk}=\rho_{j\ell}=0, the asymptotic Fisher-z statistics for edges (j,k) and (j,\ell) have correlation \rho_{k\ell}. Gaussian fourth-moment factorization supplies that simplification. General non-Gaussian data require their actual fourth moments, and nonzero shared-edge correlations require the fuller covariance formula. Dependence alone neither proves nor disproves PRDS.
An end-to-end expression analysis uses BY for a prespecified feature family, distinguishes pointwise intervals from FDR control, and estimates FDR with MCSE across independent datasets.
7 Permutation calibration
When a valid permutation scheme reproduces the relevant joint null distribution, it can account for dependence without a parametric model for every pair of tests. The same relabeling must be applied to all features: it moves a person’s whole profile together and preserves the features’ relationship. Independently shuffling each gene destroys that relationship and supplies a different reference experiment.
7.1 Westfall–Young maxT for FWER
The question is now: how large could the largest absolute statistic in the family be under the null? Each permutation contributes one maximum, rather than m separate entries to a common marginal reference.
The comparison in step 4 is with the maximum, so these are already multiplicity-adjusted p-values; do not apply Bonferroni to them again as if they were marginal p-values. The plus-one term includes the observed arrangement’s contribution. If all allowed arrangements are enumerated instead, use the exact fraction of arrangements, including the observed one, that meet the comparison.
Subset pivotality is a standard sufficient condition for extending the global-null argument: the joint distribution of statistics for any set of true nulls is unchanged by whether the other hypotheses are null or nonnull. This is a property to justify, not a consequence of running permutations. In particular, a global-null simulation alone cannot establish strong FWER control. The gene-expression permutation analysis by Rempala and Yang discusses the role and limitations of this condition.
Code
set.seed(3)
# Two groups, 200 genes in 4 correlated blocks, no true differences. Compare Bonferroni and maxT.
n1 <- 15; n2 <- 15; G <- 200; blocks <- 4
MakeData <- function() {
f <- matrix(rnorm((n1 + n2) * blocks), ncol = blocks)
f[, rep(1:blocks, each = G / blocks)] * sqrt(0.7) + matrix(rnorm((n1 + n2) * G), ncol = G) * sqrt(0.3)
}
Tstat <- function(Y, g) {
a <- Y[g == 1, , drop = FALSE]; b <- Y[g == 2, , drop = FALSE]
(colMeans(a) - colMeans(b)) / sqrt(apply(a, 2, var) / n1 + apply(b, 2, var) / n2)
}
g <- rep(1:2, c(n1, n2)); Y <- MakeData(); t_obs <- Tstat(Y, g)
B <- 499
max_perm <- replicate(B, max(abs(Tstat(Y, sample(g)))))
p_maxT <- vapply(abs(t_obs), function(t) (1 + sum(max_perm >= t)) / (B + 1), numeric(1))
c(maxT_rejections = sum(p_maxT <= 0.05),
illustrative_maxT_cutoff = unname(quantile(max_perm, 0.95)),
bonferroni_cutoff = qt(1 - 0.05 / (2 * G), df = n1 + n2 - 2)) maxT_rejections illustrative_maxT_cutoff bonferroni_cutoff
0.000000 3.948490 4.192911 This is a global-null illustration with independent people and exchangeable group labels. The maxima retain the feature dependence. A lower observed cutoff than Bonferroni in this run is an illustration, not a fixed speed or power guarantee. Strong FWER claims under mixed nulls require conditions such as subset pivotality or another justified construction.
Code
hist(max_perm, breaks = 25, col = "#2F80C866", border = "white",
main = "A joint null reference", xlab = "Maximum absolute t-statistic")
abline(v = quantile(max_perm, 0.95), lty = 2)
7.2 A pooled permutation null for FDR
Some procedures instead pool permuted statistics to estimate the number of null tail exceedances at a candidate cutoff. This supplies an estimated false count, which must be related to the observed discoveries through a specified procedure. A common null scale across tests, treatment of nonnull tests, an appropriate null-fraction estimate, and the rule choosing a cutoff all matter. Heterogeneous sample sizes or null distributions can defeat naive pooling. Dividing pooled counts by permutation count and calling the result an FDR estimate does not by itself give FDR control.
SAM (Tusher et al. 2001) is a specific method with its own modified statistic, cutoff rule, and false-discovery estimation. The paragraph above is a description of a family of ideas, not a complete implementation of SAM. Follow the actual procedure and its assumptions when using it.
Both tools inherit the requirements of the Multivariate comparisons: permutation and distance-based tests page: permute the exchangeable units, and identify the null that justifies the actual relabeling. Equal means or correlations alone do not justify unrestricted label exchangeability.
8 Permutation resolution is a separate limitation
With B random permutations and the plus-one formula, the smallest possible p-value is 1/(B+1). That floor says nothing about whether the permutation scheme is valid, and reaching a cutoff does not imply precise tail estimation.
For marginal permutation p-values followed by Bonferroni, the floor must satisfy 1/(B+1)\le\alpha/m before any rejection is possible. Thus B+1\ge m/\alpha. With 1,000 tests and \alpha=0.05, this requires at least 19,999 permutations just to reach the cutoff.
For BH the relevant threshold is kq/m, where k is the last passing rank. At m=1000, q=0.05, and B=999, the minimum p-value is 0.001. One isolated p-value at that floor cannot pass rank 1, but 20 p-values at the floor pass at rank 20, since 20(0.05)/1000=0.001. The smallest BH cutoff is therefore a limitation for a sparse discovery set, not an unconditional requirement for every possible rejection. For 10^5 tests, reaching rank 1’s 5\times10^{-7} cutoff requires about two million permutations.
The maxT adjusted p-values also have the 1/(B+1) floor, but they are compared directly with \alpha, because the maximum already accounts for the family. Pooling tests or fitting a parametric tail changes the calibration argument; neither is an automatic cure for limited permutation resolution.
9 Some practical rules
- Decide FWER or FDR before you look, and write it in the protocol.
- Report the number of tests, the procedure, and the level. “Significant after adjustment” is not a method.
- When the tests are edges, state the dependence justification. BH requires independence or a suitable positive-regression-dependence condition, not merely some positive correlations. Benjamini–Yekutieli (
p.adjust(p, "BY")) controls FDR under arbitrary dependence with valid marginal p-values, at a possible power cost. Permutation methods need their own null and design justification. - Data-dependent filtering needs a null-independence argument or a procedure accounting for selection. Variance filtering is not automatically valid for every test and design; filtering on the same group contrast that is then tested generally changes the reference distribution.
- Report effect estimates and their scale alongside adjusted p-values. Testing a zero contrast does not test whether a change exceeds a biologically relevant size; that requires its own null. Ordinary pointwise intervals also do not acquire simultaneous or post-selection coverage from an FDR adjustment.
- Multiple testing is about false positives. Nothing here protects against the bias terms in the What is your estimator actually estimating? page: a procedure that is biased for every edge produces a list that is confidently, reproducibly wrong, and no adjustment catches it.
10 Exercises
Checkpoints are completion targets for this page. Additional practice below is optional. Each checkpoint gives a completion standard, a hint, and a solution outline.
Checkpoint · MT-01 — Specify the family and error guarantee
Task type: Interpretation and computation.
You will report discoveries from 1,000 prespecified gene tests, each with a valid marginal null p-value. Define FWER, realized FDP, and FDR. Give the Bonferroni threshold for FWER 0.05. Explain why an FDR guarantee of 0.05 does not certify that at most 5% of your particular discovery list is false, and state what must be justified before choosing BH under dependence.
Then execute BH by hand on the six p-values in the worked table at q=0.05. Identify the selected rank, every rejection, and why the failed comparison at rank 2 does not stop the algorithm.
Ready to continue when: you can state the family, distinguish an expectation from a realized fraction, execute the BH step-up rule, and identify the p-value and dependence assumptions.
Transfer: If you test two contrasts for every gene, how must the family be specified before adjusting p-values?
Hint
FDR averages FDP over repeated datasets; Bonferroni uses a union bound. For BH, scan every rank and take the last passing comparison.
Solution outline — open after trying
FWER is the probability of one or more false rejections. FDP is the false discoveries divided by the number of discoveries, with denominator at least one; FDR is its expectation. The Bonferroni threshold is 0.05/1000=0.00005, requiring valid marginal p-values but no dependence restriction. A realized list can have FDP above 0.05 even under FDR control. Standard BH needs independence or an appropriate positive-regression-dependence condition; arbitrary positive pairwise correlations alone do not establish it.
In the worked list, ranks 1 and 3 pass their own comparisons; rank 3 is the last to pass. Therefore k=3 and BH rejects the hypotheses with p-values 0.004, 0.018, and 0.022. Rank 2 enters as part of that prefix. Stopping at its failed comparison would confuse BH’s step-up rule with a step-down procedure.
Checkpoint · MT-02 — Separate permutation resolution from validity and precision
Task type: Computation and interpretation.
You have 1,000 prespecified tests and 999 valid random permutations per test, using (b+1)/(B+1). At Bonferroni FWER 0.05, can any test reject? How many permutations are needed merely to reach that smallest threshold? Could BH at FDR 0.05 reject if 20 tests all attain p=0.001? Finally, state the MCSE for mean FDP over independent simulated datasets and explain why it is generally not a binomial MCSE.
Ready to continue when: you can distinguish a valid Monte Carlo p-value, its resolution, tail precision, and the error-rate summary being estimated.
Transfer: Does increasing the number of features rescue an invalid permutation scheme, even if it produces many tail values?
Hint
Compare the minimum p-value with the relevant step-specific cutoff, then use variation of entire replicate FDPs.
Solution outline — open after trying
The minimum p-value is 0.001, above the Bonferroni threshold 0.00005, so Bonferroni cannot reject. At least 19,999 permutations are needed for the minimum grid point to equal that threshold; this is a resolution requirement, not a precision guarantee. BH can reject at rank 20 because 20(0.05)/1000=0.001, assuming its p-value/dependence conditions hold. With M independent simulation replicates, MCSE of mean FDP is sample SD of replicate FDP divided by sqrt(M). A general FDP is a fraction, not a Bernoulli outcome; under a global null it reduces to the indicator of any rejection.
Additional practice (optional)
These exercises extend the ideas through further derivations, simulations, or applications. Choose those useful for your work.
- Rerun the dependence experiment with
rho = 0.3andrho = 0.95. Plot the distribution of realized FDP in each case. At what dependence does the 95th percentile of FDP exceed 2q? - The original experiment puts all 100 signals in the first two blocks. Redistribute them as five signals in each of the 20 blocks, keeping the total signal count, effect size, and within-block correlation fixed. Compare mean power, mean FDP with MCSE, and the FDP distribution. Why can signal placement matter even with the same marginal distributions?
- Estimate pooled null tail counts in the
MakeData()setting with 20 shifted genes. List the additional ingredients and assumptions required before calling the resulting cutoff a valid SAM or FDR procedure; contrast it with BH/BY using justified p-values. - Under the null, simulate expression for 30 genes in 40 samples with a block-correlation structure, compute all 435 Fisher z edge statistics, and estimate the correlation between the statistics for two edges that share a gene and two that do not. Relate what you see to the correlation structure you simulated.
- Work through the six-test BH example by hand at q=0.04 and q=0.05, then recover the same decisions from the adjusted p-values. Explain why rank 2 can be rejected even when its own comparison fails. This extends MT-01’s guarantee calculation to execution of the algorithm.
11 Reading
Primary references include the BH paper, BY paper, SAM paper, and R p.adjust documentation.
- Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B, 57(1), 289 to 300.
- Benjamini, Y. and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics, 29(4), 1165 to 1188.
- Storey, J. D. and Tibshirani, R. (2003). Statistical significance for genomewide studies. PNAS, 100(16), 9440 to 9445.
- Westfall, P. H. and Young, S. S. (1993). Resampling-Based Multiple Testing. Wiley. Chapter 2.
- Tusher, V. G., Tibshirani, R., and Chu, G. (2001). Significance analysis of microarrays applied to the ionizing radiation response. PNAS, 98(9), 5116 to 5121.
- Goeman, J. J. and Solari, A. (2014). Multiple hypothesis testing in genomics. Statistics in Medicine, 33(11), 1946 to 1978. A review of genomic multiple-testing methods and their assumptions.
- Olkin, I. and Siotani, M. (1976). Asymptotic distribution of functions of a correlation matrix. In Essays in Probability and Statistics, 235 to 251. The covariance formula for sample correlations.
- Efron, B. (2010). Large-Scale Inference. Cambridge. Chapters 2, 4, and 8 for the empirical Bayes view and the effect of correlation.