1.

Introduction

DNA microarray technology has emerged as a useful technology in genetic research. Using this technology, it is possible to quantitate nucleic acids through fluorescent intensities of signals from distinct tissue samples. The ratios of these signals can be calculated and inferences drawn as to which differences in signal intensity are meaningful.

Chen, Dougherty, and Bittner¹ studied the problem of analyzing the mRNA from two distinct tissue samples. They derived an elegant formulation for the probability density function of a signal ratio under the condition that the means of the two component signals are equal. The density function was derived with the assumption that signal expression levels across the gene population of a microarray are independent random variables. Further assumptions are that the signal expression levels follow a normal distribution and that all signal expression levels share a common coefficient of variation.

Informal discussions with biologists have suggested that the assumption of a constant coefficient of variation for the entire gene set may not always be appropriate. For example, in a recent study we utilized the microarray approach of Genome Systems, Inc., to measure the expression of RNA from female mouse liver tissues as it related to the level of serum corticosterone. Our data consisted of five pairs of samples that compared the livers of mice with high and low serum corticosterone, thus permitting the calculation of a coefficient of variation (CV) for each of the 8772 genes in the study. We found that the ranges for the CVs were substantial. For 28 of the genes, ratios from the five pairs of samples produced CVs in the range of 0.05–0.9. About 53 of the gene ratios were in the range from 1.0 to 4.9 and 19 were in the range from 5.0 to 9.9. Our experience is that wide ranges such as these are not unusual in studies involving large numbers of genes in complex biological systems. In this particular study, the 20 genes with the lowest expression intensities had an average CV of 6.5 (±2.6), a phenomenon commonly observed with low expression intensities. The 20 genes with the highest expression intensities, on the other hand, had a lower average CV of 1.0 (±1.9). Of the many genes in a large study, some may be tightly regulated while others may be inducible or repressible by factors that may not be related to experimental parameters.

Our concerns with the normality and constant coefficient of variation assumptions led us to conduct a robustness study in connection with the Chen–Dougherty–Bittner (CDB) procedure. The question of robustness is an important issue with respect to any applied statistical technique. A robust procedure is one that is affected only slightly by appreciable departures from the assumptions involved. In regard to the CDB procedure the question of robustness reduces to the following: suppose the expression levels are not normally distributed and the coefficients of variation across the entire gene set are not constant. Can we still use the CDB procedure for constructing approximately accurate confidence intervals? How likely are we to draw erroneous conclusions about observed differences in expression ratios?

If the CDB procedure is robust, it may be used regardless of whether or not theoretical assumptions of normality and constancy in the coefficient of variation are satisfied. This paper reports results of Monte Carlo studies that assessed the impact of varying sample sizes, distribution shapes, and degrees of heterogeneity in coefficients of variation. We shall demonstrate that the Chen–Dougherty–Bittner procedure is unaffected by varying sample sizes and is, indeed, quite robust to even extreme deviations from normality and to substantial inequalities in the coefficients of variation.

2.

Chen–Dougherty–Bittner Procedure

Chen et al.;¹ developed an expression for the probability density function for a signal ratio. Assuming we have n genes represented in a microarray then two signal expressions for the kth gene are denoted as X_k and Y_k, respectively. For example, these two signals might come from treatment and control conditions, where X_k represents a signal for the kth gene obtained from a sample of mRNA taken from a treated cell and Y_k represents a signal for the kth gene taken from a cell in the control condition. Let the ratio of the two signals for the kth gene be denoted by T_k. Therefore

Assuming that all genes in the array have a common coefficient of variation, c, that the distributions of X_k and Y_k are N(μ_k,c²μ_k ²), and that all signal measurements are independent random variables, the approximate probability density function for T_k is given by

The probability density function (PDF) in Eq. (3) depends only on the common coefficient of variation c, and not on the parameters of the distributions of X_k and Y_k. Thus a test of the null hypothesis in Eq. (2) can be conducted at the α level of significance by integrating Eq. (3) and finding the 100^*(α/2) and the 100^*(1−α/2) percentile points of the distribution of T_k. Performing this integration, however, requires a value for c. This value is obtained by use of a maximum likelihood estimator, obtained by calculating the likelihood function across the entire gene array and then solving for c^ as the maximizing value. Thus c can be estimated from the data as follows:

With c estimated from the data using Eq. (4), the rejection region for testing, Eq. (2), can be established. Under the assumption of Eq. (2), approximately α^*100 of the t_k values would fall outside the region of rejection, i.e. (1−α)^*100 of the t_k values would be contained in the interval between the critical values for such a region of rejection. This assumes the distribution of the ratio behaves in accordance with Eq. (3).

