# The use of median polish for feature selection

In a paper I was reading recently I came across the following bit in their data analysis section:

The data table was then split into tissues and cell lines, and the two subtables were separately median polished (the rows and columns were iteratively adjusted to have median 0) before being rejoined into a single table. We finally then selected for the subset of genes whose expression varied by at least 4-fold from the median in this sample set in at least three of the samples tested

I have to say I don't really follow the reasoning here. I was wondering if you could help me answer the following two questions:

1. Why is it desirably/helpful to adjust the median in the datasets? Why should it be done separately for different type of samples?

2. How is this not modifying the experimental data? Is this a known way of picking a number of genes/variables from a large set of data, or is it rather adhoc?

Thanks,

• Can you pls elaborate on what kind of data are you/they looking at? I think judging by what you have quoted -- to me -- the method seems very ad hoc. – suncoolsu Mar 14 '11 at 9:11
• @suncoolsu: it's microarray data, if you are familiar with the concept. If not, I could perhaps summarize it as; which genes are expressed, to what extent in the studied samples. Here's a better explanation: en.wikipedia.org/wiki/Gene_expression_profiling – posdef Mar 14 '11 at 9:50
• @suncoolsu Almost definitely Gene Expression Analysis data. – kriegar Mar 14 '11 at 9:51
• Ok - I wasn't sure, next-gen sequencing is also getting popular. – suncoolsu Mar 14 '11 at 10:38

Tukey Median Polish, algorithm is used in the RMA normalization of microarrays. As you may be aware, microarray data is quite noisy, therefore they need a more robust way of estimating the probe intensities taking into account of observations for all the probes and microarrays. This is a typical model used for normalizing intensities of probes across arrays.

$$Y_{ij} = \mu_{i} + \alpha_{j} + \epsilon_{ij}$$ $$i=1,\ldots,I \qquad j=1,\ldots, J$$

Where $Y_{ij}$ is the $log$ transformed PM intensity for the $i^{th}$probe on the $j^{th}$ array. $\epsilon_{ij}$ are background noise and they can be assumed to correspond to noise in normal linear regression. However, a distributive assumption on $\epsilon$ may be restrictive, therefore we use Tukey Median Polish to get the estimates for $\hat{\mu_i}$ and $\hat{\alpha_j}$. This is a robust way of normalizing across arrays, as we want to separate signal, the intensity due to probe, from the array effect, $\alpha$. We can obtain the signal by normalizing for the array effect $\hat{\alpha_j}$ for all the arrays. Thus, we are only left with the probe effects plus some random noise.

The link that I have quoted before uses Tukey median polish to estimate the differentially expressed genes or "interesting" genes by ranking by the probe effect. However, the paper is pretty old, and probably at that time people were still trying to figure out how to analyze microarray data. Efron's non-parametric empirical Bayesian methods paper came in 2001, but probably may not have been widely used.

However, now we understand a lot about microarrays (statistically) and are pretty sure about their statistical analysis.

Microarray data is pretty noisy and RMA (which uses Median Polish) is one of the most popular normalization methods, may be because of its simplicity. Other popular and sophisticated methods are: GCRMA, VSN. It is important to normalize as the interest is probe effect and not array effect.

As you expect, the analysis could have benefited by some methods which take advantage of information borrowing across genes. These may include, Bayesian or empirical Bayesian methods. May be the paper that you are reading is old and these techniques weren't out until then.

Regarding your second point, yes they are probably modifying the experimental data. But, I think, this modification is for a better cause, hence justifiable. The reason being

a) Microarray data are pretty noisy. When the interest is probe effect, normalizing data by RMA, GCRMA, VSN, etc. is necessary and may be taking advantage of any special structure in the data is good. But I would avoid doing the second part. This is mainly because if we don't know the structure in advance, it is better not impose a lot of assumptions.

b) Most of the microarray experiments are exploratory in their nature, that is, the researchers are trying to narrow down to a few set of "interesting" genes for further analysis or experiments. If these genes have a strong signal, modifications like normalizations should not (substantially) effect the final results.

Therefore, the modifications may be justified. But I must remark, overdoing the normalizations may lead to wrong results.

• +1 This is a much better answer than my attempt. Thanks. – kriegar Mar 14 '11 at 11:53
• @posdef. I am wondering if there was any statistician involved in the statistical analysis of the paper. – suncoolsu Mar 14 '11 at 12:08
• thanks for your thorough reply. I think the fact that this is a pre-processing step is not well explained (or just assumed to be well known) in the paper. Speaking of which, the paper is published in 2000 (in Nature) so I presume they had at least some statistician look at their methods, if not involved in writing. But of course I can only speculate.. :) – posdef Mar 14 '11 at 12:34
• @posdef. Ok- cool answers a lot of questions. 2000 was the time when people were still figuring out how to analyze microarray data. FDR wasn't fancy back then :-) – suncoolsu Mar 14 '11 at 12:56

