BACKGROUND
microarrays provide large-scale comparative gene expression profiles between biological samples by detecting differential expression for thousands of genes in parallel.
typically, systematic error  in the measurements is removed at the background correction and normalization steps, and is followed by statistical testing.
in the most common type of study, statistical testing produces a list of genes that are differentially expressed across two or more classes   <cit> .
extensive research has shown that choice of pre-processing methods designed to correct for bias in the measurements can have a substantial impact on rank ordering of gene expression fold-change  estimates for both cdna  <cit>  and oligonucleotide data  <cit>  and the affycomp website  <cit> .
one major disadvantage of fc estimates, however, is that they do not take the variance of the samples into account.
this is especially problematic because variability in gene expression measurements is partially gene-specific  <cit> , even after the variance has been stabilized by data transformation  <cit> .
there is consensus in the statistical community that statistical tests of differential expression are preferred over fc for inference  <cit> .
one advantage is that they standardize differential expression by dividing fc measurements by their associated standard error, rescaling fcs to a common metric.
moreover, associated output such as p values, effect sizes, and confidence intervals can be used for various purposes such as false positive control  <cit>  and meta-analysis  <cit> .
although microarray studies with sample sizes of five or more observations per class are becoming increasingly common, cost considerations and the need for small scale initial studies mean that many studies are conducted with smaller sample sizes.
the use of classical statistical tests such as the t-test is sub-optimal for small sample studies, however, because of low statistical power for detection of differential expression.
this has led to the development of small sample size error estimation methods which borrow information from the entire data set or from a subset of the data to improve error estimation accuracy and precision  <cit> .
there has been a number of studies analysing the performance of statistical tests applied to microarray data but few have used data where the differentially expressed genes are known in advance.
qin and kerr  <cit>  compared normalization methods and test statistics using data sets with six known differentially expressed genes and showed that the standard t-statistic performed worse than those that used variance estimates that incorporated information from all genes.
sioson et al.
<cit>  compared statistical methodologies of two software applications using qrt-pcr of a subset of genes and chen et al.
<cit>  used a consensus of differentially expressed genes across statistical methods to analyze performance.
none of these studies, however, evaluated statistical tests designed for small sample size experiments.
a number of studies have examined small n tests.
kooperberg et al.
<cit>  compared the performances of various 2-group statistical tests with empirical and simulated data.
distributions of p-values were generated based on data from the same experimental group  or data from different experimental groups which were known to differ.
the high performing tests generated p-value distributions consistent with a null distribution in the former case and produced the largest number of small p-values in the latter case.
using the limma method  <cit>  as an exemplar, they concluded that an empirical bayes approach to statistical testing provided good power while accurately controlling false positive rates for small n microarray experiments.
cui and churchill  <cit>  reviewed a number of statistical tests for use with microarray data, including one small n test, but provided no comparative analysis.
tong and wang  <cit>  used simulated data to explore the theoretical properties of shrinkage methods of estimating variance and found that they outperformed tests that use only the sample variance.
hu and wright  <cit>  and xie et al.
<cit>  both compared a number of statistical tests, including several small n tests, using a consensus list of differentially expressed genes from all methods.
hu and wright  <cit>  found that, based on false discovery rates, tests that model the variance/intensity relationship and use variance estimates generated with information from all genes performed the best.
xie et al  <cit>  showed that there was comparability of results for only a few of the methods but the lönnstedt and speed  <cit>  b-statistic, which is monotonically equivalent to the limma  <cit>  and brb  <cit>  t-statistics, had the lowest false positive rate.
jefferey et al.
<cit>  looked at the use of statistical tests, including several small n tests, in feature selection for group classification; the results varied greatly across gene list and sample size but the empirical bayes t-statistic performed well across all sample sizes.
most of these studies compared methods using comparability of results, estimated false discovery, family wise error rates, or p-value distributions to assess performance.
in this paper we focus on each method's ability to detect genes which have been spiked at different concentrations across samples .
we use simulated data, two publicly available empirical affymetrix latin-square data sets, and two cdna data sets in which the differentially expressed genes are known to assess the relative performances of fc, t-statistic, and four small sample size statistical tests which all borrow strength from other genes in different ways.
for consistency, we refer to the various tests by their software implementation.
the local pooled error method  z-statistic  <cit>  is based on pooling errors for genes with similar intensities.
the three empirical bayes methods were: the brb t-statistic  <cit>  which combines gene-specific error estimates with a common error estimate obtained from the distribution of the variances across all genes; the cybert regularized t-statistic  <cit>  which combines gene-specific error estimates with a local pooled error estimate based on genes with similar intensities; and the limma method  <cit>  which combines a fitted linear model of gene expression data with a variance estimate into a moderated t-statistic.
all methods gain degrees of freedom over the standard t-test and in theory should provide greater sensitivity with no loss in specificity.
data
latin square spike-in data
the hgu133a and hgu <dig> latin square data sets  <cit>  are based on a  <dig> ×  <dig> latin square design of "spiked-in" transcripts  with three replicates for each group.
the concentrations for the "spiked-in" transcripts were doubled for each consecutive group .
we compared performance of the statistical tests  for detecting two-fold differential expression after the data had been normalized by six popular expression summary algorithms: mas  <dig>   <cit> , dchip  <cit> , rma  <cit> , gcrma  <cit> , vsn  <cit> , and lmgene  <cit> .
data were analysed on the log scale for the first four methods and on a generalized log  scale for the latter two methods.
cdna spike-in data
two spotted cdna experiments  <cit>  were also used to compare statistical test performance.
the first experiment compared two aliquots from the same mouse liver rna sample, and the second compared one rna sample from mouse liver to a pooled sample of five different mouse tissues .
in both experiments six transcripts were "spiked-in" at a three-fold difference in concentration between the two groups.
all cdna slides measured each "spike-in" transcript at  <dig> different well locations for a total of  <dig> "spikes" in each experiment.
the data were normalized using the lowess normalization of lmgene  <cit> , with and without a glog transformation and background correction , and the vsn method  <cit> .
the lmgene lowess normalization method is equivalent to the intensity-dependent normalization used in yang et al.
<cit> .
mva plots and their corresponding lowess fits for each cdna array are shown in additional file  <dig>
simulated data
two limitations of the above spike-in data sets make them less challenging than most biologically motivated studies: variation across technical replicates is lower than that typically observed across biological replicates  <cit>  and many biological effects of interest may be smaller than two-fold.
also, there is only a small number of "spiked-in" transcripts in each latin square experiment, which can lead to increased variability in the algorithm performance metrics  <cit> .
simulated data were generated to offset these limitations: null hypothesis and "spike-in" data  were used to examine false positive rates , genes incorrectly labeled as differentially expressed, and true positive rates , genes correctly labeled as differentially expressed.
the "spike-in"data incorporated varying degrees of variability between replicate measurements and genes were differentially expressed across a range of fold changes to offset the limitations of the latin square data sets.
each simulation experiment consisted of  <dig> iterations of ten thousand genes for each of two groups  with  <dig>   <dig>  or  <dig> replicates per gene.
gene intensity values were generated by randomly selecting from a n distribution, where u is the "true" expression intensity  distribution) for a single gene in log <dig> scale and σ <dig> is the random error associated with the gene's expression measurement.
in the null hypothesis data, each gene-specific u remains constant across both groups; in the "spike-in" data, μ was upwardly regulated by a specified fold change for  <dig> genes in the treatment group.
for each "spiked-in" gene, μ was assigned a value from a set of  <dig> log <dig> values ranging from  <dig> to  <dig>  in  <dig>  increments.
the log <dig> fold change applied to the "spiked-in genes" was implemented by an additive adjustment to their u in the treatment group chosen from a range of  <dig> values from  <dig>  to  <dig>   in  <dig>  increments.
there are  <dig> "spiked-in" genes at each of the  <dig> concentration values and each of those  <dig> genes at a particular concentration value is upwardly regulated by one of the  <dig> different fold changes.
the random error of each gene is derived from one of three variance models: common, inverse gamma, and local .
in the common error model, the population variance  of the random error term  is a constant across all genes.
in the inverse gamma variance model, c is drawn randomly from an inverse gamma distribution, ig, where  <dig> is the shape parameter and b is the scale parameter.
in the local variance model, ϵ is a combination of additive  and proportional  components, as outlined in rocke and durbin  <cit> .
the expression intensity column shows how the replicate intensity values were generated for each transcript.
the random error column shows the method of calculating specific ϵ values.
the parameters column shows the values that were used to generate high, medium, and low variances, respectively.
the parameters of each variance model were varied to create data sets from low variance  to medium and high variance ; c and d show hgu133a data with the "spiked-in" genes removed.
scatter plot of intensity versus pooled log variance .
histograms of p values .
scatter plot of intensity versus log pooled variance.
histograms of p values .
orange columns in histograms correspond to the proportion of p values ≤  <dig> .
methods
see additional file  <dig> for software versions and availability.
fold-change
the fold-change  ratio is typically represented on the log <dig> scale:
fc=x¯1−x¯ <dig>
where x¯ <dig> and x¯ <dig> are the means of the two groups' logged expression values.
independent t-test
the independent t-statistic provides a standardized estimate of differential expression according to the following formula:
t=11n1+1n2x¯1−x¯2σ^pool
where
σ^pool2=σ^12+σ^22n1+n2− <dig>
σ^i <dig> and ni are the variance and number of replicates for the ith group, respectively.
σ^pool <dig> is the weighted average of the two groups gene-specific variance.
the associated probability under the null hypothesis is calculated by reference to the t distribution with n <dig> + n <dig> -  <dig> degrees of freedom.
local pooled error  z-statistic
jain et al.
's  local pooled error test is based on the model that genes with similar observed intensities have similar true variances and is calculated as follows:
z=median1−median2σpool
where
σpool2=π <dig>
the fold change is calculated using a difference of medians and the gene specific σi <dig> is obtained from a calibration curve derived from pooling the variance estimates of each gene's replicate expression measures with the variance estimates of other genes with similar expression values.
n <dig> and n <dig> are the sample sizes for the two groups.
the variance is adjusted by π/ <dig> due to the increased variability of using medians when calculating the fold change.
the associated probability of the z-statistic under the null hypothesis is calculated by reference to the standard normal distribution.
statistical tests with empirical bayes variance estimators
the following three tests use an empirical bayes method to estimate error associated with differential expression and use a statistical test identical in form to the independent t-statistic shown in equation  <dig>  the tests use a posterior variance, σ˜ <dig>  in place of the pooled variance, σpool <dig>  in the t-test.
using bayes rule, the posterior variance, σ˜ <dig>  for each gene becomes a combination of the observed gene-specific error and an estimate obtained from the prior distribution according to the following formula.
σ˜2=σ^g2+σ02dg+d <dig>
dg are the degrees of freedom and σ^g <dig> are the gene specific variances obtained from the experimental data.
d <dig> and σ <dig> are the prior degrees of freedom and variance estimates respectively  <cit> .
the three methods differ in how they estimate the parameters d <dig> and σ <dig>
limma t-statistic
smyth  <cit>  developed the hierarchical model of lönnstedt and speed  <cit>  into a general linear model that is more easily applied to microarray data.
it is implemented using the r statistical language in the limma bioconductor package  <cit> .
this method is based on a model where the variances of the residuals vary from gene to gene and are assumed to be drawn from a chi-square distribution.
the linear model is as follows:
e = xαg
where yg is the expression summary values for each gene across all arrays, x is a design matrix of full column rank and αg is the coefficient vector.
the contrasts of coefficients that are of interest for a particular gene g are defined as βg = ctαg.
although this approach is able to analyse any number of contrasts, we examine two sample comparisons only so βg can be defined as the log fold change x^1−x^ <dig> .
the contrast estimators, βg, are assumed to be normal and the residual variances, sg <dig>  are assumed to follow a scaled chi-square distribution.
under this hierarchical model the posterior variance, s˜g <dig>  takes the form:
s˜g2=d0s02+dgsg2d0+dg
where d <dig> and s <dig> are the prior degrees of freedom and variance and dg and sg <dig> are the experimental degrees of freedom and the sample variance for a particular gene g, respectively.
because we examine only two sample comparisons, dg will always be equal to n -  <dig> where n is the total number of replicates.
the two prior parameters, d <dig> and s <dig>  can be interpreted as the degrees of freedom and variance of the prior distribution respectively.
the prior parameters are estimated by fitting the logged sample variances to a scaled f distribution.
the moderated t-statistic is defined by:
t˜=11n1+1n2β^s˜g
the associated probability of the moderated t-statistic for the two sample case under the null hypothesis is calculated by reference to the t-distribution with d <dig> + dg degrees of freedom.
brb t-statistic
like the limma model, wright and simon's  <cit>  rvm t-statistic, labelled brb t-statistic after the brb software, assumes a random variance model where the variances of the residuals vary from gene to gene and are random selections from a single prior inverse gamma distribution, which is a more general form of the chi-square.
the two models differ, however, in the method of the prior parameter estimations.
in the brb model, the estimate of the gene specific error corresponds to the sample variance of the gene's replicate expression measures; the prior estimate corresponds to the mean of the inverse gamma prior distribution fitted to the sample variances.
the prior distribution is pp, where the marginal p is the scaled inverse gamma and the conditional distribution p is normal.
that is, the random error associated with each gene's measurement is assumed to be distributed normally with a mean of zero and the variance of the intensity distribution is assumed to be randomly drawn from a prior inverse gamma distribution whose parameters need to be estimated.
the posterior variance for each gene is a weighted average of its observed sample variance and the mean of the prior variance distribution.
the shape  and scale  parameters for the prior are estimated by fitting the entire collection of observed variances to an inverse gamma distribution.
rocke  <cit>  has discussed a method of moments approach to estimate the parameters while wright and simon  <cit>  have advocated a maximum likelihood approach to fit ab * σ^ <dig> to an f distribution, where n is the total number of replicates and k is the number of groups.
the maximum likelihood approach was found to have lower variability  and was the method used in this study.
the brb t-statistic for a particular gene takes the form:
t˜=11n1+1n2x¯1−x¯2σ˜
where the posterior variance is
σ˜2=σ^2+2a−1n−2+2a
σ^ <dig> is the gene specific sample variance with degrees of freedom of n -  <dig>  the prior variance estimate-1) is the mean of the fitted inverse gamma distribution with degrees of freedom of 2a.
the limma prior parameters d <dig> and s <dig> are equivalent to the brb prior parameters 2a and - <dig>  a large number of replicates will give increased weight to the observed variance while a high value for a gives increased weight to the prior variance.
a large a means that the inverse gamma is highly peaked making it more likely that the true variance for any particular gene is close to - <dig>  in the two-sample case, the associated probability of the t-statistic under the null hypothesis is calculated by reference to the t-distribution with n -  <dig> + 2a degrees of freedom.
cybert t-statistic
the baldi and long  <cit>  regularized t-statistic, labelled the cybert statistic after the cyber-t software, combines elements of both the lpe and the two previous empirical bayes methods.
like lpe, variance estimates are obtained by pooling across genes with similar intensities.
a gene specific running average of the variance is used for the prior variance estimate.
the running average variance is calculated by averaging the variance of a particular gene and the variances of the z number of genes on either side of that gene in the ordered distribution.
the total window size  for the running average is equal to 2z +  <dig>  accordingly, for pair, , ordered by xi, the running average is generated using:
σ0i2=σi−z2+…+σi2+…+σi+z2w
the degrees of freedom  for the running average is a user supplied theoretical value based on the belief in the quality of the running average variance.
the heuristic suggested by baldi and long  <cit>  of using three times the number of replicates was used for the running average degrees of freedom in the present study.
the cybert t-statistic takes the form:
t˜=11n1+1n2x¯1−x¯2σ˜
where the posterior variance is
σ˜2=σ^2+v0σ02n−2+v <dig>
x¯ <dig> and x¯ <dig> are the sample means of the two groups.
σ^ <dig> is the gene specific sample variance, n is the total number of replicates  per gene across treatment classes, vo is the degrees of freedom of the prior, and σ <dig> is the gene specific running average variance.
the associated probability of the t-statistic under the null hypothesis is calculated by reference to the t-distribution with n -  <dig> + v <dig> degrees of freedom.
RESULTS
descriptive comparisons
figures 1a and 1c show the scatter plots of log pooled variance estimates  versus average intensity across replicates of both groups for simulated and affymetrix data, respectively.
the t-statistic shows the unadjusted variance values before the small n methods are applied and can be used as a baseline for comparison.
the scatter plots for the complete set of latin square comparisons and the null hypothesis data with low and medium variances are comparable to these plots .
with two exceptions  the four small n statistical tests reduced the variability of the pooled variance estimates relative to the unadjusted estimates of the t-test for all data sets; the range of points on the y-axis in figures 1a and 1c is diminished along the smoothing curve.
lpe in particular compressed the variance because the gene-specific variance estimates for each group are generated from a calibration function, which maps a single associated variance value per group to each intensity value.
the pooled variance estimates of the brb and limma method with the latin square data normalized by mas  <dig>  showed little compression.
the degrees of freedom  assigned to the gene-specific variance is much greater than the degrees of freedom  assigned to the prior variance estimate .
consequently the adjustment applied by the brb method is negligible.
in contrast, the degrees of freedom values for the data shown in figure  <dig> for rma, gcrma, dchip, vsn, and lmgene were  <dig> ,  <dig> ,  <dig> ,  <dig>  and  <dig> for brb and  <dig> ,  <dig> ,  <dig> ,  <dig> , and  <dig>  for limma, respectively.
the brb and limma methods produced similar results for all expression summary methods with the exception of the affymetrix data normalized with vsn and lmgene which use a generalized log transformation.
although the brb and limma methods generated similar prior variance parameters, the former generated much higher prior degrees of freedom than the latter.
the brb prior degrees of freedom were on average  <dig> times higher than the limma degrees of freedom  for the vsn expression summary method and  <dig> times higher  for the lmgene expression summary method.
for these data, the degrees of freedom for the sample variance were  <dig> so the posterior variance of the brb method, which is a weighted average of the prior and sample variance, deviates only slightly from the prior variance across gene specific intensity .
the lmgene and vsn normalization methods both produced a narrow peaked sample variance distribution consisting of very small variances which will cause the brb prior degrees of freedom parameter to be excessively large.
the limma prior degrees of freedom parameter has a more appropriate value because it is calculated using the log of the sample variances which do not show such a narrow distribution .
the brb t-statistic may not be appropriate for lmgene and vsn normalized data due to this issue.
figures 1b and 1d show histograms of the p value distributions of the simulated and hgu133a data, respectively.
under the null hypothesis, the distribution of the p values should follow a uniform distribution.
deviation from a uniform distribution indicates lack of correspondence between the data and the assumptions of the statistical model.
the histograms for the complete set of latin square experiments and the null hypothesis data with low and medium variances are comparable to these plots .
the p value distributions for the t-statistic, brb t-statistic and the limma moderated t-statistic followed the theoretically expected uniform distribution for all simulation models; the cybert t-statistic deviated slightly from uniformity .
for the hgu133a data, all expression summary methods and statistical tests, with the exception of lpe, showed a rightward skew in their p value distributions and a greater than expected number of small p values .
unlike the other statistical tests, the p value distribution for the lpe z-statistic deviated substantially from uniformity for all data sets and expression summary methods.
performance
a standard method of comparing the performance of different statistical methods is to use a partial auc of an roc graph calculated using only the lower end of the ordered distributions of p values.
the pauc has a value between  <dig>  and  <dig> .
these metrics were calculated for each of the  <dig> iterations of the simulated "spike-in" data sets and for the  <dig> comparisons of each of the two latin square data sets.
a particular false positive rate that applies commonly across the different distributions is often chosen.
this has the advantage of ensuring that the same number of data points are used in the comparison but has the disadvantage that it ignores the actual p values of the data points being compared.
an alternative is to use a cutoff based on a specific p value rather than using only the rankings which is more in line with how the test results would be used in applied contexts.
we present results for the latter below and for the former in additional file  <dig>
the average pauc scores and rankings for all methods and data sets are shown in table  <dig>  the average true and false positive rates are shown in tables  <dig> and  <dig> respectively.
rankings from best to worst  of average pauc scores  of the six tests of differential expression.
the average pauc for each score is in parentheses.
the lmgene package was used to apply the loess normalization and glog transformation to the cdna data.
bc indicates that background correction was applied to the data before normalization.
rankings from best to worst  of average true positive rate  scores of the six tests of differential expression.
the average tpr for each score is in parentheses.
the lmgene package was used to apply the loess normalization and glog transformation to the cdna data.
bc indicates that background correction was applied to the data before normalization.
average false positive rates  of the six tests of differential expression.
the lmgene package was used to apply the loess normalization and glog transformation to the cdna data.
bc indicates that background correction was applied to the data before normalization.
the three empirical bayes variance estimation methods, cybert, brb, and limma, performed comparably to fc with almost all low variance data sets .
they appreciably outperformed fc, however, with the liver versus pooled cdna data set, which is characterized by higher variance than the other data sets due to the use of mrna from multiple tissues, and mas  <dig>  normalized latin square data, which also has higher variability than the other normalization methods .
greater variability randomly introduces a higher number of large fold changes between groups simply because of the greater range of expression intensities across replicate values.
fc, which ignores variability, is less able to separate the higher number of random effects from the "spiked-in" effects than those methods that take the variability into account.
for these latter data, the empirical bayes variance estimation tests of the loess/glog transformed data performed especially well .
these results show the limitations of using fc with highly variable data  <cit> .
based on average ranks, the cybert t-statistic, brb t-statistic, and the limma moderated t-statistic performed best across all data sets .
the top performers with the latin square data were able to find almost all of the differentially expressed genes and had comparable pauc scores.
the three empirical bayes variance estimation tests were in the top four performers with the low variability liver <dig> versus liver <dig> cdna data and the top three with the higher variability liver versus pooled cdna data.
for the simulated data, the limma t-statistic performed best with the local variance model, tied with the cybert t-statistic and the brb t-statistic for the inverse gamma variance model, and tied with fc and the brb t-statistic for best performance for the common variance model.
the t-statistic and lpe z-statistic were the worst performers overall.
these results are similar to those generated by the paucs with the same number of data points across tests .
in both cases, the average performances of the three empirical bayes tests were ranked in the top three, fc was ranked fourth, and the lpe and t-test were ranked last .
with the common variance model, the high fold change transcripts were the most easily identified by all methods and it was within the medium fold change range that there was the best differentiation between the methods; this was also seen with lesser effect with the local variance model .
concentration effects
we examined the data in greater detail by analyzing the true positive rates  and the false positive rates  separately, conditioned on spike-in concentration and using a p <  <dig>  threshold.
the differential expression tests found most of the latin square "spiked-in" genes at the mid to high concentrations with the hgu133a data .
there was, however, a drop in the tpr at low concentrations found across all combinations of expression summary methods and differential expression tests.
fc had a very low tpr across all concentrations with the mas  <dig>  data.
the cybert, brb, and limma tests performed as well or better than the lpe and t-test methods in limiting this effect.
for the simulated local variance model, the brb t-statistic and limma moderated t-statistic yielded substantially higher tpr values across all concentration levels, although there was a decrease in the tpr among low concentrations for all statistical tests .
for the common variance and inverse gamma simulated models, the brb t-statistic, the cybert t-statistic, the limma moderated t-statistic, and fc yielded similarly high tprs across all concentrations; the t-test and the lpe z-test yielded lower tpr values across all concentrations .
the effects of expression intensity on the fpr with the hgu133a latin square data varied across differential expression tests .
lpe showed a decreasing fpr at low concentrations, with the exception of the mas  <dig>  data.
by contrast, the fpr for the t-statistic and cybert t-statistic were more uniform across concentration.
fc showed an unstable fpr across intensity with all normalization methods.
the fc fpr curve closely follows the smoothing curve of pooled variance versus intensity from figures 1a and 1c for each normalization method.
the fc fpr increases as the average variability at a particular intensity is high and decreases where the average variability is low.
for example, with the hgu133a data normalized with dchip, there are generally many genes with high variance at low intensities and as the intensity increases the gene variances tend to decrease.
the corresponding fpr follows a similar curve, with a high fpr for genes with low intensities and a low fpr for genes with high intensities.
the fc method is unable to distinguish between random and biological differences in gene expression when the variability of the replicate measurements is high.
the cybert t-statistic and limma moderated t-statistic also showed a similar, although weaker, pattern as fc in its fpr behaviour across concentration.
this is due to the pooling of each gene-specific variance with a single prior error estimate.
genes with variances higher than the prior estimate have their variances reduced and genes with low variances have their variances increased.
this causes the denominator of the t-statistic formula to decrease where variance is generally high, making it more likely that the null hypothesis is rejected.
the opposite effect will happen when variance is relatively low.
the two tests had a stable fpr with the vsn and lmgene normalized hgu133a data as the use of a glog transformation stabilized the variance more effectively across intensity.
for the simulated local variance model, fc, the brb t-statistic and limma moderated t-statistic  showed a similar tracking of the pooled variance smoothing curve of figure 1a.
thus, the limma and brb tests' superior pauc performance  and higher true positive rate at low concentrations  for the local variance model is mitigated by its high fpr among low concentrations; for the local variance model, the t-statistic and the cybert t-statistic alone yielded the expected  <dig>  false positive rate across all concentrations.
the lpe z-statistic yielded the lowest fpr across all concentrations.
CONCLUSIONS
consistent with the consensus among analysts that "borrowing strength" across genes for estimating differential expression random error is desirable  <cit> , the three empirical bayes tests  performed best.
the t-statistic and the lpe z-statistic had reduced power in comparison.
fc was comparable to the three empirical bayes methods with data of low variability  but performed poorly with data characterized by high variability .
this is in accordance with other studies that have shown the standard t-statistic to be suboptimal  <cit>  and that fc is an inappropriate inferential test  <cit> .
lpe introduced a bias that caused it to generate fewer false positives than expected.
this bias appears to be due to a theoretical weakness and an incorrect assumption intrinsic to the test as originally proposed.
murie and nadon  <cit>  have shown that the adjustment to the lpe standard error, which depends on an asymptotic proof, is overly conservative with sample sizes smaller than  <dig>  moreover, empirical evidence suggests that the lpe's assumption of similar error variability for genes of similar expression intensity is incorrect, which leads the lpe test to overestimate p-values for low variability genes and to underestimate them for high variability genes.
recent modifications to the test have addressed this issue for microarray gene expression  <cit>  and for mass spectrometry proteomics data  <cit> .
indeed, potential differences in high and low variability genes between our studies and that of kooperberg et al.
<cit>  may explain why they concluded that the lpe test was overly liberal in contrast to our current findings and those of murie and nadon that the test is overly conservative.
an important difference in our study between the cybert test and the other two empirical bayes tests , was that only the former maintained a stable false positive error rate with data that showed an unstable variance across intensity, such as the experimental and local normal simulated data.
this suggests that these latter tests may also generate too many false positives among some biological data sets which, depending on pre-processing, often show high variance at low concentrations.
as designed in our simulations and has been argued generally, however, this variance/intensity relationship is an artifact of the inappropriate log transformation for low intensity expression values and the generalized log  transformation is able to stabilize the error variance across the entire concentration range  <cit> .
the glog transformation, as implemented with vsn and lmgene, lowers this potentially high fpr among low abundance transcripts in biological data sets while maintaining a relatively high tpr.
background correction when applied to cdna data can generate negative values which cannot be logged and requires that the negative data points be either discarded or set to a loggable value.
moreover, we found, as did qin and kerr  <cit> , that using a log transformation after applying background correction in conjunction with loess normalization either reduced performance, as with the liver vs pool data, or had a mixed effect on performance, as with the liver vs liver data.
we recommend the use of a glog transformation  when applying background correction which can transform negative values and thus avoids the loss of information.
the combination of glog and loess normalization produced comparable results for both cdna data sets with and without background correction and were amongst the top performers of all normalization methods.
an attractive alternative might be to use background correction methods which avoid negative values and stabilize the variance as a function of intensity as described in ritchie et al.
<cit> .
the performance of the statistical tests was most similar using the hgu133a data set.
with the exception of mas  <dig>  with fc, all of the statistical tests found most of the "spiked-in" genes readily, illustrating the limitations of these empirical data sets for comparing performances of statistical test as they would be applied in experiments with biological variation.
the hgu133a data set, which used technical replicates, is comparable to our simulated data with low variance, where in general the methods do not show large performance differences.
the hgu <dig> data had higher variability than the hgu133a data and also showed greater performance differences between methods.
most biological experiments will be closer to the medium and the high variance simulated data sets where differences in performance are more apparent.
the three empirical bayes methods in our study consistently performed well across normalization methods and simulated data designed to mirror the greater variability that is observed in most biological experiments.
for this reason, we expect that our results apply beyond the narrow methodological confines of our study.
nonetheless, more definitive tests of the algorithms' merits will be made possible by future spike-in experiments which are anticipated to incorporate biological variability  <cit> .
there is also a need for a similar comparative study among statistical tests which are suited to study designs that are more complex than the two-independent sample design in our study.
two of the approaches we evaluated  and at least one other that we did not examine  <cit>  have more general capabilities suited to this type of comparative study.
authors' contributions
cm conducted the computer simulations, performed the data analysis, participated in research design, and drafted the manuscript.
ow and al participated in discussion of the research, investigation of the statistical methods, and analysis of data.
rn conceived of the study, supervised research design, interpretation of results, and drafting of the manuscript.
supplementary material
additional file 1
mva plots of cdna data.
click here for file
additional file 2
source code availability.
click here for file
additional file 3
histograms of logged and unlogged variances for limma and brb statistical tests.
click here for file
additional file 4
pauc scores using a 5% false positive rate.
click here for file
additional file 5
trellis plots of performance criteria for simulated data.
click here for file
additional file 6
total average rankings of statistical tests for pauc and tpr.
click here for file
acknowledgements
we thank mathieu miron, jacek majewski, and amy norris for their helpful suggestions.

 microarrays provide large-scale comparative gene expression profiles between biological samples by detecting differential expression for thousands of genes in parallel. typically, systematic error in the measurements is removed at the background correction and normalization steps, and is followed by statistical testing. in the most common type of study, statistical testing produces a list of genes that are differentially expressed across two or more classes. extensive research has shown that choice of pre-processing methods designed to correct for bias in the measurements can have a substantial impact on rank ordering of gene expression fold-change estimates for both cdna and oligonucleotide data and the affycomp website. one major disadvantage of fc estimates, however, is that they do not take the variance of the samples into account. this is especially problematic because variability in gene expression measurements is partially gene-specific, even after the variance has been stabilized by data transformation. there is consensus in the statistical community that statistical tests of differential expression are preferred over fc for inference. one advantage is that they standardize differential expression by dividing fc measurements by their associated standard error, rescaling fcs to a common metric. moreover, associated output such as p values, effect sizes, and confidence intervals can be used for various purposes such as false positive control and meta-analysis. although microarray studies with sample sizes of five or more observations per class are becoming increasingly common, cost considerations and the need for small scale initial studies mean that many studies are conducted with smaller sample sizes. the use of classical statistical tests such as the t-test is sub-optimal for small sample studies, however, because of low statistical power for detection of differential expression. this has led to the development of small sample size error estimation methods which borrow information from the entire data set or from a subset of the data to improve error estimation accuracy and precision. there has been a number of studies analysing the performance of statistical tests applied to microarray data but few have used data where the differentially expressed genes are known in advance. qin and kerr compared normalization methods and test statistics using data sets with six known differentially expressed genes and showed that the standard t-statistic performed worse than those that used variance estimates that incorporated information from all genes. sioson et al. compared statistical methodologies of two software applications using qrt-pcr of a subset of genes and chen et al. used a consensus of differentially expressed genes across statistical methods to analyze performance. none of these studies, however, evaluated statistical tests designed for small sample size experiments. a number of studies have examined small n tests. kooperberg et al. compared the performances of various 2-group statistical tests with empirical and simulated data. distributions of p-values were generated based on data from the same experimental group or data from different experimental groups which were known to differ. the high performing tests generated p-value distributions consistent with a null distribution in the former case and produced the largest number of small p-values in the latter case. using the limma method as an exemplar, they concluded that an empirical bayes approach to statistical testing provided good power while accurately controlling false positive rates for small n microarray experiments. cui and churchill reviewed a number of statistical tests for use with microarray data, including one small n test, but provided no comparative analysis. tong and wang used simulated data to explore the theoretical properties of shrinkage methods of estimating variance and found that they outperformed tests that use only the sample variance. hu and wright and xie et al. both compared a number of statistical tests, including several small n tests, using a consensus list of differentially expressed genes from all methods. hu and wright found that, based on false discovery rates, tests that model the variance/intensity relationship and use variance estimates generated with information from all genes performed the best. xie et al showed that there was comparability of results for only a few of the methods but the lönnstedt and speed b-statistic, which is monotonically equivalent to the limma and brb t-statistics, had the lowest false positive rate. jefferey et al. looked at the use of statistical tests, including several small n tests, in feature selection for group classification; the results varied greatly across gene list and sample size but the empirical bayes t-statistic performed well across all sample sizes. most of these studies compared methods using comparability of results, estimated false discovery, family wise error rates, or p-value distributions to assess performance. in this paper we focus on each method's ability to detect genes which have been spiked at different concentrations across samples. we use simulated data, two publicly available empirical affymetrix latin-square data sets, and two cdna data sets in which the differentially expressed genes are known to assess the relative performances of fc, t-statistic, and four small sample size statistical tests which all borrow strength from other genes in different ways. for consistency, we refer to the various tests by their software implementation. the local pooled error method z-statistic is based on pooling errors for genes with similar intensities. the three empirical bayes methods were: the brb t-statistic which combines gene-specific error estimates with a common error estimate obtained from the distribution of the variances across all genes; the cybert regularized t-statistic which combines gene-specific error estimates with a local pooled error estimate based on genes with similar intensities; and the limma method which combines a fitted linear model of gene expression data with a variance estimate into a moderated t-statistic. all methods gain degrees of freedom over the standard t-test and in theory should provide greater sensitivity with no loss in specificity. data latin square spike-in data the hgu133a and hgu latin square data sets are based on a × latin square design of "spiked-in" transcripts with three replicates for each group. the concentrations for the "spiked-in" transcripts were doubled for each consecutive group. we compared performance of the statistical tests for detecting two-fold differential expression after the data had been normalized by six popular expression summary algorithms: mas, dchip, rma, gcrma, vsn, and lmgene. data were analysed on the log scale for the first four methods and on a generalized log scale for the latter two methods. cdna spike-in data two spotted cdna experiments were also used to compare statistical test performance. the first experiment compared two aliquots from the same mouse liver rna sample, and the second compared one rna sample from mouse liver to a pooled sample of five different mouse tissues. in both experiments six transcripts were "spiked-in" at a three-fold difference in concentration between the two groups. all cdna slides measured each "spike-in" transcript at different well locations for a total of "spikes" in each experiment. the data were normalized using the lowess normalization of lmgene, with and without a glog transformation and background correction, and the vsn method. the lmgene lowess normalization method is equivalent to the intensity-dependent normalization used in yang et al.. mva plots and their corresponding lowess fits for each cdna array are shown in additional file simulated data two limitations of the above spike-in data sets make them less challenging than most biologically motivated studies: variation across technical replicates is lower than that typically observed across biological replicates and many biological effects of interest may be smaller than two-fold. also, there is only a small number of "spiked-in" transcripts in each latin square experiment, which can lead to increased variability in the algorithm performance metrics. simulated data were generated to offset these limitations: null hypothesis and "spike-in" data were used to examine false positive rates, genes incorrectly labeled as differentially expressed, and true positive rates, genes correctly labeled as differentially expressed. the "spike-in"data incorporated varying degrees of variability between replicate measurements and genes were differentially expressed across a range of fold changes to offset the limitations of the latin square data sets. each simulation experiment consisted of iterations of ten thousand genes for each of two groups with or replicates per gene. gene intensity values were generated by randomly selecting from a n distribution, where u is the "true" expression intensity distribution) for a single gene in log scale and σ is the random error associated with the gene's expression measurement. in the null hypothesis data, each gene-specific u remains constant across both groups; in the "spike-in" data, μ was upwardly regulated by a specified fold change for genes in the treatment group. for each "spiked-in" gene, μ was assigned a value from a set of log values ranging from to in increments. the log fold change applied to the "spiked-in genes" was implemented by an additive adjustment to their u in the treatment group chosen from a range of values from to in increments. there are "spiked-in" genes at each of the concentration values and each of those genes at a particular concentration value is upwardly regulated by one of the different fold changes. the random error of each gene is derived from one of three variance models: common, inverse gamma, and local. in the common error model, the population variance of the random error term is a constant across all genes. in the inverse gamma variance model, c is drawn randomly from an inverse gamma distribution, ig, where is the shape parameter and b is the scale parameter. in the local variance model, ϵ is a combination of additive and proportional components, as outlined in rocke and durbin. the expression intensity column shows how the replicate intensity values were generated for each transcript. the random error column shows the method of calculating specific ϵ values. the parameters column shows the values that were used to generate high, medium, and low variances, respectively. the parameters of each variance model were varied to create data sets from low variance to medium and high variance; c and d show hgu133a data with the "spiked-in" genes removed. scatter plot of intensity versus pooled log variance. histograms of p values. scatter plot of intensity versus log pooled variance. histograms of p values. orange columns in histograms correspond to the proportion of p values ≤. methods see additional file for software versions and availability. fold-change the fold-change ratio is typically represented on the log scale: fc=x¯1−x¯ where x¯ and x¯ are the means of the two groups' logged expression values. independent t-test the independent t-statistic provides a standardized estimate of differential expression according to the following formula: t=11n1+1n2x¯1−x¯2σ^pool where σ^pool2=σ^12+σ^22n1+n2− σ^i and ni are the variance and number of replicates for the ith group, respectively. σ^pool is the weighted average of the two groups gene-specific variance. the associated probability under the null hypothesis is calculated by reference to the t distribution with n + n - degrees of freedom. local pooled error z-statistic jain et al.'s local pooled error test is based on the model that genes with similar observed intensities have similar true variances and is calculated as follows: z=median1−median2σpool where σpool2=π the fold change is calculated using a difference of medians and the gene specific σi is obtained from a calibration curve derived from pooling the variance estimates of each gene's replicate expression measures with the variance estimates of other genes with similar expression values. n and n are the sample sizes for the two groups. the variance is adjusted by π/ due to the increased variability of using medians when calculating the fold change. the associated probability of the z-statistic under the null hypothesis is calculated by reference to the standard normal distribution. statistical tests with empirical bayes variance estimators the following three tests use an empirical bayes method to estimate error associated with differential expression and use a statistical test identical in form to the independent t-statistic shown in equation the tests use a posterior variance, σ˜ in place of the pooled variance, σpool in the t-test. using bayes rule, the posterior variance, σ˜ for each gene becomes a combination of the observed gene-specific error and an estimate obtained from the prior distribution according to the following formula. σ˜2=σ^g2+σ02dg+d dg are the degrees of freedom and σ^g are the gene specific variances obtained from the experimental data. d and σ are the prior degrees of freedom and variance estimates respectively. the three methods differ in how they estimate the parameters d and σ limma t-statistic smyth developed the hierarchical model of lönnstedt and speed into a general linear model that is more easily applied to microarray data. it is implemented using the r statistical language in the limma bioconductor package. this method is based on a model where the variances of the residuals vary from gene to gene and are assumed to be drawn from a chi-square distribution. the linear model is as follows: e = xαg where yg is the expression summary values for each gene across all arrays, x is a design matrix of full column rank and αg is the coefficient vector. the contrasts of coefficients that are of interest for a particular gene g are defined as βg = ctαg. although this approach is able to analyse any number of contrasts, we examine two sample comparisons only so βg can be defined as the log fold change x^1−x^. the contrast estimators, βg, are assumed to be normal and the residual variances, sg are assumed to follow a scaled chi-square distribution. under this hierarchical model the posterior variance, s˜g takes the form: s˜g2=d0s02+dgsg2d0+dg where d and s are the prior degrees of freedom and variance and dg and sg are the experimental degrees of freedom and the sample variance for a particular gene g, respectively. because we examine only two sample comparisons, dg will always be equal to n - where n is the total number of replicates. the two prior parameters, d and s can be interpreted as the degrees of freedom and variance of the prior distribution respectively. the prior parameters are estimated by fitting the logged sample variances to a scaled f distribution. the moderated t-statistic is defined by: t˜=11n1+1n2β^s˜g the associated probability of the moderated t-statistic for the two sample case under the null hypothesis is calculated by reference to the t-distribution with d + dg degrees of freedom. brb t-statistic like the limma model, wright and simon's rvm t-statistic, labelled brb t-statistic after the brb software, assumes a random variance model where the variances of the residuals vary from gene to gene and are random selections from a single prior inverse gamma distribution, which is a more general form of the chi-square. the two models differ, however, in the method of the prior parameter estimations. in the brb model, the estimate of the gene specific error corresponds to the sample variance of the gene's replicate expression measures; the prior estimate corresponds to the mean of the inverse gamma prior distribution fitted to the sample variances. the prior distribution is pp, where the marginal p is the scaled inverse gamma and the conditional distribution p is normal. that is, the random error associated with each gene's measurement is assumed to be distributed normally with a mean of zero and the variance of the intensity distribution is assumed to be randomly drawn from a prior inverse gamma distribution whose parameters need to be estimated. the posterior variance for each gene is a weighted average of its observed sample variance and the mean of the prior variance distribution. the shape and scale parameters for the prior are estimated by fitting the entire collection of observed variances to an inverse gamma distribution. rocke has discussed a method of moments approach to estimate the parameters while wright and simon have advocated a maximum likelihood approach to fit ab * σ^ to an f distribution, where n is the total number of replicates and k is the number of groups. the maximum likelihood approach was found to have lower variability and was the method used in this study. the brb t-statistic for a particular gene takes the form: t˜=11n1+1n2x¯1−x¯2σ˜ where the posterior variance is σ˜2=σ^2+2a−1n−2+2a σ^ is the gene specific sample variance with degrees of freedom of n - the prior variance estimate-1) is the mean of the fitted inverse gamma distribution with degrees of freedom of 2a. the limma prior parameters d and s are equivalent to the brb prior parameters 2a and - a large number of replicates will give increased weight to the observed variance while a high value for a gives increased weight to the prior variance. a large a means that the inverse gamma is highly peaked making it more likely that the true variance for any particular gene is close to - in the two-sample case, the associated probability of the t-statistic under the null hypothesis is calculated by reference to the t-distribution with n - + 2a degrees of freedom. cybert t-statistic the baldi and long regularized t-statistic, labelled the cybert statistic after the cyber-t software, combines elements of both the lpe and the two previous empirical bayes methods. like lpe, variance estimates are obtained by pooling across genes with similar intensities. a gene specific running average of the variance is used for the prior variance estimate. the running average variance is calculated by averaging the variance of a particular gene and the variances of the z number of genes on either side of that gene in the ordered distribution. the total window size for the running average is equal to 2z + accordingly, for pair,, ordered by xi, the running average is generated using: σ0i2=σi−z2+…+σi2+…+σi+z2w the degrees of freedom for the running average is a user supplied theoretical value based on the belief in the quality of the running average variance. the heuristic suggested by baldi and long of using three times the number of replicates was used for the running average degrees of freedom in the present study. the cybert t-statistic takes the form: t˜=11n1+1n2x¯1−x¯2σ˜ where the posterior variance is σ˜2=σ^2+v0σ02n−2+v x¯ and x¯ are the sample means of the two groups. σ^ is the gene specific sample variance, n is the total number of replicates per gene across treatment classes, vo is the degrees of freedom of the prior, and σ is the gene specific running average variance. the associated probability of the t-statistic under the null hypothesis is calculated by reference to the t-distribution with n - + v degrees of freedom.  descriptive comparisons figures 1a and 1c show the scatter plots of log pooled variance estimates versus average intensity across replicates of both groups for simulated and affymetrix data, respectively. the t-statistic shows the unadjusted variance values before the small n methods are applied and can be used as a baseline for comparison. the scatter plots for the complete set of latin square comparisons and the null hypothesis data with low and medium variances are comparable to these plots. with two exceptions the four small n statistical tests reduced the variability of the pooled variance estimates relative to the unadjusted estimates of the t-test for all data sets; the range of points on the y-axis in figures 1a and 1c is diminished along the smoothing curve. lpe in particular compressed the variance because the gene-specific variance estimates for each group are generated from a calibration function, which maps a single associated variance value per group to each intensity value. the pooled variance estimates of the brb and limma method with the latin square data normalized by mas showed little compression. the degrees of freedom assigned to the gene-specific variance is much greater than the degrees of freedom assigned to the prior variance estimate. consequently the adjustment applied by the brb method is negligible. in contrast, the degrees of freedom values for the data shown in figure for rma, gcrma, dchip, vsn, and lmgene were,,, and for brb and,,,, and for limma, respectively. the brb and limma methods produced similar results for all expression summary methods with the exception of the affymetrix data normalized with vsn and lmgene which use a generalized log transformation. although the brb and limma methods generated similar prior variance parameters, the former generated much higher prior degrees of freedom than the latter. the brb prior degrees of freedom were on average times higher than the limma degrees of freedom for the vsn expression summary method and times higher for the lmgene expression summary method. for these data, the degrees of freedom for the sample variance were so the posterior variance of the brb method, which is a weighted average of the prior and sample variance, deviates only slightly from the prior variance across gene specific intensity. the lmgene and vsn normalization methods both produced a narrow peaked sample variance distribution consisting of very small variances which will cause the brb prior degrees of freedom parameter to be excessively large. the limma prior degrees of freedom parameter has a more appropriate value because it is calculated using the log of the sample variances which do not show such a narrow distribution. the brb t-statistic may not be appropriate for lmgene and vsn normalized data due to this issue. figures 1b and 1d show histograms of the p value distributions of the simulated and hgu133a data, respectively. under the null hypothesis, the distribution of the p values should follow a uniform distribution. deviation from a uniform distribution indicates lack of correspondence between the data and the assumptions of the statistical model. the histograms for the complete set of latin square experiments and the null hypothesis data with low and medium variances are comparable to these plots. the p value distributions for the t-statistic, brb t-statistic and the limma moderated t-statistic followed the theoretically expected uniform distribution for all simulation models; the cybert t-statistic deviated slightly from uniformity. for the hgu133a data, all expression summary methods and statistical tests, with the exception of lpe, showed a rightward skew in their p value distributions and a greater than expected number of small p values. unlike the other statistical tests, the p value distribution for the lpe z-statistic deviated substantially from uniformity for all data sets and expression summary methods. performance a standard method of comparing the performance of different statistical methods is to use a partial auc of an roc graph calculated using only the lower end of the ordered distributions of p values. the pauc has a value between and. these metrics were calculated for each of the iterations of the simulated "spike-in" data sets and for the comparisons of each of the two latin square data sets. a particular false positive rate that applies commonly across the different distributions is often chosen. this has the advantage of ensuring that the same number of data points are used in the comparison but has the disadvantage that it ignores the actual p values of the data points being compared. an alternative is to use a cutoff based on a specific p value rather than using only the rankings which is more in line with how the test results would be used in applied contexts. we present results for the latter below and for the former in additional file the average pauc scores and rankings for all methods and data sets are shown in table the average true and false positive rates are shown in tables and respectively. rankings from best to worst of average pauc scores of the six tests of differential expression. the average pauc for each score is in parentheses. the lmgene package was used to apply the loess normalization and glog transformation to the cdna data. bc indicates that background correction was applied to the data before normalization. rankings from best to worst of average true positive rate scores of the six tests of differential expression. the average tpr for each score is in parentheses. the lmgene package was used to apply the loess normalization and glog transformation to the cdna data. bc indicates that background correction was applied to the data before normalization. average false positive rates of the six tests of differential expression. the lmgene package was used to apply the loess normalization and glog transformation to the cdna data. bc indicates that background correction was applied to the data before normalization. the three empirical bayes variance estimation methods, cybert, brb, and limma, performed comparably to fc with almost all low variance data sets. they appreciably outperformed fc, however, with the liver versus pooled cdna data set, which is characterized by higher variance than the other data sets due to the use of mrna from multiple tissues, and mas normalized latin square data, which also has higher variability than the other normalization methods. greater variability randomly introduces a higher number of large fold changes between groups simply because of the greater range of expression intensities across replicate values. fc, which ignores variability, is less able to separate the higher number of random effects from the "spiked-in" effects than those methods that take the variability into account. for these latter data, the empirical bayes variance estimation tests of the loess/glog transformed data performed especially well. these results show the limitations of using fc with highly variable data. based on average ranks, the cybert t-statistic, brb t-statistic, and the limma moderated t-statistic performed best across all data sets. the top performers with the latin square data were able to find almost all of the differentially expressed genes and had comparable pauc scores. the three empirical bayes variance estimation tests were in the top four performers with the low variability liver versus liver cdna data and the top three with the higher variability liver versus pooled cdna data. for the simulated data, the limma t-statistic performed best with the local variance model, tied with the cybert t-statistic and the brb t-statistic for the inverse gamma variance model, and tied with fc and the brb t-statistic for best performance for the common variance model. the t-statistic and lpe z-statistic were the worst performers overall. these results are similar to those generated by the paucs with the same number of data points across tests. in both cases, the average performances of the three empirical bayes tests were ranked in the top three, fc was ranked fourth, and the lpe and t-test were ranked last. with the common variance model, the high fold change transcripts were the most easily identified by all methods and it was within the medium fold change range that there was the best differentiation between the methods; this was also seen with lesser effect with the local variance model. concentration effects we examined the data in greater detail by analyzing the true positive rates and the false positive rates separately, conditioned on spike-in concentration and using a p < threshold. the differential expression tests found most of the latin square "spiked-in" genes at the mid to high concentrations with the hgu133a data. there was, however, a drop in the tpr at low concentrations found across all combinations of expression summary methods and differential expression tests. fc had a very low tpr across all concentrations with the mas data. the cybert, brb, and limma tests performed as well or better than the lpe and t-test methods in limiting this effect. for the simulated local variance model, the brb t-statistic and limma moderated t-statistic yielded substantially higher tpr values across all concentration levels, although there was a decrease in the tpr among low concentrations for all statistical tests. for the common variance and inverse gamma simulated models, the brb t-statistic, the cybert t-statistic, the limma moderated t-statistic, and fc yielded similarly high tprs across all concentrations; the t-test and the lpe z-test yielded lower tpr values across all concentrations. the effects of expression intensity on the fpr with the hgu133a latin square data varied across differential expression tests. lpe showed a decreasing fpr at low concentrations, with the exception of the mas data. by contrast, the fpr for the t-statistic and cybert t-statistic were more uniform across concentration. fc showed an unstable fpr across intensity with all normalization methods. the fc fpr curve closely follows the smoothing curve of pooled variance versus intensity from figures 1a and 1c for each normalization method. the fc fpr increases as the average variability at a particular intensity is high and decreases where the average variability is low. for example, with the hgu133a data normalized with dchip, there are generally many genes with high variance at low intensities and as the intensity increases the gene variances tend to decrease. the corresponding fpr follows a similar curve, with a high fpr for genes with low intensities and a low fpr for genes with high intensities. the fc method is unable to distinguish between random and biological differences in gene expression when the variability of the replicate measurements is high. the cybert t-statistic and limma moderated t-statistic also showed a similar, although weaker, pattern as fc in its fpr behaviour across concentration. this is due to the pooling of each gene-specific variance with a single prior error estimate. genes with variances higher than the prior estimate have their variances reduced and genes with low variances have their variances increased. this causes the denominator of the t-statistic formula to decrease where variance is generally high, making it more likely that the null hypothesis is rejected. the opposite effect will happen when variance is relatively low. the two tests had a stable fpr with the vsn and lmgene normalized hgu133a data as the use of a glog transformation stabilized the variance more effectively across intensity. for the simulated local variance model, fc, the brb t-statistic and limma moderated t-statistic showed a similar tracking of the pooled variance smoothing curve of figure 1a. thus, the limma and brb tests' superior pauc performance and higher true positive rate at low concentrations for the local variance model is mitigated by its high fpr among low concentrations; for the local variance model, the t-statistic and the cybert t-statistic alone yielded the expected false positive rate across all concentrations. the lpe z-statistic yielded the lowest fpr across all concentrations.  consistent with the consensus among analysts that "borrowing strength" across genes for estimating differential expression random error is desirable, the three empirical bayes tests performed best. the t-statistic and the lpe z-statistic had reduced power in comparison. fc was comparable to the three empirical bayes methods with data of low variability but performed poorly with data characterized by high variability. this is in accordance with other studies that have shown the standard t-statistic to be suboptimal and that fc is an inappropriate inferential test. lpe introduced a bias that caused it to generate fewer false positives than expected. this bias appears to be due to a theoretical weakness and an incorrect assumption intrinsic to the test as originally proposed. murie and nadon have shown that the adjustment to the lpe standard error, which depends on an asymptotic proof, is overly conservative with sample sizes smaller than moreover, empirical evidence suggests that the lpe's assumption of similar error variability for genes of similar expression intensity is incorrect, which leads the lpe test to overestimate p-values for low variability genes and to underestimate them for high variability genes. recent modifications to the test have addressed this issue for microarray gene expression and for mass spectrometry proteomics data. indeed, potential differences in high and low variability genes between our studies and that of kooperberg et al. may explain why they concluded that the lpe test was overly liberal in contrast to our current findings and those of murie and nadon that the test is overly conservative. an important difference in our study between the cybert test and the other two empirical bayes tests, was that only the former maintained a stable false positive error rate with data that showed an unstable variance across intensity, such as the experimental and local normal simulated data. this suggests that these latter tests may also generate too many false positives among some biological data sets which, depending on pre-processing, often show high variance at low concentrations. as designed in our simulations and has been argued generally, however, this variance/intensity relationship is an artifact of the inappropriate log transformation for low intensity expression values and the generalized log transformation is able to stabilize the error variance across the entire concentration range. the glog transformation, as implemented with vsn and lmgene, lowers this potentially high fpr among low abundance transcripts in biological data sets while maintaining a relatively high tpr. background correction when applied to cdna data can generate negative values which cannot be logged and requires that the negative data points be either discarded or set to a loggable value. moreover, we found, as did qin and kerr, that using a log transformation after applying background correction in conjunction with loess normalization either reduced performance, as with the liver vs pool data, or had a mixed effect on performance, as with the liver vs liver data. we recommend the use of a glog transformation when applying background correction which can transform negative values and thus avoids the loss of information. the combination of glog and loess normalization produced comparable results for both cdna data sets with and without background correction and were amongst the top performers of all normalization methods. an attractive alternative might be to use background correction methods which avoid negative values and stabilize the variance as a function of intensity as described in ritchie et al.. the performance of the statistical tests was most similar using the hgu133a data set. with the exception of mas with fc, all of the statistical tests found most of the "spiked-in" genes readily, illustrating the limitations of these empirical data sets for comparing performances of statistical test as they would be applied in experiments with biological variation. the hgu133a data set, which used technical replicates, is comparable to our simulated data with low variance, where in general the methods do not show large performance differences. the hgu data had higher variability than the hgu133a data and also showed greater performance differences between methods. most biological experiments will be closer to the medium and the high variance simulated data sets where differences in performance are more apparent. the three empirical bayes methods in our study consistently performed well across normalization methods and simulated data designed to mirror the greater variability that is observed in most biological experiments. for this reason, we expect that our results apply beyond the narrow methodological confines of our study. nonetheless, more definitive tests of the algorithms' merits will be made possible by future spike-in experiments which are anticipated to incorporate biological variability. there is also a need for a similar comparative study among statistical tests which are suited to study designs that are more complex than the two-independent sample design in our study. two of the approaches we evaluated and at least one other that we did not examine have more general capabilities suited to this type of comparative study. authors' contributions cm conducted the computer simulations, performed the data analysis, participated in research design, and drafted the manuscript. ow and al participated in discussion of the research, investigation of the statistical methods, and analysis of data. rn conceived of the study, supervised research design, interpretation of results, and drafting of the manuscript. supplementary material additional file 1 mva plots of cdna data. click here for file additional file 2 source code availability. click here for file additional file 3 histograms of logged and unlogged variances for limma and brb statistical tests. click here for file additional file 4 pauc scores using a 5% false positive rate. click here for file additional file 5 trellis plots of performance criteria for simulated data. click here for file additional file 6 total average rankings of statistical tests for pauc and tpr. click here for file acknowledgements we thank mathieu miron, jacek majewski, and amy norris for their helpful suggestions.
