[BioC] limma - different parametrization and weights

[Gordon K Smyth]

Dear Hans-Ulrich,

See the help page for contrasts.fit(), the last paragraph of the 
details section.

Best wishes
Gordon

> Date: Tue, 04 Nov 2008 11:56:17 +0100
> From: Hans-Ulrich Klein <h.klein at uni-muenster.de>
> Subject: [BioC] limma - different parametrization and weights
> To: Bioconductor list <bioconductor at stat.math.ethz.ch>
> Message-ID: <49102A51.2030809 at uni-muenster.de>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> Dear All,
>
> I used limma with two different parametrizations. Both approaches should
> be equivalent in my opinion. However, if I use weights, the results
> differ. (Results of both approaches are equal without weights.) I
> attached an example below. Does someone know the reason for this?
>
> Regards,
> Hans-Ulrich
>
>
>
> Here is the example:
>
> library("limma")
>
> n=30
> m=100
> data = matrix(rnorm(n*m, mean=8, sd=1), ncol=n, nrow=m)
> W = matrix(rbinom(n*m, 1, p=0.8), ncol=n, nrow=m)
> W[W==0] = 1/2
> disease = factor(c(rep("D1", n/3), rep("D2", n/3), rep("D3", n/3)))
> batch = factor(sample(c("B1", "B2"), n, replace=TRUE))
>
> D1 = model.matrix(~ disease + batch)
> D2 = model.matrix(~ 0 + disease + batch)
>
> fit1 = lmFit(data, D1, weights=W)
> fit2 = lmFit(data, D2, weights=W)
>
> contrs = makeContrasts(contrasts=
>    c("diseaseD2-diseaseD1", "diseaseD3-diseaseD1"), levels=D2)
> fit2c = contrasts.fit(fit2, contrs)
>
> eFit1 = eBayes(fit1)
> eFit2 = eBayes(fit2c)
>
> topTable(eFit1, coef="diseaseD2")
>        logFC  AveExpr         t    P.Value adj.P.Val         B
> 71  1.0992458 8.414207  2.547634 0.01148965 0.6240313 -4.461615
> 79  1.0701491 8.306736  2.412824 0.01660306 0.6240313 -4.482316
> 65  0.9004096 7.956182  2.061818 0.04033472 0.6240313 -4.515367
> 73  0.8769716 7.822919  2.014672 0.04508839 0.6240313 -4.519390
> 97 -0.8838333 8.059430 -1.969606 0.05006854 0.6240313 -4.526330
> 28 -0.9103940 8.163924 -2.000665 0.04658917 0.6240313 -4.527588
> 92 -0.8426017 7.948664 -1.885870 0.06055682 0.6240313 -4.530451
> 78  0.8518031 7.921508  1.928056 0.05506380 0.6240313 -4.531817
> 96 -0.8062308 8.137124 -1.833017 0.06807629 0.6240313 -4.542318
> 76  0.7613192 8.176955  1.771061 0.07785820 0.6240313 -4.543365
>
> topTable(eFit2, coef="diseaseD2-diseaseD1")
>        logFC  AveExpr         t    P.Value adj.P.Val         B
> 71  1.0992458 8.414207  2.525839 0.01220672 0.6062981 -4.466368
> 79  1.0701491 8.306736  2.344510 0.01989314 0.6062981 -4.495402
> 28 -0.9103940 8.163924 -2.104442 0.03641183 0.6062981 -4.510339
> 65  0.9004096 7.956182  2.083941 0.03825580 0.6062981 -4.511464
> 73  0.8769716 7.822919  2.003999 0.04622820 0.6062981 -4.521191
> 78  0.8518031 7.921508  1.993785 0.04734180 0.6062981 -4.521311
> 92 -0.8426017 7.948664 -1.905554 0.05793940 0.6062981 -4.527255
> 97 -0.8838333 8.059430 -1.942207 0.05331756 0.6062981 -4.530611
> 45  0.7744936 7.745213  1.769790 0.07807041 0.6062981 -4.544968
> 8  -0.7040262 7.806587 -1.693483 0.09170026 0.6062981 -4.547372
>
>
> sessionInfo()
> R version 2.8.0 (2008-10-20)
> x86_64-pc-linux-gnu
>
> locale:
> C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] limma_2.16.2
>
>
> -- 
> Hans-Ulrich Klein
> Department of Medical Informatics and Biomathematics
> University of M?nster, Germany
> Tel.: +49 (0)251 83-58405