3.

General Method and Research Questions

The dependent variable of interest throughout this paper is the empirical containment rate (ECR) of the CDB procedure. The ECR is the proportion of times the CDB procedure correctly accepts the null hypothesis of equal signal means under specific sets of experimental conditions that are described below. The ECR is the complement of the empirical type I error rate (EER). The ECR and the EER must sum to 1.0. For example, suppose under a specific experimental setup, the proportion of times the procedure accepts the null hypothesis (the ECR) is 0.937. Then the proportion of times the procedure rejects the null hypothesis (the EER) is 0.063. Thus, in this example, ECR+EER=0.937+0.063=1.0. Since the ECR and the EER always sum exactly to 1.0, whatever factors explain the ECR also explain the EER.

The three independent variables that we consider are (1) sample size, the number of simulated pairs of signals per gene array, (2) distribution shape, the proximity of the simulated data to normality, and (3) pattern of coefficients of variation within a given gene array.

Three general questions are proposed and evaluated in this paper. (1) Is the CDB procedure robust with respect to (i) violation of the normality assumption and (ii) with respect to violation of the assumption of a common coefficient of variation across all genes in the array? (2) Can an adequate statistical model be developed to predict ECR from the independent variables in the study? (3) What accounts for the robustness or lack of robustness of this procedure?

The answer to the second question will determine to a great extent the answer to the first question. If a highly accurate predictive model can be developed for ECR, then the robustness or nonrobustness of the CDB procedure is inferable from such a model. In fact, the answer to the second question can even provide a more sensitive answer to the robustness question by informing researchers about what levels of deviation from basic assumptions do not cause substantial deviations from nominal (expected) rates of containment.

The answer to the third question involves an analysis of the dynamics of the CDB procedure to determine how it produces a decision to reject or not reject the null hypothesis. The procedure computes a ratio and then calculates a pair of critical values to be used as criteria for rejection. The latter depends on the proper estimation of the common coefficient of variation, c. Thus, non-normality and heterogeneous coefficients of variation can affect the procedure in two different ways: (1) by a possible distortion in the distribution of the ratio, or (2) by a possible effect on the estimation of c. The robustness, or lack thereof, of the CDB procedure will depend on the interplay of these two effects.

4.

Study I: Method and Results

The first simulation study investigated the effects of the three independent variables on the empirical containment rate. The three independent variables are described in turn.

5.

Study II: Method and Results

To further investigate the relationship between non-normality, CV pattern, and ECR, a more extensive study was done. Given the lack of a significant relationship between sample size and ECR, the data for all cells was simulated with a sample size of 8000. The study also took into account the fact that an interaction effect was found between the CV pattern and the distribution shape in study I.

5.3.

Results

To demonstrate that the data were generated according to specifications, the average values of the mean, variance, skewness, and excess kurtosis are reported for both the X and Y variables, for each distribution shape, at the CV pattern that involves complete homogeneity of CV. These results are presented in Table 4. The four moments of the distribution are very close on average to the intended values, although the Fleishman procedure slightly underestimates the excess kurtosis.

Table 4

Average mean, variance, skew, and excess kurtosis by distributional shape, CV pattern reflecting all CVs equal to 0.15. (Averages are based on 1000 replications.)
Statistic	Normal	Mild NN	Moderate NN	Extreme NN
Avg mean X	0.999 97	1.000 07	1.000 09	1.000 03
Avg variance X	0.022 48	0.024 90	0.022 51	0.022 49
Avg skew X	0.000 60	0.496 73	0.992 99	1.733 32
Avg xkur Xa	−0.007 67	0.238 94	0.968 38	3.633 37
Avg mean Y	1.000 01	1.000 09	0.999 95	0.999 93
Avg variance Y	0.022 50	0.022 51	0.022 48	0.022 48
Avg skew Y	−0.000 18	0.495 17	0.992 66	1.737 66
Avg xkur Ya	−0.007 14	0.232 05	0.965 50	3.639 52
a xkur=excess kurtosis.

The average estimated coefficients of variation by the CV pattern are displayed in Figure 1 for each distribution shape. The average estimated CV increases as the heterogeneity of the CV pattern increases. This relationship holds for all distribution shapes, although in the three non-normality conditions, the common CV of 0.15 is underestimated under the experimental condition of complete homogeneity of the CVs. The degree of underestimation is in direct proportion to the non-normality of the distribution.

Figure 1

Mean estimated CV by distribution shape.

