Thank you: that sent me on a better path. From a design point of view, the root of the problem was that I had too many parameters, some of which had no biological meaning. Fixing that allowed proper model-fitting, and allows me to avoid the issue of limma not (yet) detecting non-estimable coefficients. Paul > -----Original Message----- > From: Gordon Smyth [mailto:smyth@wehi.edu.au] > Sent: Wednesday, July 28, 2004 4:46 AM > To: Paul Boutros > Cc: BioConductor Mailing List > Subject: RE: [BioC] Non-conformable arguments fitting limma contrasts > > > The problem is caused by the fact that your design matrix is singular -- > the coefficient for WT is non-estimable. The limma code is supposed to > detect this eventuality, but there is still a bug which I will now fix. > > The question remains, what should the code ideally do when a user > asks for > contrasts involving non-estimable coefficients, as you have done? > The limma > code is supposed to remove from the contrast matrix any rows > corresponding > to non-estimable coefficients. In your example, you have asked > for -WT as a > contrast. This contrast is non-estimable, so the coefficient for that > contrast will be set to zero, but with NAs for t-statistics. That is the > best that can be done, but it may surprise some users. Perhaps I should > simply disallow singular design matrices? > > Gordon > > At 02:40 PM 28/07/2004, Paul Boutros wrote: > >Hi Gordon, > > > >The problem did remain after upgrading to 1.7.4 so below is the output of > >show(fit1). Let me know if the output wraps illegibly and I'll > try to send > >it as an attachment instead. > > > >Thanks, > >Paul > > > > > show(fit1); > >An object of class "MArrayLM" > >$coefficients > > HW LE LnA LnC TCDD Time WT > >WT.TCDD WT.Time > >[1,] 9.313025 9.139576 8.972385 9.087700 -0.038905622 -0.2986598 NA > >0.008103176 0.24776335 > >[2,] 9.996145 9.909185 9.915040 9.836435 0.007888667 -0.1400902 NA > >0.059780420 0.15182790 > >[3,] 9.545351 9.641814 9.689196 9.566151 0.155838277 -0.2069456 NA > >0.177727270 -0.06770628 > >[4,] 10.497738 10.503054 10.420403 10.435289 0.130853665 -0.2115233 NA > >0.009263846 0.09493219 > >[5,] 10.868702 10.825530 10.610891 10.723626 -0.036030238 -0.1316818 > >NA -0.112144889 0.21864944 > >15918 more rows ... > > > >$stdev.unscaled > > HW LE LnA LnC TCDD Time WT WT.TCDD > >WT.Time > >[1,] 0.4330127 0.4330127 0.4330127 0.4330127 0.5 0.6614378 NA 0.7071068 > >0.9354143 > >[2,] 0.4330127 0.4330127 0.4330127 0.4330127 0.5 0.6614378 NA 0.7071068 > >0.9354143 > >[3,] 0.4330127 0.4330127 0.4330127 0.4330127 0.5 0.6614378 NA 0.7071068 > >0.9354143 > >[4,] 0.4330127 0.4330127 0.4330127 0.4330127 0.5 0.6614378 NA 0.7071068 > >0.9354143 > >[5,] 0.4330127 0.4330127 0.4330127 0.4330127 0.5 0.6614378 NA 0.7071068 > >0.9354143 > >15918 more rows ... > > > >$sigma > >[1] 0.10520832 0.05827840 0.09787933 0.11734057 0.09119965 > >15918 more elements ... > > > >$df.residual > >[1] 32 32 32 32 32 > >15918 more elements ... > > > >$cov.coefficients > > [,1] [,2] [,3] [,4] > [,5] > >[,6] [,7] [,8] > >[1,] 1.875000e-01 -1.197959e-17 > > 6.250000e-02 -5.367648e-18 -1.250000e-01 -6.250000e-02 0.125 0.0625 > >[2,] -1.197959e-17 1.875000e-01 -1.601570e-17 6.250000e-02 > > 2.335778e-17 -1.945126e-17 -0.125 -0.0625 > >[3,] 6.250000e-02 -1.601570e-17 1.875000e-01 -5.668344e-18 > -1.250000e-01 > >6.250000e-02 0.125 -0.0625 > >[4,] -5.367648e-18 6.250000e-02 -5.668344e-18 1.875000e-01 > > 6.132519e-18 -4.834864e-18 -0.125 0.0625 > >[5,] -1.250000e-01 2.335778e-17 -1.250000e-01 6.132519e-18 > > 2.500000e-01 -1.250000e-01 -0.250 0.1250 > >[6,] -6.250000e-02 -1.945126e-17 6.250000e-02 -4.834864e-18 > -1.250000e-01 > >4.375000e-01 0.125 -0.4375 > >[7,] 1.250000e-01 -1.250000e-01 1.250000e-01 -1.250000e-01 > -2.500000e-01 > >1.250000e-01 0.500 -0.2500 > >[8,] 6.250000e-02 -6.250000e-02 -6.250000e-02 6.250000e-02 > > 1.250000e-01 -4.375000e-01 -0.250 0.8750 > > > >$pivot > >[1] 1 2 3 4 5 6 8 9 7 > > > >$method > >[1] "ls" > > > >$design > > HW LE LnA LnC TCDD Time WT WT.TCDD WT.Time > >RAE230A_021304_IM01T_LH.CEL 1 0 0 0 1 1 0 0 0 > >RAE230A_021304_IM02T_LH.CEL 1 0 0 0 1 1 0 0 0 > >RAE230A_021304_IM03T_LH.CEL 1 0 0 0 1 1 0 0 0 > >RAE230A_021304_IM04T_LH.CEL 1 0 0 0 1 1 0 0 0 > >RAE230A_021304_IM05T_LH.CEL 0 1 0 0 1 1 1 1 1 > >35 more rows ... > > > >$genes > >[1] "1367452_at" "1367453_at" "1367454_at" "1367455_at" "1367456_at" > >15918 more rows ... > > > >$Amean > >1367452_at 1367453_at 1367454_at 1367455_at 1367456_at > > 9.091929 9.931726 9.705881 10.519856 10.715440 > >15918 more elements ... > > > > > > > -----Original Message----- > > > From: Gordon Smyth [mailto:smyth@wehi.edu.au] > > > Sent: Tuesday, July 27, 2004 10:15 PM > > > To: paul.boutros@utoronto.ca > > > Cc: BioConductor Mailing List > > > Subject: Re: [BioC] Non-conformable arguments fitting limma contrasts > > > > > > > > > Dear Paul, > > > > > > At 11:43 AM 27/07/2004, paul.boutros@utoronto.ca wrote: > > > >Hello, > > > > > > > >I've come across an error when fitting contrasts with Affy data > > > using limma. > > > >I'm on: > > > >WinXP SP1 > > > >R1.9.1 > > > >limma 1.7.2 > > > >gcrma 1.1.0 > > > >arrays: RAE230A > > > > > > > >The error message I'm getting is: > > > > > fit2 <- contrasts.fit(fit1, contrast.matrix); > > > >Error in o %*% RUC^2 : non-conformable arguments > > > > > traceback(); > > > >1: contrasts.fit(fit1, contrast.matrix) > > > > > > I don't know what the problem is straight away. There have > however been > > > non-trivial changes to contrasts.fit() recently as limma has been > > > gradually > > > upgraded to handle more general correlation structures arising > > > from various > > > types of technical replication. Can I get you to: > > > > > > 1. upgrade to limma 1.7.4 and see if the problem goes away > > > 2. type show(fit1) and send me the output > > > > > > Note that in your example you don't strictly need to use > > > contrasts.fit() at > > > all, because all the contrasts are already available as > coefficients in > > > your design matrix. > > > > > > >Below I've appended the sequence of commands I'm using, > followed by my > > > >pData(eset), my design and contrast matrices. > > > > > > > >Briefly, what I have are the following factors: > > > >Strain: four levels (HW, LE, LnC, LnA) > > > >Time: two levels (hour3, hour19) > > > >Treatment: two levels (TCDD, COil) > > > >AHR Genotype: two levels (HW, WT) > > > > > > > >Where I believe I may have gone wrong is explicitly including the AHR > > > >genotype > > > >in the model. This factor can be determined by the strain: two > > > strains (HW, > > > >LnA) always have AHR = HW while the other two strains always > > > have AHR = WT. > > > >Including both surely feels like double-fitting to me. On the > > > other hand, > > > >the > > > >AHR genotype is *not* the only factor differing between the > > > strains which > > > >was my > > > >motivation for including it. > > > > > > model.matrix() will remove any singularities in the design > > > matrix, so there > > > should be no problem with over-parametrization. > > > > > > Gordon > > > > > > >As always, I'm open to any comments/suggestions/corrections. > > > >Paul > > > > > > > > > > > >### Start Command Sequence ### > > > >eset <- ReadAffy(filenames=cel.files, phenoData="phenodata.txt"); > > > >eset <- gcrma(eset); > > > > > > > ># create a design matrix: > > > >design <- model.matrix(~-1 + Strain + Treatment + Time + AHR + > > > >AHR:Treatment + > > > >AHR:Time, pData(eset)); > > > >colnames(design) <- c("HW", "LE", "LnA", "LnC", "TCDD", "Time", "WT", > > > >"WT.TCDD", > > > >"WT.Time"); > > > > > > > ># make a contrast matrix > > > >contrast.matrix <- makeContrasts(TCDD, Time, -WT, -WT.TCDD, -WT.Time, > > > >levels=design); > > > > > > > ># fit the model > > > >fit1 <- lmFit(eset, design); > > > >fit2 <- contrasts.fit(fit1, contrast.matrix); > > > >fit3 <- eBayes(fit2); > > > >### End Command Sequence ### > > > > > > > > > pData(eset); > > > > Strain Treatment Time AHR > > > >RAE230A_021304_IM01T_LH.CEL HW TCDD Hour3 HW > > > >RAE230A_021304_IM02T_LH.CEL HW TCDD Hour3 HW > > > >RAE230A_021304_IM03T_LH.CEL HW TCDD Hour3 HW > > > >RAE230A_021304_IM04T_LH.CEL HW TCDD Hour3 HW > > > >RAE230A_021304_IM05T_LH.CEL LE TCDD Hour3 WT > > > >RAE230A_021304_IM06T_LH.CEL LE TCDD Hour3 WT > > > >RAE230A_021304_IM07T_LH.CEL LE TCDD Hour3 WT > > > >RAE230A_021304_IM08T_LH.CEL LE TCDD Hour3 WT > > > >RAE230A_043003_IM01T_LH.CEL LE TCDD Hour19 WT > > > >RAE230A_043003_IM02T_LH.CEL LE COil Hour19 WT > > > >RAE230A_043003_IM03T_LH.CEL HW COil Hour19 HW > > > >RAE230A_043003_IM04T_LH.CEL HW TCDD Hour19 HW > > > >RAE230A_060403_IM01T_LH.CEL LE COil Hour19 WT > > > >RAE230A_060403_IM02T_LH.CEL LE COil Hour19 WT > > > >RAE230A_060403_IM03T_LH.CEL LE COil Hour19 WT > > > >RAE230A_060403_IM04T_LH.CEL HW COil Hour19 HW > > > >RAE230A_060403_IM05T_LH.CEL HW COil Hour19 HW > > > >RAE230A_060403_IM06T_LH.CEL LE TCDD Hour19 WT > > > >RAE230A_060403_IM07T_LH.CEL LE TCDD Hour19 WT > > > >RAE230A_060403_IM08T_LH.CEL LE TCDD Hour19 WT > > > >RAE230A_060403_IM09T_LH.CEL HW TCDD Hour19 HW > > > >RAE230A_060403_IM10T_LH.CEL HW TCDD Hour19 HW > > > >RAE230A_070303_IM01T_LH.CEL HW TCDD Hour19 HW > > > >RAE230A_070303_IM02T_LH.CEL HW COil Hour19 HW > > > >RAE230A_102803_IM01T_LH.CEL LnA COil Hour19 HW > > > >RAE230A_102803_IM02T_LH.CEL LnA COil Hour19 HW > > > >RAE230A_102803_IM03T_LH.CEL LnA COil Hour19 HW > > > >RAE230A_102803_IM04T_LH.CEL LnA COil Hour19 HW > > > >RAE230A_102803_IM05T_LH.CEL LnA TCDD Hour19 HW > > > >RAE230A_102803_IM06T_LH.CEL LnA TCDD Hour19 HW > > > >RAE230A_102803_IM07T_LH.CEL LnA TCDD Hour19 HW > > > >RAE230A_102803_IM08T_LH.CEL LnC COil Hour19 WT > > > >RAE230A_102803_IM09T_LH.CEL LnC COil Hour19 WT > > > >RAE230A_102803_IM10T_LH.CEL LnC COil Hour19 WT > > > >RAE230A_102803_IM11T_LH.CEL LnC COil Hour19 WT > > > >RAE230A_102803_IM12T_LH.CEL LnC TCDD Hour19 WT > > > >RAE230A_102803_IM13T_LH.CEL LnC TCDD Hour19 WT > > > >RAE230A_102803_IM14T_LH.CEL LnC TCDD Hour19 WT > > > >RAE230A_102803_IM15T_LH.CEL LnC TCDD Hour19 WT > > > >RAE230A_102803_IM16T_LH.CEL LnA TCDD Hour19 HW > > > > > > > > > design; > > > > > > > > HW LE LnA LnC TCDD Time WT > WT.TCDD WT.Time > > > >RAE230A_021304_IM01T_LH.CEL 1 0 0 0 1 1 0 > 0 0 > > > >RAE230A_021304_IM02T_LH.CEL 1 0 0 0 1 1 0 > 0 0 > > > >RAE230A_021304_IM03T_LH.CEL 1 0 0 0 1 1 0 > 0 0 > > > >RAE230A_021304_IM04T_LH.CEL 1 0 0 0 1 1 0 > 0 0 > > > >RAE230A_021304_IM05T_LH.CEL 0 1 0 0 1 1 1 > 1 1 > > > >RAE230A_021304_IM06T_LH.CEL 0 1 0 0 1 1 1 > 1 1 > > > >RAE230A_021304_IM07T_LH.CEL 0 1 0 0 1 1 1 > 1 1 > > > >RAE230A_021304_IM08T_LH.CEL 0 1 0 0 1 1 1 > 1 1 > > > >RAE230A_043003_IM01T_LH.CEL 0 1 0 0 1 0 1 > 1 0 > > > >RAE230A_043003_IM02T_LH.CEL 0 1 0 0 0 0 1 > 0 0 > > > >RAE230A_043003_IM03T_LH.CEL 1 0 0 0 0 0 0 > 0 0 > > > >RAE230A_043003_IM04T_LH.CEL 1 0 0 0 1 0 0 > 0 0 > > > >RAE230A_060403_IM01T_LH.CEL 0 1 0 0 0 0 1 > 0 0 > > > >RAE230A_060403_IM02T_LH.CEL 0 1 0 0 0 0 1 > 0 0 > > > >RAE230A_060403_IM03T_LH.CEL 0 1 0 0 0 0 1 > 0 0 > > > >RAE230A_060403_IM04T_LH.CEL 1 0 0 0 0 0 0 > 0 0 > > > >RAE230A_060403_IM05T_LH.CEL 1 0 0 0 0 0 0 > 0 0 > > > >RAE230A_060403_IM06T_LH.CEL 0 1 0 0 1 0 1 > 1 0 > > > >RAE230A_060403_IM07T_LH.CEL 0 1 0 0 1 0 1 > 1 0 > > > >RAE230A_060403_IM08T_LH.CEL 0 1 0 0 1 0 1 > 1 0 > > > >RAE230A_060403_IM09T_LH.CEL 1 0 0 0 1 0 0 > 0 0 > > > >RAE230A_060403_IM10T_LH.CEL 1 0 0 0 1 0 0 > 0 0 > > > >RAE230A_070303_IM01T_LH.CEL 1 0 0 0 1 0 0 > 0 0 > > > >RAE230A_070303_IM02T_LH.CEL 1 0 0 0 0 0 0 > 0 0 > > > >RAE230A_102803_IM01T_LH.CEL 0 0 1 0 0 0 0 > 0 0 > > > >RAE230A_102803_IM02T_LH.CEL 0 0 1 0 0 0 0 > 0 0 > > > >RAE230A_102803_IM03T_LH.CEL 0 0 1 0 0 0 0 > 0 0 > > > >RAE230A_102803_IM04T_LH.CEL 0 0 1 0 0 0 0 > 0 0 > > > >RAE230A_102803_IM05T_LH.CEL 0 0 1 0 1 0 0 > 0 0 > > > >RAE230A_102803_IM06T_LH.CEL 0 0 1 0 1 0 0 > 0 0 > > > >RAE230A_102803_IM07T_LH.CEL 0 0 1 0 1 0 0 > 0 0 > > > >RAE230A_102803_IM08T_LH.CEL 0 0 0 1 0 0 1 > 0 0 > > > >RAE230A_102803_IM09T_LH.CEL 0 0 0 1 0 0 1 > 0 0 > > > >RAE230A_102803_IM10T_LH.CEL 0 0 0 1 0 0 1 > 0 0 > > > >RAE230A_102803_IM11T_LH.CEL 0 0 0 1 0 0 1 > 0 0 > > > >RAE230A_102803_IM12T_LH.CEL 0 0 0 1 1 0 1 > 1 0 > > > >RAE230A_102803_IM13T_LH.CEL 0 0 0 1 1 0 1 > 1 0 > > > >RAE230A_102803_IM14T_LH.CEL 0 0 0 1 1 0 1 > 1 0 > > > >RAE230A_102803_IM15T_LH.CEL 0 0 0 1 1 0 1 > 1 0 > > > >RAE230A_102803_IM16T_LH.CEL 0 0 1 0 1 0 0 > 0 0 > > > >attr(,"assign") > > > >[1] 1 1 1 1 2 3 4 5 6 > > > >attr(,"contrasts") > > > >attr(,"contrasts")$Strain > > > >[1] "contr.treatment" > > > > > > > >attr(,"contrasts")$Treatment > > > >[1] "contr.treatment" > > > > > > > >attr(,"contrasts")$Time > > > >[1] "contr.treatment" > > > > > > > >attr(,"contrasts")$AHR > > > >[1] "contr.treatment" > > > > > > > > > contrast.matrix; > > > > TCDD Time -WT -WT.TCDD -WT.Time > > > >HW 0 0 0 0 0 > > > >LE 0 0 0 0 0 > > > >LnA 0 0 0 0 0 > > > >LnC 0 0 0 0 0 > > > >TCDD 1 0 0 0 0 > > > >Time 0 1 0 0 0 > > > >WT 0 0 -1 0 0 > > > >WT.TCDD 0 0 0 -1 0 > > > >WT.Time 0 0 0 0 -1 > > > >