Generalized shrinkage Flike statistics for testing an interaction term in gene expression analysis in the presence of heteroscedasticity
 Jie Yang^{1}Email author,
 George Casella^{2, 4} and
 Lauren M McIntyre^{2, 3, 4}
DOI: 10.1186/1471210512427
© Yang et al; licensee BioMed Central Ltd. 2011
Received: 10 June 2010
Accepted: 1 November 2011
Published: 1 November 2011
Abstract
Background
Many analyses of gene expression data involve hypothesis tests of an interaction term between two fixed effects, typically tested using a residual variance. In expression studies, the issue of variance heteroscedasticity has received much attention, and previous work has focused on either betweengene or withingene heteroscedasticity. However, in a single experiment, heteroscedasticity may exist both within and between genes. Here we develop flexible shrinkage error estimators considering both betweengene and withingene heteroscedasticity and use them to construct Flike test statistics for testing interactions, with cutoff values obtained by permutation. These permutation tests are complicated, and several permutation tests are investigated here.
Results
Our proposed test statistics are compared with other existing shrinkagetype test statistics through extensive simulation studies and a real data example. The results show that the choice of permutation procedures has dramatically more influence on detection power than the choice of F or Flike test statistics. When both types of gene heteroscedasticity exist, our proposed test statistics can control preselected typeI errors and are more powerful. Raw data permutation is not valid in this setting. Whether unrestricted or restricted residual permutation should be used depends on the specific type of test statistic.
Conclusions
The Flike test statistic that uses the proposed flexible shrinkage error estimator considering both types of gene heteroscedasticity and unrestricted residual permutation can provide a statistically valid and powerful test. Therefore, we recommended that it should always applied in the analysis of real gene expression data analysis to test an interaction term.
Background
The regulation of gene expression starts when a cell's DNA is transcribed into mRNA. The simultaneous expression profiles of many genes under different circumstances can provide insight into physiological processes. Using modern technologies in gene expression experiments such as oligonucleotide arrays [1], and cDNA spotted arrays [2], many scientists have made novel discoveries about complex biological processes of yeast [3, 4], drosophila [5], mice [6], humans [7], and other species. Recently one such study also included RNAseq [8]. Statistical methodologies and issues involved in microarray data analysis have been widely reviewed [9–12], and it is expected that many of the same issues will need to be addressed with RNAseq.
The analysis of variance (ANOVA) model is a popular statistical modeling method for the analysis of microarrays. Since its introduction by Kerr et al. [13], it has been extensively examined for use in this setting [14–21]. Kerr et al. constructed an ANOVA model that included the gene effect as a fixed effect. This model assumes identically and independently distributed residual errors across genes. The advantage of this model is that the large number of genes involved in a microarray experiment results in huge degrees of freedom for the error estimate, which can lead to a very powerful test. However, the common assumption of homoscedasticity may not hold true in this setting [22]. One alternative is to use an ANOVA model for each gene, but the resulting test statistics from genespecific models may have limited power because the biological sample size for each gene in a microarray experiment is usually small.
To address this problem of limited power, researchers have proposed other methods for obtaining more information across genes, ranging from a simple equalweighted average of a genespecific error estimate and the global average of all genespecific error estimates (F_{2} statistic proposed by Wu et al. [19] to empirical Bayesian modeling of all genespecific errors [23–26]. Other variations [27–29] used different variance modeling strategies to address the heteroscedasticity problem, but no clear winner has emerged [30]. Huang and Liu [31] extended the test statistics proposed by Cui et al. [28] by assuming a normal distribution on the mean and then deriving an empirical Bayes likelihood ratio test. The resulting test statistic shrinks both the mean and variances.
In addition to the problem of betweengene heteroscedasticity, we must also be concerned with withingene heteroscedasticity. For example, in the study of simple differential gene expression between a treatment group and a control group, the variance in the treatment arm may differ from that in the control arm. Some approaches to this problem include a general Bayesian framework to model heteroscedastic error in a single generalized linear mixed model setting [32] and a structural model placed on the error variances specific to each gene and treatment combination [33].
As gene expression studies become more popular, the complexity of the experiment increases. Instead of only simple treatment and control experiments, two or more factor experiments are being conducted. This increase in experiment complexity has led to many scientific questions involving the hypothesis testing of an interaction between two factors. For example, testing a probe by genotype interaction can result in inferences about polymorphism in the probe, such as single nucleotide polymorphism (SNP) and insertiondeletion (indel) [34–37]; testing a probe by sex can imply that alternative splicing occurs between male and female subjects [38]; and in pharmacogenomic studies, testing the genotypedrug/treatment or genotypedisease interaction may be of interest [39]. Thus far, all the development of ANOVA methods for microarray studies has focused on tests of main effects.
Here, a generalized shrinkage estimator incorporating both within and betweengene heteroscedasticities is developed (see Lehmann and Cesella [40] for a review of shrinkage estimation). In any given experiment, both withingene and betweengene heteroscedasticity may exist; thus, taking these possibilities into account should lead to an improved test statistic. Moreover, given the increasing complexity of recent studies and the burgeoning interest in hypotheses that involve interactions, we focus on an improved shrinkagebased Ftest for interaction terms.
Methods
Here we develop new shrinkage estimates for the error term and show how to use these estimates to construct Flike statistics. We then estimate the null distribution of these statistics by using permutation tests.
Shrinkage error estimators
i = 1,...,I, k = 1,..., K, and let ${\sigma}_{i,k}^{2}$ be the true variance of gene i in group k. When the experimental design is balanced, ${\widehat{\sigma}}_{i,k}^{2}$ is the residual mean square for gene i in group k and $\nu {\widehat{\sigma}}_{i,k}^{2}\u2215{\sigma}_{i,k}^{2}~{\chi}_{\nu}^{2}$, where ν represents the degrees of freedom for the error estimates.
 1.
Specific values for each genegroup combination
 2.
Genespecific values that are the same across all other groups
 3.
Groupspecific values that are the same across genes but different across groups
 4.
A single point representing the underlying common error
Correspondingly, these targets are correct when (1) there are both withingene and betweengene heteroscedasticity; (2) there is only betweengene heteroscedasticity; (3) there is only withingene heteroscedasticity; and (4) all error variances are identical. We now develop a generalized shrinkage error estimator using these four shrinkage targets.
where $\stackrel{\u0303}{\theta}=\left({\theta}_{1,1},...,{\theta}_{1,K,}...,{\theta}_{I,1},...,{\theta}_{I,K}\right),\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{0.3em}{0ex}}\stackrel{\u0303}{\alpha}=\left({\alpha}_{1},...,{\alpha}_{I}\right)$ represents the genespecific mean differences, and $\stackrel{\u0303}{\beta}=\left({\beta}_{1},...,{\beta}_{K}\right)$ models different means with respect to different classes of the subgroups.
the empirical Bayes estimator for τ^{2} is $\Sigma {\left({X}_{i,k}\widehat{\mu}{\widehat{\alpha}}_{i}{\widehat{\beta}}_{k}\right)}^{2}\u2215\left[IK\left(I+K1\right)2\right]{\sigma}^{2}$.
where X_{ i }is the residual variance estimate from a genespecific model, and m and σ^{2} are the mean and variance of log $\frac{{\chi}^{2}K\nu}{K\nu}$. The underlying assumption for this estimator is that there is no betweengene heteroscedasticity, as this estimator shrinks every genespecific error estimator toward one target. Therefore, it will overshrink the genespecific error estimates when gene heteroscedasticity exists. In comparison, generalized shrinkage error estimators are flexible in terms of incorporating a different type of heteroscedasticity. Some degrees of freedom are used for incorporating the heteroscedasticity. However, the gain is that the error estimator is then closer to the underlying distribution and should lead to better performance of the resultant Flike test statistics as shown in the results section.
Results from raw data permutation
Restricted?  Data set  F _{1}  F _{2}  F _{3}  F _{ Cui }  F _{ Gen }  F _{ Gengene }  F _{ Gengrp } 

