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) 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. 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

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

Gordon Smyth <smyth@wehi.edu.au> writes: > 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? I think you need to be careful here. Many users of limma (as evident from this list) are not too familiar with estimable contrasts. However, _I_ think it would be unresonable to disallow singular matrices as advanced users sometimes may have need for them (I have used them frequently, allthough not in micro array analysis). Suggestion: add an argument to limma such as singular.ok = FALSE Hence singular matrices would be disallowed per default. This is how "lm" deals with the problem - however in lm's case the default is TRUE which I think would be dangerous for limma given the users. /Kasper > > 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 >> > > > _______________________________________________ > Bioconductor mailing list > Bioconductor@stat.math.ethz.ch > https://www.stat.math.ethz.ch/mailman/listinfo/bioconductor > -- Kasper Daniel Hansen, Research Assistant Department of Biostatistics, University of Copenhagen