[BioC] Non-conformable arguments fitting limma contrasts
Gordon Smyth
smyth at wehi.edu.au
Wed Jul 28 10:45:49 CEST 2004
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 at wehi.edu.au]
> > Sent: Tuesday, July 27, 2004 10:15 PM
> > To: paul.boutros at 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 at 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
> >
More information about the Bioconductor
mailing list