I have measurements from 12 mice, grouped in two conditions. I each mouse I have measurements from 4 tissues. The design is not balanced, 5 mice in condition1 and 7 in condition2.

After reading the bioconductor edgeR manual I have set up the following model:

design <- model.matrix(~ Condition + Tissue + Condition:Mouse + Condition:Tissue)

I have then manually removed terms with no observations.

design <- design[, -which(colnames(design) %in% c("Condition1:Mouse6", "Condition1:Mouse7"))]

I can fit the model, but when I set up contrasts I get confused. I want to get the main effect of being in condition2, the main effect of each tissue and the interaction effect for each tissue.

 [1] "(Intercept)"  "Condition2"    "Tissue2"              
 [4] "Tissue3"  "Tissue4"   "Condition1:Mouse2"            
 [7] "Condition2:Mouse2"    "Condition1:Mouse3" "Condition2:Mouse3"         
 [10] "Condition1:Mouse4"    "Condition2:Mouse4" "Condition1:Mouse5"            
 [13] "Condition2:Mouse5"    "Condition2:Mouse6" "Condition2:Mouse7"         
 [16] "Condition2:Tissue2"   "Condition2:Tissue3"  "Condition2:Tissue4" 

To get the main effect of condition I use the contrast


Comparing e.g. Tissue 1 and 2 or 2 and 3 I do like this


But say I want to get the difference between condition 1 and 2 in tissue 1. Given how the model matrix is built tissue 1 is not present. Can I do a comparison like:


And how about an anova like test for changes across all the tissue:condition interaction terms. The contrast


Would compare only tissue2, 3 and 4, and with what?

Many questions here. Quick recap. 1. Is it OK to remove empty terms in an unbalanced design? 2. How do I find the effect of the level in the tissues category that is missing in the model matrix. 3. How do I find changes happening across all levels of the interaction terms (related to the missing tissue1 problem).

A link to a page with a thorough R model matrix/linear model guide would also be appreciated.

  • 1
    $\begingroup$ This seems like a tedious approach. Is there a reason you're not using a higher-level mixed-model-fitting function such as lmer in the lme4 package? $\endgroup$ – Russ Lenth Aug 17 '14 at 23:44
  • $\begingroup$ The data points are really vectors with sequence counts from 25.000 genes. The edgeR package was designed for this particular case. $\endgroup$ – TIM Aug 18 '14 at 9:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.