[BioC] limma - different parametrization and weights
Hans-Ulrich Klein
h.klein at uni-muenster.de
Tue Nov 4 11:56:17 CET 2008
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
More information about the Bioconductor
mailing list