The impact of measurement errors in the identification of regulatory networks

Background There are several studies in the literature depicting measurement error in gene expression data and also, several others about regulatory network models. However, only a little fraction describes a combination of measurement error in mathematical regulatory networks and shows how to identify these networks under different rates of noise. Results This article investigates the effects of measurement error on the estimation of the parameters in regulatory networks. Simulation studies indicate that, in both time series (dependent) and non-time series (independent) data, the measurement error strongly affects the estimated parameters of the regulatory network models, biasing them as predicted by the theory. Moreover, when testing the parameters of the regulatory network models, p-values computed by ignoring the measurement error are not reliable, since the rate of false positives are not controlled under the null hypothesis. In order to overcome these problems, we present an improved version of the Ordinary Least Square estimator in independent (regression models) and dependent (autoregressive models) data when the variables are subject to noises. Moreover, measurement error estimation procedures for microarrays are also described. Simulation results also show that both corrected methods perform better than the standard ones (i.e., ignoring measurement error). The proposed methodologies are illustrated using microarray data from lung cancer patients and mouse liver time series data. Conclusions Measurement error dangerously affects the identification of regulatory network models, thus, they must be reduced or taken into account in order to avoid erroneous conclusions. This could be one of the reasons for high biological false positive rates identified in actual regulatory network models.


Background
There has been an increasing interest among bioinformaticians in the problem of quantifying correctly gene expression in a given sample. It is well accepted that the observed gene expression value is a combination of the "true" gene expression signal with intrinsic biological variation (natural fluctuation) and a variation caused by the measuring process, also known as measurement error. Studies have documented the presence of sizable measurement error in data collected mainly from Open Access microarrays and also by other approaches such as Real Time RT-PCR, Northern blot, CAGE, SAGE, etc [1,2]. This measurement error can be easily observed when two technical replicates are plotted in a MA (M is the logarithm of the intensity ratio and A is the mean of the logged intensities for a dot in the plot) or scatter plots. Frequently, a considerable dispersion can be observed. This dispersion is due to the measurement error, since, in theory, technical replicates (same samples) must present the same quantifications. In general, these fluctuations are derived from probe sequence, hybridization problems, high background fluorescence, signal quantification procedures (image analysis), etc [3,4]. In the last few years, a considerable number of reports on the problem of quantifying and separating "true" gene expression signal from noise [5][6][7] has been published with the main aim to find differentially expressed genes [8,9]. Despite these results in gene expression analysis and a large amount of research performed in modeling regulatory networks (Bayesian networks [10,11], Boolean networks [12,13], Relevance networks [14], Graphical Gaussian models [15], Differential equations [16], etc), only a fraction of the statistical studies use procedures designed for modeling networks taking into account measurement error.
Frequently, Ordinary Least Squares (OLS) and methods related to it, such as Pearson and Spearman correlations [17], ridge, lasso and elastic net regressions [18] are widely used as estimators to quantify the strength of association between gene expression signals and model regulatory network structures. In the time series context, estimation process of Autoregressive (AR) [19][20][21][22] models also use OLS to identify which gene is or is not Granger causing another gene. Generally, a regression is carried out between the target gene and its potential predictors in order to test which predictor gene has, at a gene expression level, association with the target gene.
It is well known in the statistical literature that, when the measurement errors are ignored in the estimation process, OLS and its variants become inconsistent (i.e., even increasing the sample size the estimates do not converge to the true values). More precisely, the estimation of the slope parameters is attenuated [23] and consequently, regulatory network models become biased. Moreover, there is no control of type I error since standard OLS was not designed to treat measurement error. In this context, an adequate inference treatment must be considered for the model parameters in order to avoid inconsistent estimators. Usually, measurement equations are added to the model to capture the measurement errors effect, therefore, producing consistent estimators, efficient and asymptotically normally distributed. A careful and deep exposition on the inferential process is presented in [23] and the references therein. Although there are studies referring to problems caused by measurement errors in the statistical literature, there is a gap in the network modeling theory which must be filled in to avoid misinterpretation and distort conclusions from the inferential process. Here, we focus on the development and present some important statistical tools to be applied in OLS-based and VAR network models taking into account the measurement errors effect. We also conduct simulation studies in order to evaluate the impact of the measurement error in the identification of gene regulatory networks using the standard OLS in both conditions, time series and non time series data. Surprisingly, both the simulations and theory described that, in the presence of measurement error, the estimated coefficients are biased even increasing the amount of observations, and the statistical tests are not controlling the rate of false positives properly. These results were also observed in time series context, where the autoregressive coefficients were strongly affected. Thus, a corrected version of the OLS estimator for independent (in the regression context) and dependent (in the autoregressive context) data containing measurement error were developed. Results in both, simulated and actual biological data are illustrated. Moreover, two procedures to estimate measurement error in microarrays are presented.

