# R model.matrix and makeContrast. Understanding model and possible contrast

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.

colnames(design)
[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

c(0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)


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

c(0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
c(0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0)


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:

c(0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,-1,-1)


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

c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1)


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.

• 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? – rvl Aug 17 '14 at 23:44
• The data points are really vectors with sequence counts from 25.000 genes. The edgeR package was designed for this particular case. – TIM Aug 18 '14 at 9:06