The average, minimum, and maximum values of ECR for each CV pattern by distribution shape are given in Table 5. A quick summary of the findings shows that of the 120 000 replications in the entire study, the lowest ECR observed was 0.916 (minimum for CV pattern 30 with Extreme NN) and the highest ECR observed was 0.959 (maximum for CV pattern 1 with Mild NN). For the condition combining the most serious deviation from normality (Extreme NN) and the most extreme heterogeneity of coefficients of variation (CV pattern 30), the average ECR was 0.922. The ECRs are represented graphically in Figure 2. Also represented in Figure 2 are the ECR15s. Note that for the CV pattern reflecting complete homogeneity [standard deviation (SD) of CV=0 ] and with a normal distribution of the signals, both the ECR and the ECR15 are precisely equal to 0.95. The two sets of containment rates exhibit the same decreasing pattern for all distribution shapes.

Figure 2

Average ECR and ECR15 by SD of CV.

Table 5

Average, minimum, and maximum ECR by CV pattern by distribution shape (based on 1000 replications).
CV Pattern	Normal			Mild NN			Moderate NN			Extreme NN
	Avg	Min	Max	Avg	Min	Max	Avg	Min	Max	Avg	Min	Max
1	0.950	0.945	0.955	0.953	0.948	0.959	0.951	0.946	0.956	0.937	0.932	0.942
2	0.950	0.946	0.955	0.953	0.946	0.958	0.951	0.945	0.956	0.937	0.932	0.943
3	0.950	0.944	0.955	0.953	0.947	0.959	0.951	0.945	0.956	0.937	0.932	0.944
4	0.950	0.945	0.955	0.953	0.947	0.958	0.950	0.945	0.955	0.937	0.932	0.943
5	0.949	0.944	0.954	0.952	0.947	0.957	0.950	0.944	0.954	0.937	0.931	0.943
6	0.949	0.943	0.954	0.951	0.945	0.957	0.949	0.944	0.954	0.937	0.931	0.943
7	0.948	0.943	0.953	0.951	0.946	0.956	0.949	0.944	0.954	0.936	0.931	0.941
8	0.948	0.941	0.953	0.950	0.944	0.956	0.948	0.942	0.953	0.936	0.931	0.942
9	0.947	0.940	0.952	0.949	0.944	0.954	0.947	0.940	0.952	0.936	0.930	0.941
10	0.946	0.940	0.952	0.948	0.942	0.953	0.946	0.941	0.951	0.936	0.930	0.942
11	0.945	0.941	0.951	0.946	0.941	0.951	0.945	0.939	0.952	0.935	0.930	0.940
12	0.944	0.938	0.950	0.945	0.939	0.950	0.944	0.937	0.949	0.935	0.930	0.941
13	0.943	0.937	0.949	0.944	0.938	0.949	0.943	0.937	0.948	0.934	0.928	0.940
14	0.942	0.937	0.947	0.943	0.936	0.948	0.941	0.935	0.947	0.933	0.928	0.938
15	0.941	0.936	0.946	0.941	0.936	0.946	0.940	0.935	0.945	0.933	0.927	0.939
16	0.940	0.934	0.945	0.940	0.935	0.944	0.939	0.933	0.944	0.932	0.927	0.938
17	0.939	0.932	0.944	0.938	0.932	0.943	0.939	0.932	0.943	0.931	0.926	0.937
18	0.938	0.933	0.944	0.937	0.931	0.941	0.936	0.931	0.943	0.931	0.925	0.937
19	0.937	0.931	0.942	0.936	0.931	0.943	0.935	0.930	0.941	0.930	0.923	0.934
20	0.936	0.930	0.941	0.935	0.930	0.940	0.934	0.929	0.939	0.929	0.924	0.934
21	0.935	0.930	0.940	0.933	0.928	0.938	0.933	0.927	0.938	0.928	0.923	0.933
22	0.934	0.929	0.940	0.932	0.927	0.938	0.931	0.926	0.937	0.927	0.922	0.933
23	0.933	0.928	0.939	0.931	0.924	0.937	0.930	0.924	0.936	0.927	0.921	0.932
24	0.932	0.928	0.938	0.930	0.924	0.936	0.929	0.924	0.934	0.926	0.919	0.931
25	0.931	0.926	0.937	0.929	0.924	0.934	0.928	0.922	0.934	0.925	0.920	0.930
26	0.931	0.926	0.936	0.928	0.923	0.934	0.927	0.923	0.932	0.924	0.919	0.930
27	0.930	0.925	0.936	0.928	0.923	0.932	0.926	0.922	0.931	0.924	0.916	0.929
28	0.929	0.924	0.936	0.927	0.921	0.932	0.926	0.921	0.931	0.922	0.917	0.929
29	0.929	0.923	0.934	0.926	0.920	0.931	0.925	0.919	0.931	0.922	0.917	0.927
30	0.928	0.923	0.934	0.925	0.920	0.930	0.924	0.919	0.930	0.922	0.916	0.927
Overall	0.940	0.923	0.955	0.940	0.920	0.959	0.939	0.919	0.956	0.931	0.916	0.944

