[BioC] Two interaction coefficients not estimable in limma's lmFit() in a 2^3 factorial design - does this change my contrast matrix now?
Massimo Pinto
pintarello at gmail.com
Mon Aug 10 12:07:00 CEST 2009
James and all,
I thought of reporting here a few final remarks on this issue, in case
anyone will be every taking advantage from this thread.
Here I was trying to estimate 8 parameters when I only had 6 data
points (with four replicates each). At maximum, I could estimate just
6 parameters.
In my 2^3 factorial design, model.matrix() was neatly preparing for me
a design matrix with eight parameters, though the last two could not
be estimated by lmFit(). It turned out that I needed those last two
parameters and that I had to drop other two 'in bewteen', which could
not be estimated in my experiment as they made no real sense. They
were not needed in my model. Therefore, I proceeded from my design
matrix by cutting two of the columns that were created by
model.matrix().
> Dose <- factor(targets$Dose, levels=c("Cn", "1Gy")) # It's very important that you define the levels now, with their order, as otherwise he will do it for you and this may not be the way you want it.
> Lab <- factor(targets$Lab, levels=c("ISS", "LNGS"))
> Ageing <- factor(targets$Ageing, levels=c("t0", "6mo"))
> disegno <- model.matrix(~Dose*Ageing*Lab) # makes a design matrix with several interaction parameters.
> designo <- disegno[,1:3]
> designo1 <- disegno[,5] # I don't need column 4
> designo <- cbind(designo, designo1)
> colnames(designo[4])="Dose1Gy:Ageing6mo"
> designo <- cbind(designo, disegno[,7:8]) # I don't need column 6 either
> is.fullrank(designo)
[1] TRUE
After this, my lmFit() went well and it did not take too long before I
could define all the contrasts of my interest.
Massimo
Massimo Pinto
Post Doctoral Research Fellow
Enrico Fermi Centre and Italian Public Health Research Institute (ISS), Rome
http://claimid.com/massimopinto
On Fri, Aug 7, 2009 at 4:34 PM, Massimo Pinto<pintarello at gmail.com> wrote:
> thank you Jim,
>
> you are indeed correct: I did not specify how I created my matrix.
>
> Here it comes
> ====
> Dose <- factor(targets$Dose, levels=c("Cn", "1Gy")) # It's very
> important that you define the levels now, with their order, as
> otherwise he will do it for you and this may not be the way you want
> it.
> Lab <- factor(targets$Lab, levels=c("ISS", "LNGS"))
> Ageing <- factor(targets$Ageing, levels=c("t0", "6mo"))
> disegno <- model.matrix(~Dose*Ageing*Lab) # makes a design matrix with
> several interaction parameters and a total of eight param.
> ====
>
> given that my data set includes six independent samples (with four
> biological replicates each), I suppose I cannot estimate more than 6
> parameters. But model.matrix() went on anyway...
>
> I suppose I could erase the last two columns of my 'disegno'
>
> Yours
> Massimo
>
> On Wed, Aug 5, 2009 at 6:30 PM, James W. MacDonald<jmacdon at med.umich.edu> wrote:
>> Hi Massimo,
>>
>> Massimo Pinto wrote:
>>>
>>> Greetings all,
>>>
>>> I have noticed that several users have presented the issue of
>>> parameters not estimated in some cases of linear model fittings with
>>> the limma function lmFit().
>>>
>>> In trying to implement my 2^3 factorial design, I have followed
>>> Bioinformatics and Computational Biology Solutions using R and
>>> Bioconductor's examples as at Chapter 14 (Multifactor experiments).
>>> I have encountered my own problem with two such parameters:
>>>
>>>> fit <- lmFit(esetSub, disegno)
>>>
>>> Coefficients not estimable: Ageing6mo:LabLNGS Dose1Gy:Ageing6mo:LabLNGS
>>> Warning message:
>>> Partial NA coefficients for 2176 probe(s)
>>
>> Well, you don't show how you built this design matrix, but I would bet you
>> didn't use model.matrix because the matrix isn't of full rank, so you can't
>> solve for all the coefficients you are trying to estimate here.
>>
>> The lmFit function is nice enough to let you know that, but you could have
>> checked using either is.fullrank() to see if the matrix is of full rank, or
>> nonEstimable() to see which coefficients aren't going to be estimated.
>>
>> Best,
>>
>> Jim
>>
>>
>>>
>>> whereby
>>>
>>>> disegno
>>>
>>> (Intercept) Dose1Gy Ageing6mo LabLNGS Dose1Gy:Ageing6mo
>>> Dose1Gy:LabLNGS Ageing6mo:LabLNGS Dose1Gy:Ageing6mo:LabLNGS
>>> 1 1 0 0 0 0
>>> 0 0 0
>>> 2 1 0 0 0 0
>>> 0 0 0
>>> 3 1 0 0 0 0
>>> 0 0 0
>>> 4 1 0 0 0 0
>>> 0 0 0
>>> 5 1 1 0 0 0
>>> 0 0 0
>>> 6 1 1 0 0 0
>>> 0 0 0
>>> 7 1 1 0 0 0
>>> 0 0 0
>>> 8 1 1 0 0 0
>>> 0 0 0
>>> 9 1 0 1 0 0
>>> 0 0 0
>>> 10 1 0 1 0 0
>>> 0 0 0
>>> 11 1 0 1 0 0
>>> 0 0 0
>>> 12 1 0 1 0 0
>>> 0 0 0
>>> 13 1 1 1 0 1
>>> 0 0 0
>>> 14 1 1 1 0 1
>>> 0 0 0
>>> 15 1 1 1 0 1
>>> 0 0 0
>>> 16 1 1 1 0 1
>>> 0 0 0
>>> 17 1 0 1 1 0
>>> 0 1 0
>>> 18 1 0 1 1 0
>>> 0 1 0
>>> 19 1 0 1 1 0
>>> 0 1 0
>>> 20 1 0 1 1 0
>>> 0 1 0
>>> 21 1 1 1 1 1
>>> 1 1 1
>>> 22 1 1 1 1 1
>>> 1 1 1
>>> 23 1 1 1 1 1
>>> 1 1 1
>>> 24 1 1 1 1 1
>>> 1 1 1
>>> attr(,"assign")
>>> [1] 0 1 2 3 4 5 6 7
>>> attr(,"contrasts")
>>> attr(,"contrasts")$Dose
>>> [1] "contr.treatment"
>>>
>>> attr(,"contrasts")$Ageing
>>> [1] "contr.treatment"
>>>
>>> attr(,"contrasts")$Lab
>>> [1] "contr.treatment"
>>>
>>> I have checked my design and it does make good sense to me. But, I
>>> believe I was asking too much from my parameter estimation.
>>> I am wondering how is this going to change the way I design my
>>> contrast matrix? I am trying to define my contrast matrix following
>>> the formulation of null hypotheses as at the 'yellow' book, p244.
>>>
>>> Thanking you very much,
>>>
>>> Massimo
>>>
>>> Massimo Pinto
>>> Post Doctoral Research Fellow
>>> Enrico Fermi Centre and Italian Public Health Research Institute (ISS),
>>> Rome
>>> http://claimid.com/massimopinto
>>>
>>> _______________________________________________
>>> Bioconductor mailing list
>>> Bioconductor at stat.math.ethz.ch
>>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>>> Search the archives:
>>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>>
>> --
>> James W. MacDonald, M.S.
>> Biostatistician
>> Douglas Lab
>> University of Michigan
>> Department of Human Genetics
>> 5912 Buhl
>> 1241 E. Catherine St.
>> Ann Arbor MI 48109-5618
>> 734-615-7826
>>
>
More information about the Bioconductor
mailing list