MicroArray, tesing if a sample is the same with high variance data.

Question

I'll explain the problem but what I am looking for is a few suggested methods to approach this problem. You don't need to know what a microarray but if you are interested look here link text The info below is simplified, not addressed to a micro-biologist (I am not one)

Data:

There are about 250,000 probes on an array that each measure a part of a dna sample. An array is specific to an individual(germ, person....) If the array is designed for my dna and we put yours on it some probes will match (we are both human I think) and some will not. If they match the measurement is high(the probe glows) as the dna at a probe matched less it glows less. There is a very high variance at each probe across the same dna. So it is very difficult to know if a individual sample matches (same dna) at the probe level.

If the data was nice we could just do a simple Hypothesis test to determine if the sample probe was the same as the dna the array was made for. Because of the high variance this will not work.

Question:

We are looking for a way to consider the measurements at other probes to improve the test at a probe. We kinda know this should work (based on the DNA) but we are not sure how to choose the probes that will best assist for each of the set of 250,000 probes or what statistical method would be best to use the information about the other probes to test an individual probe.

Different example, maybe easier to understand:

A person has 6 characteristics, lets say weight, height, age, income, date of birth (db), education, race. I measure characteristics for myself but my measurement tools are really bad so I make several measurements but I have a very large variance in my measurements. Now I measure another person and want to know if they have a different income than myself. But my measurement tools are so bad almost any other person’s measured income would be within my confidence interval. (I can only measure the other person once, It kills them ☺ ) But I know that other characteristics should be able to help me determine if the income is the same. I don’t know which ones (weight, height, age, db, education, race). I need a method to choose the most useful/helpful characteristics. Now lets say I know that height and age can help determine if our income are the same or different. What would be the statistical method used?

I have what might be called a training set. That is in know what is and is not the same.

@sheldon-cooper answer, "Note: neither of this will give a "significance score"", A confidence value is good enough. I don't need a statistical "Test"

As an MD I hope that this is an acceptable question for MathOverFlow.net. and I hope you can understand the question as abridged as it is. If some would like to help me rewrite or has suggestion on how I can rewrite the question without all the bio background and as a more specific stats question, let me know.

Answer 1

In general, the approach of using additional measurements of other values (not the one you are interested in directly) was useful in many problems in the past, so it sounds like a good idea in your case as well. Here are a few things to look at:

Binary classification

If the task is to determine whether a given probe is a 'match', you can just pose it as a binary classification problem. The input features to a classifier will be ALL the probes in the array. The classifier will determine automatically which of these are useful. If the variance at the target probe is actually small, it may decide to just look at that probe only. If other sites (probes) are helpful too, it may decide to look at them as well. Classifiers can automatically construct complex rules that take into account one or many probes, as needed.

This can be simple to try since there are packages readily available and there are not many parameters to tune. I'd try using a classifier called SVM, it works well in many cases, can deal with limited training data (looks like that's an issue for you). Again, there are SVM packages available that you can download and just run (libsvm is one example, but there are many others).

The disadvantage of this is that it may become resource-consuming if you want a classifier for each probe. Your input dimensionality is 250,000. A simple linear classifier needs roughly one number per dimension, so 250,000 numbers. If you need such a classifier for every probe, this will require 250,000^2 numbers total, which is a lot. There are ways to go around this, but that's one concern.

Collaborative filtering

Collaborative filtering basically allows you to predict the value of one item from the values of related items. Say you know that I liked two movies M1, M2, and disliked two others, M3 and M4. Using collaborative filtering, you can predict whether or not I will like a new movie M5. In your case, you are not doing a completely blind prediction for M5; you do have some (although very noisy) observation. So you would combine that observation with the prediction from other features.

The advantage of this is that there is only one model (as opposed to one classifier per probe). So this is easier to scale.

One model would work something like this. It would group probes that are either all "same" or all "different" most of the time. (These groups would be learned by the model from training data.) Then, given a new array, it would use these groups and data from other probes in each group to determine whether a group is, again, "same" or "different". This question had a link to what looked like an easy intro to collaborative filtering.

The main disadvantage is, I'm not sure you can find an implementation of the model you need for microarrays online (or anywhere). You'd probably need to derive something suitable for your data, which is not difficult, but requires knowing what to do.