TEs, DEs, and EEs were computed and are represented graphically in Figure 3 as a function of the heterogeneity of the CV pattern. The EE shows a decreasing pattern for all distribution shapes. This reflects the fact that the estimated coefficient of variation increases with the heterogeneity in the CV pattern. Likewise, the displacement in the graphs of EE by the distribution corresponds to the displacement in the graphs of the estimated coefficients of variation by distribution shape. The pattern in the DE shows an opposing pattern to that in the EE. DE increases for every distribution shape as the heterogeneity of the CV pattern increases. Furthermore, under the homogeneity of the CV condition, there was no observed distortion effect for the normal distribution shape but there were negative distortion effects for each of the non-normal shapes. As a general rule, the DE is more positive as the distribution shape gets closer to normality.

Figure 3

Average distortion effect, estimation effect, and total effect by SD of CV.

The moments of the signal ratio distributions are displayed in Table 6. The mean, variance, skewness, and kurtosis all increase as a function of CV heterogeneity. Both Pearson and Spearman correlations between these moments and distortion effects are also presented in Table 6. Distortion effects are strongly positively correlated with the four distribution indices. In the case of the ratios under the normality condition, the clear effects of extreme outliers are observed both in the case of the variance of the ratio and in the case of the Pearson correlation with DE. The nonparametric Spearman correlation was therefore included to provide an indicator of the relationship of the ratio variance and the distortion effect that is not sensitive to extreme values.

Table 6

Average mean, average variance, average skew, and average kurtosis of ratio by selected-CV patterns. (All statistics are based on 1000 replications, and correlations are based on a sample size of 30 000.)
Pattern	Normality				Mild non-normality
	Avg. mean	Avg. variance	Avg. skew	Avg. kurtosis	Avg. mean	Avg. variance	Avg. skew	Avg. kurtosis
1	1.0242	0.052	0.791	4.723	1.0223	0.047	0.612	3.547
4	1.0244	0.052	0.828	4.965	1.0225	0.047	0.629	3.643
7	1.0246	0.054	0.943	5.937	1.0228	0.048	0.689	3.978
10	1.0257	0.056	1.230	13.628	1.0233	0.050	0.778	4.487
13	1.0269	0.060	1.511	14.867	1.0241	0.052	0.906	5.235
16	1.0282	0.069	2.454	60.083	1.0252	0.055	1.054	6.115
19	1.0303	0.084	3.687	111.489	1.0263	0.059	1.225	7.199
22	1.0365	97.865	6.918	313.816	1.0282	0.063	1.411	8.462
25	1.0368	1.187	10.438	535.867	1.0300	0.068	1.612	9.898
28	1.0434	22.511	14.630	821.201	1.0318	0.075	1.820	11.562
30	1.0444	5.148	16.541	927.335	1.0338	0.080	1.993	13.094
ra	0.161c	0.013d	0.410c	0.295c	0.788c	0.977c	0.966c	0.919c
rhob	0.841c	0.969c	0.921c	0.899c	0.770c	0.989c	0.973c	0.968c
Pattern	Moderate non-normality				Extreme non-normality
	Avg. mean	Avg. variance	Avg. skew	Avg. kurtosis	Avg. mean	Avg. variance	Avg. skew	Avg. kurtosis
1	1.0207	0.043	0.609	3.596	1.0188	0.040	0.816	4.850
4	1.0208	0.044	0.626	3.683	1.0190	0.040	0.834	4.961
7	1.0209	0.044	0.672	3.938	1.0191	0.041	0.883	5.281
10	1.0213	0.046	0.746	4.348	1.0194	0.042	0.958	5.762
13	1.0221	0.047	0.843	4.892	1.0200	0.043	1.064	6.450
16	1.0228	0.049	0.957	5.529	1.0206	0.045	1.183	7.207
19	1.0238	0.052	1.077	6.195	1.0213	0.047	1.311	8.049
22	1.0252	0.055	1.217	6.993	1.0223	0.049	1.450	8.930
25	1.0264	0.059	1.343	7.715	1.0233	0.052	1.582	9.794
28	1.0282	0.063	1.476	8.502	1.0245	0.056	1.710	10.647
30	1.0291	0.067	1.563	9.030	1.0253	0.058	1.799	11.259
ra	0.723c	0.982c	0.971c	0.958c	0.651c	0.979c	0.948c	0.925c
rhob	0.705c	0.987c	0.971c	0.968c	0.631c	0.980c	0.948c	0.939c
a Pearson’s correlation between the distortion effect and the variable.
b Spearman’s correlation between the distortion effect and the variable.
cp<0.01.
dp<0.05.

