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