I am trying to analyse a gene expression dataset of about 240 samples (Illumina microarray).

There are multiple questions I need to ask (healthy vs diseased, effect of treatment in diseased, differences between responders and non responders) and I need to correct for multiple covariates including batch (samples were run in 5 batches).

My initial strategy was to select different samples for every question (for example when I analyse the effect of treatment I only select diseased samples) and always include batch and other covariates like sex and age in the model.

The problem is that for some questions I can only use 15 samples so I am not sure I am estimating the effect of batches correctly.

Would it be better to use all the samples to estimate the effect of batch and other covariates and then work on residuals to answer specific questions?

Are there other strategies I could use?

Any help would be greatly appreciated.

Here is a schema of the variables I have:

study_group: healthy, disease1, disease2 treatment (for disease 2 only): drug A, drug B response_status (for disease 2 only): RESPONDER, NON_RESPONDER visit_number: 1 (all study_groups), 2 (disease 1 and 2), 3 (disease 2 only), 4 (disease 2 only), 5 (disease 2 only) age, ethnicity, batch, sex as confounders

If for example I want to ask the question "does gene expression differ between responders and non responders at visit 1" should I select only samples for which I have response_status or should I use all samples and assign the category "none" to treatment and response_status for healthy and disease 1 samples?


1 Answer 1


If you estimate each contrast using only a subset of samples, you throw away information. E.g. if you use only 24 observations out of 240, then only 24 residuals are used to estimate the variance of error term in limma, and, all other things being equal, your inference is less precise.

I think the way to do it is to use all of the observations to estimate all of the contrasts. Remember that, if you drop a factor that explains some of the response variance from the model, the variance explained by that factor is added to the error term variance. As a result, p-values for the factor(s) remaining in the model can go up. Therefore, it's wrong to think that if factor A is of interest and factor B is not related to the asked question, then it's only A that has to be in the model.

If you have microarray data, then you should have continuous intensity measurements for each of 240 cases – that's your response. Use all of the observations, just make sure to set the contrasts right. E.g. if you want to test drug A vs B then you can't just include “drug” main effect into the model, because it's not defined for Disease1 and for Healthy subjects. Instead, you should create a new factor with the levels like “Disease2_drugA” and “Disease2_DrugB”, and the corresponding contrast will be “ Disease2_drugA vs Disease2_DrugB”. A similar problem is described here. In your case, you may even need to “merge” three or more factors in that manner.

  • $\begingroup$ The reason why I thought I had to subset is that I have no values for the variable response for healthy people for example because they have not been treated. If I use NA in R the corresponding samples get automatically eliminated. Should I use a fake variable like None for those samples? $\endgroup$
    – Rossella
    Jul 19, 2016 at 7:32
  • $\begingroup$ Could you post a an approximate snapshot of your dataset and specify what the response and predictors are. $\endgroup$
    – Nik Tuzov
    Jul 19, 2016 at 14:28
  • $\begingroup$ I have edited the post to include a description of the variables in the dataset. $\endgroup$
    – Rossella
    Jul 20, 2016 at 16:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.