The increases in the moments of the ratio distribution as a function of CV heterogeneity are moderated by the distribution shape of the constituent signals. If we take the moments of the ratio distribution for the normal condition and CV homogeneity as normative, it is clear that all of the non-normal shapes show a smaller dispersion in the signal ratio at low levels of CV heterogeneity than this normative level. This accounts for the negative distortion effects in the non-normal conditions at low levels of heterogeneity of CV. The signal ratio dispersion increases as the CV heterogeneity increases, and the differences in ratio dispersion among distribution shapes of the constituent signals are consistent with the differences in the graphs of DE by distribution shape of the signals.

6.

Discussion

This study was designed to answer three questions. First we asked, “Is the CDB procedure robust with respect to (1) deviations from normality in the shapes of the signal distributions and to (2) deviations from homogeneity of the coefficients of variation of the signal distributions?” The answers to these questions are “Yes.” The worst case of a deviation from 0.950 was 0.916 in 120 000 replications. Each of the distribution shapes had an average ECR between 0.922 and 0.928 under even the worst pattern of CV heterogeneity. For most cases of CV heterogeneity the resulting ECRs were even closer to the nominal containment rate of 0.950.

In practice this means that even if expression levels in microarrays are not normally distributed and even if the coefficients of variation across the gene set are not constant, the Chen–Dougherty–Bittner procedure may still be used for constructing highly accurate confidence intervals. Even under extreme deviations from the assumptions of Chen, Dougherty, and Bittner, an error in rejecting the null hypothesis that we would expect to occur only 5 of the time, would, in reality, occur no more often than 8 of the time.

Another question posed was, “Can an adequate explanatory model be developed to predict the ECR from the independent variables in the study?” The answer is, again, “Yes.” Sample size played no role in explaining empirical containment rates. The ECR decreased as the CV heterogeneity increased. Distribution shape moderated this relationship. In general, the ECR was closest to the nominal level of 0.95 for the normality condition and farthest away from 0.95 for the extreme non-normality condition. The results for the two intermediate non-normality conditions showed similarity to each other. Each had ECR values that were too high under CV homogeneity. However, ECR values dipped below those for the normality condition as the CV heterogeneity increased.

A four-parameter logistic model was used to describe the data. It showed an outstanding fit to the data for all cases of non-normality. The model permitted an estimate of what the ECR would be if the CV heterogeneity were to become infinitely large. The results were consistent with the previously expressed finding that the CDB procedure is robust against violations of its relevant assumptions.

A final question involved an assessment of why the CDB procedure is robust. The decision as to whether or not to reject the null hypothesis of equality of means has two parts. A ratio is computed and a confidence interval is obtained. If the ratio is in the interval, the null hypothesis is not rejected. Otherwise it is rejected. Thus the rate of containment of the procedure is a function of the distribution of the signal ratio and the size of the confidence interval. The size of the confidence interval is strictly a function of the estimate of the common coefficient of variation.

The TE was shown to be the sum of the EE and the DE. We found that, for each distribution shape, the two effects offset each other. For normal data, large positive distortion effects are greatly offset by large negative estimation effects, thus stabilizing the total effects and yielding ECRs that are very close to the nominal level. Other distribution shapes show a similar pattern of rising distortion effects accompanied by falling estimation effects. The non-normality conditions have negative distortion effects offset by positive estimation effects at low CV heterogeneity levels. However, the trends in the graphs are the same for each distribution shape. For the extreme non-normality condition, the actual level of the EE remains positive over all CV patterns. The two effects of CV heterogeneity have opposite directions. They serve to mitigate each other. Increased dispersion in the ratio distribution is checked by increases in the estimated common coefficient of variation, leading to concomitantly larger confidence intervals. The net effect is relative stability in the ECR.

In summary, we have demonstrated that the Chen–Dougherty–Bittner procedure is robust to deviations of assumptions regarding normality and constancy in the coefficients of variation in the signal distributions. We have developed a statistical model that explains how the empirical containment rate would be affected if heterogeneity in the coefficient of variation were to become infinitely large. Finally, we have shown that distortion effects and estimation effects offset each other in contributing to total effects on empirical containment rates, further explaining the reasons for the robustness of the procedure.