microarrays are used to measure the ( relative ) expression levels of thousands of genes ( or expressed sequence tags ) . a comparison of gene expression in cells or tissues from two conditions may provide useful information on important biological processes or functions . it is now known that simply using fold changes , as in the earlier days , is unreliable and inefficient . many proposals have appeared in the literature . in particular , it has been noticed that it may be necessary to design an experiment that uses multiple arrays ( or multiple spots on each array ) containing multiple measurements for each gene under each condition . one reason is that because of a high noise - to - signal ratio , a single array may not provide enough information that can be reliably extracted . more important , multiple measurements from each gene make it possible to assess the potentially different variability of genes . the problem then seems to fall within the traditional two - sample comparison in statistics . two of the best known two - sample statistical tests are the two - sample t - test and the wilcoxon test ( or equivalently , mann - whitney test ) . the t - test is parametric and is based on the assumption that the gene - expression levels have normal distributions . in contrast , the wilcoxon test is nonparametric and is based on the ranks of observed gene - expression levels . although the t - test is robust to departures from normality and the wilcoxon test does not depend on the normality assumption , the problem is that under non - normal situations the t - test may be too conservative , and hence , as with the wilcoxon test , may have too low power , especially when the sample size is small , which is the case for most microarray experiments . these points have been verified in two case studies using real data . in a class of nonparametric approaches , a version of the two - sample t - statistic is used but its null distribution is estimated nonparametrically , rather than directly assumed to be a t - distribution . in addition , some earlier studies have suggested that the variability of gene expression may be related to the mean expression . therefore , it implies that the t - statistic being used should be based on unequal variances for the two samples . an important and natural question often asked by biologists is how many replicates are required . for microarray experiments , unlike many other experimental contexts , this issue has rarely been discussed in the literature . to our knowledge , the only exception is the work by black and doerge , which , however , is for the situation where parametric statistical methods are applied to detect expression changes . in this paper , we discuss the problem when a nonparametric method , the normal mixture model approach , is used to detect differential expression . but to facilitate calculations of sample size , the formulation is slightly changed from their original one . nonparametric methods of microarray data analysis have been pioneered by efron and tibshirani and co - workers . they take advantage of the presence of replicates and thus can impose much weaker modeling assumptions . for instance , the parametric methods of black and doerge   depend on the assumption on the log - normal or gamma distribution of gene - expression levels , whereas the mixture model approach does not have such a distributional assumption and directly estimates distributions related to random errors . note that modeling the distribution of random errors has advantages over direct modeling of expression levels , and is a common practice in applied statistics . for example , gene - expression levels may be correlated ( for example , as a result of coexpression of some genes ) whereas random errors can be more reasonably assumed to be independent . this is similar to modeling longitudinal data using a linear mixed - effects model : the responses from each subject ( corresponding to a group of coregulated genes here ) are in general correlated , but the measurement errors from the same subject can be considered to be independent after incorporating a random - subject effect in the model . note that the random effect will be canceled out from the t - statistic for each gene . our proposal here also shows an attractive feature of the mixture model approach , as compared to the other two nonparametric approaches , because it is still unclear how the sample size / power calculation can be done in the other two approaches . the problem of calculating the number of replicates required in a microarray experiment is similar to that of sample size / power calculations in clinical trials and other experiment designs ; the ( to - be - determined ) sample size in microarray experiments refers to the number of replicates , whereas the number of genes is not an issue here . as usual , we assume that the replicates are ( approximately ) independent with each other , whether they are drawn from the same individual or multiple individuals . in general , the required sample size depends on several factors : the true magnitude of the change of gene expression ( say , d ) , the desired statistical power ( that is , probability ) ( ) to detect the change , and the specified type i error rate ( ) . the problem of how to calculate the number of replicates for any given triplet ( d , , ) is equivalent to that of how the power   depends on the pair ( d , ) and the number of replicates , which we consider in the paper . , the expression level can refer to a summary measure of relative red - to - green - channel intensities in a fluorescence - labeled cdna array , a radioactive intensity of a radiolabeled cdna array ( as used in the example later ) , or a summary difference of the perfect match ( pm ) and mismatch ( mm ) scores from an oligonucleotide array . the gene - expression levels may have been suitably preprocessed , including dimension reduction , data normalization and data transformation . we consider a generic situation that , for each gene i , i = 1,2 , ... , n , we have ( relative ) expression levels x1i , ... , xmi from m microarrays under condition 1 , and y1i , ... , ymi from m arrays under condition 2 . a general statistical model is assumed for gene expression data : xji = ( 1),i + ji ,  yli = ( 2),i + eli , where ( 1),i and ( 2),i are the mean expression levels for gene i under the two conditions respectively , and ji and eli are independent random errors with means and variances for any j = 1 , ... , m and i = 1 , ... , n. it is assumed that random errors ji/(1),i and   are randomly taken respectively from one of two ( not necessarily equal ) distributions that are symmetric about their mean 0 . note that the above assumption on the distributions of random errors , not on that of gene expression levels ( that is , xji and yli ) , is often reasonable , and similar assumptions are common in other statistical applications . in addition , we do not assume that the expression levels of all the genes have an equal variance , because some previous studies   have found that the variance ( c),i ( for c = 1,2 ) of gene - expression levels may depend on the mean expression ( c),i . also , we do not even need to assume that ( 1),i = ( 2),i unless ( 1),i = ( 2),i . a goal is to detect all genes with ( 1),i  ( 2),i . this can be accomplished through statistical hypothesis testing . to test the null hypothesis h0 : ( 1),i = ( 2),i , we use a t - type test statistic or score note that the mean and variance of zi are whereas the mean e(zi ) = 0 under h0 . hence , it can be seen that a large absolute value of zi , |zi| , gives evidence against h0 . as the number of arrays ( that is , m ) increases , the variance of the test statistic zi decreases . hence , it is possible to reject h0 ( that is , detect differential expression for gene i ) with any e(zi )  0 if m is large enough . in other words , if the type i error rate and other parameters are fixed , then the statistical power of the test will increase as m increases . this is the key point that motivates the discussion on sample size calculations . to determine the cut - off point for |zi| to reject h0 , we need to know or estimate the distribution of zi under h0 , the null distribution f0 . in a parametric approach , based on some full distributional assumptions for xji and yji , one may derive the null distribution f0 , such as in a two - sample t - test . however , the validity of such a parametric method critically depends on the correctness of assumed distributions , which of course is not guaranteed . here , we consider a nonparametric approach : a finite normal mixture model is used to estimate f0 nonparametrically . for instance , using expression levels of some housekeeping genes that are known to have non - differential expression , one can construct their zi scores and then estimate f0 using the obtained zi scores . in practice , however , there may be only a small number of or no housekeeping genes in a given experiment . here , following the basic idea in a class of nonparametric methods , we construct a null score zi for each gene and then use these null scores to estimate f0 nonparametrically . the null score is constructed from the same observed gene expression data as used in zi : under the assumption that ji and eji have symmetric distributions , then ji and -ji have the same distribution , and eji and -eji have the same distribution . thus , by comparing the form of zi with that of zi , we know that the distribution of zi is exactly f0 , the null distribution for zi ( under h0 ) . note that under h0 , ( 1),i = ( 2),i , and hence ( 1),i = ( 2),i ( since we assume that ( c),i only depends on ( c),i ) , then thus zi and zi have the same distribution f0 under h0 . we use all zi values across all genes to estimate f0 . in practice , ( c),i ( for c = 1 , 2 ) are unknown , and can be estimated using the sample standard deviations ( sds ) s(c),i . although the sample sd s(c),i is asymptotically unbiased , if m and n are small , s(c),i may not be stable , and some modifications may be necessary . in any case , substituting ( c),i by any suitable estimates , we can calculate the scores zi values and zi values , on the basis of which we can estimate f0 and f respectively . by comparing f0 and f , we can gain insight about genes with altered expression ( that is , ( 1),i  ( 2),i ) . we assume that all the zi values for i = 1 , ... , n are a random sample from f0 ; thus we can use the observed zi values to estimate f0 . specifically , it is assumed that where   ( z ; ar , vr ) denotes the density function of a normal distribution n(ar , vr ) with mean ar and variance vr , and r values are mixing proportions . g0 represents all unknown parameters { r , ar , vr ) : r = 1, ... g0 } in a g0-component mixture model . among others , a normal mixture is essentially nonparametric and flexible , and easy to use with stable tail probabilities . a mixture model can be fitted by maximum likelihood using the expectation - maximization ( em ) algorithm . the number of components can be selected adaptively using the akaike information criterion ( aic )   or the bayesian information criterion ( bic ) . in using the aic or bic , one first fits a series of models with various values of g0 , then picks up the g0 corresponding to the first local minimum of aic or bic . once we obtain an estimate of the null distribution f0 , we can determine the cut - off point of the rejection region for testing h0 . in general , as for a two - sample test , the rejection region can be selected in the tails of f0 because , under the null hypothesis , zi should be close to the center of f0 , whereas if there is differential expression for gene i , zi is likely to be in one of the two tails of f0 . for example , if we are only interested in detecting upregulated genes , we can choose the rejection region at the right - tail of f0 . our proposed method works for any specified way of determining the rejection region . as f0 should be symmetric about its mean 0 , and often we are interested in both up- and downregulated genes , we propose to take the rejection region at the two tails of f0 , { z : f0(z ) < c } , where the constant c > 0 is the cut - off point and depends on the specified ( gene - specific ) type i error rate . as usual , c > 0 is chosen such that the rejection rate under h0 is exactly : where   ( . ; a , v ) is the corresponding cumulative distribution function for   ( . ; a , v ) . using a numerical algorithm , such as the bisection method , we can solve the above equation to obtain c for any given . for microarray data , because we are testing h0 for each gene , the multiple test problem arises and some control on it is necessary . for instance , if we want to maintain the genome - wide type i error rate at the usual 5% level , then the bonferroni - adjusted gene - specific ( that is , test - specific ) type i error rate is   = 0.05/n , where n is the total number of genes to be tested . once c is determined , we can calculate the power as a function of d , the magnitude of the expression change targeted to be detected . if ( 1),i = ( 2),i , d can be interpreted as the change of the mean expression levels from condition 1 to condition 2 . otherwise , it can be regarded as the difference of ( variation ) standardized mean expression levels . specifically , we have the power function unsurprisingly , we can see that ( d , ) will increase as |d| increases . the effects of having more replicates will reduce the variability of f0 , leading to larger ( d , ) for any given d. now we describe how to calculate replicate numbers based on some pilot data taken from earlier studies . we use zm , i to explicitly denote the zi scores in ( 2 ) with m replicates . based on the data we can estimate the density function f0,m ( z;g0 ) of zm , ivalues as a normal mixture from now on , we treat f0,m as known in equation ( 5 ) . with estimated f0,m , we want to estimate the density function f0,mk for zmk , i , the zi scores based on mk replicates ( with k > 1 ) . if we can have an estimate of f0,mk , then we can obtain the corresponding power function ( d , ) for mk replicates in the same way as described earlier for m replicates . of course , we assume that our pilot data are drawn from only m arrays under each of the two experimental conditions , and thus we do not observe any zmk , i based on mk arrays . however , we show next that it is possible to generate zmk , i values from zm , i values . note that we can draw random realizations of zm , i from the estimated f0,m ( see pan et al .   or the example below ) . suppose zm , i values ( for j = 1,2 , ... , k ) are k independent realizations of zm , i , then it is easy to show that have the distribution f0,mk . thus , the density function for zmk , i values is for example , if we triple the number of replicates , the resulting density function is the number of components of f0,mk may be too large . for example , if the number of components is g0 = 3 for m = n = 2 , the corresponding numbers of components for m = n = 4 , m = n = 6 and m = n = 8 are , respectively , g0 = 9 , g0 = 27 and g0 = 81 . in fact , some of these components maybe very similar or have a negligible role , hence the form of f0,mk , may be simplified . in the extreme situation , as mk   , by the central limit theorem , the mixture model will reduce to a single - component normal distribution . hence , we propose a simulation - based method to select a more parsimonious model for f0,mk . on the basis of the mixture model f0,m in equation ( 5 ) , we can generate a random sample of zm , i values , from which we can calculate zmk , i values using equation ( 6 ) . using zmk , i values we can fit a normal mixture model for f0,mk . as we shall show later , we find such a fitted mixture model often contains a smaller number of components than g0 , as dictated in equation ( 7 ) , leading to a simplified form of f0,mk . in summary , suppose that we have pilot gene expression data xji and yij from m arrays under each condition . use zi , m and the normal mixture model ( 5 ) to estimate f0,m . determine the cutoff point c for the rejection region using formula ( 3 ) , in which f0 is replaced with the estimated f0,m . calculate the power function ( d , ) using formula ( 4 ) , in which f0 is replaced with the estimated f0,m . step 5 . for any given k > 1 , use formula ( 7 ) or ( 6 ) to estimate f0,mk . determine the cutoff point c for the rejection region using formula ( 3 ) , in which f0 is replaced with the estimated f0,mk . calculate the power function   ( d , ) using formulae ( 4 ) , in which f0 is replaced with the estimated f0,mk . repeat steps 5 to 7 until all k > 1 of interest have been tried . after the power functions for many possible mk replicates have been obtained , we can determine an appropriate number of replicates by considering all the factors involved , the desired power and type i error rate , the targeted expression changes and other experimental constraints . to understand the pathogenesis of otitis media , a study was conducted to identify genes involved in response to pneumococcal middle - ear infection and to study their roles in otitis media . radioactively labeled dna microarrays were applied to the mrna analysis of 1,176 genes in middle - ear mucosa of rats with and without subacute pneumococcal middle - ear infection . the data are available for the control group and for the pneumococcal middle - ear infection group . a more detailed description of how the data were collected and their public availability was provided in pan et al . for the purpose of sample size calculations and to mimic many practical situations with only a small number of replicates , we only use m = n = 2 arrays from each group . we first take a natural logarithm transformation for all the observed gene - expression levels ( that is , radioactive intensities ) so that the resulting distributions are less skewed ( which will reduce the number of components of a fitted mixture model ) . then , for each microarray , we standardize the transformed gene - expression levels by subtracting their median . because of the small m = 2 one way is to add a small constant as suggested by efron et al . . here we follow the idea of lin et al .   and use a loess smoother   to nonparametrically model the sample sds in terms of the mean expression levels ( figure 1 ) . note that an alternative use of sd or its modification in calculating z2,i values will not change the basic idea and the following steps in sample size calculations . however , for the purposes of general illustration , we choose g0 = 2 as the fitted model : f0,2(z ) = 0.76 ( z;-0.0415 , 1.3117 ) + 0.24 ( z;0.0700 , 2.6970 ) . figure 2a presents the histogram of zi values and the fitted f0 with g0 = 1 and 2 . there is not much difference between the two fitted f0,2 , both of which fit the data well . in particular , f0,2 does not look like a t - distribution with small degrees of freedom , as predicted from the t - test . first , we draw a random number pi from { 1 , 2 } with probability 0.76 and 0.24 respectively . second , if the drawn pi = 1 , zi is randomly drawn from a normal distribution ( z ; -0.0415 , 1.3117 ) ; otherwise , it is drawn from ( z ; 0.0700 , 2.6970 ) . from the generated z2,i values , following expression ( 6 ) we generated three simulated data sets : z2k , i values , i = 1 , ... , 1,176 for k = 2 , 3 and 4 . then a normal mixture model was fitted to each data set . from table 1 , it can be seen that a single - component normal distribution was selected in each case . in figure 2 , each of the fitted normal distributions , n(-0.0494 , 0.8226 ) , n(-0.0644 , 0.5383 ) and n(-0.0438 , 0.4206 ) , is compared with its theoretically derived mixture model in equation ( 7 ) ; they are all very close . here we see that using simulated data to fit a mixture model results in a much - simplified model . for example , for k = 4 , it is a fitted single - component model versus a 2 = 16-component model in equation ( 7 ) . note that , as predicted , all the means of the fitted models are all essentially 0 , and their variances decrease as k increases . if we want to have only one expected false - positive result from testing each of 1,176 non - differentially expressed genes , the gene - specific ( or test - specific ) type i error rate is   = 1/1176 = 0.09% . using formula ( 3 ) and fitted - mixture model f0,2k , the cut - off points c are determined . then the power functions   ( d , ) are drawn in figure 3 , which may help make a decision on the required number of replicates . for instance , if we want to detect an expression change d = 3 with probability at least 80% and with   = 0.09% , then six replicates are needed . also , with just two replicates , the power to detect a change as high as 4 is very low , smaller than 30% . for instance , based on the pilot data , we can estimate the d values for some selected genes ( with the sample means and sample sds substituting the true means and sds in the formula for d ) , from which one can determine a range of d values of interest . figures 4,5,6 give the results for testing n = 1,000 , 5,000 and 10,000 genes , respectively , while controlling the genome - wide type i error rate at the usual 5% level . it can be seen that as n increases , we also need a larger number of arrays to maintain the power of the statistical test when other parameters are fixed . for instance , for n = 10,000 ( figure 6 ) , even eight replicates can not detect a change as large as d = 3 with 80% power , but six replicates can detect a change d = 4 with 80% power . we consider a generic situation that , for each gene i , i = 1,2 , ... , n , we have ( relative ) expression levels x1i , ... , xmi from m microarrays under condition 1 , and y1i , ... , ymi from m arrays under condition 2 . a general statistical model is assumed for gene expression data : xji = ( 1),i + ji ,  yli = ( 2),i + eli , where ( 1),i and ( 2),i are the mean expression levels for gene i under the two conditions respectively , and ji and eli are independent random errors with means and variances for any j = 1 , ... , m and i = 1 , ... , n. it is assumed that random errors ji/(1),i and   are randomly taken respectively from one of two ( not necessarily equal ) distributions that are symmetric about their mean 0 . note that the above assumption on the distributions of random errors , not on that of gene expression levels ( that is , xji and yli ) , is often reasonable , and similar assumptions are common in other statistical applications . in addition , we do not assume that the expression levels of all the genes have an equal variance , because some previous studies   have found that the variance ( c),i ( for c = 1,2 ) of gene - expression levels may depend on the mean expression ( c),i . also , we do not even need to assume that ( 1),i = ( 2),i unless ( 1),i = ( 2),i . a goal is to detect all genes with ( 1),i  ( 2),i . this can be accomplished through statistical hypothesis testing . to test the null hypothesis h0 : ( 1),i = ( 2),i , we use a t - type test statistic or score note that the mean and variance of zi are whereas the mean e(zi ) = 0 under h0 . hence , it can be seen that a large absolute value of zi , |zi| , gives evidence against h0 . as the number of arrays ( that is , m ) increases , the variance of the test statistic zi decreases . hence , it is possible to reject h0 ( that is , detect differential expression for gene i ) with any e(zi )  0 if m is large enough . in other words , if the type i error rate and other parameters are fixed , then the statistical power of the test will increase as m increases . this is the key point that motivates the discussion on sample size calculations . to determine the cut - off point for |zi| to reject h0 , we need to know or estimate the distribution of zi under h0 , the null distribution f0 . in a parametric approach , based on some full distributional assumptions for xji and yji , one may derive the null distribution f0 , such as in a two - sample t - test . however , the validity of such a parametric method critically depends on the correctness of assumed distributions , which of course is not guaranteed . here , we consider a nonparametric approach : a finite normal mixture model is used to estimate f0 nonparametrically . for instance , using expression levels of some housekeeping genes that are known to have non - differential expression , one can construct their zi scores and then estimate f0 using the obtained zi scores . in practice , however , there may be only a small number of or no housekeeping genes in a given experiment . here , following the basic idea in a class of nonparametric methods , we construct a null score zi for each gene and then use these null scores to estimate f0 nonparametrically . the null score is constructed from the same observed gene expression data as used in zi : under the assumption that ji and eji have symmetric distributions , then ji and -ji have the same distribution , and eji and -eji have the same distribution . thus , by comparing the form of zi with that of zi , we know that the distribution of zi is exactly f0 , the null distribution for zi ( under h0 ) . note that under h0 , ( 1),i = ( 2),i , and hence ( 1),i = ( 2),i ( since we assume that ( c),i only depends on ( c),i ) , then thus zi and zi have the same distribution f0 under h0 . we use all zi values across all genes to estimate f0 . in practice , ( c),i ( for c = 1 , 2 ) are unknown , and can be estimated using the sample standard deviations ( sds ) s(c),i . although the sample sd s(c),i is asymptotically unbiased , if m and n are small , s(c),i may not be stable , and some modifications may be necessary . in any case , substituting ( c),i by any suitable estimates , we can calculate the scores zi values and zi values , on the basis of which we can estimate f0 and f respectively . by comparing f0 and f , we can gain insight about genes with altered expression ( that is , ( 1),i  ( 2),i ) . we assume that all the zi values for i = 1 , ... , n are a random sample from f0 ; thus we can use the observed zi values to estimate f0 . specifically , it is assumed that where   ( z ; ar , vr ) denotes the density function of a normal distribution n(ar , vr ) with mean ar and variance vr , and r values are mixing proportions . g0 represents all unknown parameters { r , ar , vr ) : r = 1, ... g0 } in a g0-component mixture model . among others , a normal mixture is essentially nonparametric and flexible , and easy to use with stable tail probabilities . a mixture model can be fitted by maximum likelihood using the expectation - maximization ( em ) algorithm . the number of components can be selected adaptively using the akaike information criterion ( aic )   or the bayesian information criterion ( bic ) . in using the aic or bic , one first fits a series of models with various values of g0 , then picks up the g0 corresponding to the first local minimum of aic or bic . once we obtain an estimate of the null distribution f0 , we can determine the cut - off point of the rejection region for testing h0 . in general , as for a two - sample test , the rejection region can be selected in the tails of f0 because , under the null hypothesis , zi should be close to the center of f0 , whereas if there is differential expression for gene i , zi is likely to be in one of the two tails of f0 . for example , if we are only interested in detecting upregulated genes , we can choose the rejection region at the right - tail of f0 . as f0 should be symmetric about its mean 0 , and often we are interested in both up- and downregulated genes , we propose to take the rejection region at the two tails of f0 , { z : f0(z ) < c } , where the constant c > 0 is the cut - off point and depends on the specified ( gene - specific ) type i error rate . as usual , c > 0 is chosen such that the rejection rate under h0 is exactly : where   ( . ; a , v ) is the corresponding cumulative distribution function for   ( . ; a , v ) . using a numerical algorithm , such as the bisection method for microarray data , because we are testing h0 for each gene , the multiple test problem arises and some control on it is necessary . for instance , if we want to maintain the genome - wide type i error rate at the usual 5% level , then the bonferroni - adjusted gene - specific ( that is , test - specific ) type i error rate is   = 0.05/n , where n is the total number of genes to be tested . once c is determined , we can calculate the power as a function of d , the magnitude of the expression change targeted to be detected . if ( 1),i = ( 2),i , d can be interpreted as the change of the mean expression levels from condition 1 to condition 2 . otherwise , it can be regarded as the difference of ( variation ) standardized mean expression levels . specifically , we have the power function unsurprisingly , we can see that ( d , ) will increase as |d| increases . the effects of having more replicates will reduce the variability of f0 , leading to larger ( d , ) for any given d. now we describe how to calculate replicate numbers based on some pilot data taken from earlier studies . we use zm , i to explicitly denote the zi scores in ( 2 ) with m replicates . based on the data we can estimate the density function f0,m ( z;g0 ) of zm , ivalues as a normal mixture from now on , we treat f0,m as known in equation ( 5 ) . with estimated f0,m , we want to estimate the density function f0,mk for zmk , i , the zi scores based on mk replicates ( with k > 1 ) . if we can have an estimate of f0,mk , then we can obtain the corresponding power function ( d , ) for mk replicates in the same way as described earlier for m replicates . of course , we assume that our pilot data are drawn from only m arrays under each of the two experimental conditions , and thus we do not observe any zmk , i based on mk arrays . however , we show next that it is possible to generate zmk , i values from zm , i values . note that we can draw random realizations of zm , i from the estimated f0,m ( see pan et al .   or the example below ) . suppose zm , i values ( for j = 1,2 , ... , k ) are k independent realizations of zm , i , then it is easy to show that have the distribution f0,mk . thus , the density function for zmk , i values is for example , if we triple the number of replicates , the resulting density function is the number of components of f0,mk may be too large . for example , if the number of components is g0 = 3 for m = n = 2 , the corresponding numbers of components for m = n = 4 , m = n = 6 and m = n = 8 are , respectively , g0 = 9 , g0 = 27 and g0 = 81 . in fact , some of these components maybe very similar or have a negligible role , hence the form of f0,mk , may be simplified . in the extreme situation , as mk   , by the central limit theorem , the mixture model will reduce to a single - component normal distribution . hence , we propose a simulation - based method to select a more parsimonious model for f0,mk . on the basis of the mixture model f0,m in equation ( 5 ) , we can generate a random sample of zm , i values , from which we can calculate zmk , i values using equation ( 6 ) . using zmk , i values we can fit a normal mixture model for f0,mk . as we shall show later , we find such a fitted mixture model often contains a smaller number of components than g0 , as dictated in equation ( 7 ) , leading to a simplified form of f0,mk . in summary , our proposed method of calculating the required replicate number works in the following steps . suppose that we have pilot gene expression data xji and yij from m arrays under each condition . use zi , m and the normal mixture model ( 5 ) to estimate f0,m . determine the cutoff point c for the rejection region using formula ( 3 ) , in which f0 is replaced with the estimated f0,m . calculate the power function ( d , ) using formula ( 4 ) , in which f0 is replaced with the estimated f0,m . any given k > 1 , use formula ( 7 ) or ( 6 ) to estimate f0,mk . determine the cutoff point c for the rejection region using formula ( 3 ) , in which f0 is replaced with the estimated f0,mk . calculate the power function   ( d , ) using formulae ( 4 ) , in which f0 is replaced with the estimated f0,mk . repeat steps 5 to 7 until all k > 1 of interest have been tried . after the power functions for many possible mk replicates have been obtained , we can determine an appropriate number of replicates by considering all the factors involved , the desired power and type i error rate , the targeted expression changes and other experimental constraints . to understand the pathogenesis of otitis media , a study was conducted to identify genes involved in response to pneumococcal middle - ear infection and to study their roles in otitis media . radioactively labeled dna microarrays were applied to the mrna analysis of 1,176 genes in middle - ear mucosa of rats with and without subacute pneumococcal middle - ear infection . the data are available for the control group and for the pneumococcal middle - ear infection group . a more detailed description of how the data were collected and their public availability was provided in pan et al . for the purpose of sample size calculations and to mimic many practical situations with only a small number of replicates , we only use m = n = 2 arrays from each group . we first take a natural logarithm transformation for all the observed gene - expression levels ( that is , radioactive intensities ) so that the resulting distributions are less skewed ( which will reduce the number of components of a fitted mixture model ) . then , for each microarray , we standardize the transformed gene - expression levels by subtracting their median . because of the small m = 2 one way is to add a small constant as suggested by efron et al . . here we follow the idea of lin et al .   and use a loess smoother   to nonparametrically model the sample sds in terms of the mean expression levels ( figure 1 ) . note that an alternative use of sd or its modification in calculating z2,i values will not change the basic idea and the following steps in sample size calculations . however , for the purposes of general illustration , we choose g0 = 2 as the fitted model : f0,2(z ) = 0.76 ( z;-0.0415 , 1.3117 ) + 0.24 ( z;0.0700 , 2.6970 ) . figure 2a presents the histogram of zi values and the fitted f0 with g0 = 1 and 2 . there is not much difference between the two fitted f0,2 , both of which fit the data well . in particular , f0,2 does not look like a t - distribution with small degrees of freedom , as predicted from the t - test . first , we draw a random number pi from { 1 , 2 } with probability 0.76 and 0.24 respectively . second , if the drawn pi = 1 , zi is randomly drawn from a normal distribution ( z ; -0.0415 , 1.3117 ) ; otherwise , it is drawn from ( z ; 0.0700 , 2.6970 ) . from the generated z2,i values , following expression ( 6 ) we generated three simulated data sets : z2k , i values , i = 1 , ... , 1,176 for k = 2 , 3 and 4 . then a normal mixture model was fitted to each data set . from table 1 , it can be seen that a single - component normal distribution was selected in each case . in figure 2 , each of the fitted normal distributions , n(-0.0494 , 0.8226 ) , n(-0.0644 , 0.5383 ) and n(-0.0438 , 0.4206 ) , is compared with its theoretically derived mixture model in equation ( 7 ) ; they are all very close . here we see that using simulated data to fit a mixture model results in a much - simplified model . for example , for k = 4 , it is a fitted single - component model versus a 2 = 16-component model in equation ( 7 ) . note that , as predicted , all the means of the fitted models are all essentially 0 , and their variances decrease as k increases . if we want to have only one expected false - positive result from testing each of 1,176 non - differentially expressed genes , the gene - specific ( or test - specific ) type i error rate is   = 1/1176 = 0.09% . using formula ( 3 ) and fitted - mixture model f0,2k , the cut - off points c are determined . then the power functions   ( d , ) are drawn in figure 3 , which may help make a decision on the required number of replicates . for instance , if we want to detect an expression change d = 3 with probability at least 80% and with   = 0.09% , then six replicates are needed . also , with just two replicates , the power to detect a change as high as 4 is very low , smaller than 30% . for instance , based on the pilot data , we can estimate the d values for some selected genes ( with the sample means and sample sds substituting the true means and sds in the formula for d ) , from which one can determine a range of d values of interest . figures 4,5,6 give the results for testing n = 1,000 , 5,000 and 10,000 genes , respectively , while controlling the genome - wide type i error rate at the usual 5% level . it can be seen that as n increases , we also need a larger number of arrays to maintain the power of the statistical test when other parameters are fixed . for instance , for n = 10,000 ( figure 6 ) , even eight replicates can not detect a change as large as d = 3 with 80% power , but six replicates can detect a change d = 4 with 80% power . this method is designed for the situation where the mixture approach is going to be taken to analyze the data . note that any method for sample size / power calculations has to depend on a specific statistical test to be used in data analysis ; this explains why there is a huge literature on the topic for clinical trials . however , because of the close relation between the mixture approach and the other two recently proposed nonparametric approaches - the empirical bayes method   and the statistical analysis of microarray ( sam ) method   - our proposed method can be also applied to provide some useful guideline for designing microarray experiments even when one of the latter two approaches ( or other approaches ) is planned to be used for data analysis in a later stage . for instance , even though the null distribution f0 is estimated using the null scores zi in our proposal , there maybe alternative ways of estimating f0 , such as using an alternative nonparametric method ( for example , kernel or local likelihood ) , rather than the finite normal mixture model , to estimate f0 or using the test statistics , zi , of a large number of housekeeping genes to estimate f0 . some modifications to the test statistic zi and the null statistic zi are also possible , especially when we consider differential gene expression across more than two conditions . these are all interesting topics we are investigating now . in most sample size / power calculations , some pilot data are needed to provide reasonable estimates of some parameters needed for subsequent calculations . however , because of the rapid development of microarray technology , the latter is not likely and we expect a researcher will have to do his or her own pilot study . a particular challenge is how to obtain good estimates of the variances of gene expression levels from a small number of replicates . in our example the proposed method is straightforward to statisticians and can be implemented in many existing statistical packages . sample standard deviations of expression levels and their loess smoothers as a function of the average expression levels for the two conditions respectively . histograms and estimated distribution density functions . ( a - d ) two , four , six and eight replicates ( z2 - z8 ) , respectively . in ( a ) , the solid and dotted lines are the fitted one- and two - component mixtures . in ( b - d ) , the solid and dotted lines are the fitted and the theoretically derived mixtures . power   ( d , ) as a function of the magnitude of expression changes d and the number of replicates , with the gene - specific type i error rate   = 0.09% for the middle - ear data . power ( d , ) as a function of the magnitude of expression changes d and the number of replicates , with the gene - specific type i error rate   = 0.05/1,000 for the middle - ear data . power   ( d , ) as a function of the magnitude of expression changes d and the number of replicates , with the gene - specific type i error rate   = 0.05/5,000 for the middle - ear data . power   ( d , ) as a function of the magnitude of expression changes d and the number of replicates , with the gene - specific type i error rate   = 0.05/10,000 for the middle ear data .
backgroundit has been recognized that replicates of arrays ( or spots ) may be necessary for reliably detecting differentially expressed genes in microarray experiments . 
 however , the often - asked question of how many replicates are required has barely been addressed in the literature . in general , the answer depends on several factors : a given magnitude of expression change , a desired statistical power ( that is , probability ) to detect it , a specified type i error rate , and the statistical method being used to detect the change . here , we discuss how to calculate the number of replicates in the context of applying a nonparametric statistical method , the normal mixture model approach , to detect changes in gene expression.resultsthe methodology is applied to a data set containing expression levels of 1,176 genes in rats with and without pneumococcal middle - ear infection . 
 we illustrate how to calculate the power functions for 2 , 4 , 6 and 8 replicates.conclusionsthe proposed method is potentially useful in designing microarray experiments to discover differentially expressed genes . 
 the same idea can be applied to other statistical methods .