YES  nullce  5.05(0.07)  5.06(0.08)  5.12(0.17)  5.09(0.10)  5.09(0.10)  5.05(0.08)  5.11(0.10) 
nullgh  5.02(0.07)  5.13(0.16)  5.26(0.20)  5.03(0.07)  5.07(0.12)  5.03(0.07)  5.11(0.16)  
nullwgh  4.97(0.07)  4.96(0.09)  4.93(0.18)  4.99(0.08)  4.99(0.12)  4.96(0.09)  5.01(0.16)  
nullbgh  5.02(0.07)  4.99(0.17)  5.03(0.21)  5.02(0.07)  5.02(0.15)  5.01(0.09)  5.03(0.18)  
NO  nullce  5.10(0.07)  5.06(0.08)  5.06(0.08)  7.4(0.12)  5.15(0.09)  5.12(0.08)  5.08(0.08) 
nullgh  5.08(0.07)  5.12(0.16)  5.12(0.12)  7.4(0.09)  5.10(0.11)  5.07(0.10)  5.12(0.09)  
nullwgh  12.31(0.10)  7.56(0.10)  4.61(0.10)  17.37(0.14)  5.32(0.11)  5.07(0.09)  5.87(0.11)  
nullbgh  12.31(0.11)  6.63(0.17)  5.55(0.19)  15.68(0.12)  6.30(0.12)  6.10(0.11)  6.30(0.11) 
Shrinkage Flike statistics
The error variance estimator in this test uses data from only gene i. In oligonucleotide microarray models, the degrees of freedom for the error estimate can be small because the sample size of RNA is usually small, and hence the power of F_{1} can be limited.
Permutation tests
For the proposed generalized shrinkage Flike test statistics, the null distributions are not known named distributions. Therefore, an empirical approach such as a permutation test can be used to estimate the null distributions. The permutation test for interaction is complicated, because there is no exact permutation test for such a purpose [43]. We therefore must consider an approximate permutation method for testing an interaction term in a crossed fixed/mixed model [44, 45].
Permutation approaches developed previously focused on a single ANOVA model. In the typical gene expression study, thousands of ANOVA models are considered simultaneously. The additional complexity of the shrinkage Flike statistics indicates that Monte Carlo studies are needed to investigate the performance of residual permutation and raw data permutation, with restrictions or not, in a geneexpression analysis. The choice of permutation procedures is critical for assessing the performance of a test statistic.
For all the modified Flike statistics presented in the previous section, the null distributions can only be approximated empirically, but permutation procedures can be used to find the approximate null distribution of all the F and Flike statistics. The important issues in performing a permutation analysis include the choice of the exchangeable units under the null hypothesis, the choice of using restricted permutation or not, and the choice of residual permutation or raw data permutation. These choices influence the power of a test statistic.
Residual permutation using residuals from a reduced model and unrestricted raw data permutation can be used to approximate the null distribution of a statistic for testing an interaction term [44]. When using F_{1} to test an interaction term in a single ANOVA model, the residual permutation leads to a more powerful test than unrestricted raw data permutation [44]. However, in gene expression analysis, thousands of genespecific ANOVA models are simultaneously considered, and for a particular genespecific ANOVA model, information from other genespecific ANOVA models is used to construct the shrinkage error estimate. Hence, both residual permutation and raw data permutation were investigated. Furthermore, both restricted and unrestricted permutations were studied, because the permutation units are exchangeable only within each particular group when withingene heteroscedasticity is present across those subgroups.
Results
The properties of this shrinkage estimator are compared with those of other existing F and Flike statistics that have been proposed and described in the "Shrinkage Flike statistics" section.
Simulation studies
The purpose of these simulation studies was to compare the performances of F_{1}, F_{2}, F_{3}, F_{ Cui }, F_{ Gen }, F_{ Gengene }, and F_{ Gengrp }in terms of type I error and power and to compare the results of a particular Flike statistic using four different permutation strategies: restricted/unrestricted residual permutation and restricted/unrestricted raw data permutation.
In these simulation studies, 100 genes with two probes for each gene and three replicates from each of two lines were simulated to mimic a splitplot design in a general oligonucleotide microarray experiment. The genespecific ANOVA model in which data were generated from the model, y_{ plr }= P_{ p }+ L_{ l }+ RL_{ rl }+ PL_{ pl }+ ϵ_{ plr }, wp = 1, 2, l = 1,2, r = 1,2,3, where P, L, RL, and PL represent probe, line, replicates from a particular line, and the interaction between probe and line, respectively.
Replicates were nested within each line, and RL is usually treated as a random effect during the modelfitting procedure, which results in a correlation between probes from the same biological sample. In the simulated data sets, the correlation between genes was 0. As many as 900 simulation runs were carried out to compare the performances of F_{1}, F_{2}, F_{3}, F_{ Cui }, F_{ Gen }, F_{ Gengene }, and F_{ Gengrp }based on different permutation procedures. The four permutations tested were unrestricted residual permutation, restricted residual permutation with respect to each line, unrestricted raw data permutation, and restricted raw data permutation with respect to each line. The residuals permuted were from a reduced fixed model with fixed effects for only line and probe.
Results from residual permutation
Restricted?  Data set  F _{1}  F _{2}  F _{3}  F _{ Cui }  F _{ Gen }  F _{ Gengene }  F _{ Gengrp } 