Results and discussions
In order to evaluate the performance of conventional OLS and VAR methods in practice, simulations were carried out in artificial data with absence and presence of measurement error. Noise was added at different rates, and sample size was increased in order to evaluate the consistence of conventional and proposed approaches.
In the following we give a brief explanation about the usual and proposed methods. Let x and y be variables (gene expression values) with the following relationship y = a + bx + ε, where ε is the random error (intrinsic biological variation) of the model with zero mean and finite variance. In general, we are interested in estimating the parameters a and b to make inferences about them. In practice, we take a sample x i , y i for i = 1,..., n and use these quantities to obtain estimates for the parameters of interest. However, it is not always possible to observe directly the values of x and y because sometimes they are latent values, i.e., they are masked by measurement errors derived by the measurement process in microarrays, for example. Then, instead of observing the true variables, we observe surrogate variables X and Y which carry an error, that is X = x + 1 and Y = y + 2 , where 1 and 2 are measurement errors. Generally, what is done in practice is a naive solution, since it simply replaces x with X and y with Y in the regression equation and uses the OLS approach to estimate the parameters. That is, based on the equation Y = a + bX + ε, estimators are built. On the other hand, the proposed approach is slightly different. The latter considers three equations, namely: y = a + bx + ε, X = x + 1 and Y = y + 2 and uses them to estimate the model parameters. This little difference can result great impact in the estimators properties of each approach. Notice that the former produces inconsistent estimators and the latter produces consistent estimators when the data contains measurement error. The same idea can be applied in the time series context.
The corrected versions of the OLS estimators in both independent and dependent data were compared to their conventional forms in order to evaluate the performance under gene expression data containing measurement error. The standard OLS and VAR models are particular cases of the proposed models in the case when the measurement error is absent. Firstly, simulations were performed in regression models. Table 1 illustrates average coefficients estimated by standard OLS in 10,000 Monte Carlo simulations. Notice that increasing the rate of measurement error, more attenuated become the estimated coefficients, i.e., the estimates are shifted towards zero. Table 2 illustrates the percentage of rejected hypotheses in 10,000 Monte Carlo simulations. Analyzing when b 1 = 0, i.e., when there is no association between the corresponding covariate and the response variable, Table 2 shows that the OLS approach does not control, at a 5% nominal level, the rate of false positives. The larger the sample size, the worst the OLS performance, as it was expected to be. On the other hand, the coefficients of the corrected OLS are unbiased (Table 1 values between brackets) and converge to "true" value when sample's size becomes larger. Moreover, the rate of false positives are actually controlled under the null hypothesis ( Table 2 -values between brackets).
Analyzing Figure 1A, we conclude that, the standard OLS is not controlling the rate of false positives in the presence of measurement error for any significance level (p-value threshold). On the other hand, Figure 1B describes the consistency of the test performed by the corrected OLS, i.e., the uniform distribution of p-values illustrates that the rate of false positives is actually controlled under any considered threshold, since the    Average OLS estimated coefficients and corrected OLS (between brackets) in 10,000 simulations. The model is described in (Simulations section, simulation I -independent data)."-" means that it was not possible to calculate due to high measurement error in comparison to sample's size. EM: Standard deviation of the Error of Measure. n: Number of samples. Notice that, in the presence of measurement error, the coefficients (b) estimated by the corrected OLS (between brackets) converge to the "true" values, while the estimated by standard OLS do not.
uniform distribution emerges for p-values when the distribution of the statistic is correctly specified (otherwise, the p-value distribution may not be uniform).
In the time series case, similar results were observed. The standard VAR estimates produce biased coefficients in the presence of measurement error (Table 3). Moreover, there is no control of the type I error in both, autoregressive and cross-autoregressive coefficients (in all the text, in order to simplify the notation, autoregressive coefficient will denote the auto-loop, i.e., the coefficient related to z i, t-r → z i, t and cross-autoregressive coefficient will represent the coefficient for z j, t-r → j and r < t) ( Table 4). Analyzing the results produced by the proposed VAR model (Table 3 -values between brackets), it is possible to observe that the estimated coefficients converge to the true value as time series length increases. Notice that, the results produced by the standard VAR model indicate that, increasing the sample size does not imply in convergence of the estimates to the true values (Table 3). By observing Table 4, we see that the corrected VAR approach is actually controlling the rate of false positives in the set significance level (p < 0.05). Figure 2 emphasizes this result. Figure 2A and 2B describe the p-value distributions of autoregressive and cross-autoregressive coefficients of standard VAR under the null hypothesis (b 0 = 0 (autoregressive) and b 1 = 0 (cross-autoregressive)). Notice that when b 0 = b 1 = 0, the p-value distributions should be uniform in the interval [0,1]. However, there is a high concentration around zero, demonstrating that the rates of false positives are inflated (and consequently not controlled) in both autoregressive or cross-autoregressive cases. In Figures 2C and 2D, the p-value distributions are uniform, i.e., the test under the null hypothesis (b 0 = 0 and b 1 = 0) using the corrected VAR model is actually controlling the type I error in autoregressive and cross-autoregressive coefficients (uniform distribution). Figures 3 and 4 illustrates the corrected power curves for both, OLS and VAR. The corrected power curve P c (a) can be defined as where a is the adopted type I error nominal level, P(a(a)) is the power using the true probability of the type I error, namely a(a). Notice that the corrected     power is just the power penalized by the distance between a(a) and a. This correction in the power is necessary because under the null hypothesis, the power has to be the nominal level, and for comparing powers from different statistics it must be done using the same nominal level.
For a good statistic, notice that under an alternative hypothesis and when n → ∞, the corrected power P c (a) converges to one because a n ( ) α α → →∞ and P a n ( ( )) α → →∞ 1 .
On the one hand, for a statistic that does not control the rate of false positives, for example, when a is set to 5% and the true probability of the type I error is a(a) = 0.08, since a(a)/a is greater than one, P c (a) will not converge to one. On the other hand, for a good statistic, the rate a (a)/a converges to one when n → ∞, then P c (a) will converge to one. Analyzing Figures 3 and 4, it is possible to verify that, for standard OLS and VAR approaches (dashed lines), the ratio a(a)/a increases faster than the corresponding powers P(a(a)), i.e., the dashed lines is decreasing as n increases. Notice on Tables 2 and 4 that the rates of false positives (a(a)) increase as n increases, and consequently, in our specific case, the ratio a(a)/a increases and the corrected power P c (a(a)) converges to zero. On the other hand, the proposed methods (full lines) keep the false positive rates controlled while the corrected power increases as n increases. It can be observed by the full lines converging to one (Figures 2 and 4) and also on Tables 2 and 4. The variations present in the curves are probably due to variations in Monte Carlo simulations, since these variations decreased (become smoother) when the number of simulations was increased from 5,000 to 10,000 and from 10,000 to 15,000. In order to illustrate the performance of standard and corrected OLS and VAR approaches in actual biological data, firstly, the measurement error was estimated using the method described in the Measurement error estimation section (No technical replicates subsection). Then, the TP53 network was constructed using a dataset composed of 400 microarrays. Table 5 illustrates the results of a multivariate regression using OLS. Four genes known to be direct targets of TP53 were selected, namely, MDM2, FAS, BAX and MAP4, and a multivariate network was constructed using OLS. In fact, these four genes were actually identified as targets of TP53 (high t-statistics). Notice that comparing the standard and corrected OLS estimators, it is possible to conclude that the t-statistics are different probably due to the biased standard OLS estimator in the presence of measurement error. Table 6 shows the application of both, standard and proposed VAR models in a set of well known five genes related to circadian rhythm, namely, CLOCK, CRY2, PER2, PER3 and DBP. The genes CRY2, PER2, PER3 and DBP are known to be regulated by the complex BMAL1-CLOCK in mammals [24]. A VAR process of order one was adjusted and applied in a multivariate manner. Notice that also in the time series data, the estimators presented different results due to measurement error.

