I have p values from a lot of tests and would like to know whether there is actually something significant after correcting for multiple testing. The complication: my tests are not independent. The method I am thinking about (a variant of Fisher's Product Method, Zaykin et al., Genet Epidemiol, 2002) needs the correlation between the p values.

In order to estimate this correlation, I am currently thinking about bootstrapping cases, running the analyses and correlating the resulting vectors of p values. Does anyone have a better idea? Or even a better idea for my original problem (correcting for multiple testing in correlated tests)?

Background: I am logistically regressing whether or not my subjects are suffering from a particular disease on the interaction between their genotype (AA, Aa or aa) and a covariate. However, the genotype is actually a lot (30-250) of Single Nucleotide Polymorphisms (SNPs), which are certainly not independent but in Linkage Disequilibrium.


This is actually a hot topic in Genomewide analysis studies (GWAS)! I am not sure the method you are thinking of is the most appropriate in this context. Pooling of p-values was described by some authors, but in a different context (replication studies or meta-analysis, see e.g. (1) for a recent review). Combining SNP p-values by Fisher's method is generally used when one wants to derive an unique p-value for a given gene; this allows to work at the gene level, and reduce the amount of dimensionality of subsequent testing, but as you said the non-independence between markers (arising from spatial colocation or linkage disiquilibrium, LD) introduce a bias. More powerful alternatives rely on resampling procedures, for example the use of maxT statistics for combining p-value and working at the gene level or when one is interested in pathway-based approaches, see e.g. (2) (§2.4 p. 93 provides details on their approach).

My main concerns with bootstraping (with replacement) would be that you are introducing an artificial form of relatedness, or in other words you create virtual twins, hence altering Hardy-Weinberg equilibrium (but also minimum allele frequency and call rate). This would not be the case with a permutation approach where you permute individual labels and keep the genotyping data as is. Usually, the plink software can give you raw and permuted p-values, although it uses (by default) an adaptive testing strategy with a sliding window that allows to stop running all permutations (say 1000 per SNP) if it appears that the SNP under consideration is not "interesting"; it also has option for computing maxT, see the online help.

But given the low number of SNPs you are considering, I would suggest relying on FDR-based or maxT tests as implemented in the multtest R package (see mt.maxT), but the definitive guide to resampling strategies for genomic application is Multiple Testing Procedures with Applications to Genomics, from Dudoit & van der Laan (Springer, 2008). See also Andrea Foulkes's book on genetics with R, which is reviewed in the JSS. She has great material on multiple testing procedures.

Further Notes

Many authors have pointed to the fact that simple multiple testing correcting methods such as the Bonferroni or Sidak are too stringent for adjusting the results for the individual SNPs. Moreover, neither of these methods take into account the correlation that exists between SNPs due to LD which tags the genetic variation across gene regions. Other alternative have bee proposed, like a derivative of Holm's method for multiple comparison (3), Hidden Markov Model (4), conditional or positive FDR (5) or derivative thereof (6), to name a few. So-called gap statistics or sliding window have been proved successful in some case, but you'll find a good review in (7) and (8).

I've also heard of methods that make effective use of the haplotype structure or LD, e.g. (9), but I never used them. They seem, however, more related to estimating the correlation between markers, not p-value as you meant. But in fact, you might better think in terms of the dependency structure between successive test statistics, than between correlated p-values.


  1. Cantor, RM, Lange, K and Sinsheimer, JS. Prioritizing GWAS Results: A Review of Statistical Methods and Recommendations for Their Application. Am J Hum Genet. 2010 86(1): 6–22.
  2. Corley, RP, Zeiger, JS, Crowley, T et al. Association of candidate genes with antisocial drug dependence in adolescents. Drug and Alcohol Dependence 2008 96: 90–98.
  3. Dalmasso, C, Génin, E and Trégouet DA. A Weighted-Holm Procedure Accounting for Allele Frequencies in Genomewide Association Studies. Genetics 2008 180(1): 697–702.
  4. Wei, Z, Sun, W, Wang, K, and Hakonarson, H. Multiple Testing in Genome-Wide Association Studies via Hidden Markov Models. Bioinformatics 2009 25(21): 2802-2808.
  5. Broberg, P. A comparative review of estimates of the proportion unchanged genes and the false discovery rate. BMC Bioinformatics 2005 6: 199.
  6. Need, AC, Ge, D, Weale, ME, et a. A Genome-Wide Investigation of SNPs and CNVs in Schizophrenia. PLoS Genet. 2009 5(2): e1000373.
  7. Han, B, Kang, HM, and Eskin, E. Rapid and Accurate Multiple Testing Correction and Power Estimation for Millions of Correlated Markers. PLoS Genetics 2009
  8. Liang, Y and Kelemen, A. Statistical advances and challenges for analyzing correlated high dimensional snp data in genomic study for complex diseases. Statistics Surveys 2008 2 :43–60. -- the best recent review ever
  9. Nyholt, DR. A Simple Correction for Multiple Testing for Single-Nucleotide Polymorphisms in Linkage Disequilibrium with Each Other. Am J Hum Genet. 2004 74(4): 765–769.
  10. Nicodemus, KK, Liu, W, Chase, GA, Tsai, Y-Y, and Fallin, MD. Comparison of type I error for multiple test corrections in large single-nucleotide polymorphism studies using principal components versus haplotype blocking algorithms. BMC Genetics 2005; 6(Suppl 1): S78.
  11. Peng, Q, Zhao, J, and Xue, F. PCA-based bootstrap confidence interval tests for gene-disease association involving multiple SNPs. BMC Genetics 2010, 11:6
  12. Li, M, Romero, R, Fu, WJ, and Cui, Y (2010). Mapping Haplotype-haplotype Interactions with Adaptive LASSO. BMC Genetics 2010, 11:79 -- although not directly related to the question, it covers haplotype-based analysis/epistatic effect
  • 1
    $\begingroup$ Wow, thanks for going to all this trouble! I understand your qualms about bootstrapping, and I'm almost convinced. I think my main complication is the numerical covariate I have that will certainly be necessary (either by itself or in interaction with genotype), and that seems to rule out mt.maxT and plink, although I may need to look into plink again. But I will certainly dig through the references you provided! $\endgroup$ – Stephan Kolassa Sep 18 '10 at 15:56
  • $\begingroup$ You can always work with the residuals of your GLM to get ride of your covariates, though you lost some Df that will be difficult to account or reintroduce afterwards (e.g. for computing p-value). $\endgroup$ – chl Sep 18 '10 at 16:08
  • $\begingroup$ Hm, residuals from my logistic regression? Would that be legitimate? $\endgroup$ – Stephan Kolassa Sep 18 '10 at 18:09
  • $\begingroup$ Yes, why not? It is not uncommon to remove the variance accounted for by other covariates and then move on 2nd-level analysis with your residualized data. It is often faster (for instance, plink is pretty slow with categorical covariates, while it's ok with continuous ones; snpMatrix or simply glm() performs quite better on this point but you cannot embed a lot of SNPs within glm()...); the problem is that getting the corrected p-value at the end of your 2nd analysis is rather tricky (because you have to account for the parameters already estimated). $\endgroup$ – chl Sep 18 '10 at 18:27
  • $\begingroup$ For an illustration of how people are working with residuals, see for example p. 466 of Heck et al. Investigation of 17 candidate genes for personality traits confirms effects of the HTR2A gene on novelty seeking. Genes, brain, and behavior (2009) vol. 8 (4) pp. 464-72 $\endgroup$ – chl Sep 18 '10 at 18:31

Using a method like bonferroni is fine, the problem is that if you have many tests you are not likely to find many "discoveries".

You can go with the FDR approach for dependent tests (see here for details) the problem is that I am not sure if you can say upfront if your correlations are all positive ones.

In R you can do simple FDR with p.adjust. For more complex things I would take a look at multcomp, but I didn't go through it to see for solutions in cases of dependencies.

Good luck.

  • 1
    $\begingroup$ Hi Tal, thanks! Bonferroni does not seem appropriate to me - if one of my SNPs is causal and others are associated with it, there should be a signal, and Bonferroni has always looked too conservative to me (I usually prefer Holm's stepwise correction). The FDR you link to and p.adjust do not consider combined evidence (and the FDR still requires me to understand the correlation of my tests, the original question). multcomp may help, though at first glance it seems like it deals more with multiple tests within a single model, whereas I have multiple models. I'll dig deeper... $\endgroup$ – Stephan Kolassa Sep 18 '10 at 14:12
  • $\begingroup$ Hello Stephan. I understand, sorry for not helping more. Good luck! Tal $\endgroup$ – Tal Galili Sep 18 '10 at 20:14
  • $\begingroup$ Hello Stephan, I still think you can still use the method=BY (for Benjamini Hochberg Yekuteli Procedure) in p.adjust in R, as pointed by Tal. Definitely, using Bonferroni can be conservative. $\endgroup$ – suncoolsu Sep 22 '10 at 3:44
  • $\begingroup$ suncoolsu, I think that this method only works when the correlation is positive (not negative) between the variables. Cheers. $\endgroup$ – Tal Galili Sep 22 '10 at 13:57

I think Multivariate Normal Models are being used to model the correlated p-values and to get the right type of multiple testing corrections. Rapid and Accurate Multiple Testing Correction and Power Estimation for Millions of Correlated Markers. PLoS Genet 2009 talks about them and also gives other references. It sounds similar to what you were talking about, but I think other than getting a more accurate global p-value correction, LD structure knowledge should also be use to remove spurious positives arising from markers correlated with causal markers.


I am looking for a working solution for exactly the same problem. The best I found is the Null Unrestricted Bootstrap introduced by Foulkes Andrea in his book Applied Statistical Genetics with R(2009). Contrary to all bunch of other articles and books he considers specifically the regressions. Besides other methods he advises the Null Unrestricted Bootstrap, which is suitable where one cannot easily compute residuals (as in my case, where I model many independent regressions (basically simple correlations), each with the same response variable and different snip). I found this method to be also called the maxT method.

> attach(fms)
> Actn3Bin <- > data.frame(actn3_r577x!="TT",actn3_rs540874!="AA",actn3_rs1815739!="TT",actn3_1671064!="GG")
> Mod <- summary(lm(NDRM.CH~.,data=Actn3Bin))
> CoefObs <- as.vector(Mod$coefficients[-1,1]) 
> B <-1000
> TestStatBoot <- matrix(nrow=B,ncol=NSnps)
> for (i in 1:B){
+    SampID <- sample(1:Nobs,size=Nobs, replace=T)
+    Ynew <- NDRM.CH[!MissDat][SampID]
+    Xnew <- Actn3BinC[SampID,]
+    CoefBoot <- summary(lm(Ynew~.,data=Xnew))$coefficients[-1,1]
+    SEBoot <- summary(lm(Ynew~.,data=Xnew))$coefficients[-1,2]
+    if (length(CoefBoot)==length(CoefObs)){
+       TestStatBoot[i,] <- (CoefBoot-CoefObs)/SEBoot
+    }
+ }

Once we have the all the TestStatBoot matrix (in rows we have bootstrap replications, and in columns we have bootstrapped $\hat{\vec{T}^*}$ statistics) we find for which $T_{\text{crit.}}$ we observe exactly $\alpha=0.05$ percent of more significant $\hat{\vec{T}^*}$ statistics (more significant means that with bigger absolute value than $T_{\text{crit.}}$).

We report $i$-th model component significant, if its $\hat{\vec{T}_i} > T_{\text{crit.}}$

The last step can be accomplished with this code

p.value<-0.05 # The target alpha threshold
library(gtools) # for binsearch

p.level<-(1-pnorm(q=ans$where[[1]]/digits))*2 #two-sided.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.