YES  nullce  4.59(0.07)  4.15(0.08)  3.63(0.14)  4.09(0.09)  4.55(0.07)  3.23(0.06)  4.44(0.08) 
nullgh  4.57(0.07)  4.1(0.13)  3.95(0.16)  4.49(0.07)  4.6(0.07)  4.61(0.07)  4.38(0.07)  
nullwgh  6.74(0.08)  4.33(0.08)  3.51(0.14)  6.49(0.09)  4.38(0.07)  4.2(0.07)  2.78(0.09)  
nullbgh  6.74(0.08)  4.35(0.16)  4.07(0.19)  6.58(0.08)  4.36(0.07)  4.16(0.07)  3.64(0.07)  
NO  nullce  5.1(0.07)  4.99(0.08)  4.5(0.08)  4.59(0.08)  4.99(0.07)  4.1(0.07)  4.68(0.07) 
nullgh  5.1(0.07)  4.83(0.1)  4.59(0.11)  5.08(0.07)  4.99(0.07)  5.01(0.07)  4.95(0.07)  
nullwgh  10.75(0.09)  8.46(0.09)  7.6(0.09)  12.37(0.11)  5.03(0.08)  6.43(0.08)  4.93(0.08)  
nullbgh  10.75(0.1)  8.38(0.17)  8.07(0.19)  10.79(0.1)  5.02(0.08)  6.38(0.08)  6.73(0.08) 
Further simulations to compare the rejection rates were conducted. Only results from residual permutation are shown because it was found that raw data permutation was less powerful than residual permutation. This is consistent with the findings of Anderson and Ter Braak [44]. Figure 1 shows the estimated average null hypothesis rejection rate curves from all Flike statistics and both restricted and unrestricted residual permutation procedures. The xaxis represents the average null hypothesis rejection rate using F_{1} and the tabulated pvalues. The solid line shows that the corresponding statistic controls the prespecified CWER, and the dashed line shows that the corresponding CWER was inflated. In general, restricted residual permutation is less powerful than unrestricted residual permutation. For example, the power of all statistics from unrestricted residual permutation almost doubled in some cases where heteroscedasticity existed.
Drosophila data
The data used in this study are from a gene expression comparison study between D. melanogaster and D. simulans [46]. Expression of 10 genotypes of each species was measured in male flies. In D. simulans, each genotype was measured separately, and in D. melanogaster, a pool of 10 genotypes was measured. All genotypes (individual or pooled) were independently isolated and hybridized three times. The goal of the original study was to provide a genomewide approach to identifying candidate genes potentially responsible for adaptation and speciation in D. simulans and D. melanogaster. In this study, we focus on identifying sequence differences between genotypes in D. simulans based on hybridization profiles. Withingene heteroscedasticity is expected because the genotypes come from different lines. The proposed generalized shrinkage Flike test statistics F_{ Gen }, F_{ Gengene }, and F_{ Gengrp }were compared with F_{2}, F_{3} with restricted residual permutation, which could control prespecified CWER for any variance structure in simulation studies. Furthermore, Smyth's moderated Ftest statistic [25] without multiple testing adjustment and controlling the false discovery rate (FDR) at 5% were used for comparison. As the main interest is in sequence difference, the focus is on the test of interaction between line and probe. The split plot model described above is used. SAS program codes are included in the additional files (additional file 1 and additional file 2).
Probe sets with significant line*probe terms found by Flike test statistics and appropriate residual permutation procedures and Smyth's moderated Ftest statistic
Test statistic  Restricted permutation?  Number of probe sets found  True false discovery rate  Power 

