[BioC] Non-conformable arguments fitting limma contrasts
Kasper Daniel Hansen
k.hansen at biostat.ku.dk
Wed Jul 28 22:09:19 CEST 2004
Gordon Smyth <smyth at 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 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
>> >
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://www.stat.math.ethz.ch/mailman/listinfo/bioconductor
>
--
Kasper Daniel Hansen, Research Assistant
Department of Biostatistics, University of Copenhagen
More information about the Bioconductor
mailing list