\documentclass[article]{ajs} \usepackage{amssymb,amsmath} %\newcommand{\E}{\mbox{E}} \newcommand{\var}{\mbox{var}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% almost as usual \author{Rohmatul Fajriyah\\ Graz University of Technology and Universitas Islam Indonesia} %\And %Second Author\\Plus Affiliation} \title{A Study of Convolution Models for Background Correction of BeadArrays} %% for pretty printing and a nice hypersummary also set: \Plainauthor{Rohmatul Fajriyah} %% comma-separated \Plaintitle{A Study of Convolution Models for Background Correction of BeadArrays} %% without formatting \Shorttitle{Convolution Models for Background Correction of BeadArrays} %% a short title (if necessary) %% an abstract and keywords \Abstract{ The robust multi-array average (RMA), since its introduction in \cite{Iri03a, Iri03b, Iri06}, has gained popularity among bioinformaticians. It has evolved from the exponential-normal convolution to the gamma-normal convolution, from single to two channels and from the Affymetrix to the Illumina platform. The Illumina design provides two probe types: the regular and the control probes. This design is very suitable for studying the probability distribution of both and one can apply a convolution model to compute the true intensity estimator. In this paper, we study the existing convolution models for background correction of Illumina BeadArrays in the literature and give a new estimator for the true intensity, assuming that the intensity value is exponentially or gamma distributed and the noise has lognormal distribution. Our study shows that one of our proposed models, the gamma-lognormal with the method of moments for parameters estimation, is the optimal one for the benchmarking data set with benchmarking criteria, while the gamma-normal model has the best performance for the benchmarking data set with simulation criteria. For the publicly available data sets, the gamma-normal and the exponential-gamma models with maximum likelihood estimation method can not be used and our proposed models exponential-lognormal and gamma-lognormal have the best performance, showing a moderate error in background correction and in the parametrization. } \Keywords{convolution, background correction, BeadArrays} \Plainkeywords{convolution, background correction, BeadArrays} %% without formatting %% at least one keyword must be supplied %% publication information %% NOTE: Typically, this can be left commented and will be filled out by the technical editor %% \Volume{50} %% \Issue{9} %% \Month{June} %% \Year{2012} %% \Submitdate{2012-06-04} %% \Acceptdate{2012-06-04} %% \setcounter{page}{1} \Pages{1--xx} %% The address of (at least) one author should be given %% in the following format: \Address{ Rohmatul Fajriyah\\ Institute of Statistics\\ Graz University of Technology\\ 8010 Graz, Austria\\ E-mail: \email{fajriyah@student.tugraz.at}\\ Department of Statistics\\ Universitas Islam Indonesia\\ Jl. Kaliurang Km. 14.4 \\ 55584 Jogjakarta, Indonesia\\ E-mail: \email{rfajriyah@fmipa.uii.ac.id}\\ URL: \url{http://rfajriyah.staff.uii.ac.id/} %URL: \url{http://www.statistik.tuwien.ac.at/public/templ} } %% It is also possible to add a telephone and fax number %% before the e-mail in the following format: %% Telephone: +43/512/507-7103 %% Fax: +43/512/507-2851 %% for those who use Sweave please include the following line (with % symbols): %% need no \usepackage{Sweave.sty} %% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} %% include your article here, just as usual \section{Introduction} \label{sec1} There are various processes in producing data from microarray experiments and each process contributes noise to the data. The noise can be of two types, biological and non-biological. Non-biological noise should be avoided or at least minimized. Sources for the non-biological noise are, for example, the chip itself, the scanner, or fluctuations in the electric network. Therefore, the data needs to be adjusted. The pre-processing will adjust the intensity value \citep{Hub05b,Hub05a} and provides an estimate of the true intensity. To estimate the intensity value, researchers proposed additive and multiplicative models and also {\it additive-muliplicative error models}, see e.g.\ \cite{Hub04}. In case of additive models, the underlying distribution is generally chosen as normal (log-normal), exponential, or a Gamma-$t$ mixture in the parametric approach (\cite{All09}, \cite{Bol03}, \cite{Che11}, \cite{Hoc06}, \cite{Iri03a,Iri03b,Iri06}, and \cite{Pla11,Pla12}). \cite{Iri03a,Iri03b,Iri06} and \cite{Bol03}, on the Affymetrix platform, have estimated the true intensity values based on a convolution model in the background correction step of their robust multi-array average (RMA) pre-processing method. They assumed that the true intensity is exponentially distributed and the background noise is normally distributed. \cite{Pla11,Pla12} showed that the RMA model (in \cite{Bol03} and \cite{Iri03a,Iri03b,Iri06}) does not fit Illumina BeadArrays: using the exponential-normal convolution leads to a large distance between the observed and the modeled intensities. They proposed, instead, the implementation of a gamma distribution for the intensity value and normal distribution for the noise. The simulation study of \cite{Pla11,Pla12} showed that the gamma-normal model performs better than the existing exponential-normal convolution model, giving a more accurate and correct fit for the observed intensities in Illumina BeadArray. Using a Gamma distribution for the intensity values in Illumina BeadArrays has been first suggested by \cite{Xie09}. The studies of \cite{Bae07} (on the background correction of the image processing) and \cite{Che11} show that the noise distribution is usually skewed in different degrees. In their studies, based on simulated and real data sets, \cite{Bae07} conclude that the gamma distribution is well suited for the noise. It accounts for the intensities with a positive lower bound and is very flexible in its shape, including asymmetric exponential type and symmetric normal type. The proposed convolution of exponential-gamma distribution by \cite{Che11} improves the intensity estimation and the detection of differentially expressed genes in the case when the intensity to noise ratio is large and the noise has a skewed distribution. In view of the remarks above, it is natural to model both the true intensity and the background noise in Illumina BeadArrays as gamma distributed. In an earlier version of this paper we have developed an estimator for the true intensity based on the gamma-gamma convolution model of RMA. However, this model does not fit very well the Illumina benchmarking data set. Independently, \cite{Tri13} proposed and applied the gamma-gamma model to pre-process Illumina methylation arrays. In this paper we introduce a new model for background correction in Illumina BeadArrays where the true intensity value is exponentially or gamma distributed and the noise has a lognormal distribution. As we will see, this model avoids the difficulties with the gamma-gamma model and has an overall satisfactory performance. We note that a new method reducing the bias of the maximum likelihood estimator of the shape parameter of the gamma distribution was proposed by \cite{Zha13}. But since our samples are very large, bias is not a problem in our studies. We compare the performance of the models on the Illumina spike-in data set, based on various criteria: root and mean square error (RMSE), $L_1$ error, Kullback-Leibler (K-L) coefficient, and some adapted criteria from Affycomp \citep{Cop04}. These criteria are measuring the reproducibility, accuracy, precision, specificity, and sensitivity of the expression measure of each model. We then provide a simulation study to measure the consistency of the error of background correction and the parametrization. The description and some details on these criteria and simulation can be found at \url{http://rfajriyah.staff.uii.ac.id/category/supplementary/} Our paper is organized as follows. In Section~\ref{sec3} we review previous work related to the background correction for Illumina BeadArrays. Our proposed models are described in Sections~\ref{sec41} and \ref{sec42}. Section~\ref{sec43} explains the benchmarking and the simulation studies on the benchmarking data set and Section~\ref{sec44} compares the performance of all models in the public data sets. Finally, Section~\ref{sec5} states the conclusions and indications of future work. \section{Previous work}\label{sec3} Affymetrix is the pioneer and most widely used platform for microarray gene expression experiments. The tools and algorithms to handle the data are numerous, both free and commercial. Some methods for pre-processing are available. Examples for the background correction step are: MAS5.0 by Affymetrix, multiplicative model based expression index (MMBE) by \cite{LiW01}, RMA in \cite{Iri03a,Iri03b,Iri06} and \cite{Bol03}, GC-RMA by \cite{Wu04} and maximum likelihood estimation based on the normal-exponential convolution model by \cite{Sil09}. Illumina is one of the alternative platforms and is increasingly popular. A few statistical methods have been developed for BeadArray data and there is no consensus yet for the pre-processing steps \citep{Shi10}. \cite{Xie09} mention that for the background correction step, Illumina bead studio gives two options (no background correction and background substraction) and the packages for BeadArrays in R provide three options (no background correction, background substraction and RMA background correction). \cite{Din08} extended the RMA model by proposing the model-based background correction method (MBCB) and showed that their model leads to a more precise determination of the gene expression and a better biological interpretation of Illumina BeadArray data. The studies of \cite{Che11} and \cite{Pla11,Pla12} show that their background correction models are made by adapting the RMA Affymetrix model. As \cite{For12}, pointed out, most preprocessing methods for Illumina BeadArrays are taken from the Affymetrix microarray platform. In general, the background correction is applied toward each array, where in each array there are probes, probesets and genes (terminology for the Affymetrix platform) or bead and bead-type level probes (terminology for the Illumina platform). At the Illumina platform, each gene is only targeted by one bead-type, which has been represented by about 30 time replications. If we can have a raw benchmarking data set, then it is possible to have all bead-type level probes of the raw data intensities. The current publicly available benchmarking data set for the Illumina platform is the raw data from the bead studio, which is the average of the bead-type level probes, not background corrected and of unnormalized intensity. Therefore, the background correction in this paper is applied to the gene intensity in each array. Suppose we have $J$ arrays and for each array there are $I$ regular genes. Our convolution model is as follows: \begin{equation} \label{eq31} P_{ij} = S_{ij} + B_{ij}\,, \qquad i = 1, \dots,I, \quad j = 1, \dots,J\,, \end{equation} where $P_{ij}, S_{ij}$, and $B_{ij}$ are the observed signal, true signal, and noise intensity, respectively. The $S_{ij}$ are i.i.d.\ and so are the $B_{ij}$; the $S_{ij}$ are independent of the $B_{ij}$. The $P_{ij}$ are observable. Our task is to recover the unknown signals $S_{ij}$ from the $P_{ij}$. To do this, for each array $j$ we also have $M$ observable noise intensities $B_{0mj}$, $m = 1, \dots, M$; the $B_{0mj}$ are i.i.d.\ with the same distribution as the $B_{ij}$ and are independent of all of the $S_{ij}, B_{ij}$. The $S_{ij}$ and $B_{ij}$ have a known type of distribution (exponential, gamma, normal, lognormal) with unknown parameters. \subsection{Background correction by RMA} \label{sec31} The RMA method was developed for the Affymetrix platform, where the design of arrays is different from the Illumina one. In this platform, one has perfect match and mismatch probe design. The RMA uses only the perfect match probe, which is the targeted probe of the intended gene. For those not familiar with the Affymetrix platform, refer to \cite{Bol04}, \cite{Bol03}, and \cite{Iri03a,Iri03b,Iri06} for further information. In modeling the intensity values, the RMA model (\cite{Bol03} and \cite{Iri03a,Iri03b,Iri06}) assumes that the intensity values are affected by the noise of the chip. By referring to Equation~(\ref{eq31}), in the RMA model $P_{ij}$ is the observed bead-type level probe intensity, $S_{ij}$ is the true signal with $S_{ij} \sim f_1(s_{ij}; \theta_j) = \textrm{Exp}(\theta_j)$, $\theta_j > 0$, and $B_{ij}$ is the background noise of the chip with $B_{ij} \sim f_2(b_{ij}; \mu_j, \sigma_j^2) = \mathcal{N}(\mu_j, \sigma_j^2)$, $\mu_j \in \mathbb{R}$, $\sigma_j^2 > 0$. To avoid negative intensity values, we truncate $B_{ij}$ at $0$ from below, i.e.\ we replace $B_{ij}$ by $\max\{B_{ij}, 0\}$; this will not change its density function $f_2(b_{ij}; \mu_j, \sigma_j^2)$ for $b_{ij} > 0$. Assuming independence, the joint density of the two-dimensional random variable $(S_{ij}, B_{ij})$ is \begin{displaymath} f_{S_{ij}, B_{ij}}(s_{ij}, b_{ij}; \mu_j, \sigma_j^2, \theta_j) = \theta_j\exp\Big\{-\theta_j s_{ij}\Big\}f_2(b_{ij}; \mu_j, \sigma_j^2), \qquad b_{ij}, s_{ij} > 0\,. \end{displaymath} Furthermore, the transformation formula for two-dimensional densities gives that joint density of $S$ and $P$ is \begin{multline} \label{eq33} f_{S_{ij}, P_{ij}}(s_{ij}, p_{ij}; \mu_j, \sigma_j^2, \theta_j) \\ = \theta_j\exp\left\{\frac{\theta_j^2\sigma_j^2}2 - \theta_j(p_{ij} - \mu_j)\right\} f_2(s_{ij}; \mu_{S.P,j}, \sigma_j^2), \qquad 0 < s_{ij} < p_{ij}, \end{multline} where $\mu_{S.P,j} = p_{ij} - \mu_j - \theta_j\sigma_j^2$. From (\ref{eq33}) we get the marginal density of $P_{ij}$ and the conditional density of $S_{ij}$ given $P_{ij}$ as \begin{gather*} f_{P_{ij}}(p_{ij}; \theta_j, \mu_j, \sigma_j^2) = \theta_j\exp\left\{\frac{\theta_j^2\sigma_j^2}2 - \theta_j(p_{ij} - \mu_j)\right\} \left(\Phi\left(\frac{\mu_{S.P,j}}{\sigma_j}\right) + \Phi\left(\frac{p_{ij}-\mu_{S.P,j}}{\sigma_j} \right) - 1\right) \\ f_{S_{ij} | P_{ij}}(s_{ij} | p_{ij}; \mu_{S.P,j}, \sigma_j^2) = \frac{f_2(s_{ij}; \mu_{S.P,j}, \sigma_j^2)} {\Phi\left(\frac{\mu_{S.P,j}}{\sigma_j}\right) + \Phi\left(\frac{p_{ij}-\mu_{S.P,j}}{\sigma_j}\right) - 1}. \end{gather*} The background adjusted intensity is computed as the conditional expectation of the true signal given the observed intensity, i.e. \begin{displaymath} \E(S_{ij}|P_{ij} = p_{ij}) = \frac 1{\Phi\left(\frac{\mu_{S.P,j}}{\sigma_j}\right) + \Phi\left(\frac{p_{ij} - \mu_{S.P,j}}{\sigma_j}\right) - 1} \int_0^{p_{ij}} s_{ij}f_2(s_{ij}; \mu_{S.P,j}, \sigma_j^2) ds_{ij}. \end{displaymath} The substitution $s_{ij} = \mu_{S.P,j} + \sigma_j t$ yields \begin{multline*} \int_0^{p_{ij}} s_{ij}f_2(s_{ij}; \mu_{S.P,j}, \sigma_j^2) ds_{ij}\\ = \mu_{S.P,j} \left(\Phi\left(\frac{p_{ij} - \mu_{S.P,j}}{\sigma_j}\right) + \Phi\left(\frac{\mu_{S.P,j}}{\sigma_j}\right) - 1\right) + \sigma_j \left(\phi\left(\frac{\mu_{S.P,j}}{\sigma_j}\right) - \phi\left(\frac{p_{ij} - \mu_{S.P,j}}{\sigma_j}\right)\right) \end{multline*} and thus (see \cite{Bol03}) \begin{equation} \label{eq36} \E(S_{ij}|P_{ij}=p_{ij}) = \mu_{S.P,j} + \sigma_j \frac{\phi(\frac{\mu_{S.P,j}}{\sigma_j}) - \phi(\frac{p_{ij}-\mu_{S.P,j}}{\sigma_j})} {\Phi(\frac{\mu_{S.P,j}}{\sigma_j}) + \Phi(\frac{p-\mu_{S.P,j}}{\sigma_j}) - 1}. \end{equation} Note that modelling the noise as a truncated normal variable has the consequence that the noise equals 0 with a positive probability $p_0$, a rather unpleasant feature of the model. As pointed out in \cite{Xie09}, however, in practical cases $p_0$ is rather small, so this problem can be disregarded. To avoid this difficulty, one can model the noise as the absolute value of a $\mathcal{N}(\mu, \sigma^2)$ variable, which changes the calculations above. However, since in this paper we will provide a background correction model fitting the reality considerably better, we do not give the details here. \subsection{Exponential-normal MBCB} \label{sec32} \cite{Xie09} use the same underlying distributions in (\ref{eq31}) for background correction. The difference with the RMA (\cite{Bol03} and \cite{Iri03a,Iri03b,Iri06}) are \begin{enumerate} \item \cite{Xie09} use $+\infty$ as the upper bound of the integral to compute the marginal density function and the conditional expectation of the true intensity value. On the other hand, RMA uses $p$ as the upper bound of the integration. The background corrected intensity value of \cite{Xie09} is \begin{displaymath} \E(S_{ij}|P_{ij} = p_{ij}) = \mu_{S.P,j} + \sigma_j\frac{\phi(\frac{\mu_{S.P,j}}{\sigma_j})}{\Phi(\frac{\mu_{S.P,j}}{\sigma_j})}. \end{displaymath} \item Under the convolution model (\ref{eq31}), where the true intensity value is assumed exponentially distributed and the noise is normally distributed, we need to estimate the parameters $\theta_j$, $\mu_j$, and $\sigma_j^2$. \cite{Xie09} offer three estimation methods: the method of moments, maximum likelihood estimation, and a Bayesian approach. On the other hand, RMA applies the \emph{ad-hoc} method. \end{enumerate} \cite{Din08} use the exponential-normal convolution model to correct the background of the Illumina platform by using Markov chain Monte Carlo simulation. \subsection{Gamma-normal convolution} \label{sec33} \cite{Pla11,Pla12} introduced gamma-normal convolution to model the background correction of Illumina BeadArray. The model is based on the RMA background correction of Affymetrix GeneChip. \cite{Pla11,Pla12} assume that the intensity value is gamma distributed and the noise is normally distributed. Under model (\ref{eq31}), $f_{P_{ij}}$ is the convolution of $f_{S_{ij}}$ and $f_{B_{ij}}$. The background corrected intensity is computed as the conditional expectation of $S_{ij}$ given $P_{ij} = p_{ij}$, i.e. \begin{equation} \label{eq38} \E(S_{ij}|P_{ij}=p_{ij}) = \frac{\int s_{ij} f^{\textrm{gam}}_{\alpha_j, \beta_j}(s_{ij}) f^{\textrm{norm}}_{\mu_j, \sigma_j}(p_{ij}-s_{ij}) ds_{ij}} {\int f^{\textrm{gam}}_{\alpha_j, \beta_j}(s_{ij}) f^{\textrm{norm}}_{\mu_j, \sigma_j}(p_{ij}-s_{ij}) ds_{ij}}, \end{equation} where \begin{displaymath} f^{\textrm{gam}}_{\alpha_j, \beta_j}(s) = \frac{\beta_j^{\alpha_j} s^{\alpha_j - 1}\exp(-\beta_j s)}{\Gamma(\alpha_j)}, \qquad \alpha_j, \beta_j, s > 0 \end{displaymath} is the gamma density. When $S_{ij}$ is gamma distributed and $B_{ij}$ is normally distributed, then (\ref{eq38}) has no analytic expression like (\ref{eq36}). \cite{Pla11,Pla12} implemented the Fast Fourier Transform to estimate the parameters and to correct the background. For the background correction with Fast Fourier Transform, the corrected intensity (\ref{eq38}) is rewritten as \begin{displaymath} \E(S_{ij}|P_{ij}=p_{ij}) = \frac{\alpha_j\beta_j\int f^{\textrm{gam}}_{\alpha_j + 1, \beta_j}(s_{ij}) f^{\textrm{norm}}_{\mu_j, \sigma_j}(p_{ij}-s_{ij}) ds_{ij}} {\int f^{\textrm{gam}}_{\alpha_j, \beta_j}(s_{ij}) f^{\textrm{norm}}_{\mu_j, \sigma_j}(p_{ij}-s_{ij}) ds_{ij}}, \end{displaymath} since $s_{ij}f^\textrm{gam}_{\alpha_j, \beta_j}(s_{ij}) = \alpha_j\beta_jf^\textrm{gam}_{\alpha_j + 1, \beta_j}(s_{ij})$ is valid for every $s_{ij} > 0$. \newpage \subsection{Exponential-gamma convolution}\label{sec34} Under model (\ref{eq31}), \cite{Che11} proposed for the distribution of the true intensity and its noise the exponential and gamma distribution, respectively. Therefore, $S_{ij} \sim f_1(s_{ij}; \theta_j) = \textrm{Exp}(\theta_j)$ and $B_{ij} \sim f_2(b_{ij}; \alpha_j, \beta_j) = \textrm{Gamma}(\alpha_j, \beta_j)$, where $s_{ij}, b_{ij}, \theta_j, \alpha_j, \beta_j > 0$. The corrected background intensity of \cite{Che11} is \begin{displaymath} \E(S_{ij}|P_{ij} = p_{ij}) = p_{ij} - \frac{\int_0^{p_{ij}}b_{ij}^{\alpha_j}\exp(-(\beta_j - \theta_j)b_{ij}) db_{ij}} {\int_0^{p_{ij}}b_{ij}^{\alpha_j-1}\exp(-(\beta_j - \theta_j)b_{ij}) db_{ij}}. \end{displaymath} \section{Results}\label{sec4} Now we present the results of the two proposed convolution models: the exponential-lognormal convolution in Section \ref{sec41} and the gamma-lognormal convolution in Section \ref{sec42}. In each section, the formula for the background corrected intensity value is derived and methods to estimate the parameters are explained. Section \ref{sec43} present the benchmarking results, i.e.\ a performance comparison of all models at the Illumina spike-in data set. Section \ref{sec44} present the performance comparison of all models for the public data sets. The simulation study results are presented in Sections \ref{sec43} and \ref{sec44}. For further details on the benchmarking criteria, supplemental plots and simulations see {Supplementary\_Materials} at \url{http://rfajriyah.staff.uii.ac.id/category/supplementary/}. \subsection{Exponential-lognormal convolution}\label{sec41} \paragraph{Background correction} Consider model (\ref{eq31}) when the true intensity $S_{ij}$ is exponentially distributed, i.e.\ \begin{displaymath} S_{ij} \sim f_1(s_{ij}; \theta_j) = \theta_j\exp(-\theta_j s_{ij}), \qquad \theta_j, s_{ij} > 0, \end{displaymath} and the background noise $B_{ij}$ is lognormally distributed, i.e.\ \begin{displaymath} B_{ij} \sim f_2(b_{ij}; \mu_j, \sigma_j^2) = \frac 1{b_{ij}\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(\log b_{ij}-\mu_{j})^2}{2\sigma_j^2}\right), \qquad \mu_j \in \mathbb{R},\ \sigma_j^2, b_{ij} > 0. \end{displaymath} Then the joint density function of $S_{ij}$ and $B_{ij}$ equals \begin{displaymath} f_{S_{ij}, B_{ij}}(s_{ij}, b_{ij}) = \frac{\theta_j\exp(-\theta_j s_{ij})}{b_{ij}\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(\log b_{ij} -\mu_j)^2}{2\sigma_j^2}\right), \end{displaymath} and thus the joint density function of $S_{ij}$ and $P_{ij}$ is \begin{displaymath} f_{S_{ij}, P_{ij}}(s_{ij}, p_{ij}) = \frac{\theta_j\exp(-\theta_j s_{ij})}{(p_{ij}-s_{ij})\sigma_j\sqrt{2\pi}} \exp\left(\displaystyle -\frac{(\log (p_{ij}-s_{ij})-\mu_j)^2}{2\sigma_j^2}\right), \qquad s_{ij} < p_{ij}. \end{displaymath} Consequently, the marginal density function of $P_{ij}$ equals \begin{displaymath} f_{P_{ij}}(p_{ij}) = \int_0^{p_{ij}}\frac{\theta_j\exp(-\theta_j s_{ij})} {(p_{ij}-s_{ij})\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(\log(p_{ij}-s_{ij})-\mu_j)^2}{2\sigma_j^2}\right) ds_{ij}. \end{displaymath} Using the substitution $\log(p_{ij}-s_{ij}) = z_{ij}$ we get \begin{align*} f_{P_{ij}}(p_{ij}) &= \int_{-\infty}^{\log p_{ij}} \frac{\theta_j\exp(-\theta_j(p_{ij}-e^{z_{ij}}))}{\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(z_{ij}-\mu_j)^2}{2\sigma_j^2}\right) dz_{ij} \\ &= \frac{\theta_j \exp(-\theta_j p_{ij})}{\sigma_j\sqrt{2\pi}} \int_{-\infty}^{\log p_{ij}} \exp\left(-\frac{(z_{ij}-\mu_j)^2}{2\sigma_j^2}\right) \sum_{k=0}^{\infty}\frac{\theta_j^k e^{kz_{ij}}}{k!} dz_{ij} \\ &= \frac{\theta_j \exp(-\theta_j p_{ij})}{\sigma_j\sqrt{2\pi}} \sum_{k=0}^{\infty}\frac{\theta_j^k}{k!} \int_{-\infty}^{\log p_{ij}} \exp\left(-\frac{(z_{ij}-\mu_j)^2}{2\sigma_j^2} + kz_{ij}\right) dz_{ij} \\ &= \frac{\theta_j \exp(-\theta_j p_{ij})}{\sigma_j\sqrt{2\pi}} \sum_{k=0}^{\infty} \frac{\theta_j^k}{k!}\exp\left(k\left(\mu_j + \frac k2\sigma_j^2\right)\right) \int_{-\infty}^{\log p_{ij}}\!\!\exp\left(\!-\frac{(z_{ij}-(\mu_j+k\sigma_j^2))^2}{2\sigma_j^2}\!\right) dz_{ij} \\ &= \theta_j\exp(-\theta_j p_{ij}) C_{a,j}, \end{align*} where \begin{displaymath} C_{a,j} = \sum_{k=0}^{\infty} \frac{\theta_j^k}{k!}\exp\left(k\left(\mu_j + \frac k2\sigma_j^2\right)\right) \Phi\left(\frac{\log p_{ij}-(\mu_j+k\sigma_j^2)}{\sigma_j}\right). \end{displaymath} The conditional density function of $S_{ij}$ given $P_{ij}=p_{ij}$ is now obtained as \begin{displaymath} f_{S_{ij}|P_{ij}}(s_{ij}|p_{ij}) = \frac{\exp(\theta_j(p_{ij}-s_{ij}))} {(p_{ij}-s_{ij})\sigma_j\sqrt{2\pi}C_{a,j}} \exp\left(\displaystyle -\frac{(\log (p_{ij}-s_{ij})-\mu_j)^2}{2\sigma_j^2}\right) \end{displaymath} with conditional mean \begin{displaymath} \E(S_{ij}|P_{ij}=p_{ij}) = \frac{\exp(\theta_j p_{ij})}{C_{a,j}} \int_0^{p_{ij}}\frac{s_{ij}\exp(-\theta_j s_{ij})}{(p_{ij}-s_{ij})\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(\log(p_{ij}-s_{ij})-\mu_j)^2}{2\sigma_j^2}\right) ds_{ij}. \end{displaymath} Using the substitution $\log(p_{ij}-s_{ij}) = z_{ij}$, this equals \begin{align*} \E(S_{ij}|P_{ij}=p_{ij}) &= \frac{p_{ij}}{C_{a,j}} \int_{-\infty}^{\log p_{ij}} \frac{\left(1 - \frac{e^{z_{ij}}}{p_{ij}}\right) \exp(\theta_j e^{z_{ij}})}{\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(z_{ij}-\mu_j)^2}{2\sigma_j^2}\right)dz_{ij} \\ &= \frac{p_{ij}}{C_{a,j}} \left[C_{a,j} - \int_{-\infty}^{\log p_{ij}} \frac{\frac{e^{z_{ij}}}{p_{ij}}\exp(\theta_je^{z_{ij}})}{\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(z_{ij} -\mu_j)^2}{2\sigma_j^2}\right) dz_{ij}\right] \\ &= p_{ij} - \frac{\exp\left(\mu_j+\frac{\sigma_j^2}2\right)}{C_{a,j}} \int_{-\infty}^{\log p_{ij}} \frac{\exp(\theta_j e^{z_{ij}})}{\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(z_{ij}-(\mu_j+\sigma_j^2))^2}{2\sigma_j^2}\right) dz_{ij} \\ &= p_{ij} - \frac{\exp\left(\mu_j+\frac{\sigma_j^2}2\right)}{C_{a,j}} \sum_{k=0}^\infty \frac{\theta_j^k}{k!}\exp\left(k\left(\mu_j+\frac{k+2}2\sigma_j^2\right)\right) \\ & \qquad\qquad\times \int_{-\infty}^{\log p_{ij}} \frac 1{\sigma_{j}\sqrt{2\pi}} \exp\left(-\frac{(z_{ij}-(\mu_j+(k+1)\sigma_j^2))^2}{2\sigma_j^2}\right) dz_{ij} \\ &= p_{ij} - \frac{C_{b,j}}{C_{a,j}}\exp\left(\mu_j+\frac{\sigma_j^2}2\right), \end{align*} where \begin{displaymath} C_{b,j} = \sum_{k=0}^\infty \frac{\theta_j^k}{k!}\exp\left(k\left(\mu_j+\frac{k+2}2\sigma_j^2\right)\right) \Phi\left(\frac{\log p_{ij} - (\mu_j+(k+1)\sigma_j^2)}{\sigma_j}\right). \end{displaymath} \paragraph{Parameter estimation} To estimate the parameters $\theta_j, \mu_j$, and $\sigma_j^2$, $j = 1, \dots,J$ in the exponential-lognormal model we can use various methods. \begin{enumerate} \item Maximum likelihood estimation (MLE): \\ This is implemented by applying the \emph{optim} function in R to maximize the log-likelihood function of the $j$th array \begin{eqnarray*} \sum_{i=1}^I (\log \theta_j - \theta_jp_{ij} + \log C_{a,j;K}) + \sum_{m=1}^M \left(-\frac{(\log b_{0mj}-\mu_j)^2}{2\sigma_j^2} -\frac 12\log(2\pi\sigma^2_jb^2_{0mj})\right), \end{eqnarray*} where $p_{ij}$ and $b_{0mj}$ are the observed values of $P_{ij}$ and $B_{0mj}$. Note that $C_{a,j}$ in the log-likelihood function is defined by the infinite series at the end of the previous section. However, the terms of this infinite series decrease very rapidly and thus we can cut off the series at a proper index $K$ giving $C_{a,j;K}$, making it suitable for its computation in R. The index $K$ is chosen as the smallest integer for which $|(C_{a,j;K} - C_{a,j;K-1})/C_{a,j;K-1}| < 0.001$ holds. \item Method of moments: \\ Note that the method of moments estimator of the parameter $\theta$ in an exponential distribution is the reciprocal of the sample mean. Since the $S_{ij}$ are not observable, but the $B_{0mj}$ are, we consider Equation (\ref{eq31}) and estimate $1/\theta_j$ by the difference $\textrm{mean}(p_{ij}) - \textrm{mean}(b_{0mj})$. Further, the parameters $\mu_j$ and $\sigma_j^2$ in the lognormal part are estimated by the sample mean and variance of the observed $\log b_{0mj}$ values. \item Plug-in estimator: \\ \begin{enumerate} \item We calculate the MLE of the parameter in the model of $S_{ij}$ by utilizing the observed sample $p_{ij}$ instead of $s_{ij}$. This is justified by the fact that, in some sense, the distribution of $S_{ij}$ is similar to the distribution of $P_{ij}$. \item We estimate $\mu_j$ and $\sigma_j^2$ through MLE based on $B_{0mj}$ as described above. \end{enumerate} \end{enumerate} \subsection{Gamma-lognormal convolution} \label{sec42} \paragraph{Background correction} Consider now model (\ref{eq31}) and assume that the true intensity $S_{ij}$ is gamma distributed, i.e.\ \begin{displaymath} S_{ij} \sim f_1(s_{ij}; \alpha_j, \beta_j) = \frac{\beta_j^{\alpha_j} s_{ij}^{\alpha_j-1} \exp(-s_{ij}\beta_j)}{\Gamma(\alpha_j)}, \qquad \alpha_j, \beta_j, s_{ij} > 0, \end{displaymath} and that the background noise $B_{ij}$ is lognormally distributed, i.e.\ \begin{displaymath} B_{ij} \sim f_2(b_{ij}; \mu_j, \sigma_j^2) = \frac 1{b_{ij}\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(\log b_{ij} - \mu_j)^2}{2\sigma_j^2}\right), \qquad \mu_j \in \mathbb{R},\ \sigma_j^2 > 0. \end{displaymath} The joint density function of $S_{ij}$ and $B_{ij}$ is \begin{displaymath} f_{S_{ij}, B_{ij}}(s_{ij}, b_{ij}) = \frac{\beta_j^{\alpha_j} s_{ij}^{\alpha_j-1}\exp(-s_{ij}\beta_j)}{\Gamma(\alpha_j)} \frac 1{b_{ij}\sigma_j\sqrt{2\pi}}\exp\left(-\frac{(\log b_{ij} - \mu_j)^2}{2\sigma_j^2}\right) \end{displaymath} and thus the joint density function of $S_{ij}$ and $P_{ij}$ is \begin{displaymath} f_{S_{ij}, P_{ij}}(s_{ij}, p_{ij}) = \frac{\beta_j^{\alpha_j}s_{ij}^{\alpha_j-1}\exp(-s_{ij}\beta_j)}{\Gamma(\alpha_j)} \frac 1{(p_{ij} - s_{ij})\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(\log (p_{ij} - s_{ij}) - \mu_j)^2}{2\sigma_j^2}\right). \end{displaymath} Hence, the marginal density function of $P_{ij}$ is obtained as \begin{displaymath} f_{P_{ij}}(p_{ij}) = \int_0^{p_{ij}} \frac{\beta_j^{\alpha_j}s_{ij}^{\alpha_j-1}\exp(-s_{ij}\beta_j)}{\Gamma(\alpha_j)} \frac 1{(p_{ij} - s_{ij})\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(\log (p_{ij} - s_{ij}) - \mu_j)^2}{2\sigma_j^2}\right)ds_{ij}. \end{displaymath} Using the substitution $\log(p_{ij} - s_{ij}) = z_{ij}$ we get \begin{align*} f_{P_{ij}}(p_{ij}) &= \int_{-\infty}^{\log p_{ij}} \frac{\beta_j^{\alpha_j}p_{ij}^{\alpha_j-1}\left(1-\frac{e^{z_{ij}}}{p_{ij}}\right)^{\alpha_j-1} \exp(-p_{ij}\beta_j + e^{z_{ij}\beta_j})} {\Gamma(\alpha_j)\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(z_{ij} - \mu_j)^2}{2\sigma_j^2}\right) dz_{ij} \\ &= \frac{\beta_j^{\alpha_j} p_{ij}^{\alpha_j-1}\exp(-p_{ij}\beta_j)}{\Gamma(\alpha_j)\sigma_j\sqrt{2\pi}} \int_{-\infty}^{\log p_{ij}} \left(1 - \frac{e^{z_{ij}}}{p_{ij}}\right)^{\alpha_j-1}\exp(e^{z_{ij}\beta_j}) \exp\left(-\frac{(z_{ij} - \mu_j)^2}{2\sigma_j^2}\right) dz_{ij} \\ &= \frac{\beta_j^{\alpha_j} p_{ij}^{\alpha_j-1}\exp(-p_{ij}\beta_j)}{\Gamma(\alpha_j)}C_{c,j}, \end{align*} where \begin{displaymath} C_{c,j} = \sum_{k=0}^\infty\sum_{n=0}^\infty \frac{(-1)^k {\alpha_j-1 \choose k}}{p_{ij}^kn!} \beta_j^n \exp\left((k+n)\left(\mu_j + (k+n)\frac{\sigma_j^2}2\right)\right) \Phi\left(\frac{\log p_{ij}-(\mu_j+(k+n)\sigma_j^2)}{\sigma_j}\right). \end{displaymath} The conditional density function is now obtained as \begin{displaymath} f_{S_{ij}|P_{ij}}(s_{ij}|p_{ij}) = \frac{\exp((p_{ij}-s_{ij})\beta_j) s_{ij}^{\alpha_j-1}} {C_{c,j}p_{ij}^{\alpha_{j}-1}(p_{ij}-s_{ij}) \sigma_{j}\sqrt{2\pi}} \exp\left(-\frac{(\log(p_{ij} - s_{ij}) - \mu_{j})^2}{2\sigma_j^2}\right) \end{displaymath} with respective conditional mean \begin{displaymath} \E(S_{ij}|P_{ij} = p_{ij}) = \frac{e^{p_{ij}\beta_j}}{C_{c,j}p_{ij}^{\alpha_j-1}} \int_0^{p_{ij}}\frac{s_{ij}^{\alpha_j}e^{-s_{ij}\beta_j}}{(p_{ij}-s_{ij})\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(\log(p_{ij} - s_{ij}) - \mu_{j})^2}{2\sigma_j^2}\right) ds_{ij}. \end{displaymath} Substituting $\log(p_{ij} - s_{ij}) = z_{ij}$ this becomes \begin{eqnarray} \label{eq413} \E(S_{ij}|P_{ij} = p_{ij}) &=& \frac{p_{ij}}{C_{c,j}} \int_{-\infty}^{\log p_{ij}} \frac{\left(1-\frac{e^{z_{ij}}}{p_{ij}}\right)^{\alpha_j}\exp\left(e^{z_{ij}\beta_j}\right)} {\sigma_j\sqrt{2\pi}} \exp\left(-\frac{(z_{ij} - \mu_{j})^2}{2\sigma_j^2}\right) dz_{ij}\nonumber \\ &=& \frac{p_{ij}C_{d,j}}{C_{c,j}}, \end{eqnarray} where \begin{displaymath} C_{d,j} = \sum_{k=0}^\infty\sum_{n=0}^\infty \frac{(-1)^k{\alpha_j \choose k}}{p_{ij}^kn!} \beta_j^n\exp\left((k+n)\left(\mu_j + (k+n)\frac{\sigma_j^2}2\right)\right) \Phi\left(\frac{\log p_{ij}-(\mu_j+(k+n)\sigma_j^2)}{\sigma_j}\right). \end{displaymath} \paragraph{Parameter estimation} To estimate the parameters $\alpha_j, \beta_j, \mu_j$, and $\sigma_j^2$ in (\ref{eq413}), we can use either of the following methods. \begin{enumerate} \item Maximum likelihood estimation:\\ This is implemented by applying the \emph{optim} function in R to maximize the log-likelihood function of the $j$th array \begin{eqnarray*} && \sum_{i=1}^I\big(\log(C_{c,j;K}) +(\alpha_j-1)\log(p_{ij}) - p_{ij}\beta_j + \alpha_j\log(\beta_j) - \log(\Gamma(\alpha_j))\big) \\ && + \sum_{m=1}^M \left(-\frac{(\log b_{0mj}-\mu_j)^2}{2\sigma_j^2} -\frac 12\log(2\pi\sigma^2_jb^2_{0mj})\right). \end{eqnarray*} Similarly to the exponential-lognormal model, in the computation of $C_{c,j}$ the cutoff index $K$ is chosen according to the criteria $|(C_{c,j;K} - C_{c,j;K-1})/C_{c,j;K-1}| < 0.002$. \item Method of moments:\\ The implementation of this method for the $j$th array is done by recalling that in case of a gamma distribution with parameters $\alpha_j$ and $\beta_j$, the method of moments estimator for $\alpha_j/\beta_j$ and $\alpha_j/\beta_j^2$ is the sample mean and the sample variance, respectively. Thus, considering Equation (\ref{eq31}) we get $$ \hat\beta_j = \frac{\textrm{mean}(s_{ij})}{\textrm{var}(s_{ij})} = \frac{\textrm{mean}(p_{ij})-\textrm{mean}(b_{0mj})}{\textrm{var}(p_{ij})-\textrm{var}(b_{0mj})}, $$ as also $$ \hat\alpha_j = \hat\beta_j \ \textrm{mean}(s_{ij}) = \frac{(\textrm{mean}(s_{ij}))^2}{\textrm{var}(s_{ij})} = \frac{(\textrm{mean}(p_{ij})-\textrm{mean}(b_{0mj}))^2}{\textrm{var}(p_{ij})-\textrm{var}(b_{0mj})}. $$ Furthermore, $\mu_j$ and $\sigma_j^2$ are estimated by the mean and variance of $\log b_{0mj}$. \item Plug-in estimator: \\ In equation (\ref{eq31}), $P_{ij}$ and $B_{0mj}$ are observable intensities. Therefore, the plug-in estimator is implemented by \begin{enumerate} \item estimating $\alpha_j$ and $\beta_j$ through MLE based on the $p_{ij}$ values, and \item estimating $\mu_j$ and $\sigma_j^2$ through MLE based on the $b_{0mj}$ values. \end{enumerate} \end{enumerate} \subsection{Benchmarking} \label{sec43} \paragraph{Benchmarking data set} Illumina platform has provided a benchmarking data set, the Illumina spike-in \cite{Dun08a}. These spike-in probes are targeting bacterial and viral genes absent from the mouse genome. These were added at specific concentrations on each sample. Therefore the change in expression level of a particular spike between samples is known a priori. The expression levels of the non-spikes should not change between samples. There are twelve different concentrations of spike: 1000 picomolar (pM), 300, 100, 30, 10, 3, 1, 0.3, 0.1, 0.03, 0.01, and 0.00 pM. It was replicated four times. Therefore, there are 48 samples and each sample has regular and control bead-type level probes. There are approximately about 48000 bead-type level probes for each sample and in addition the 33 spike-in bead-type level probes are added into it. For the control, there are 1616 bead-type level probes. These control experiments are the benchmarking data sets of Illumina and are used to compare low-level analysis methods such as in Affymetrix platform. \paragraph{Performance studies} We compare all convolution models: \cite{Iri03a,Iri03b,Iri06} and \cite{Bol03}: RMA (Exponential-Normal), \cite{Pla11,Pla12}: Gamma-Normal, \cite{Che11}: Exponential-Gamma, \cite{Xie09}: Exponential-Normal adjusted for Illumina BeadArrays with MLE for the parameters, Bayesian approach and the method of moments, and the proposed models: exponential-lognormal and gamma-lognormal. We will call the methods above, respectively, as follows: ENr, GN, EG, ENm, ENm, ENn, ELNn, ELNm, ELNp, GLNn, GLNm, and GLNp. We use the MBCB package (\cite{All09} and \cite{Xie09}) to adjust the intensity values of these existing models ENr, ENm, ENmc and ENn. Except that, the GN uses the NormalGamma package (\cite{Pla11}). \begin{table}[!t] \begin{center} \caption{Reproducibility of each method toward the Illumina spike-in concentration}\label{Tab1} \begin{tabular}{lrr} \hline Model & RMSE & K-L\\ \hline ENr & 1.346 & 51310 \\ ENn & 1.407 & 41010 \\ ENm & 1.483 & 23170 \\ ENmc & 1.483 & 23170 \\ EG & 1.470 & 20660 \\ GN & 1.521 & 58480 \\ ELNn & 1.411 & 41200 \\ ELNm & 1.489 & 21280 \\ ELNp & 1.423 & 37800 \\ GLNn & \textbf{1.323} & \textbf{4333} \\ GLNm & 1.510 & 29630 \\ GLNp & 10.700 & $-$115400 \\ \hline \end{tabular} \end{center} \end{table} \begin{table}[!h] \begin{center} \caption{Reproducibility of each method toward the Illumina spike-in based on the experiment data}\label{Tab2} \begin{tabular}{lrr} \hline Model & RMSE & K-L\\ \hline ENr & 7.251 & 1141000 \\ ENn & 7.127 & 1062000 \\ ENm & 6.927 & 926500 \\ ENmc & 6.927 & 926200 \\ EG & 6.919 & 907900 \\ GN & 7.100 & 1183000 \\ ELNn & 7.124 & 1062000 \\ ELNm & 6.904 & 911600 \\ ELNp & 7.092 & 1035000 \\ GLNn & \textbf{6.825} & \textbf{793400} \\ GLNm & 6.937 & 968400 \\ \hline \end{tabular} \end{center} \end{table} Table~\ref{Tab1} shows that the GLNn reproduces the Illumina concentration better than others. The ENr shows the closest performance toward the GLNn. Note that the computation of the Kullback-Leibler coefficient is implemented in each array $j$ based on the nominal concentrations $O$ in Table~\ref{Tab1} and the observed intensities $P$ in Table \ref{Tab2}, and the value in each table is the median Kullback-Leibler coefficient based on $J = 42$ arrays. The Kullback-Leibler coefficient for two positive sequences $(X_{1j}, \dots, X_{Ij}), (S_{1j}, \dots, S_{Ij})$ is computed as $\textrm{K-L}_j = \sum_{i=1}^I X_{ij}\log (X_{ij}/S_{ij})$, which can also be negative if $S_{ij} > X_{ij}$ for all or for most $i$. This is a sign that the $S$ are overestimating $X$, where $X$ could be $O$ (Table~\ref{Tab1}) or $P$ (Table~\ref{Tab2}). Therefore we exclude the GLNp model from further comparisons. The behavior of GLNp which is different from other models, also shown at the supplemental plots. Table~\ref{Tab2} shows how each method reproduces the data from the experiment. We see that GLNn can be considered to reproduce it better than others, based on the RMSE, and the Kullback-Leibler coefficient. Tables~\ref{Tab1} and \ref{Tab2} provide insight about how the performance comparison among the models would be conducted further. In the first part, we compute the adopted Affycomp benchmarking criteria, based on the data after background correction and their log transformation. In the second part, in the simulation, the MSE$_{bc}$ and the $L_1$ error will be computed based on the log transformation of the experiment and the nominal concentration data. The log transformations that we use here, respectively, for the benchmarking and the FFPE data sets are as follows \begin{displaymath} y = \log_2\left(x + \sqrt{x^2 + 1}\right) \qquad\textrm{and}\qquad y = \log_2\left(x + 1 + \sqrt{x^2 + 1}\right), \end{displaymath} where $x$ is the nominal concentration $O$ or the observed intensity value $P$. \paragraph{\textit{First part}} \label{sec4321} In Table~\ref{Tab3} it is shown that the ENr method provides the smallest variation and IQR and the GLNn model provides the smallest 99.9\% percentiles of log fold change for the non spike-in between replicates. The largest variation, IQR, and 99.9\% percentiles, respectively, are observed for the GLNm, the ELNm, and the GN method. \begin{table}[!t] \begin{center} \caption{Median SD, IQR, and 99.9\% percentiles of log fold change for non spike-in between replicates for each model.}\label{Tab3} \begin{tabular}{lccc} \hline Model & Median SD & IQR & 99.9\% \\ \hline ENr & \textbf{0.027} & \textbf{0.062} & 0.415 \\ ENn & 0.043 & 0.089 & 0.441 \\ ENm & 0.069 & 0.139 & 0.486 \\ ENmc & 0.069 & 0.139 & 0.486 \\ EG & 0.065 & 0.134 & 0.477 \\ GN & 0.051 & 0.098 & 0.520 \\ ELNn & 0.045 & 0.093 & 0.442 \\ ELNm & 0.071 & 0.145 & 0.489 \\ ELNp & 0.049 & 0.100 & 0.449 \\ GLNn & 0.038 & 0.075 & \textbf{0.398} \\ GLNm & 0.076 & 0.080 & 0.507 \\ \hline \end{tabular} % \end{center} %\end{table} % %\begin{table}[!h] % \begin{center} \caption{The signal detect $R^2$ by regressing the nominal and observed value for each model for the Illumina spike-in.}\label{Tab5} \begin{tabular}{lcrrr} \hline Model & $R^2$ & Low.$R^2$ & Med.$R^2$ & High.$R^2$ \\ \hline ENr & 0.959 & 0.618 & \textbf{0.698} & \textbf{0.559}\\ ENn & 0.958 & 0.622 & 0.695 & 0.557 \\ ENm & 0.957 & 0.635 & 0.695 & 0.558 \\ ENmc & 0.957 & 0.635 & 0.695 & 0.558 \\ EG & 0.957 & 0.633 & 0.695 & 0.558 \\ GN & 0.956 & \textbf{0.650} & 0.697 & 0.555 \\ ELNn & 0.958 & 0.624 & 0.695 & 0.557 \\ ELNm & 0.957 & 0.636 & 0.694 & 0.558 \\ ELNp & 0.958 & 0.627 & 0.695 & 0.557 \\ GLNn &\textbf{0.960} & 0.609 & 0.696 & 0.558 \\ GLNm & 0.956 & 0.637 & 0.694 & 0.558 \\ \hline \end{tabular} \end{center} \end{table} In Table~\ref{Tab5} it is shown that, in general, all methods perform similar to each other. The GLNn models have the highest signal detect $R^2$. The GN model has the highest $R^2$ at low concentration but has the lowest $R^2$ at high concentration. This means that the GN model works better at low concentration. On the other hand the ENr shows that it works better at medium and high concentrations, which is followed closely by GLNn model. If we divide the concentrations into two categories, where high concentration means that the nominal concentration is at least 3~pM and low concentration means that the nominal concentration is at most 1~pM, the GLNn model has the highest $R^2$ (the data is not shown here). It means, in general and at high concentrations, the GLNn offers a better fit than other models. As in Table~\ref{Tab5}, Table~\ref{Tab6} shows that all models have similar performance, although the GLNn model has the highest $R^2$ of nominal concentration against observed log-fold-change. \begin{table}[!t] \begin{center} \caption{The $R^2$ observed log-fold-change against nominal log-fold-changes for the spike-in genes.}\label{Tab6} \begin{tabular}{lcc} \hline Model & Obs-intended-fc.$R^2$ & Obs-(low)-int-fc.$R^2$ \\ \hline ENr & 0.976 & 0.989 \\ ENn & 0.974 & 0.990 \\ ENm & 0.972 & 0.985 \\ ENmc & 0.972 & 0.985 \\ EG & 0.972 & 0.986 \\ GN & 0.970 & 0.987 \\ ELNn & 0.974 & 0.990 \\ ELNm & 0.972 & 0.985 \\ ELNp & 0.973 & 0.990 \\ GLNn & \textbf{0.978} & \textbf{0.991} \\ GLNm & 0.971 & 0.984 \\ \hline \end{tabular} % \end{center} %\end{table} % %\begin{table}[!h] % \begin{center} \caption{AUC value for each model.}\label{Tab7} \begin{tabular}{lccccc} \hline Model & Low & Medium & High & Average & All \\ & concentration AUC & concentration AUC & concentration AUC & AUC & \\ \hline ENr & 0.450 & 0.987 & \textbf{0.785} & 0.585 & 0.886 \\ ENn & 0.518 & 0.987 & 0.764 & 0.631 & 0.899 \\ ENm & 0.573 & 0.987 & 0.741 & 0.667 & 0.911 \\ ENmc & 0.573 & 0.987 & 0.741 & 0.667 & 0.911 \\ EG & 0.567 & 0.987 & 0.746 & 0.664 & 0.910 \\ GN & 0.552 & 0.987 & 0.723 & 0.651 & 0.904 \\ ELNn & 0.524 & 0.987 & 0.763 & 0.635 & 0.900 \\ ELNm & 0.574 & 0.987 & 0.741 & 0.668 & 0.912 \\ ELNp & 0.534 & 0.987 & 0.761 & 0.642 & 0.902 \\ GLNn & 0.498 & \textbf{0.987} & 0.784 & 0.619 & 0.896 \\ GLNm & \textbf{0.579} & 0.987 & 0.730 & \textbf{0.671} & \textbf{0.913} \\ \hline \end{tabular} \end{center} \end{table} Table~\ref{Tab7} provides the results from the computation of the AUC value. The table shows that all models have a better accuracy at medium concentrations than at low and high concentrations. The ENr performs very poor at the low concentrations, but the GLNm performs best. At high concentrations, the ENr performs the best and it is followed by the GLNn. But in general, the highest AUC is achieved by all the model with the MLE parameter estimation methods: the GLNm, ELNm, and ENm. The computation, which is based on all arrays, provides the results where all models have the AUC greater than 0.9. According to \cite{Zhu10}, the AUC between 0.9 and 1.0 is classified as excellent in measuring the accuracy. Therefore, based on Table~\ref{Tab7}, we can identify that there are some models excellency accurate in predicting the gene expression. \paragraph{\textit{Second part}} \label{sec4322} We do $N=100$ simulations to assess the performance of each model. The bias of the background correction is assessed by the MSE$_{bc}$, and the bias of the parameterization is assessed by the $L_1$ error. \begin{table}[!t] \begin{center} \caption{The simulation results on the spike-in data set.}\label{Tab8} \begin{tabular}{lcccccc} \hline Model & MSE$_{bc}$ & \multicolumn{4}{c}{$L_1$ error} \\ \cline{3-6} & & $\alpha$ & $\beta$ & $\mu$ & $\sigma$ \\ \hline ENr & 0.045 & & 0.664 & 46.58 & 11.44 \\ ENn & 0.049 & & 0.625 & 41.61 & 2.806 \\ ENm & 0.038 & & 0.610 & 58.92 & 2.040 \\ ENmc & 0.036 & & 0.610 & 62.77 & 2.039 \\ GN & \textbf{0.030} & \textbf{0.0003} & \textbf{0.007} & 0.013 & \textbf{0.015} \\ ELNn & 0.048 & & 0.009 & 0.0004 & 0.018 \\ ELNm & 0.039 & & 0.840 & 0.0004 & 0.018 \\ ELNp & 0.061 & & 0.472 & 0.0004 & 0.018 \\ GLNn & 0.216 & 0.052 & 0.055 & 0.0004 & 0.018 \\ GLNm & 84.37 & 38.86 & 0.851 & \textbf{0.0003} & 0.017 \\ \hline \end{tabular} \end{center} \end{table} From the simulation results in Table~\ref{Tab8} we can see that results for the EG model are not available, because the MBCB package did not work at the log transformation that we have chosen. The GN model performs best, by providing the smallest bias for the background correction and the parameters. A close performance is achieved by the ELN, particularly by the ELNn. The GLNn does not have an optimal performance on the MSE$_{bc}$, but we still can consider its performance good, concerning that the bias of the parameters are similar to other proposed models and GN. One of the proposed models, GLNm has the highest bias on the MSE$_{bc}$ and the parameter $\alpha$. In our view this happens because we use an approximation in estimating the true intensity value. The EN models (ENr, ENm, ENn and ENmc) have considerably better performance at MSE$_{bc}$, but are not good at the parametrization. The bias on the parametrization of the noise is higher than in other models. \subsection{The public data sets} \label{sec44} Based on the results from Section~\ref{sec43}, we compare the performance of these models on some public data sets. We would like to know how good these models are in real data samples. Here, we choose to use the formalin-fixed, paraffin-embedded (FFPE) data sets from \cite{Wal12}: the FFPE of tumors from colorectal cancer patients (GSE32651, 1003 samples), breast cancer metastases of the lymph node and autopsy tissues (GSE32490: GSE32489, 120 samples). Each sample has 24526 bead-type level probes. The links for the data set are \url{http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32651} and \url{http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32490}. Currently the FFPE archival samples are widely available in million and it is a great source of information in medical studies about some diseases, for example cancer. This data type is suffering from the RNA degradation, which leads to poor performance in array-based studies. However, the Illumina's DASL assays could provide high-quality data from this degraded RNA samples. Comparing the performance of these background correction models certainly would help researchers to choose the appropriate background correction for their data, particularly if their data is the FFPE type. The background correction for the FFPE data set is implemented in three steps: \begin{description} \item[step 1] Do the quality control (QC) to the raw FFPE data. In this paper, we used the \emph{ffpe} package in R (\cite{Wal14}). \item[step 2] Do the data transformation $\log_2(P_{ij} + 1 + \sqrt{P_{ij}^2 + 1}$ to the raw FFPE data after QC and estimate the background correction parameters based on it. The estimators of true intensity value and the background correction are based on the regular and negative control probe intensity data, respectively. \item[step 3] Compute the true intensity value (the adjusted intensity estimator) based on the BC parameters at step 2. \end{description} The results of our computation are in Tables~\ref{Tab9} and \ref{Tab10}. From these tables we can see that there are no EG and GN models. Neither of these models can work on these data sets. For some samples in the data set, both models fail to compute the parameters which has the consequence that the true intensity value cannot be provided. We decided to remove the EG and GN models from further comparisons in both FFPE data sets. Here we provide the results of the rest of the models only. \begin{table}[!t] \begin{center} \caption{Simulation results on the GSE32651 data set.} \label{Tab9} \begin{tabular}{lrrrrrr} \hline Model & MSE$_{bc}$ & \multicolumn{4}{c}{$L_1$ error} \\ \cline{3-6} & & $\alpha$ & $\beta$ & $\mu$ & $\sigma$ \\ \hline ENr & 0.058 & & 0.657 & 297.67 & 22.61 \\ ENn & 0.281 & & 0.572 & 1.984 & 2.627 \\ ENm & 0.086 & & 0.591 & 28.05 & 1.715 \\ ENmc & \textbf{0.036} & & 0.610 & 62.77 & 2.039 \\ ELNn & 0.275 & & \textbf{0.025} & 0.001 & 0.018 \\ ELNm & 0.059 & & 0.826 & \textbf{0.0004} & 0.018 \\ ELNp & 0.672 & & 0.487 & 0.001 & 0.018 \\ GLNn & 0.838 & \textbf{0.340} & 0.527 & 0.001 & 0.018 \\ GLNm & 84.15 & 71.87 & 0.887 & 0.0005 & \textbf{0.018} \\ \hline \end{tabular} % \end{center} %\end{table} % %\begin{table}[!h] % \begin{center} \caption{Simulation results on the GSE32489 data set.} \label{Tab10} \begin{tabular}{lrrrrrr} \hline Model & MSE$_{bc}$ & \multicolumn{4}{c}{$L_1$ error} \\ \cline{3-6} & & $\alpha$ & $\beta$ & $\mu$ & $\sigma$ \\ \hline ENr & \textbf{0.093} & & 0.665 & 67.17 & 9.511 \\ ENn & 0.863 & & 0.509 & 0.936 & 2.039 \\ ENm & 0.182 & & 0.558 & 14.20 & 1.712 \\ ENmc & 0.184 & & 0.556 & 14.18 & 1.049 \\ ELNn & 1.055 & & 0.857 & 0.002 & 0.018 \\ ELNm & 0.116 & & 0.781 & 0.001 & 0.018 \\ ELNp & 1.247 & & 0.461 & 0.002 & 0.018 \\ GLNn & 1.348 & \textbf{0.332} & \textbf{0.497} & 0.002 & \textbf{0.018} \\ GLNm & 164.98 & 22.24 & 0.805 & \textbf{0.0004} & \textbf{0.018} \\ \hline \end{tabular} \end{center} \end{table} Tables~\ref{Tab9} and \ref{Tab10} consistently show that the bias of the parameters of noise in the EN models are higher than the proposed models. For the parameter $\beta$, the ELNn has the smallest bias and it is followed by the ELNp and the GLNn. With regard to the bias of the background correction, the EN models show the smallest bias in both of FFPE data sets. \section{Conclusions and indication of future work}\label{sec5} We have studied additive models of background correction for BeadArrays and proposed some new models where the true intensity is assumed to have exponential or gamma distribution and the noise is lognormally distributed. We have derived the estimator of the true intensity value of the proposed models. Further, we compared the performance of all models, based on the benchmarking and public data sets. In the benchmarking data set we adopted the criteria from the Affycomp \citep{Cop04} and for the simulation study we used the criteria which have been used in \cite{Xie09}, \cite{Che11} and \cite{Pla11,Pla12}. For the public data sets, we only used the criteria for the simulation study. We have seen in Sections~\ref{sec4321} and \ref{sec4322} that EN, EG, GN and GLN perform rather similar. However, the GLNn model has provided the highest reproducibility in comparison to other models. From the Affycomp criteria we can provide the following points: \begin{enumerate} \item The ENr and GLNn provide the lowest variation between replicates and all models using the MLE estimation method have a higher variation than others. \item The GLNn model has the highest signal detect $R^2$ in general and in high concentration. This means the GLNn model is the best fitted for the gene expression. \item The GLNn model, based on the MvA plot, produces the least number of genes which should not be expressed but are nevertheless expressed. On the other hand, the GN model provides the largest number of such genes. \item All models with the MLE estimation method have a higher average AUC value, which means that they provide a better accuracy in predicting the gene expression. \item The ENr and GLNn have the lowest IQR of log fold-change between replicates. \item Points 1 and 2 show that the GLNn and ENr are more accurate and precise in modelling the gene expression and points 3 and 5 show that the specificity and sensitivity of the GLNn and ENr model are better than others. \end{enumerate} In the simulation study, the best performance in estimating the signal by measuring its background correction and parametrization errors is achieved by the GN model. It is followed by our proposed ELN models. It has been shown that the GLNn does not perform optimally at the MSE$_{bc}$ criterion, but for the parametrization this model still can be considered good. In the FFPE public data set, the GN and EG models cannot be implemented. This is in strong contrast with the fact that in the simulation study of benchmarking data set, the GN model has the best performance. The EN models show the highest bias in the parametrization in both public data sets and the lowest bias in the background correction. Our proposed models, except the GLNm, show the lowest bias in the parametrization in both data sets and a moderate bias in the background correction. Based on the results from the benchmarking data and the public data sets, we would suggest researchers the following: \begin{enumerate} \item if the GN model works properly at the data set at hand (i.e. the estimated signals in all arrays can be computed by this model and the simulation criteria for this data at this model are low) then use the GN model to correct the background. \item if the GN model fails, then use our proposed models, particularly the GLNn model. The reason for not choosing the ELN models is that the value of the parameter $\alpha$ from the benchmarking data set is less than 1, around 0.2. Therefore, the gamma model is more appropriate to model the true intensity distribution than the exponential one. We believe that the right approximative computation of the GLN models will lead to a better performance than the current approximation. The ELN models perform better than the original EN models, due to the fact that not only the regular probes, but also the control probes are skew-distributed \citep{Che11}. Therefore, these models could be the second choice after the GLN, when the GN model does not work. \item With regard to the computation time, at the benchmarking data set the EN models are working faster than the others. They are followed by the ELNp, ELNn, and EG. The GLNn and the ENmc are the third fastest, then come the GN and the ELNm, which are followed by the GLNm as the slowest one. \end{enumerate} One of the purposes of using microarray technology is finding the genes which are expressed differentially due to some disease or condition. Therefore, it is important to investigate the effect of bias of the background correction and the parametrization toward the differentially expressed genes. This will be our future work. \subsection*{Acknowledgements} This paper is part of the author's Ph.D dissertation written under the direction of Prof.\ Istv\'an Berkes. We would like to thank Herwig Friedl, Ernst Stadlober, Levi Waldron, and an anonymous referee for their comments leading to a substantial improvement of the paper. We thank Florian Endel for discussion and suggestion how to make the R programming more efficient. We thank Yevgeniy Grigoryev and Ken Laing for their generous permission to use their pictures in the Supplementary material at \url{http://rfajriyah.staff.uii.ac.id/category/supplementary/}. Financial support from the Austrian Science Fund (FWF), Project P24302-N18 is gratefully acknowledged. \bibliography{emaRef} \end{document}