[BioC] limma, voom, weights: factorial versus nested interaction results

[Reiner Schulz]

hi Ryan,

the voom call is:
counts.voom= voom( counts, design= design.matrix, plot=T,
lib.size=nf*colSums( counts), normalize.method= 'none')
where nf= calcNormFactors( counts) and design.matrix as below.

fitting is done using:
fit.limma= lmFit( counts.voom$E, design= design.matrix, weights=
counts.voom$weights)
and:
fit.limma.alt= lmFit( counts.voom$E, design= design.matrix.alt, weights=
counts.voom$weights)
where design.matrix.alt as below. so, the only difference is in the
design matrix; expression values and weights are the same. yet, i
observe the unexpected inequalities below. without 'weights=
counts.voom$weights', everything is as expected.

also, relevant bits from sessionInfo():
R version 2.15.2 (2012-10-26)
Platform: x86_64-pc-linux-gnu (64-bit)
limma_3.14.3
edgeR_3.0.7
DESeq_1.10.1

best, r

On 02/01/13 19:47, Ryan C. Thompson wrote:
> It would probably help to have the exact code that you used to call
> voom. voom can take a lot of arguments, and I'm not sure which ones
> would make a difference here.
> 
> On 01/02/2013 07:31 AM, Reiner Schulz wrote:
>> Dear all,
>>
>> Am analysing an RNA-seq experiment with a two-way factorial design, and
>> am observing different results for two contrasts of interest depending
>> on whether I use a factorial model matrix or a nested interaction model
>> matrix. I only observe this if the data are pre-processed w/ voom:
>>> counts.voom= voom( ...
>> i.e., when there are observational level weights passed to lmFit. If I
>> only pass counts.voom$E (no observational level weights) to lmFit, the
>> results for the two contrasts of interest are identical for the
>> factorial and nested interaction model matrices.
>> Is this to be expected, and if so, what is the explanation?
>>
>> Much appreciated, Reiner
>>
>> PS: some more details ...
>>
>>> strain
>>   [1] MUT MUT MUT MUT MUT MUT MUT MUT MUT MUT MUT MUT MUT WT  WT  WT  WT
>>   WT  WT
>> [20] WT  WT  WT  WT  WT  WT  WT
>> Levels: WT MUT
>>> treat
>>   [1] U  T2 T1 U  T2 T1 U  U  T2 T1 U  T2 T1 U  T2 T1 U  T2 T1 U  U  T2
>> T1 U  T2
>> [26] T1
>> Levels: U T1 T2
>>
>>> cond= factor( paste( sep='.', strain, treat), lev=c('WT.U', 'MUT.U',
>> 'WT.T1', 'MUT.T1', 'WT.T2', 'MUT.T2'))
>>> design.matrix= model.matrix( ~0 + cond)
>>> design.matrix
>>     WT.U MUT.U WT.T1 MUT.T1 WT.T2 MUT.T2
>> 1  0    1     0     0      0      0
>> 2  0    0     0     0      0      1
>> 3  0    0     0     1      0      0
>> 4  0    1     0     0      0      0
>> 5  0    0     0     0      0      1
>> 6  0    0     0     1      0      0
>> 7  0    1     0     0      0      0
>> 8  0    1     0     0      0      0
>> 9  0    0     0     0      0      1
>> 10 0    0     0     1      0      0
>> 11 0    1     0     0      0      0
>> 12 0    0     0     0      0      1
>> 13 0    0     0     1      0      0
>> 14 1    0     0     0      0      0
>> 15 0    0     0     0      1      0
>> 16 0    0     1     0      0      0
>> 17 1    0     0     0      0      0
>> 18 0    0     0     0      1      0
>> 19 0    0     1     0      0      0
>> 20 1    0     0     0      0      0
>> 21 1    0     0     0      0      0
>> 22 0    0     0     0      1      0
>> 23 0    0     1     0      0      0
>> 24 1    0     0     0      0      0
>> 25 0    0     0     0      1      0
>> 26 0    0     1     0      0      0
>>
>>> design.matrix.alt= model.matrix( ~strain + strain:treat)
>>> design.matrix.alt
>>     Intercept MUT WT:T1 MUT:T1 WT:T2 MUT:T2
>> 1  1         1   0     0      0     0
>> 2  1         1   0     0      0     1
>> 3  1         1   0     1      0     0
>> 4  1         1   0     0      0     0
>> 5  1         1   0     0      0     1
>> 6  1         1   0     1      0     0
>> 7  1         1   0     0      0     0
>> 8  1         1   0     0      0     0
>> 9  1         1   0     0      0     1
>> 10 1         1   0     1      0     0
>> 11 1         1   0     0      0     0
>> 12 1         1   0     0      0     1
>> 13 1         1   0     1      0     0
>> 14 1         0   0     0      0     0
>> 15 1         0   0     0      1     0
>> 16 1         0   1     0      0     0
>> 17 1         0   0     0      0     0
>> 18 1         0   0     0      1     0
>> 19 1         0   1     0      0     0
>> 20 1         0   0     0      0     0
>> 21 1         0   0     0      0     0
>> 22 1         0   0     0      1     0
>> 23 1         0   1     0      0     0
>> 24 1         0   0     0      0     0
>> 25 1         0   0     0      1     0
>> 26 1         0   1     0      0     0
>>
>> Observed equalities (as expected):
>> MUT == MUT.U - WT.U
>> WT:T1 == WT.T1 - WT.U
>> MUT:T1 == MUT.T1 - MUT.U
>> WT:T2 == WT.T2 - WT.U
>> MUT:T2 == MUT.T2 - MUT.U
>> MUT:T1 - WT:T1 == (MUT.T1 - MUT.U) - (WT.T1 - WT.U)
>> MUT:T2 - WT:T2 == (MUT.T2 - MUT.U) - (WT.T2 - WT.U)
>>
>> Observed inequalities (only w/ voom pre-processed data; unexpected):
>> WT:T1 - WT:T2 != (WT.T1 - WT.U) - (WT.T2 - WT.U)
>> MUT:T1 - MUT:T2 != (MUT.T1 - MUT.U) - (MUT.T2 - MUT.U)
>>
>> -- 
>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives:
>> http://news.gmane.org/gmane.science.biology.informatics.conductor
> 
> 

-- 
Reiner Schulz, Ph.D.
Senior Lecturer in Bioinformatics and Epigenomics
Dept. of Medical and Molecular Genetics
King's College London
Guy's Hospital, 8th floor Tower Wing
London SE1 9RT
https://atlas.genetics.kcl.ac.uk/~rschulz