Another disadvantage is, you'll probably need more training data for collaborative filtering to work.

Note: neither of this will give a "significance score" (as in t-test, for example). You'll get a confidence value, but you won't be able to say that the findings are statistically significant with a certain confidence.

Answer 2

I'll try to approach the question from a very naiive point of view (you've already got two high-tech answers and, given my lack of special knowledge, it would be foolish for me to try to compete with them using some other high-tech). So, I assume that the model is the following. You have $m=36$ samples with $n=250,000$ positions and one $n$-position array against which you want to match. You believe (this is the most shaky part because I'm still at loss about what exactly that belief is) that the matching in the first position occurs if and only if it also occurs in some other $k$ unknown positions (or, at least, that matching in those $k$ positions makes the matching in the first position more likely and vice versa). You have the measuring tool that outputs some numerical value $X_i$' for the position $i$ that has mean $1$ when the match occurs and mean $0$ when it doesn't but you also have a lot of random noise in it, the noises for different positions being independent and having pretty much the same variance of size comparable to $1$. Your question is which $k$ extra positions to use to support the match in the first position and how to do it.

It seems to me that it would be reasonable to try the weighted sum rule $X_1+w\sum_{j=1}^k X_{i_j}>T$ for making the decision (with some $w\in(0,1)$). If you agree to that, the choice of $k$ unknown positions $i_j$ is very simple: just do all your 36 measurements and in each position $i$ look at the difference of the average value of $X_i$ over the good samples (those that match in the first position) and bad ones (those that do not). Choose the $k$ positions with the largest values of this difference. Now, play with $w$ and $T$ to make as few errors on your 36 samples as possible. Then just cross your fingers and hope for the best.

This is more or less equivalent to what an intelligent binary classifier would do except you can have more control over what is going on and can incorporate some a priori knowledge into this scheme (for instance, you may know what $k$ is, which positions are certainly irrelevant, and so on). You may want to assign different weights to different supporting probes but this gives you too much freedom and may be more harmful than useful on a small sample size. The crucial decision is what $k$ to use. I surmise you have some idea coming from your knowledge of the DNA structure but if not, keep in mind that to make it useful, you should have $(1+kw)^2\gg 1+kw^2$ (the mean difference should be noticeably larger than the standard deviation).

Will that work? I still suspect that 36 is far too small a number to do the correct choice of a few supporting positions out of 250,000 based solely on the differences in averages so you should use your knowledge of biology to reduce the number of possible supporting positions as much as you can (and then you choose $k$ largest differences out of those).

Answer 3

What you have is a classic case of a high dimensional, low signal to noise ratio signal. There are a lot of ways to proceed, but ultimately you will want to know about three different effects:

1. Bayesian estimation.
2. Dimension reduction.
3. Superefficient estimation.

These three ideas are frequently conflated by people with an informal understanding of them. So let's clear that up right away.

1. The Bayesian estimation can always be applied. To estimate the probability of success p of independent flips of a possibly unfair coin you can always have a belief that the coin is fair or loaded, and this can lead to a different posterior distribution. You cannot reduce the dimension of the parameter space - unless you decide that you have what passes for 'revealed wisdom' as to that probability. It is also impossible to perform superefficient estimation in this case because the parameter space does not admit it.

Note that in high dimensional cases, the Bayesian prior may reflect beliefs about noise, parameter spaces, and high dimensional estimation in general, as opposed to the traditional understanding that the Bayesian prior reflects belief about the parameter value. The farther you go into the high dimensional / low signal to noise regime, the less you expect to find anything passing for "expert opinion" that relates to the parameter and the more you expect to find all sorts of "expert opinion" about how hard it is to estimate the parameter. This is an important point: You can throw out the bathwater and keep the baby by knowing a lot about babies. But, and this is the one that most statisticians and engineers don't keep in the front window: You can ALSO throw out the bathwater and keep the baby by knowing a lot about bathwater. One of the strategies for exploiting the high availability of noise is that you can study the noise pretty easily. Another important point is that your willingness to believe that you do not have a good idea how to parameterize the probability density means that you will be considering things like Jeffreys' rule for a prior - the unique prior which provides expectations invariant under change of parameterization. Jeffreys' rule is very thin soup as Bayesian priors go - but in the high dimensional low signal to noise ratio it can be significant. It represents another very important principle of this sort of work: "Don't know? Don't Care". There are a lot of things you are not going to know in your situation; you should arrange as much as possible to line up what you are not going to know with what you will not care about. Don't know a good parameterization? Then appeal to a prior (e.g. Jeffreys' rule) which does not depend on the choice of parameterization.

As an example, in the parameterization of poles and zeros, a finite dimensional linear time invariant system has the Jeffreys' prior given by the hyperbolic transfinite diameter of the set of poles and zeros. It turns out that poles and zeros is a provably exponentially bad general parameterization for the finite dimensional linear time invariant system in high dimensions. But you can use the Jeffreys' prior and know that the expectations you would compute with this bad parameterization will be the same (at least on the blackboard) as if you had computed them in some unknown 'good' parameterization.

1. Dimension reduction. A high dimensional model is by definition capable of dimension reduction. One can, by various means, map the high dimensional parameter space to a lower dimensional parameter space. For example one can locally project onto the eigenspaces of the Fisher information matrix which have large eigenvalues. A lot of naive information theory is along the lines that "fewer parameters is better". It turns out that is false in general, but sometimes true. There are many "Information criteria" which seek to choose the dimension of the parameter space based on how much the likelihood function increases with the parameter. Be skeptical of these. In actual fact, every finite dimensional parameterization has some sort of bias. It is normally difficult to reduce the dimension intelligently unless you have some sort of side information. For example, if you have a translation invariant system, then dimension reduction becomes very feasible, although still very technical. Dimension reduction always interacts with choice of parameterization. Practical model reduction is normally an ad hoc approach. You want to cultivate a deep understanding of Bayesian and superefficient estimation before you choose a form of dimension reduction. However, lots of people skip that step. There is a wide world of ad hoc dimension reduction. Typically, one gets these by adding a small ad hoc penalty to the likelihood function. For example an L1 norm penalty tends to produce parameter estimates with many components equal to zero - because the L1 ball is a cross-polytope (e.g. octahedron) and the vertices of the cross polytope are the standard basis vectors. This is called compressed sensing, and it is a very active area. Needless to say, the estimates you get from this sort of approach depend critically on the coordinate system - it is a good idea to think through the coordinate system BEFORE applying compressed sensing. We will see an echo of this idea a bit later in superefficient estimation. However it is important to avoid confounding dimension reduction with superefficient estimation; so to prove that I give the example of a parameter with two real unconstrained components. You can do compressed sensing in two dimensions. You cannot do superefficient estimation in the plane. Another very important aspect here is "Who's Asking?". If you are estimating a high dimensional parameter, but the only use that will be made of that parameter is to examine the first component? Stop doing that, OK? It is very worthwhile to parameterize the DECISIONS that will be made from the estimation and then look at the preimage of the decisions in the parameter space. Essentially you want to compose the likelihood function (which is usually how you certify your belief about the observations that have been parameterized) with the decision function, and then maximize that (in the presence of your dimension reduction). You can consider the decision function another piece of the baby/bathwater separation, or you can also look at maximal application of don't care to relieve the pressure on what you have to know.

   1. Superefficient estimation. In the 1950s, following the discovery of the information inequality (formerly known as Cramer-Rao bound), statisticians thought that the best estimates would be unbiased estimates. To their surprise, embarassment, and dismay, they were able to show that in three dimensions and higher, this was not true. The basic example is James-Stein estimation. Probably because the word "biased" sounds bad, people adhere much much longer to "unbiased" estimation than they should have. Plus, the other main flavor of biased estimation, Bayesian estimation, was embroiled in an internecine philosophical war among the statisticians. It is appropriate for academic statisticians to attend to the foundations of their subject, and fiercely adhere to there convictions. That's what the 'academy' is for. However, you are faced with the high dimensional low signal to noise ratio case, which dissolves all philosophy. You are in a situation where unbiased estimation will do badly, and an expert with an informative Bayesian prior is not to be found. But you can prove that superefficient estimation will be good in some important senses (and certainly better than unbiased estimation). So you will do it. In the end, the easiest thing to do is to transform your parameter so that the local Fisher information is the identity (we call this the Fisher Unitary parameterization) and then apply a mild modification of James-Stein estimation (which you can find on Wikipedia). Yes, there are other things that you can do, but this one is as good as any if you can do it. There are some more ad hoc methods, mostly called "shrinkage" estimation. There is a large literature, and things like Beran's REACT theory are worth using, as well as the big back yard of wavelet shrinkage. Don't get too excited about the wavelet in wavelet shrinkage - it's just another coordinate transformation in this business (sorry, Harmonic Analysts). None of these methods can beat a Fisher unitary coordinate transformation IF you can find one. Oddly enough, a lot of the work that goes into having a Fisher unitary coordinate transformation is choosing a parameterization which affords you one. The global geometry of the parameter space and superefficient estimation interact very strongly. Go read Komaki's paper:

KOMAKI, F. (2006) Shrinkage priors for Bayesian prediction. to appear in Annals of Statistics. (http://arxiv.org/pdf/math/0607021)

which makes that clear. It is thought that if the Brownian diffusion (with Fisher information as Riemannian pseudo-metric) on your parameter space is transient, then you can do superefficient estimation, and not if it is recurrent. This corresponds to the heat kernel on the parameter space, etc. This is very well known in differential geometry to be global information about the manifold. Note that the Bayesian prior and most forms of model reduction are entirely local information. This is a huge and purely mathematical distinction. Do not be confused by the fact that you can show that after you do shrinkage, that there EXISTS a Bayesian prior which agrees with the shrinkage estimate; that just means that shrinkage estimates have some properties in common with Bayesian estimates (from the point of view of decisions and loss functions, etc.) but it does not give you that prior until you construct the shrinkage estimate. Superefficient estimation is GLOBAL.

One other excruciatingly scary aspect of superefficient estimation is that it is extremely "non-unique". There is an inexplicable arbitrary choice; typically of a "point of superefficiency", but it is more like a measure of superefficiency that started out as an atomic measure. You have to choose this thing. There is no reason whatever for you to make one or another available choice. You might as well derive your choice from your social security number. And the parameter estimate that you get, as well as the performance of that estimate, depends on that choice. This is a very important reason that statisticians hated this kind of estimation. But they also hate the situation where you have tons of variance and grams of data, and that is where you are. You can prove (e.g. Komaki's paper gives an example of such proof) that your estimation will be better if you make this choice, so you're going to do it. Just don't expect to ever understand much about that choice. Apply the dont' know/don't care postulate - you will not know, so you're better off not caring. The defense of your estimation is the theorem that proves it's better.

It should be very clear now that these three "nonclassical" effects in estimation theory are really distinct. I think most people don't really understand that. And to some extent it's easy to see why.

Suppose you have an overparameterized generalized linear model (GLM), so your Fisher information is singular, but you do something like iteratively reweighted least squares for the Fisher scoring (because that's what you do with a GLM) and it turns out that the software you use solves with say Householder QR with column pivoting. It's down there under the covers enough that many statisticians performing this estimation would not necessarily know that is what was happening. But because the QR with column pivoting regularizes the system, effectively it is estimating the parameter in a reduced dimension system, where the reduction was done in a Fisher unitary coordinate system (because of the R in the QR factorization). It's really hard for people to understand what they are really doing when they are not aware of the effect of each step. We used to call this "idiot regularization" but I think "sleepwalker" is more accurate than "idiot". But what if the software package used modified Cholesky to solve the system? Well that actually amounts to a form of shrinkage (again in a Fisher unitary coordinate system), it can also be considered a form ("maximum a posteriori") Bayesian prior.

So in order to sort out what these effects do, and what that means you should do, you need a reasonably deep understanding of the custody of your digits all the way from the data being observed through the final decision being taken (treatment, grant proposal, etc.).

At this point (if you're still with me) you might want to know why didn't I just write a one line recommendation of some method suitable for beginners.

Well there isn't one. If you really are in the high dimensional, low signal to noise case, then a high quality result is only available to someone who understands a pretty big part of estimation theory. This is because all three of these distinct effects can improve your result; but in each situation it is difficult to predict what combination will be the most successful. In fact the reason you cannot prescribe a method in advance is precisely because you are in the regime you are in. You can bail out and do something out of a canned program, but that has a good risk of leaving a good deal of the value on the table.