Correcting for noise in RT-qPCR gene expression data

Question

I have a training set of RT-qPCR gene expression data (not run in triplicate) for a batch of samples with two phenotypes $A$ and $B$ on which I've trained a "logistic regression classifier".

I also have another smaller set of samples which have been run twice as technical replicates on the RT-qPCR machine in order to study the effect of qPCR noise on my classifier.

I have a particular cut-off probability $P$ for my logistic regression model, such that we must be at least $P$ confident to classify a test sample as $A$ (similarly for $B$), otherwise no classification is returned.

I would like to incorporate the effect of RT-qPCR noise into my logistic classifier in order to return a more accurate classification probability for a particular test sample.

My model for the RT-qPCR noise for a particular sample $i$ is given by

$$X_i\sim N(\mu_i, \sigma^2I).$$

The logit is then given by

$$t_i = \beta_0 + \beta^TX_i,$$

with distribution

$$t_i\sim N(\beta_0+\beta^T\mu_i, \sigma_t^2),$$

where $\sigma_t^2=\|\beta\|_2^2\sigma^2$.

Since $\mu_i$ is a nuisance parameter, we eliminate it by taking the difference of the two replicates for a particular sample

$$Z_i=t_i^{(1)} - t_i^{(2)} \sim N(0, 2\sigma_t^2).$$

We can then form the MLE of $\sigma_t^2$ by computing

$$\hat{\sigma}_t^2=\frac{1}{2N}\sum_{i=1}^NZ_i^2.$$

Therefore, for a particular test sample $x_i$, we compute the expected value of the logit-normal distribution with location $t_i$ and squared scale $\hat{\sigma}_t^2$

$$p_i=\int_0^1x\frac{1}{\sqrt{2\pi\hat{\sigma}_t^2}}\frac{1}{x(1-x)}e^{-\frac{(\text{logit}(x)-t_i)^2}{2\hat{\sigma}_t^2}}dx.$$

The resulting $p_i$ is then our corrected probability of phenotype $A$.

Is this a sensible way to achieve my stated goal?

Here is an example plot of the necessary correction when the squared scale is $\sigma^2=2$,

• The black line corresponds to a qPCR assay with no noise.

• The violet line indicates the corrected probability.

For example, if the regression model says that a sample has a 75% chance of being phenotype $A$, then in the corrected model, this actually corresponds to only having about a 65% chance of being phenotype $A$.

Answer 1

You're actually not correct about the functional form of the logit. It's nonlinear, but you've written a linear function.

That matters here because once you understand the functional form of logistic regression, you can write it in a way that relates it to your specific generative model for your data. That can tell you what features you ought to use if your model assumptions are correct. You can find the details in PRML section 4.2 (pdf). The essence is:

• If your class labels are generated from a mixture of two classes with prior probabilities $\pi_1$ and $\pi_2$
• If your covariates are generated from exponential-family distributions with densities proportional to $\exp(\beta_1^T\Phi(x))$ and $\exp(\beta_2^T\Phi(x))$

then the posterior probability that an observation $x$ is in class 1 is

$$\frac{\pi_1\exp(\beta_1^T\Phi(x)) }{\pi_1\exp(\beta_1^T\Phi(x)) + \pi_2\exp(\beta_2^T\Phi(x))$}$$. If you plug in known class labels and mangle this for a while with algebra, you end up with a logistic regression using features $\Phi(x)$. The intercept term of the logistic regression turns into a function of $\pi_1/\pi_2$.

If your data are Gaussian, then $\Phi(x)$ would generally include $x$ and $x^2$ and all the cross-terms if $x$ is multivariate. If there are special assumptions about equal covariance between classes, then the squared terms and cross terms may not be needed.