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

[Ryan C. Thompson]

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