Monte Carlo simulations were used to assess the FWE and power of the permutational multiple testing procedure, and compare it with use of the Bonferroni procedure. A genotype relative risk model was assumed at each SNP j,where ψ1 represents the risk of disease with one copy of the allele of interest divided by the risk of disease with no copies. Similarly, ψ2 represents the risk of disease with two copies of the allele of interest divided by the risk of disease with no copies. Correlated multivariate genotype data for triads was generated by first obtaining haplotype data in linkage disequilibrium (LD) for the parents and then applying an inheritance model.
To generate haplotype data in LD for the parents of cases, SNPs on the same strand were assumed to be in linked blocks of length nb
, with SNPs between blocks independent. For the purpose of the simulations the proportion of the allele of interest in the population, f
, is assumed to be the same for each SNP. Based on the probabilities of each mating type given a diseased case (given in of Schaid & Sommer, 1993
), the proportion of the allele of interest for parents of cases, p
, is calculated. A multivariate haplotype is formed one SNP at a time, starting by determining that the first SNP contains the allele of interest with probability p
. Consecutive SNPs in the same block are assumed to have LD parameter, D00
− p0· p·0
, where pij
represents (for i
= 0, 1, or ·) the proportion of haplotypes that have i
alleles of interest for the first SNP and j
alleles of interest for the second SNP and where if i
equals · then the proportion is evaluated for the corresponding SNP to have either allele. Subsequent SNPs in the block are then evaluated conditionally on the previous SNP, based on the probabilities
Each parental chromosome is generated independently.
Once haplotypes for the parents are generated, genotypes for the children are generated using a random crossover model. In this model, the child inherits from the same strand until a random crossover event occurs with probability pc, where inheritance is then from the other strand.
A total of 500 triads were simulated for each replication of the simulation experiment, with 100 SNPs in blocks of size nb = 5. The LD parameter was D00 = 0.5 for the parental haplotypes within a block. The LD for the case children was controlled by the value of pc.
shows the results of the null simulations for testing at level 0.05. Each simulation consisted of 100,000 replications. Both procedures control the FWE at the desired level of 5%. However, we note that the Bonferroni procedure becomes overly conservative when the correlation between SNPs within blocks is increased (correlation increases as pc decreases).
Familywise error (%) of the tests under complete null configurationsa
shows the results of simulations under three different nonnull models. The nonnull models correspond to dominant, recessive, and multiplicative disease inheritance patterns. Each simulation consisted of 10,000 replications. In each simulation, there were 5 nonnull SNPs out of 100. The average power over the nonnull SNPs is reported. In the cases of low correlation, there is very little difference in power between the methods. However, in the cases of high correlation, the power is substantially higher in the permutational method. This is in agreement with similar comparisons between the Bonferroni and two-sample permutational procedures for controlling the FWE.
Power (%) of the tests under alternative configurationsa
Neural Tube Defects and Clefts in Ireland
As part of a candidate gene study, we analyzed 1339 SNPs from 93 genes on 277 complete NTD triads from the Republic of Ireland. Birth defects are ideal candidates for triad studies because the parents are usually easily identified and likely to be willing to agree to participate. Using multiplicity correction is extremely harsh and the Bonferroni corrected p-values are all 1.0 after truncation, despite the smallest unadjusted p-value being 0.002233. The Bonferroni multiplier is 1339, so the Bonferroni adjusted p-value for the most significant SNP is not even close to being below 1.0 as 1339 × 0.002233 = 3.0. In contrast the permutational procedure gives adjusted p-values below 1.0 (smallest permutational adjusted p-value was 0.74), much smaller than the Bonferroni. Thus, although none of the adjusted p-values are close to being significant, it is clear that the adjusted p-values from the Bonferroni procedure are extremely conservative when compared to those from the permutational procedure.
As an example to see how the permutational adjustment compares to the Bonferroni when some of the adjusted p-values are relatively small, we present a subanalysis of the above experiment. This is presented as an example to examine the relative size of the adjusted p-values, and not to represent what we consider appropriate control for multiplicity. We consider 18 SNPs from a single gene on 277 complete NTD triads. The p-values adjusted only for the 18 SNPs are presented in . One sees again that the permutational procedure gives smaller adjusted p-values than the Bonferroni. Moreover, the improvement is quite large when considered as a proportion of the Bonferroni adjusted p-value. For the smallest unadjusted p-value, the permutational procedure yields an adjusted p-value more than 40% smaller than the corresponding Bonferroni adjusted p-value.
Unadjusted and adjusted p-values for SNPs from a single gene on NTD triads
A final example is given from an independent study. As part of an analysis of oral clefts in Ireland (Carter et al., 2010
), 31 SNPs on 250 complete cleft palate only case triads were analyzed. Again, this is presented as an example to compare the adjusted p-values of the procedures, and not to represent what we consider appropriate control for multiplicity, or to represent a complete analysis of cleft cases on these SNPs. The p-values corresponding to the 10 most significant SNPs, adjusted for all 31 SNPs are presented in . In this case, the smallest adjusted p-value of the permutational procedure is almost 50% smaller than the corresponding Bonferroni adjusted p-value.
Unadjusted and adjusted p-values for SNPs from a single gene on cleft triads