F _{2}  Yes  124  22.6%  12.4% 
F _{3}  Yes  187  29.4%  17.1% 
F _{ Gen }  No  453  29.5%  41.1% 
F _{ Gengene }  No  455  28.8%  41.7% 
F _{ Gengrp }  No  474  28.9%  43.4% 
F _{ Gen }  Yes  136  24.3%  13.3% 
F _{ Gengene }  Yes  122  22.1%  12.2% 
F _{ Gen grp }  Yes  116  21.5%  11.7% 
moderatedF  1  N/A  535  34.1%  75.5% 
moderatedF  2  N/A  813  34.4%  68.8% 
Discussion
For gene expression analysis, ANOVA models have been a popular modeling technique. Based on ANOVA models, flexible shrinkage Flike test statistics were developed to account for both the withingene and betweengene heteroscedasticities. The emphasis here is on testing an interaction term, as this case is of increasing interest to biologists, and there is no clear existing theory on the most powerful, valid approach for such statistics. For all Flike statistics studied here, their null distributions were approximated empirically through permutations. Four different permutation procedures were investigated for eight different Flike statistical tests of the interaction term.
As expected, we found that when an error estimator overshrinks, the resulting Flike statistic cannot control the prespecified CWER. For example, F_{ Gengene }is an overshrinkage error estimator when there is withingene heteroscedasticity. As a result, compared with generalized shrinkage Flike statistics, it is not valid when withingene heteroscedasticity exists. Undershrinkage is also important, as it will lead to a conservative test and lower power. This is clearly demonstrated when the common error can be assumed and the most powerful valid test is F_{ Gengrp }.
The most striking result was the impact of the permutation procedures. Although this was not completely unexpected [43–45], the effect of the permutation procedures is dramatic and worthy of special attention. Unrestricted raw data permutation could not control prespecified CWER when there was withingene heteroscedasticity. Restricted raw data permutation could be used, but it was less powerful than residual permutation. Also consistent with findings from Anderson and Ter Braak [44], restricted permutations are less powerful than unrestricted permutations. However, unrestricted permutations are valid only for a common error and when betweengene heteroscedasticity exists for our proposed shrinkage statistics; they are not valid in combination with F_{2}, F_{3}, or F_{ Cui }. For F_{ Gengrp }, the unrestricted permutation can also be used in cases having withingene heteroscedasticity, while only F_{ Gen }is valid with unrestricted permutation in all cases in terms of controlling prespecified CWER. Interestingly, the power gain from using the correct shrinkage target F_{ Gengrp }rather than F_{ Gen }is far less than that of using unrestricted permutation. The result is that F_{3} is never the most powerful choice when testing an interaction term.
The correct shrinkage target can lead to the most powerful test statistic. As one of the reviewers suggested, a statistical test may be applied to help pick the best shrinkage target before obtaining shrinkage error estimates. However, this extra testing step may inflate the CWER of the test statistic when there is gene heteroscedasticity. For example, when there are both types of gene heteroscedasticities, it is possible that the above test suggests only withingene heteroscedasticities exist, and F_{ Gengrp }is shown to inflate the CWER. There is minimal penalty to using the shrinkage estimator we propose, so we recommend setting the shrinkage target in the full space spanned by group and gene and using unrestricted permutation to compensate for the possible power loss in fewer degrees of freedom left for estimating the errors.
Conclusions
The proposed generalized shrinkage Flike statistic with shrinkage targets located in a space spanned by gene and another group, F_{ Gen }, with unrestricted residual permutation is always valid in terms of having a prespecified CWER. This statistic has reasonable power in most cases; thus, it is generally recommended to be applied to test an interaction term in the analysis of real gene expression data.
List of abbreviations
 CWER:

comparisonwise error rate
 FDR:

false discovery rate
 indel:

insertion and deletion
 SNP:

single nucleotide polymorphism
Declarations
Acknowledgements
We thank Brandon Walts for identifying true SNP positions; Angela J. McArthur and David R. Galloway for their help in scientific editing; associate editor and three reviewers for their constructive comments that much improved this manuscript. This research was supported by NIH 1R01GM077618 (McIntyre), NIH 1R01GM081704 (Casella).
Authors’ Affiliations
References
 Fodor SPA: Massively parallel genomics. Science 1997, 277: 393–395. 10.1126/science.277.5324.393View ArticleGoogle Scholar
 Schena M, Shalon D, Heller R, Chai A, Brown PO, Davis RW: Parallel human genome analysis: microarraybased expression monitoring of 1000 genes. Proceedings of National Academy Science 1996, 93: 10614–10619. 10.1073/pnas.93.20.10614View ArticleGoogle Scholar
 Galitski T, Saldanha AJ, Styles CA, Lander ES, Fink GR: Ploidy regulation of gene expression in yeast. Science 1999, 285: 251–254. 10.1126/science.285.5425.251View ArticlePubMedGoogle Scholar
 Tu BP, Kudlicki A, Rowicka M, McKnight SL: Logic of the yeast metabolic cycle: Temporal compartmentalization of cellular processes. Science 2005, 310: 1152–1158. 10.1126/science.1120499View ArticlePubMedGoogle Scholar
 White KP, Rifkin SA, Hurban P, Hogness DS: Microarray analysis of Drosophila development during metamorphosis. Science 1999, 286: 2179–2184. 10.1126/science.286.5447.2179View ArticlePubMedGoogle Scholar
 Chabas D, Baranzini SE, Mitchell D, Bernard CCA, Rittling SR, Denhardt DT, Sobel RA, Lock C, Karpuj M, Pedotti R, Heller R, Oksenberg JR, Steinman L: The influence of the proinflammatory cytokine, Osteopontin, on autoimmune demyelinating disease. Science 2001, 294: 1731–1735. 10.1126/science.1062960View ArticlePubMedGoogle Scholar
 Sebat J, Lakshmi B, Troge J, Alexander J, Young J, Lundin P, Maner S, Massa H, Walker M, Chi MY, Navin N, Lucito R, Healy J, Hicks J, Ye K, Reiner A, Gilliam TC, Trask B, Patterson N, Zetterberg A, Wigler M: Largescale copy number polymorphism in the human genome. Science 2005, 305: 525–528.View ArticleGoogle Scholar
 Blekhman R, Marioni JC, Zumbo P, Stephens M, Gilad Y: Sexspecific and lineagespecific alternative splicing in primates. Genome Research 2010, 20(2):180–189. 10.1101/gr.099226.109PubMed CentralView ArticlePubMedGoogle Scholar
 Butte A: The use and analysis of microarray data. Nature Reviews 2002, 1: 951–960. 10.1038/nrd961PubMedGoogle Scholar
 Churchill GA: Fundamentals of experimental design for cDNA microarrays. Nature Genetics 2002, 32: 490–495. 10.1038/ng1031View ArticlePubMedGoogle Scholar
 Craig BA, Black MA, Doerge RW: Gene expression data: The technology and statistical analysis. Journal of Agricultural, Biological, and Environmental Statistics 2003, 8(1):1–28. 10.1198/1085711031256View ArticleGoogle Scholar
 Allison DA, Cui X, Page GP, Sabripour M: Microarray data analysis: from disarray to consolidation and consensus. Nature Reviews Genetics 2006, 7: 55–65. 10.1038/nrg1749View ArticlePubMedGoogle Scholar
 Kerr MK, Martin M, Churchill GA: Analysis of variance for gene expression microarray data. Journal of Computational Biology 2000, 7: 819–837. 10.1089/10665270050514954View ArticlePubMedGoogle Scholar
 Kerr MK, Churchill GA: Statistical design and the analysis of gene expression microarrays. Genetical Research 2001, 77: 123–128.PubMedGoogle Scholar
 Kerr MK, Churchill GA: Experimental design for gene expression microarrays. Biostatistics 2001, 2: 183–201. 10.1093/biostatistics/2.2.183View ArticlePubMedGoogle Scholar
 Pritchard CC, Hsu L, Delrow J, Nelson PS: Project normal: Defining normal variation in mouse gene expression. Proceedings of the National Academy of SciencesUSA 2001, 98: 13266–13271. 10.1073/pnas.221465998View ArticleGoogle Scholar
 Wolfinger RD, Gibson G, Wolfinger ED, Bennett L, Hamadeh H, Bushel P, Ashfari C, Paules RS: Assessing gene significance from cDNA microarray expression data via mixed models. Journal of Computational Biology 2001, 8(6):625–637. 10.1089/106652701753307520View ArticlePubMedGoogle Scholar
 Kerr MK, Afshari CA, Bennett L, Bushel P, Martinez J, Walker NJ, Churchill GA: Statistical analysis of a gene expression microarray experiment with replication. Statistica Sinica 2002, 12: 203–217.Google Scholar
 Wu H, Kerr MK, Cui XQ, Churchill GA: MAANOVA: A Software package for the analysis of spotted cDNA microarray experiments, In. In The analysis of gene expression data: methods and software. Springer; 2002:313–341.Google Scholar
 Chu T, Weir B, Wolfinger R: A systematic statistical linear modeling approach to oligonucleotide array experiments. Mathematical Biosciences 2002, 176: 35–51. 10.1016/S00255564(01)001079View ArticlePubMedGoogle Scholar
 Wayne ML, Pan YJ, Nuzhdin SV, McIntyre LM: Additivity and transacting effects on gene expression in male Drosophila simulans . Genetics 2004, 168: 1413–1420. 10.1534/genetics.104.030973PubMed CentralView ArticlePubMedGoogle Scholar
 Cui X, Churchill GA: Statistical tests for differential expression in cDNA microarray experiments. Genome Biology 2003, 4(4):201.View ArticleGoogle Scholar
 Baldi P, Long AD: A Bayesian framework for the analysis of microarray expression data: Regularized ttest and statistical inferences of gene changes. Bioinformatics 2001, 17: 509–519. 10.1093/bioinformatics/17.6.509View ArticlePubMedGoogle Scholar
 Lönnstedt I, Speed T: Replicated microarray data. Statistica Sinca 2002, 12: 31–46.Google Scholar
 Smyth GK: Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 2004., 3: No. 1, Article 3 No. 1, Article 3Google Scholar
 Tong TJ, Wang YD: Optimal shrinkage estimation of variances with applications to microarray data analysis. Journal of the American Statistical Association 2007, 102: 113–122. 10.1198/016214506000001266View ArticleGoogle Scholar
 Tusher VG, Tibshirani R, Chu G: Significance analysis of microarrays applied to ionizing radiation response. The Preceedings of National Academy Science 2001, 989: 5116–5121.View ArticleGoogle Scholar
 Cui X, Hwang JTG, Qiu J, Blades NJ, Churchill GA: Improved statistical tests for differential gene expression by shrinking variance components estimates. Biostatistics 2005, 6: 59–75. 10.1093/biostatistics/kxh018View ArticleGoogle Scholar
 Feng S, Wolfinger RD, Chu TM, Gibson GC, McGraw LA: Empirical Bayes analysis of variance component models for microarray data. Journal of Agricultural, Biological & Environmental Statistics 2006, 1113: 197–190.View ArticleGoogle Scholar
 Kim SY, Lee JW, Sohn IS: Comparison of various statistical methods for identifying differential gene expression in replicated microarray data. Statistical Methods in Medical Research 2006, 15: 3–20. 10.1191/0962280206sm423oaView ArticlePubMedGoogle Scholar
 Hwang JTG, Liu P: Optimal tests shrinking both means and variances applicable to microarray data analysis. In preprint 2007–03. Department of Statistics, Iowa State University, IA; 2007.Google Scholar
 Kizilkaya K, Tempelman RJ: A general approach to mixed effects modeling of residual variances in generalized linear mixed models. Genetics Selection Evolution 2005, 37: 31–56. 10.1186/1297968637131View ArticleGoogle Scholar
 Jaffrezic F, Marot G, Degrelle S, Isabelle H, Foulley JL: A structural mixed model for variances in differential gene expression studies. Genetical Research 2007, 89(1):19–25. 10.1017/S0016672307008646View ArticlePubMedGoogle Scholar
 Rostoks N, Borevitz JO, Hedley PE, Russell J, Mudie S, Morris J, Cardle L, Marshall DF, Waugh R: Singlefeature polymorphism discovery in the Barley transcriptome. Genome Biology 2005, 6: R54. 10.1186/gb200566r54PubMed CentralView ArticlePubMedGoogle Scholar
 Kirst M, Caldo R, Casati P, Tanimoto G, Walbot V, Wise RP, Buckler ES: Genetic iversity contribution to errors in short oligonucleotide microarray analysis. Plant Biotechnology Journal 2006, 4: 489–498.PubMedGoogle Scholar
 Zhang X, Shiu SH, Cal A, Borevitz JO: Global analysis of genetic, epigenetic and transcriptional polymorphisms in Arabidopsis Thaliana Using whole genome tiling arrays. PLoS Genetics 2008, 4(3):e1000032. 10.1371/journal.pgen.1000032PubMed CentralView ArticlePubMedGoogle Scholar
 Zhang X, Borevitz JO: Global Analysis of Allelespecific Expression in Arabidopsis Thaliana. Genetics 2009, 182(4):943–954. 10.1534/genetics.109.103499PubMed CentralView ArticlePubMedGoogle Scholar
 McIntyre LM, Bono LM, Genissel A, Westerman R, Junk D, TelonisScott M, Harshman L, Wayne ML, Kopp A, Nuzhdin SV: Sex specific expression of alternative transcripts in Drosophila. Genome Biology 2006, 7: R79. 10.1186/gb200678r79PubMed CentralView ArticlePubMedGoogle Scholar
 Kelly P, Zhou YH, Whitehead J, Stallard N, Bowman C: Sequentially testing for a genedrug interaction in a genomewide analysis. Statistics in Medicine 2008, 27: 2022–2034. 10.1002/sim.3059View ArticlePubMedGoogle Scholar
 Lehmann EL, Casella G: Theory of Point Estimation. 2nd edition. New York: SpringerVerlag; 1998.Google Scholar
 Pounds S: Computational enhancement of a shrinkagebased analysis of variance Ftest proposed for differential gene expression analysis. Biostatistics 2007, 83: 505–506.View ArticleGoogle Scholar
 Neter J, Wasserman W, Kutner MH: Applied Linear Statistical Models: Regression, Analysis of Variance, and Experimental Designs. 3rd edition. Irwin, Inc; 1990.Google Scholar
 Edgington ES: Randomization Tests. 3rd edition. Marcel Dekker, New York; 1995. (1995) (1995)Google Scholar
 Anderson MJ, Ter Braak CJF: Permutation tests for multifactorial analysis of anova. Journal of Statistical Computation and Simulation 2003, 732: 85–113.View ArticleGoogle Scholar
 Churchill GA, Doerge RW: Naive application of permutation testing leads to nflated type I error rates. Genetics 2008, 178: 609–610. 10.1534/genetics.107.074609PubMed CentralView ArticlePubMedGoogle Scholar
 Nuzhdin SV, Wayne ML, Harmon KL, McIntyre LM: Common pattern of evolution of gene expression level and protein sequence in drosophila. Molecular Biology and Evolution 2004, 21: 1308–1317. 10.1093/molbev/msh128View ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.