Biostatistics Advance Access originally published online on March 24, 2006
Biostatistics 2006 7(4):630-641; doi:10.1093/biostatistics/kxj032
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
A Laplace mixture model for identification of differential expression in microarray experiments
Ecole Polytechnique Fédérale de Lausanne, Institute of Mathematics, EPFL-FSB-IMA, Station 8, CH-1015 Lausanne, Switzerland anthony.davison{at}epfl.ch
* To whom correspondence should be addressed.
| SUMMARY |
|---|
|
|
|---|
Microarrays have become an important tool for studying the molecular basis of complex disease traits and fundamental biological processes. A common purpose of microarray experiments is the detection of genes that are differentially expressed under two conditions, such as treatment versus control or wild type versus knockout. We introduce a Laplace mixture model as a long-tailed alternative to the normal distribution when identifying differentially expressed genes in microarray experiments, and provide an extension to asymmetric over- or underexpression. This model permits greater flexibility than models in current use as it has the potential, at least with sufficient data, to accommodate both whole genome and restricted coverage arrays. We also propose likelihood approaches to hyperparameter estimation which are equally applicable in the Normal mixture case. The Laplace model appears to give some improvement in fit to data, though simulation studies show that our method performs similarly to several other statistical approaches to the problem of identification of differential expression.
Keywords: Laplace distribution; Marginal likelihood; Microarray experiment; Mixture model; REML
| 1. INTRODUCTION |
|---|
|
|
|---|
Microarrays have become an important tool for studying the molecular basis of complex disease traits and fundamental biological processes. Two-channel microarrays, such as spotted cDNA or long oligonucleotide arrays, measure relative gene expression in two samples. Once preprocessed, the data from such arrays take the form of normalized base 2 logarithm of the expression ratios. A common purpose of microarray experiments is the detection of genes that are differentially expressed under two conditions, such as treatment versus control or wild type versus knockout. Numerous statistical methods have been proposed for identification of differential expression, with new ones continuing to be introduced.
In early analyses such as Schena and others (1995
, 1996
), fold change between conditions exceeding a constant was used to identify differentially expressed genes. This method performs poorly, however, because it ignores the different variability of expression across genes. Such variability can be taken into account by using a t-test on the average log-fold change, but variation in gene expression is poorly estimated with small numbers of replicates, so genes with artificially low variance may be selected even if they are not truly differentially expressed. Numerous refinements have been proposed in order to reduce the numbers of false-positive and false-negative results. Tusher and others (2001)
and Efron and others (2000)
suggested adding a constant to the t denominator so that it does not become too small. Lönnstedt and Speed (2002)
proposed a normal mixture model for the gene expression data and defined a log posterior odds statistic, their B-statistic, for ranking genes. Their approach was extended to linear models for more general designs by Smyth (2004)
; these methods are implemented in the R package "limma" as part of the BioConductor project (www.bioconductor.org). Gottardo and others (2003)
not only use a similar approach but also include a heuristic, iterative method for estimating the proportion of differentially expressed genes.
The present paper makes three main contributions: we consider robust and asymmetric variants of the mixture modeling approach, we propose two new approaches to parameter estimation, and we perform a comparative study intended to elucidate key features of such methods. The methods studied include those listed above, along with our own proposal of a Laplace mixture as a long-tailed alternative to the normal distribution for mean gene expression across replicate arrays proposed by Lönnstedt and Speed (2002)
.
This paper is organized as follows: Section 2 describes the model, with estimation of the hyperparameters discussed in Section 3. Simulation studies outlined in Section 4 show that this method performs similarly to several other statistical approaches to the problem of identification of differential expression. In Section 5, we apply our method to a published data set and compare its results with other methods. The paper ends with a brief discussion.
| 2. LAPLACE MIXTURE MODEL |
|---|
|
|
|---|
Lönnstedt and Speed (2002)
propose the use of a mixture of a point mass at zero and a normal distribution for mean gene expression, and of an inverse gamma distribution for the single-gene variances. Gaussian variation is the exception rather than the rule in practice, however, so we consider using a Laplace distribution as a potentially more realistic longer tailed alternative to the normal law. We assume that the data are normalized base 2 logarithms of fold changes, as would be the case for two-channel arrays. It is not difficult to apply our approach to single-channel technologies, by modeling the single-channel log intensities, and using a contrast matrix to obtain corresponding log ratios.
We suppose initially that it is reasonable to expect approximately symmetric over- and underexpression, as when an array encompasses an entire genome. Arrays with selected genome coverage are considered in Section 2.2.
We model normalized relative expression measures on each of G genes independently, but in this section, we ease the notation by confining our attention to data from a single gene. Although in reality genes do interact with each other, the independence assumption is a useful simplification. The measured relative expression levels for this gene in n replicates,
, are taken to be independent and normally distributed with mean relative expression level
and variance
, i.e. 
We denote the average and sample variance of the observations y by
and
.
We let
represent the probability that this gene is differentially expressed, and assume that in this case the mean expression level
has a Laplace distribution with mean zero and variance 2
, so that the density is given by
![]() |
If the gene is not differentially expressed, its mean expression level is zero. Thus,
may be regarded as a random draw from the mixture distribution
, where
represents the Laplace distribution for differential expression and
denotes the distribution which places unit mass at
. Conditional on
, we assume that
, where
represents a type of generalized signal to noise ratio. Thus, if the variability
of expression levels across replicates is large, the variation of
will also be large. This assumption, which simplifies subsequent mathematical developments, is checkable from the data; in the cases we have examined it seems fairly reasonable. Variability of measured expression levels is not constant across genes, and we follow Lönnstedt and Speed (2002)
in assuming an inverse gamma distribution for
, that is
, where
and
are shape and scale parameters, respectively. We use a similar notation for the gamma distribution
.
We treat
, V,
, and
as hyperparameters, and in Section 3.1 discuss their estimation.
A straightforward but involved calculation shows that the marginal posterior distribution of
, like the prior for
, is a mixture of continuous and discrete components. It may be expressed as
![]() | (2.1) |
where
is a Dirac delta function and provided
,
![]() | (2.2) |
Here
represents the gamma function,
,
![]() |
and
denotes the cumulative distribution function of a Student t variable with
degrees of freedom. Thus h, the posterior density of
given that a gene is differentially expressed, comprises two off-centered Student t densities placed back-to-back at the origin.
The posterior odds of differential expression may be written as
![]() | (2.3) |
The second term on the right of (2.3) is the Bayes factor for differential expression
relative to
.
Most methods for the detection of differential expression assume equal over- and underexpression, but this may be unrealistic when arrays contain genes chosen for their special interest to the investigator, or in other cases where genome coverage is restricted. We therefore consider an asymmetric extension of our mixture model.
The density of an asymmetric Laplace variable, which may be obtained as the difference of two independent exponential variables with different means, may be expressed as (Kotz and others, 2001
, Chapter 3)
![]() |
where
. Taking
yields the symmetric Laplace distribution, while
gives increasing weight to the right tail of the density, and conversely when
. We write the corresponding distribution as
, in terms of which the asymmetric Laplace mixture model for
may be written as
. A calculation generalizing that in the symmetric case shows that when
the posterior density of
is again given by (2.1) and (2.2), but with
![]() |
The posterior odds of differential expression are again given by (2.3), and the hyperparameters
, V,
,
, and
may be estimated by an empirical Bayes procedure, as in the symmetric case.
The discussion above relates to data from a single gene, but in applications many genes will be considered, making it useful to reduce the entire marginal posterior density of
to a single value. One way to do this is to consider the posterior median
of
. If the posterior probability that
exceeds one-half, then
, whereas if the posterior probability that
exceeds one-half, then
. If neither of these events occurs, then the point mass at
ensures that
. Thus, the use of a posterior median to estimate
yields a thresholding rule according to which
, if the posterior distribution of
is too closely concentrated around zero. This method has been used for selection of non-zero components in Bayesian implementations of wavelet smoothing (Johnstone and Silverman, 2005
).
| 3. ESTIMATION PROCEDURES |
|---|
|
|
|---|
We now suppose that measurements
are available on gene g, for
, and that the genes behave independently. The posterior density of the expression level for gene g,
, depends on the data y and on hyperparameters
in the symmetric case and
in the asymmetric case. Consider the symmetric model under which
|
|
for
. A simple and in principle an efficient approach to estimating the hyperparameters is to maximize the log marginal likelihood
![]() | (3.1) |
based on the marginal densities
![]() |
Rather than adopt such a procedure, Lönnstedt and Speed (2002)
suggested fixing the proportion of differentially expressed genes
at an arbitrary value and then estimating other parameters using the method of moments, while Gottardo and others (2003)
use a heuristic approach to reduce instability in the estimates. Below we describe a two-stage likelihood procedure which can yield fairly stable estimates.
The basis of the two-stage approach is the elementary fact that under a model with Gaussian errors, the average and sample variance
and
are jointly sufficient statistics. That is,
![]() |
where the terms on the right are normal and chi-square density functions. As the marginal density of
depends only on
, hyperparameters
of the density for
can be estimated by maximizing the marginal likelihood obtained as the product of
![]() | (3.2) |
The resulting estimates do not take into account any information about
contained in
, but have the useful property that they do not depend on parameters
of the prior density for
. Estimates of
obtained from the marginal likelihood based on (3.2) can be substituted into (3.1), which is then maximized with respect to the remaining hyperparameters. The result is a two-stage procedure, which is compared to an overall optimization of (3.1) in Section 3.2.
An important special case is where the density of
is inverse gamma, that is,
is gamma distributed with shape and scale parameters
and
, respectively. Then the marginal distribution of
is
, and estimates
,
are readily obtained by maximum likelihood (ML) estimation based on this distribution. Moreover, the suitability of the inverse gamma distribution can be indirectly assessed by plotting the ordered
against quantiles of the fitted F distribution. We have found that with this model, standard routines can be applied without difficulty in both optimization stages.
The above approach to estimation of the hyperparameters associated with the variances
applies also to designed experiments, provided they are replicated, because it is independent of any superstructure for the
. It also extends to cases where the error structure is equivariant, that is, where we may write
, with the
having a known density w. Let
represent the ordered
. The differences of order statistics
depend only on
and w, and so
may be estimated using a likelihood based on the marginal densities
![]() |
Typically, these integrals must be obtained numerically or approximated, unlike when w is the normal density and the prior density of
is inverse gamma.
In this section, we compare our two-stage and certain other approaches to hyperparameter estimation using simulated data. The simulation study is based on three scenarios. Initially, we generated 100 data sets each with
genes and
replicates from the Laplace mixture model with
and
. An example simulated dataset is represented in Figure 1.
|
For Scenario I, we simulate data sets with parameter values chosen to match the simulations of Lönnstedt and Speed (2002)
are generated from the inverse gamma distribution with parameter values
and
estimated from the scavenger receptor BI transgenic mice data set (Callow and others, 2000
and
. For Scenario II, we generate data sets using parameters estimated from a data set on 150 defence-related genes on the plant Arabidopsis thaliana (Reymond and others, 2000
and
. Scenario III uses parameter values intermediate between the other two, taking the mean and variance of the error variances to be 5 and 50, respectively; this yields
and
.
For each simulated dataset, we estimated the parameters by our two-stage marginal likelihood procedure (MML) and by ordinary maximum likelihood (ML). Table 1 gives the bias and root mean square error for the three scenarios, based on 100 datasets. Both procedures converged readily and appear to perform well, but as one would expect on theoretical grounds, ordinary maximum likelihood performs at least as well as marginal likelihood estimation in most cases. One perhaps surprising aspect is the high biases observed for estimation of
, using either procedure.
|
In a second study, we generated data from the asymmetric Laplace mixture model. For each of various scenarios, we generated 100 data sets each with 2000 genes. As this model is likely to be considered for subsets of genes selected because they are thought to be important, the probability of differential expression
was taken to be appreciably higher than the first study (
or 0.75). The asymmetry parameter
was varied (
), with the remaining parameters taken to be similar to the values for the symmetric cases described above. We also considered the effect of taking
replicates as well as setting
.
Unfortunately, with the low numbers of replicates considered here there is difficulty in estimating most parameters under full ML, especially for values of
substantially different from 0. Even with REML estimation, typically, only the parameters
and
are well estimated unless
is close to 0. For some of the parameter sets, we generated 100 replicates, and in these cases estimation appeared to be improved. Apparently, large numbers of replicates are required for acceptable parameter estimates in this more complex mixture model.
| 4. COMPARATIVE STUDY |
|---|
|
|
|---|
A number of methods have been proposed for identifying differentially expressed genes, with new ones continuing to be introduced. We compare the performance of several of these after brief descriptions of those we consider.
In order to generalize and compare our work with that of Lönnstedt and Speed (2002)
, we include the statistics they used, both versions of our statistic, and another recently introduced method. These seven statistics and their justification are
L -statLaplace mixture model
AL -statAsymmetric Laplace mixture model
B -statNormal mixture model (Lönnstedt and Speed, 2002
)
-statNormal mixture (Gottardo and others, 2003
)
M -statAverage
fold change
t -statt-statistic based on individual gene standard deviation
Pen-t -statPenalized t-statistic (Efron and others, 2000
).
Note that when, as is the case here, all genes are equally replicated, the B-statistic reduces to the moderated t-statistic (Smyth, 2004
).
Simulated data sets were generated under different mixture models. For each scenario, we generated 100 data sets with 3000 genes each and
, 4, or 8 replicates, as in Section3.2. Parameter values for the simulations were chosen based on our analysis of a publicly available Arabidopsis thaliana data set, presented in Section 5.
The parameter corresponding to the proportion of differentially expressed genes (
or p) is estimated from the data for the L- and AL-statistics, and set to 0.1 for the rest of the statistics.
Parameters for L-stat and AL-stat are estimated by the two-step procedure described in Section 3. The B-stat is computed using the "eBayes" function of the R package "limma" (Smyth, 2004
). This package is available as part of the open source BioConductor project (Gentleman and others, 2004
). We computed the
-stat with default values of the function "B1.bayes" in the package "amd" that works with older versions of R; similar functionality is now available through the BioConductor packages "rama" and "bridge."
We generated data from the Laplace mixture model, as described in Section3.2, with parameter values estimated by applying this model to the Arabidopsis data set.
Table 2 gives the generating parameter values, with summaries of estimates from the L and AL models for sample sizes
, 4, and 8.
|
We also simulated data sets under the normal mixture model as for the Laplace mixture model. Because the findings are very close to those for data generated under the Laplace model, we do not show detailed results. The B-statistic and L-statistic perform quite similarly, although with data generated from the normal mixture model the B-statistic performs slightly better.
For each gene, the statistics are calculated for each of the 100 data sets. For a range of cutoff values for all the methods, the numbers of false-positive and false-negative results are calculated in each data set. For each cutoff value of a given statistic, the observed numbers of false-positive and false-negative genes are then averaged over the 100 data sets. For each method, a sufficient range of cutoffs was chosen to be able to observe the behavior of the receiver operating characteristic (ROC) curves over a large range. The average numbers of false-positive genes are plotted against the average numbers of false-negative genes for the statistics in variants of ROC plots. Curves falling more steeply toward the lower left corner correspond to more powerful test procedures.
A few of these plots are presented in Figure 2. The AL-statistic performed practically identical to the L-statistic, so only the curve corresponding to the L-statistic is shown, along with those of the other six statistics. We summarize our general findings from examining several such curves.
|
Several of the methods perform broadly similarly across simulation conditions: L-stat, AL-stat, B-stat,
-stat, and Pen-t-stat. In these examples, mean
fold change M also shows reasonably good performance as the great majority of gene-specific variances are small. Unsurprisingly, the t-statistic has poor performance for these sample sizes and cannot be recommended.
To reduce computing time in these simulations, we generated arrays consisting of a few thousand genes rather than the tens of thousands of genes more typically encountered in practice. We also examined our scenarios using 20 000 genes, varying
from 0.01 to 0.10, in order to assess the sensitivity method performance in this more realistic situation. Within each generating scenario, relative method performance was broadly similar across values of
.
| 5. EXAMPLE: Arabidopsis thaliana DATA SET |
|---|
|
|
|---|
We present an example using data publicly available from the Stanford Microarray Database (SMD) (http://genome-www5.stanford.edu/). The experiments were carried out to compare the general effect on disease resistance RNA transcript levels of Arabidopsis thaliana infected by rhizobacterium Pseudomonas thivervalensis (strain MLG45) to axenic control plants (Cartieaux and others, 2003
). Here we consider the experiment on plant leaves, which has slides containing 16 416 spots for four biological replicates hybridized to a common reference (SMD Experiment ID numbers 27084, 27000, 26995, and 26718). Further details are available from Cartieaux and others (2003)
.
From the raw channel intensities, we computed print-tip loess normalized log ratios for each spot (Yang and others, 2002
). To avoid inflated variance for low-intensity spots, we did not carry out any background correction; we also removed control spots. From these normalized log ratios, we estimated the parameters for the Laplace and Asymmetric Laplace models (Table 3). These values were used as a starting point for part of the simulation study of Section 4.
|
Here we examine the fit of the Laplace mixture model to the data by comparing gene-specific means and variances of the true and simulated data.
The values of the log ratios seen in the simulated data sets tend to be larger than in the original data, although quantile-quantile (QQ) plots comparing genewise data means to simulated means indicate a reasonably good fit (data not shown).
Under the proposed mixture models, the marginal distribution of the sample variances follows an
distribution, where n is the number of replicates. We can thus check the fit of the model to the data set variances by plotting the sample variances against the quantiles of the
distribution. Figure 3(a) shows that the agreement seems very good overall, with some deviation for the highest few variances.
|
To further investigate the fit of the Laplace mixture model, and compare it with that of the normal mixture model (Lönnstedt and Speed, 2002
and the sample variance
.
Under the Laplace model, we find that
![]() | (5.1) |
which is quadratic in
, with coefficients that can be determined from n and the ML estimates of
and V.
Under the normal mixture model (Lönnstedt and Speed, 2002
), the mean
for differentially expressed genes is distributed as normal with mean zero and variance
, where
. Thus, the sample mean square,
, is linear in the sample variance (
):
![]() | (5.2) |
Here the coefficient can be determined by assuming
, the default value for the B-statistic, and with V estimated by the method of moments (Lönnstedt and Speed, 2002
).
We further investigate the relationship between the sample variances and the sample means by fitting generalized linear models. The models account for the squared means in terms of linear (normal) or linear and quadratic (Laplace) functions of the variances, using the gamma family with identity link. The quadratic effect is significant at level around 0.001, suggesting a strong preference for the Laplace model over the normal one.
| 6. DISCUSSION |
|---|
|
|
|---|
This study makes three main contributions to statistical methodology for analysis of microarray experiments. First, we have introduced the Laplace mixture model. This model permits greater flexibility than models in current use as it has the potential, at least with sufficient data, to accommodate both whole genome and restricted coverage arrays. In addition, it also appears to provide a somewhat better fit to data. We have also proposed fundamentally sound approaches for estimating the proportion of differentially expressed genes that also appear to perform well in practice. Last, we have further characterized and compared several methods of identifying differential expression.
Our original motivation for proposing the Laplace mixture model was that the normal mixture model of Lönnstedt and Speed (2002)
did not fit our data well. We were initially surprised to find that this lack of fit only marginally affects performance: the Laplace model fits better but the performance of both is similar, so gene rankings from these models are anticipated to be alike across a variety of data sets. Thus, although there may be some room for "fine tuning," we do not argue that the Laplace mixture should supplant the commonly used normal mixture in practice. Rather, we believe that efforts aimed at further model refinement might more profitably be focused on solving new and more challenging problems.
Our proposed procedures for estimating the proportion of differentially expressed genes seem to perform quite well in the simulations. Unlike the heuristic procedure of Gottardo and others (2003)
, our likelihood approach does not require an arbitrary criterion for calling a gene differentially expressed. This approach could also be readily incorporated into existing normal mixture model-based software. However, our approach might in some cases give poor results, in which case the obvious work-around is to simply fix the value of
, as suggested by Lönnstedt and Speed (2002)
.
In the comparison study, several proceduresthose which in essence use some type of global shrinkage to modify the denominatorappear to perform quite similarly across simulation conditions. Some forms of local shrinkage give worse performance than global shrinkage (Kooperberg and others, 2005
); we did not consider them here. These results underline the futility of endeavors to create substantially better procedures by further tinkering. The procedure that will be optimal in a given experiment depends on many factors which in practice are typically unknown, such as the form of the data generating distribution and the true parameter values, including the proportion of differentially expressed genes. The key aspect seems to be how the smallest and largest single-gene variances are treated.
Finally, we would like to draw attention to the issue of software. Carrying out this type of study and, more importantly, data analyses in general, depends heavily on the availability of readily usable software. Here, we mostly relied upon the open source R-based BioConductor project (Gentleman and others, 2004
), which provides an ideal framework for software distribution as well as a large number of quality software tools for statistical genomics computing. We are in the process of preparing a more user-friendly version of our software as an R package for BioConductor. We urge those who create software intended for wide dissemination also to consider implementation in R and contribution to BioConductor.
| ACKNOWLEDGMENTS |
|---|
This work was funded by the Swiss National Science Foundation through the National Centre for Competence in Research in Plant Survival. Conflict of Interest: None declared.
| REFERENCES |
|---|
|
|
|---|
-
Callow MJ, Dudoit S, Gong EL, Speed TS, Rubin EM. (2000) Microarray expression profiling identifies genes with altered expression in HDL-deficient mice. Genome Research 10:20229.
Cartieaux F, Thibaud M-C, Zimmerli L, Lessard P, Sarrobert C, David P, Gherbaud A, Robaglia C, Somerville S, Nussaume L. (2003) Transcriptome analysis of Arabidopsis colonized by a plant-growth promoting rhizobacterium reveals a general effect on disease resistance. The Plant Journal 36:17790.[CrossRef][Web of Science][Medline]
Efron B, Tibshirani RJ, Goss V, Chu G. (2000) Microarrays and their use in a comparative experiment. Technical report. (Department of Health Research and Policy, Stanford University, Stanford, CA.).
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J. (2004) Bioconductor: open software development for computational biology and bioinformatics. Genome Biology 5:R80 http://genomebiology.com/2004/5/10/R80.[CrossRef][Medline]
Gottardo R, Pannucci JA, Kuske CR, Brettin T. (2003) Statistical analysis of microarray data: a Bayesian approach. Biostatistics 4:597620.[Abstract]
Johnstone IM and Silverman BW. (2005) Empirical Bayes selection of wavelet thresholds. Annals of Statistics 33:170052.[CrossRef]
Kooperberg C, Aragaki A, Strand AD, Olson JM. (2005) Significance testing for small microarray experiments. Statistics in Medicine 24:228198.[CrossRef][Web of Science][Medline]
Kotz S, Kozubowski TJ, Podgórski K. (2001) The Laplace Distribution and Generalizations(Birkhäuser., Boston).
Lönnstedt I and Speed TP. (2002) Replicated microarray data. Statistica Sinica 12:3146.
Reymond P, Weber H, Damond H, Farmer EE. (2000) Differential gene expression in response to mechanical wounding and insect feeding in Arabidopsis. Plant Cell 12:70719.
Schena M, Shalon D, Davis RW, Brown PO. (1995) Quantitative monitoring of gene expression patterns with a complementary DNA microarray. Science 270:46770.
Schena M, Shalon D, Heller R, Chai A, Brown PO, Davis RW. (1996) Parallel human genome analysis: microarray-based expression monitoring of 1000 genes. Proceedings of the National Academy of Sciences of the United States of America 93:106149.
Smyth GK. (2004) Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 3:.
Tusher VG, Tibshirani RJ, Chu G. (2001) Significance analysis of microarrays applied to the ionising radiation response. Proceedings of the National Academy of Sciences of the United States of America 98:511624.
Yang YH, Dudoit S, Luu P, Speed TP. (2002) Normalisation issues for cDNA microarray data. In Bittner ML, Chen Y, Dorsel AN, Dougherty ER (Eds.). Microarrays: Optical Technologies and Informatics. Proceedings of SPIE(SPIE, Bellingham, WA) pp. 14152.
Received December 15, 2005; revised February 18, 2006; revised March 10, 2006; accepted for publication March 23, 2006.
![]()
CiteULike
Connotea
Del.icio.us What's this?
This article has been cited by other articles:
![]() |
A. C. Davison and D. Mastropietro Saddlepoint approximation for mixture models Biometrika, June 1, 2009; 96(2): 479 - 486. [Abstract] [PDF] |
||||
![]() |
K. Lo and R. Gottardo Flexible empirical Bayes models for differential gene expression Bioinformatics, February 1, 2007; 23(3): 328 - 335. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||













= 0.1 and the number of replicates is n = 2.
,
,
,
; (b)
,
,
.



