[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