[BioC] Non-conformable arguments fitting limma contrasts
Paul Boutros
Paul.Boutros at utoronto.ca
Wed Aug 11 06:06:11 CEST 2004
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 at 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 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