[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)
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"
More information about the Bioconductor
mailing list