# [BioC] blocking in limma gives an estimability problem

Pedro López Romero plopez at cnic.es
Mon Jun 12 18:49:05 CEST 2006

```Dear List,

I have a problem of estimability using LIMMA that I do not know how to
solve. I think that I am using the right model, but I can not solve the
problem so I need someone to help me out to figure out what is going on.

I have an experiment with paired samples within two different strains. The
data is as follows:

Array		mouse		strain	treatment
1		1		wt		c
2		2		wt		c
3		3		wt		c
4		1		wt		t
5		2		wt		t
6		3		wt		t

7		4		ko		c
8		5		ko		c
9		6		ko		c
10		4		ko		t
11		5		ko		t
12		6		ko		t

It is exactly the same as the example 8.3 in the limma user guide, except
that I have two strains, WT and KO with the same paired structure.
I am using a linear model with the animal effect as a block effect to
consider the pairing of the data,

My code is the following:

block=factor(c("A1","A2","A3","A1","A2","A3",
"A4","A5","A6","A4","A5","A6"),
levels=c("A1","A2","A3","A4","A5","A6") )

contrasts(block)=contr.sum(6)

ttoEFF=factor(c("wtC","wtC","wtC","wtM","wtM","wtM",
"koC","koC","koC","koM","koM","koM"),
levels=c("wtC","wtM","koC","koM"))

contrasts(ttoEFF)=contr.sum(4)

design=model.matrix(~ -1 + ttoEFF + block)
colnames(design)
=c("wtC","wtM","koC","koM","A1","A2","A3","A4","A5")

CM= makeContrasts(wtM-koM,wtC-koC,levels=design)

fit=lmFit(eset,design)
fit2=contrasts.fit(fit,CM)
fit2=eBayes(fit2)

>From " fit=lmFit(eset,design)" I get  that the coefficient for the animal 5
is not estimable.

I am considering the nested nature of the animals to its strain by the
codification that I am using for the animals, and I have enough degrees of
freedom to estimate the whole set of parameters of the model, so I do not
understand why parameter A5 can not be estimated.

If I work only in the arrays 1 to 6, considering only the strain WT
separately from each other, the things work fine. For arrays 7 to 12
corresponding to the KO strain the things work ok as well , but if I extend
the model to include the whole data,  the parameter A5 corresponding to the
animal 5 is not estimable and I don’t understand why.

I want to use the whole data, because I am also interested in the
comparisons between lines:

CM= makeContrasts(wtC - wtM,koC-koM,levels=design)

I will appreciate any help on this.
Thank you very much.

Pedro.-

My session info is:

> sessionInfo()
Version 2.3.0 (2006-04-24)
i386-pc-mingw32

attached base packages:
[1] "splines"   "tools"     "methods"   "stats"     "graphics"  "grDevices"
"utils"     "datasets"
[9] "base"

other attached packages:
statmod        qvalue          MASS           sma   geneplotter
affycoretools       GOstats
"1.2.4"       "1.5.0"      "7.2-27"      "0.5.15"      "1.10.0"
"1.4.0"       "1.6.0"
Category      hgu95av2          KEGG      multtest    genefilter
survival        xtable
"1.4.0"      "1.12.0"      "1.12.0"      "1.10.0"      "1.10.0"
"2.24"       "1.3-1"
RBGL            GO         graph         Ruuid       cluster