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

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.

• Do you have a model for how the results should look if you know how the sample DNA compares with the array's DNA? – Douglas Zare Feb 17 '10 at 3:52
• "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". What exactly do these sentences mean? Do you mean that even if there is a match, you may easily fail to detect it and come up with a wrong result on each particular probe or you mean something else? – fedja Feb 17 '10 at 4:02
• I have several samples that have been exposed to the dna that matches the array. But again the variance is large. In addition you can't say that with the right dna the probe should have a reading of say 100 or some specific measurement. – Lisa Feb 17 '10 at 4:04
• @Fedja, yes, Let me try an example, lets say I have 10 measurements (1-10) for the for the base DNA I have stats for each measurement location. Now I take an unknown sample and find the measurements. The variance in the base data is such that almost any sample at the individual measurement level will be in the 95% confidante int at the measurement point. I know that some of the other measurement points can help but I need a method to identify the right ones and a method to use them once identified. Maybe I will try to edit the original post with a no bio example. – Lisa Feb 17 '10 at 4:12
• All right. So, if we subtract $M_i$ and divide by $\sigma_i$, the question boils down to whether a given set of measurements $Y_i$ can come from the independent samples of the standard normal law. The difficulty is that we want to detect the (possibly present) non-zero means that are of size less than $1$, so each individual measurement doesn't tell us much. Also, we do not know if the means are positive or negative (Or do we? Is the mean measurement for a non-matching sample always below the one for a matching one? If so, this may help a lot) – fedja Feb 17 '10 at 5:04

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.

• I have considered this. Lots for me to learn. I have about 12 samples that match the array do derive the mean and sd at each of the probes. And I have a small training set of about 40 sample that don't and again there are about 250,000 probes. Do you have an classifier method/algorithm you would recommended I look at? – Lisa Feb 17 '10 at 5:33
• "collaborative filtering", How does a single model work, how do you choose the helper probes? I guess I should read about this more. That said I had assumed I would need a "test" for each probe. I can do the programing so that's not a problem. – Lisa Feb 17 '10 at 5:48
• The array is specific to a individual and those are the only arrays I have. I have 15 or so other individuals and they have been sampled 1-3 times each. (these are parasites and I have 15 different strains/hybridization, I have sampled each strain 1-3 times. The array is specific to one of the strains) – Lisa Feb 17 '10 at 5:54
• Is the original individual also a parasite, similar to the other 15 parasites? If yes, you could use other parasites as examples too. A training example is really a pair of microarrays with a label "same" or "different". You have 12 arrays for the main subject; you can build 66 positive pairs from them. But if you have a parasite you sampled 3 times, you can build 3 more positive pairs from it too. That can be useful if you think they are all similar. If they are very different (e.g. one is a human and others are worms), this will not work. – user3035 Feb 17 '10 at 6:02
• @sheldon-cooper, Determining if the sample is a match to the array (the same parasite) is a solved problem. I am interested in finding unique parasites. That is they are different in unique was at specific probes. So yes is a unknown parasite different at a specific probe? – Lisa Feb 17 '10 at 6:11

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).

• ok I like "simple". I might not (I am sure I don't) fully get how you suggest I choose the supporting positions. My primary question is why having large difference in say, position 3 is likely to help me make a decision about position 1. I guess my question is that this does not seem to use a "correlation" to choose the supporting positions. I have thought of calculating a correlation and then using this to apply an approach similar to yours. In fact I think I will start there. – Lisa Feb 17 '10 at 15:39
• Under my assumptions, the correlation with the random variable that takes two values $0$ and $1$ is just the difference of the means I mentioned times some factor almost independent of $i$, so what I said and what you are saying now are almost the same thing in this model. If 111000 is what you want to correlate with, my means test will not decide which of 111011 and 110100 should be used as a better support, but your correlation test will prefer the first. I leave it to you to decide if it is good or bad. – fedja Feb 17 '10 at 18:01

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.

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.

• Mathoverflow insists on changing the display of numbers in my sections - I actually have typed in 1. 2. and 3. I do not know why they display as 1. 1. and 1. in the final. – Andrew Mullhaupt Feb 19 '10 at 17:38