[BioC] blocking in limma gives an estimability problem

[Pedro López Romero]

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)

	
DE=decideTests(fit2,method="nestedF",adjust.method="BY",p.value=0.05)

>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
codelink      annotate 
      "1.8.0"      "1.12.0"      "1.10.0"      "1.10.0"      "1.10.5"
"1.0.0"      "1.10.0" 
        limma          affy        affyio       Biobase 
      "2.6.0"      "1.10.0"       "1.0.0"      "1.10.0"