Gordon, A combination of upgrading limma (to 1.7.1) and fixing the design matrix got things working. I haven't converted to R formula-language, but I'll look into that as well. As always, thanks for the outstanding help/support, Paul > -----Original Message----- > From: Gordon Smyth [mailto:smyth@wehi.edu.au] > Sent: Thursday, June 24, 2004 10:46 PM > To: paul.boutros@utoronto.ca > Cc: BioConductor Mailing List; LHollins@PICR.man.ac.uk > Subject: Re: [BioC] Error in La.chol(t(A/s)/s) fitting limma models > > > Paul, > > There are several things going on here, a couple of which are my fault. > > 1. The lmFit() function does handle singular design matrices, in the same > way that the lm() function in the base package does. Any > coefficients which > are inestimable will simply be set to NA. So my earlier post > https://stat.ethz.ch/pipermail/bioconductor/2004-May/004833.html claiming > that singular matrices produce an error messages from lmFit() was > not correct. > > 2. Prior to BioC release 1.4, the contrasts.fit() function was only > intended to handle orthogonal design matrices arising from > one-way layouts > (as described in the help entry). From BioC release 1.4, > contrasts.fit() is > intended to handle general design and contrast matrices. But I had > overlooked the need to allow for singular design matrices. Hence > the error > message you have got. I will correct this because I would like > contrasts.fit() to handle without errors any MArrayLM object produced by > lmFit(). > > 3. Your design matrix below is not correct. Notice that column 7 is equal > to the sum of columns 3 and 4, hence the singularity. You don't give > reproducible code to compute the design matrix, rather you have produced > the matrix by hand editing, and you don't give pData(eset), so we can't > tell you what the correct code should have been. Since you have 4 factors > in your experiment, it might well be worthwhile to use the formula > representations of factorial designs in R and the resulting > parametrizations. (I know these can be unintuitive for many people.) > > 4. You might not even need a contrast matrix if the coefficients of > interest to you are already in your linear model. > > Gordon > > At 11:36 AM 25/06/2004, paul.boutros@utoronto.ca wrote: > >Hello, > > > >I'm having some problems fitting a linear model: > > > >### BEGIN CODE #################### > ># libraries > >library(gcrma); > >library(limma); > > > ># Normalize via GCRMA > >eset <- justGCRMA(); > > > ># Save normalized data to disk > >write.exprs(eset, "gcrma.txt"); > > > ># look at phenotypic data > >pData(eset); > > > ># create a dummy design matrix: > >design <- model.matrix(~ -1 + factor(c > >(1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,6,6, > 6,6,6,6,6,6,7,8 > >))); > >colnames(design) <- c > >("basal", "LnC", "HW", "LnA", "drug", "time", "resistance", > >"drug_resistance"); > > > ># then make it real manually > >design <- edit(design) > > > ># make a contrast matrix > >contrast.matrix <- makeContrasts(drug, time, resistance, drug_resistance, > >levels=design); > > > ># fit the model > >fit1 <- lmFit(eset, design); > >fit2 <- contrasts.fit(fit1, contrast.matrix); > >### END CODE ###################### > > > >The error: > >Error in La.chol(t(A/s)/s) : the leading minor of order 4 is not positive > >definite > > > > > design; > > basal LnC HW LnA drug time resistance drug_resistance > >1 1 0 1 0 1 1 1 1 > >2 1 0 1 0 1 1 1 1 > >3 1 0 1 0 1 1 1 1 > >4 1 0 1 0 1 1 1 1 > >5 1 0 0 0 1 1 0 0 > >6 1 0 0 0 1 1 0 0 > >7 1 0 0 0 1 1 0 0 > >8 1 0 0 0 1 1 0 0 > >9 1 0 0 0 1 0 0 0 > >10 1 0 0 0 0 0 0 0 > >11 1 0 1 0 0 0 1 0 > >12 1 0 1 0 1 0 1 1 > >13 1 0 0 0 0 0 0 0 > >14 1 0 0 0 0 0 0 0 > >15 1 0 0 0 0 0 0 0 > >16 1 0 1 0 0 0 1 0 > >17 1 0 1 0 0 0 1 0 > >18 1 0 0 0 1 0 0 0 > >19 1 0 0 0 1 0 0 0 > >20 1 0 0 0 1 0 0 0 > >21 1 0 1 0 1 0 1 1 > >22 1 0 1 0 1 0 1 1 > >23 1 0 1 0 1 0 1 1 > >24 1 0 1 0 0 0 1 0 > >25 1 0 0 1 0 0 1 0 > >26 1 0 0 1 0 0 1 0 > >27 1 0 0 1 0 0 1 0 > >28 1 0 0 1 0 0 1 0 > >29 1 0 0 1 1 0 1 1 > >30 1 0 0 1 1 0 1 1 > >31 1 0 0 1 1 0 1 1 > >32 1 1 0 0 0 0 0 0 > >33 1 1 0 0 0 0 0 0 > >34 1 1 0 0 0 0 0 0 > >35 1 1 0 0 0 0 0 0 > >36 1 1 0 0 1 0 0 0 > >37 1 1 0 0 1 0 0 0 > >38 1 1 0 0 1 0 0 0 > >39 1 1 0 0 1 0 0 0 > >40 1 0 0 1 1 0 1 1 > > > > > contrast.matrix; > > drug time resistance drug_resistance > >basal 0 0 0 0 > >LnC 0 0 0 0 > >HW 0 0 0 0 > >LnA 0 0 0 0 > >drug 1 0 0 0 > >time 0 1 0 0 > >resistance 0 0 1 0 > >drug_resistance 0 0 0 1 > > > >This problem has come up previously: > >https://stat.ethz.ch/pipermail/bioconductor/2004-May/004802.html > > > >And Gordon responded that it was a design-matrix not of full-rank, and > >that an > >error should have been thrown by lmFit. > >https://stat.ethz.ch/pipermail/bioconductor/2004-May/004833.html > > > >So two things: > >1. That error doesn't seem to be thrown in my case: not sure if that's a > >feature or a bug, but it does go against what Gordon indicated > the expected > >behaviour would be. > > > >2. When I examine my fitted object: > > > > > fit1; > >An object of class "MArrayLM" > >$coefficients > > basal LnC HW LnA drug > time > >resistance drug_resistance > >[1,] 9.610458 -0.09169877 0.15457049 -0.14573919 0.02652003 - > >0.3049501 NA -0.10234624 > >[2,] 9.187651 -0.12767969 0.08154156 0.05242885 0.14016087 - > >0.1239801 NA -0.17493971 > >[3,] 9.153906 -0.05092235 -0.10669856 0.09775670 0.46374663 - > >0.3449803 NA -0.20873256 > >[4,] 11.083109 -0.08410092 -0.04708751 -0.11681116 0.20109553 - > >0.3318239 NA -0.05279189 > >[5,] 11.708261 -0.14465592 -0.03052111 -0.25919099 -0.11353194 - > >0.1448190 NA 0.05919881 > >15918 more rows ... > > > >And it's apparent that resistance is not being fitted at all, > indicating (I > >think) that my matrix isn't of full rank. I'm not catching how I > >mis-specified > >the matrix, though. > > > >The experiment has four factors: > >strain: four levels > >treatment: two levels > >time: two levels > >resistance: two levels > > > >I wasn't sure how to fit the four-level strain, so I fitted > three factors as > >separate levels. I also explicitly fitted the drug-resistance > interaction > >into > >the design-matrix. I'd guess the latter is the source of my > problems, but > >any > >pointers on how I ought to have handled both these issues are > very, very much > >appreciated. > > > >Paul >