You may find some clues in pages 4 and 5 of this

It is a method of calculating residuals for the model $$y_{i,j} = m + a_i + b_j + e_{i,j}$$ by calculating values for $m$, $a_i$ and $b_j$ so that if the $e_{i,j}$ are tabulated, the median of each row and of each column is 0.

The more conventional approach amounts to calculating values for $m$, $a_i$ and $b_j$ so that the mean (or sum) of each row and each column of residuals is 0.

The advantage of using the median is robustness to a small number of outliers; the disadvantage is that you are throwing away potentially useful information if there are no outliers.

• thanks for the answer, and the reference link. However I can't see how this model applies to the problem at hand. given that the data is comparative expression values (read: abundance) how can one define $a_i$, $b_j$ and $e_{i,j}$?? – posdef Mar 14 '11 at 8:55
• If instead you take an abundance model like $n_{i,j} = n_i \,q_j + e_{i,j}$ or one like $log(n_{i,j}) = log(n) + log(p_i) + log(q_j) + e_{i,j}$ then you can do essentially the same thing, making the the median of each row and of each column of the residuals table equal to 0. – Henry Mar 14 '11 at 12:10
• @Henry What information is "thrown out" with median polish when there are no "outliers" (and what exactly do you mean by "outlier" anyway)? After all, you can reconstruct the data exactly by means of the grand median, the row and column medians, and the residuals, all of which constitute the output of median polish. If you mean the residuals are discarded, then in what sense is "mean polish" (equivalent to OLS) any different in this regard? – whuber Mar 14 '11 at 15:21
• @whuber: The residuals are kept in both cases. The mean polish takes account of how far away observations are from the centre (in a sense, it balances the weights of the residuals) while the median polish only looks at whether they are above or below the centre (in a sense, it balances the numbers of residuals). So the weight information is unused when using the median as the centre; this can be good when some of the substantial weights/residuals are so dubious that the result for the centre cannot be trusted, but involves not using information if not. – Henry Mar 14 '11 at 15:56
• @Henry If you can recover all the original data from the polish, then how is "information" not "used"? BTW, median polish does not behave as you seem to describe: its residuals are the differences in values, not in ranks, of the data. – whuber Mar 14 '11 at 15:59

Looks like you are reading a paper that has some gene differential expression analysis. Having done some research involving microarray chips, I can share what little knowledge (hopefully correct) I have about using median polish.

Using median polish during the summarization step of microarray preprocessing is somewhat of a standard way to rid data of outliers with perfect match probe only chips (at least for RMA).

Median polish for microarray data is where you have the chip effect and probe effect as your rows and columns:

for each probe set (composed of n number of the same probe) on x chips:

chip1    chip2    chip3   ...  chipx
probe1      iv       iv       iv   ...     iv
probe2      iv       iv       iv   ...     iv
probe3      iv       iv       iv   ...     iv
...
proben      iv       iv       iv   ...     iv

where iv are intensity values

Because of the variability of the probe intensities, almost all analysis of microarray data is preprocessed using some sort of background correction and normalization before summarization.

here are some links to the bioC mailing list threads that talk about using median polish vs other methods:

https://stat.ethz.ch/pipermail/bioconductor/2004-May/004752.html

https://stat.ethz.ch/pipermail/bioconductor/2004-May/004734.html

Data from tissues and cell lines are usually analysed separately because when cells are cultured their expression profiles change dramatically from collected tissue samples. Without having more of the paper it is difficult to say whether or not processing the samples separately was appropriate.

Normalization, background correction, and summarization steps in the analysis pipeline are all modifications of experimental data, but in it's unprocessed state, the chip effects, batch effects, processing effects would overshadow any signal for analysis. These microarray experiments generate lists of genes that are candidates for follow up experiments (qPCR, etc) to confirm the results.

As far as being ad hoc, ask 5 people what fold difference is required for a gene to be considered differentially expressed and you will come up with at least 3 different answers.

• Thanks for the updates on your answer, I think I am starting to get an idea now. So if I understand correctly, the median polishing is used to assess the technical variability with regards to the probe and the chip? ... before the experiment is summed up to 1 matrix holding expression values for genes under different conditions? – posdef Mar 14 '11 at 11:06
• @posdef from my understanding yes. For each probeset on a chip (probes of the same sequence) there are probes scattered throughout. plmimagegallery.bmbolstad.com for some pseudo images of chips. In addition to the variability within a single chip, there is variability between chips. Because of the technical variability, algorithms are run on the raw intensity values to obtain a single "expression value" for the probeset. The matrix of these values are then fit to determine if the genes are differentially expressed under different conditions. – kriegar Mar 14 '11 at 11:46