[BioC] Non-conformable arguments fitting limma contrasts
Gordon Smyth
smyth at wehi.edu.au
Wed Jul 28 04:14:34 CEST 2004
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