OLS model
PER3, PER3, PER3 Comparison of the usual and proposed methods in actual biological data is a difficult task since no one knows the "true" values. However, as observed in the simulation results, it is possible to conclude that the   Standard VAR average estimated coefficients and corrected VAR (between brackets) in 10,000 simulations. The first line contains the strength of association between predictors and response variables as described in simulation II. "-" means that it was not possible to calculate due to high measurement error when compared to sample's size. EM: Standard deviation of the Error of Measure. n: Number of samples. Notice that, in the presence of measurement error, the coefficients (b) estimated by the corrected VAR (between brackets) converge to the "true" values, while the estimated by standard VAR do not.      corrected approaches provide more reasonable results than biased standard methods.
In order to uncover more details about the performance of both, OLS and VAR, other experiments were conducted. These experiments consist in adding correlation in the residues and testing other null hypothesis (data not shown). The results obtained ignoring the errors by these methods can be compiled as: Therefore, these results demonstrate that improved methods to construct regulatory networks become necessary, since it is known that genes belong to an intrincate network, i.e., the covariates may be correlated and, moreover, gene expression quantification processes such as microarray technology measure with considerable error. If these conditions are ignored, one may obtain distort results and, consequently, conclude that there is a relationship between gene expressions where there is no association.
Construction of large networks is a challenge in bioinformatics. The methods proposed here do not allow the identification of networks when the number of variables is larger than the number of observations. Increasing the number of variables, the estimates become imprecise and the chances of obtaining multicollinearity problems also increases. In the presence of multicollinearity, one may use a feature selection procedure such as a stepwise (forward or backward, for example) in order to choose the optimum set of predictors.
Analyzing Pearson correlation coefficient, one can observe that it is simply a normalized linear regression coefficient (OLS) between -1 and 1. Therefore, Pearson correlation-based methods such as Relevance networks [14] or Graphical Gaussian models [26] need further studies in order to evaluate if they are also superestimating the rate of false positives and attenuating the coefficients like OLS. Moreover, Pearson correlation is widely used in order to test linear correlation between a certain gene expression signal and another characteristic such as prognostic, phenotype, tumor grade etc. Since these covariates may be measured with error, it is also crucial to develop a corrected Pearson correlation.
In order to develop a corrected Pearson correlation for measurement error, verify that it is possible to use the improved OLS presented in this report. The corrected Pearson correlation (r) between two random variables X and Y, both measured with error is given by where b should be estimated by using the corrected OLS (i.e., by simultaneously considering the three equations: y = a + bx + ε, X = x + 1 and Y = y + 2 ), s X and s Y are the standard deviations of the observed variables X and Y, respectively, and σ 1 and σ 2 are the standard deviations of the error of measure 1 and 2 , respectively. In this way, the estimate of the corrected version of the Pearson correlation is consistent (the larger the sample size, the smaller estimation error tends to be). Notice that, the difference between the corrected and  Although the examples provided here only treat regulatory network models, the proposed approaches can be applied in a straightforward manner also to estimate linear relationships between random variables measured with error.