three empirical bayes methods were the most effective statistical tests across simulated and both 2-colour cdna and affymetrix experimental data. all three empirical bayes methods assessed in our study are good choices for statistical tests for small n microarray studies for both affymetrix and cdna data. this can lead to poor random error estimates and inaccurate statistical tests of differential expression. we compare the performance of the standard t-test, fold change, and four small n statistical test methods designed to circumvent these problems. microarray studies often have small samples sizes, however, due to cost constraints or specimen availability. the generalized log transformation was shown to be beneficial in improving results with certain data sets, in particular high variance cdna data. the local pooled error test introduced a bias that lowered false positive rates below theoretically expected values and had lower power relative to the top performers. we report results of various normalization methods for empirical microarray data and of various random error models for simulated data. pre-processing of data influences performance and the proper combination of pre-processing and statistical testing is necessary for obtaining the best results. the cybert regularized t-statistic in particular was able to maintain expected false positive rates with simulated data showing high variances at low gene intensities, although at the cost of low true positive rates.
 dna microarrays provide data for genome wide patterns of expression between observation classes. microarray studies often have small samples sizes, however, due to cost constraints or specimen availability. this can lead to poor random error estimates and inaccurate statistical tests of differential expression. we compare the performance of the standard t-test, fold change, and four small n statistical test methods designed to circumvent these problems. we report results of various normalization methods for empirical microarray data and of various random error models for simulated data.  three empirical bayes methods were the most effective statistical tests across simulated and both 2-colour cdna and affymetrix experimental data. the cybert regularized t-statistic in particular was able to maintain expected false positive rates with simulated data showing high variances at low gene intensities, although at the cost of low true positive rates. the local pooled error test introduced a bias that lowered false positive rates below theoretically expected values and had lower power relative to the top performers. the standard two-sample t-test and fold change were also found to be sub-optimal for detecting differentially expressed genes. the generalized log transformation was shown to be beneficial in improving results with certain data sets, in particular high variance cdna data.  pre-processing of data influences performance and the proper combination of pre-processing and statistical testing is necessary for obtaining the best results. all three empirical bayes methods assessed in our study are good choices for statistical tests for small n microarray studies for both affymetrix and cdna data. choice of method for a particular study will depend on software and normalization preferences. 