Conclusions
Unfortunately, avoiding measurement error in a complete manner is a very difficult task, however, it can be minimized in the measuring (experimental) process and treated during the data analysis step. Here, we have shown evidence that presence of the measurement errors has a high impact in regulatory network models. In order to overcome this problem, approaches in both major data conditions, independent and time series data were proposed in addition to measurement error estimation procedures. Further studies are necessary in order to verify how is the performance of other regulatory networks (Bayesian networks, Structural Equation models, Graphical Gaussian models, Relevance networks, etc) in the presence of measurement error.

Methods
In this section, standard Ordinary Least Squares and Vector Autoregressive models will be described. Furthermore, corrected methods for measurement error will also be presented. Finally, the model used in the simulations will be detailed.

Ordinary least squares
In a multivariate regression model, let x 1 , x 2 ,..., x p be p predictor variables (genes) possibly being related to a response variable y (gene). The conventional linear regression model states that gene y is composed of an intercept or constant a which is the basal expression level of y, the predictors or gene expressions x j 's (j = 1,..., p) which relationship with y is represented by b = (b 1 ,..., b p ) ⊤ (the sign of b j represents the relationship between y and x j , i.e., positive or negative association), and a random error ε, which accounts for an intrinsic biological variation (this is not the measurement error).
for i = 1,..., n. In matrix notation, it is described as where The entire vector of error terms, ε i = (ε i1 ,..., ε iq ) ⊤ , are assumed to be independent and identically distributed as a q-variate normal distribution with zero vector mean and positive definite covariance matrix Σ ε for all i = 1,..., n, where q is the number of response variables. Notice that the proposed method is considering the homoscedastic case, i.e., the covariance matrix Σ ε does not change with i. Let Σ yx , Σ xx and Σ yy be the covariances of (y, x), (x, x) and (y, y), respectively. These covariance matrices could be estimated by:ˆ( Then, the intercept a is estimated aŝα β = − y x (14) and the estimator for the model's coefficient is given bŷˆβ The asymptotic variance-covariance matrix of vec(β ) and its estimate are given, respectively, by w h e r e ⊗ i s t h e K r o n e c k e r p r o d u c t a n d  The asymptotic distribution of vec(β ) is given by 19) and the test is described by: This test may be performed using the Wald-type statistic expressed as where C is a matrix of contrasts (usually, C = I). For more details about the matrix of contrasts, see [27]. Under the null hypothesis, (20) has a limit c 2 (d) distribution, where d = rank(C) gives the number of linear restrictions.
Ordinary least squares with measurement error Now, we shall study models of the regression type where one is unable to observe expression values of genes x and y (as described before) directly. Instead of observing x and y, one observes the sum and Y y = + 2 (22) with where 1~N (0, ∑ ∑ 1 ) independent of 2~N (0, ∑ ∑ 2 ) with ∑ ∑ 1 and ∑ ∑ 2 known are called as measurement errors, i.e., the variation derived by the measurement process (for example, the measurement error introduced when analyzing microarrays), ε~N(0, Σ ε ) is the random error (intrinsic biological variation) and x~N (μ x , Σ xx ), y~N (μ y , Σ yy ) with μ y = a + bμ x and Σ yy = bΣ xx b ⊤ + Σ ε .
The matrices ∑ ∑ 1 and ∑ ∑ 2 are given by ∑ ∑ i.e., the measurement errors may be different for each variable. Notice that the components of the measurement error's vector may be correlated but the entire vectors are independent.
Let Σ YX , Σ XX and Σ YY be the sample covariances of (Y, X), (X, X) and (Y, Y), respectively. These covariance matrices could be estimated by substituting x and y by X and Y in equations (8)(9)(10)(11)(12). Then, the intercept a is estimated aŝα and the estimator for the model's coefficient is given bŷˆβ Notice that∑ XX is estimated using equation (10) and ∑ 1 must be known a priori (it can be estimated using the procedures described in the section "Measurement errors estimation").
The asymptotic variance-covariance matrix of vec(β ) and its estimate are given, respectively, by (the proof is in the Appendix) andˆˆ(ˆ) where I q denotes the q × q identity matrix and Notice that, in the absence of measurement error, i e., ∑ ∑ ∑ ∑ 1 2 0 = = the corrected OLS is exactly equal to standard OLS. Furthermore, it is noteworthy that this asymptotic variance is similar to the one presented by [23] but in a multivariate manner with no correlation in the errors.

Hypothesis testing
Similar to the OLS with no measurement error, the interest in a simple regression model (y i = a + bx i + ε i ) lies in testing the strength of the relationship between the predictor gene x and the response gene y. The asymptotic distribution of vec(β ) is given by and the test is similar to the previous case (standard OLS) described by: This test may be performed using the Wald-type statistic expressed as where C is a matrix of contrasts. Under the null hypothesis, (33) follows a c 2 distribution with rank(C) degrees of freedom.

Vector autoregressive model
Here we define the usual VAR model as defined in Lütkepohl (2006) [28]. Vector autoregressive model with measurement error Now, the VAR model with measurement error will be presented.
Let z t be the "true" variables that are not directly observed. Let Z t be the observed surrogate variables which have an additive structure given by where Z t = (Z 1t , Z 2t ,..., Z pt ) ⊤ is the surrogate vector and t = ( 1t , 2t ,..., pt ) ⊤ is the measurement error vector. In most cases, if the usual conditional ML estimator is adopted for the observations subject to errors, i.e., replacing z t with Z t in the equation (34), the estimator of b will be biased as well as its asymptotic variance. Therefore, in order to overcome this limitation the measurement errors should be included in the estimation procedure. Nevertheless, the model (34) plus the equation (48) is not identifiable, since the covariance matrices of ε t and t are confounded. This problem can be avoided considering known the variance of t .
Let t~N (0, Σ ) be the measurement error with Σ known (refer to section Measurement error estimation for details about how to estimate Σ ). Then, the parameters of the model (34) of zeros with one's in the |l| th diagonal above (below) the main diagonal if l > 0 (l < 0) and J 0 is a (r × r) matrix of zeros.
Notice that, if r = 1 we have the VAR(1) model and the asymptotic covariance simplifies to The i th element of vec(β ), is asymptotically normally distributed with standard error given by the square root of i th diagonal element of ∑β . Thus, we can construct hypotheses testing on the individual coefficients, or in more general form of contrasts where C is a matrix of contrasts (C = I, for instance) and m is usually a (p × 1) vector or zeros.
Under the null hypothesis, (59) has a limiting c 2 (d) distribution where d = rank(C) gives the number of linear restrictions. This test is useful to identify, in a statistical sense (controlling the rate of false positives), which gene (predictor variable) is Granger causing another gene (response variable).

Measurement error estimation
Here, two methods to estimate measurement error are proposed. One when technical replicates are available and another one in the case when they are not available.

Technical replicates
When technical replicates are available, measurement error estimation may be performed by applying a strategy extending the methods described by Dahlberg (1940) [30] (more details about Dahlberg's method in the Appendix). For microarray data, it is known that the variance varies along the spots (heteroscedasticity) due to variations in experimental conditions (efficiency of dye incorporation, washing process, etc) [31]. Moreover, it is known that Dahlberg's approach is not suitable in the presence of systematic errors. Therefore, the application of the Dahlberg's formula is not straightforward. In order to overcome this problem, we suggest the following algorithm [32].
Let W and W' be two microarrays, where W' is the technical replicate of W.
1. Perform a non-linear regression such as splines smoothing between log(W) and log(W'), i.e., log(W') = f(log(W)) + ε 1 . Notice that the logarithm was calculated as a variance stabilizer (due to the high variance observed in microarray data). This is a common practice in microarray analysis; 2. Apply again the splines smoothing between ε 1 2 and log(W), i.e., ε 1 2 = g(log(W)) + ε 2 ; 3. Calculateˆ( ( )) δ = g log W 2 . This is a possible estimate for the standard deviation of the measurement error. Notice that with this process, we obtain oneδ i for each spot i = 1,..., m, where m is the number of spots in the microarray, also in the presence of heteroscedasticity.

No technical replicates
However, unfortunately, technical replicates is not always available. To this case, we have developed a strategy based on negative control probes and housekeeping genes frequently provided in commercial microarrays. Technically, housekeeping genes and negative controls should not change their expression levels [33]. Therefore, any variation measured by them can be understood as measurement error. In order to overcome the problem of heteroscedasticity in microarrays, we present a method based on splines smoothing. The main idea of this method consists in estimating how much of the total variance (intrinsic biological variation + measurement error) is due to measurement error. The method is as follows: 1. Let S be the set of all probes in the microarray and H be the set of housekeeping genes and negative controls. Calculate the mean and variance for each probe of S and H; 2. Perform a splines smoothing in both sets of probes separately, i.e., a splines smoothing var(H) = f(mean (H)) + ε 1 and var(S\{H}) = g(mean(S\{H})) + ε 2 , where H is a matrix containing the expression values of each housekeeping gene and negative controls in each row and S\{H} is a matrix containing the expression values of the remaining set of probes in each row. The functions f and g may be represented by a linear combination of spline functions j j (·), i.e., they may be written as where d is the number of knots used in the spline expansion (d may be obtained by selecting the value that minimizes the Generalized Cross Validation). mean(H) and var(H) (or mean(S\{H}) and var(S\ {H})) are vectors containing the mean and variance values of each row of H (or S\{H}), respectively. In this step, the smoothed curvesf andĝ represent the estimated variance for each probe. Notice that the smoothed curve in housekeeping genes and negative controlsf represents the estimated measurement error for each gene expression level. Moreover, the smoothed curve in the remaining set of probesĝ represents the total variance (intrinsic biological variance + measurement error) for each gene expression level; 3. Divide the smoothed curvef (obtained in step 2) by the other smoothed curveĝ . Notice that this ratio (ˆ/f g ) is the estimation of measurement error in percentage of the total variance for each probe. With this percentage, it is possible to estimate the variance of the measurement error for each probe.

Simulations
In order to evaluate the behavior of both, standard and proposed methods, we have conducted two simulations in small, moderate and large samples sizes (50, 100, 200 and 400). Computations were performed on the R software (a free software environment for statistical computing and graphics) [34]. For each group of simulation, 10,000 Monte Carlo samples were generated. Simulation I is for independent data and Simulation II for time series data.

Simulation I -independent data
In order to evaluate the performance of both, usual and corrected OLS methods, a controlled structure was defined. Let x and y be gene expression values where one is interested in examining if a certain gene x i (i = 1,..., p; p = 9) is linearly correlated to gene y (q = 1) partialized by other genes. This situation can be represented by the following structure where where x i~N (0, 1), ε~N (0, Σ ε ) is the intrinsic biological random variation and 1~N (0, ∑ ∑ 1 ) independent of 2~N (0, ∑ ∑ 2 ) are the measurement errors, with ∑ ∑ 1 = ∑ ∑ 2 varying from 0 to 0.8. The standard deviation Σ ε is defined by . . . . . 9 9 1 0 2 0 2 0 2 1 In order to become the simulation more realistic (since actual biological gene expression signals are generally quite correlated), notice that Σ ε is not a diagonal matrix, i.e., there are little correlations between the predictors. The sample's size varied from 50 to 400.

Simulation II -time series data
In the time series case, the data has some peculiarities which are not present in the independent data. Time series data are known to be autocorrelated (past values associated with future values) and also contemporaneously correlated (contemporaneous correlation between time series). Considering these characteristics, a similar structure described in the previous section was designed. Let X t and Y t being gene expression time series data and one is interested in verifying if certain gene x i, t (i = 1,..., p; p = 9) Granger causes gene y t (q = 1). This problem can be modeled by a VAR process of order one as described below: 0  1  1 1 1  2 2 1  3 3 1  4 4 1  5 5 1   6 , , ,  . , The observed variables X t and Y t are defined by The time series length varied from 50 to 400.
Notice that b 0 is the autoregressive coefficient and all time series X i for i = 1,...,9 are autocorrelated and also contemporaneously correlated (Σ ε is not a diagonal matrix).

Actual biological data
The standard and proposed OLS methods were applied to lung cancer gene expression data collected by [35]. This dataset is composed of 400 microarrays, each of which constructed using a different cDNA obtained from a different patient. Standard and corrected VAR approaches were applied to mouse liver time series data collected by [36]. This data is composed by 48 time points distributed at intervals of 1 hour.
Proof of the asymptotic variance of b -equation (29) Here, we proof equation (29), i.e., the asymptotic variance of b in the multivariate case with no correlated errors.
Consider the following model: where X i and Y i are the observed vectors with dimensions p × 1 and q × 1, respectively, a is the model intercept (q × 1), b is a (q × p) matrix of slope parameters, ε i is a white noise vector with mean zero and covariance matrix Σ ε . The joint distribution of 1i , 2i , ε i and x i is given by In this section, we investigate the asymptotic distribution ofˆˆ,

X X X X
The proof idea, similar to presented in [23], has two main steps. The first step consists in showing that vec (β ) -vec(b ⊤ ) can be written as linear combinations of a vectorial mean. In the second one, we must demonstrate that this vectorial mean has an asymptotic normal distribution. Therefore, we need some auxiliar results for proving the asymptotic result, which are exposed in two propositions below.  That is, x 1 ..., x n is an iid sequence and we can use the central limit theory, which says that ) . As nδ W is asymptotically normally distributed for all δ ≠ 0 r then, by the Cramer-Wold device [37], we have that Then, by the Propositions (1) and (2), we have that Dahlberg's error Consider the following model: where Z ij is the measure obtained in one experiment (microarray), i is the sample index i = 1,..., m, m is the number of spost in the microarray, j is the replicate number (j = 1, 2 in the case of duplicates), μ i is the unknown true value of the measure and ij is the error of measure.
Then, assume that E( ij ) = 0 and Var( ij ) = δ 2 . Thus, one quantification of the quality of measure is the standard deviation of ij , i.e., δ . Notice that the lower is the standard deviation of the error of measure (δ ), the lower is the measurement error.
Assuming that there is no bias (systematic error), one intuitive estimator for 2 2 δ is is exactly the Dahlberg's formula proposed in [30].