[BioC] Can I input ordinal variables into a model in Limma?
Scott Robinson
Scott.Robinson at glasgow.ac.uk
Wed Mar 12 16:04:11 CET 2014
Dear Gordon,
I implemented this and am wondering how you interpret logFC from topTable with regards the linear trend? I have tried replicating it, but no joy.
Thanks,
Scott
PS code below:
design <- model.matrix(~copd + asthma + smoking + sex)
> design
(Intercept) copd.L copd.Q asthmaSEVERE smokingEX-SMOKER
1 1 -7.850462e-17 -0.8164966 0 0
2 1 -7.850462e-17 -0.8164966 0 0
3 1 -7.850462e-17 -0.8164966 0 0
4 1 -7.850462e-17 -0.8164966 0 0
5 1 -7.850462e-17 -0.8164966 0 0
6 1 -7.850462e-17 -0.8164966 0 0
7 1 -7.850462e-17 -0.8164966 0 0
8 1 -7.850462e-17 -0.8164966 0 0
9 1 7.071068e-01 0.4082483 0 0
10 1 7.071068e-01 0.4082483 0 0
11 1 -7.850462e-17 -0.8164966 0 0
12 1 7.071068e-01 0.4082483 0 0
13 1 7.071068e-01 0.4082483 0 0
14 1 -7.850462e-17 -0.8164966 0 0
15 1 7.071068e-01 0.4082483 0 0
16 1 -7.071068e-01 0.4082483 0 0
17 1 -7.071068e-01 0.4082483 0 0
18 1 -7.071068e-01 0.4082483 0 0
19 1 -7.071068e-01 0.4082483 0 0
20 1 -7.071068e-01 0.4082483 0 0
21 1 -7.071068e-01 0.4082483 0 0
22 1 -7.071068e-01 0.4082483 0 0
23 1 -7.071068e-01 0.4082483 0 0
24 1 -7.071068e-01 0.4082483 0 0
25 1 -7.071068e-01 0.4082483 0 0
26 1 -7.071068e-01 0.4082483 0 0
27 1 -7.071068e-01 0.4082483 0 0
28 1 -7.071068e-01 0.4082483 0 0
29 1 -7.071068e-01 0.4082483 0 0
30 1 -7.071068e-01 0.4082483 0 0
31 1 -7.071068e-01 0.4082483 1 0
32 1 -7.071068e-01 0.4082483 1 0
33 1 -7.071068e-01 0.4082483 1 0
34 1 -7.071068e-01 0.4082483 1 0
35 1 -7.071068e-01 0.4082483 1 0
36 1 -7.071068e-01 0.4082483 1 0
37 1 -7.071068e-01 0.4082483 1 0
38 1 -7.071068e-01 0.4082483 1 0
39 1 -7.071068e-01 0.4082483 1 0
40 1 -7.071068e-01 0.4082483 1 0
41 1 -7.071068e-01 0.4082483 1 0
42 1 -7.071068e-01 0.4082483 1 0
43 1 -7.071068e-01 0.4082483 1 0
44 1 -7.071068e-01 0.4082483 1 0
45 1 -7.071068e-01 0.4082483 1 0
46 1 -7.071068e-01 0.4082483 0 0
47 1 -7.071068e-01 0.4082483 0 0
48 1 -7.071068e-01 0.4082483 0 0
49 1 -7.071068e-01 0.4082483 0 0
50 1 -7.071068e-01 0.4082483 0 0
51 1 -7.071068e-01 0.4082483 0 0
52 1 -7.071068e-01 0.4082483 0 0
53 1 -7.071068e-01 0.4082483 0 0
54 1 -7.071068e-01 0.4082483 0 0
55 1 -7.071068e-01 0.4082483 0 0
56 1 -7.071068e-01 0.4082483 0 0
57 1 -7.071068e-01 0.4082483 0 0
58 1 -7.071068e-01 0.4082483 0 0
59 1 -7.071068e-01 0.4082483 0 0
60 1 -7.071068e-01 0.4082483 0 0
61 1 -7.071068e-01 0.4082483 1 0
62 1 -7.071068e-01 0.4082483 1 0
63 1 -7.071068e-01 0.4082483 1 0
64 1 -7.071068e-01 0.4082483 1 0
65 1 -7.071068e-01 0.4082483 1 0
66 1 -7.071068e-01 0.4082483 1 0
67 1 -7.071068e-01 0.4082483 1 0
68 1 -7.071068e-01 0.4082483 1 0
69 1 -7.071068e-01 0.4082483 1 0
70 1 -7.071068e-01 0.4082483 1 0
71 1 -7.071068e-01 0.4082483 1 0
72 1 -7.071068e-01 0.4082483 1 0
73 1 -7.071068e-01 0.4082483 1 0
74 1 -7.071068e-01 0.4082483 1 0
75 1 -7.071068e-01 0.4082483 1 0
76 1 7.071068e-01 0.4082483 0 1
77 1 -7.850462e-17 -0.8164966 0 1
78 1 -7.850462e-17 -0.8164966 0 1
79 1 7.071068e-01 0.4082483 0 1
80 1 7.071068e-01 0.4082483 0 1
81 1 7.071068e-01 0.4082483 0 1
82 1 -7.850462e-17 -0.8164966 0 1
83 1 7.071068e-01 0.4082483 0 1
84 1 7.071068e-01 0.4082483 0 1
85 1 -7.850462e-17 -0.8164966 0 1
86 1 7.071068e-01 0.4082483 0 1
87 1 -7.850462e-17 -0.8164966 0 1
88 1 7.071068e-01 0.4082483 0 1
89 1 -7.850462e-17 -0.8164966 0 1
90 1 -7.850462e-17 -0.8164966 0 1
91 1 -7.850462e-17 -0.8164966 0 0
92 1 -7.850462e-17 -0.8164966 0 0
93 1 7.071068e-01 0.4082483 0 0
94 1 7.071068e-01 0.4082483 0 0
95 1 -7.071068e-01 0.4082483 0 0
96 1 -7.071068e-01 0.4082483 0 0
97 1 -7.071068e-01 0.4082483 0 0
98 1 -7.071068e-01 0.4082483 0 0
99 1 -7.071068e-01 0.4082483 1 0
100 1 -7.071068e-01 0.4082483 1 0
101 1 -7.071068e-01 0.4082483 1 0
102 1 -7.071068e-01 0.4082483 1 0
103 1 -7.071068e-01 0.4082483 0 0
104 1 -7.071068e-01 0.4082483 0 0
105 1 -7.071068e-01 0.4082483 0 0
106 1 -7.071068e-01 0.4082483 0 0
smokingSMOKER sexM
1 1 1
2 1 0
3 1 0
4 1 1
5 1 0
6 1 0
7 1 1
8 1 0
9 1 1
10 1 1
11 1 1
12 1 1
13 1 0
14 1 1
15 1 0
16 0 0
17 0 1
18 0 0
19 0 0
20 0 1
21 0 0
22 0 0
23 0 0
24 0 0
25 0 1
26 0 1
27 0 1
28 0 0
29 0 0
30 0 0
31 1 0
32 1 1
33 1 0
34 1 0
35 1 1
36 1 0
37 1 1
38 1 1
39 1 0
40 1 0
41 1 1
42 1 0
43 1 0
44 1 1
45 1 1
46 1 0
47 1 0
48 1 1
49 1 0
50 1 1
51 1 1
52 1 1
53 1 0
54 1 0
55 1 1
56 1 0
57 1 1
58 1 0
59 1 1
60 1 0
61 0 0
62 0 0
63 0 0
64 0 1
65 0 1
66 0 1
67 0 1
68 0 0
69 0 1
70 0 1
71 0 1
72 0 1
73 0 1
74 0 0
75 0 0
76 0 1
77 0 0
78 0 1
79 0 1
80 0 0
81 0 1
82 0 0
83 0 1
84 0 0
85 0 1
86 0 0
87 0 0
88 0 0
89 0 1
90 0 1
91 1 0
92 1 1
93 1 1
94 1 0
95 0 0
96 0 1
97 0 1
98 0 0
99 1 1
100 1 0
101 1 0
102 1 0
103 1 1
104 1 1
105 1 1
106 1 0
attr(,"assign")
[1] 0 1 1 2 3 3 4
attr(,"contrasts")
attr(,"contrasts")$copd
[1] "contr.poly"
attr(,"contrasts")$asthma
[1] "contr.treatment"
attr(,"contrasts")$smoking
[1] "contr.treatment"
attr(,"contrasts")$sex
[1] "contr.treatment"
>
> corfit <- duplicateCorrelation(data, design, ndups = 1, block = as.factor(techRep))
> fit <- lmFit(data, design, block = as.factor(techRep), cor = corfit$consensus)
> fit <- eBayes(fit)
>
> topTable(fit, coef="copd.L", adjust="BH", number=1)
ID logFC AveExpr t P.Value adj.P.Val B
14244 219133_at -1.06972 7.116668 -5.540948 2.366908e-07 0.006463079 6.574808
>
> #all copd vs all non-copd
> mean(data[14244,which(copd=="GOLD_III" | copd=="GOLD_II")]) - mean(data[14244,which(copd=="control")])
[1] -0.5029798
> #all only-copd vs all healthy
> mean(data[14244,which(copd=="GOLD_III" | copd=="GOLD_II")]) - mean(data[14244,which(copd=="control" & asthma=="control")])
[1] -0.5752098
> #severe vs moderate
> mean(data[14244,which(copd=="GOLD_III")]) - mean(data[14244,which(copd=="GOLD_II")])
[1] -0.3566912
> #severe vs normal
> mean(data[14244,which(copd=="GOLD_III")]) - mean(data[14244,which(copd=="control" & asthma=="control")])
[1] -0.7745372
________________________________________
From: Scott Robinson
Sent: 01 September 2013 16:28
To: Gordon K Smyth
Cc: Bioconductor mailing list
Subject: RE: Can I input ordinal variables into a model in Limma?
Thanks very much Gordon.
________________________________________
From: Gordon K Smyth [smyth at wehi.EDU.AU]
Sent: 01 September 2013 01:51
To: Scott Robinson
Cc: Bioconductor mailing list
Subject: Can I input ordinal variables into a model in Limma?
Dear Scott,
An ordinal variable is just a special case of a categorical variable, and
you include it in the limma design matrix just as you would any other
categorical variable. For example:
DiseaseState <- ordered(DiseaseState,
levels=c("mild", "moderate", "severe"))
design <- model.matrix(~DiseaseState)
etc. By default, the first coefficient in this design will be the overall
mean, the second will be linear trend (monotonic increase or decrease),
and the third will be quadratic trend ("moderate" more extreme than the
other two levels rather than intermediate between them).
Best wishes
Gordon
> Date: Fri, 30 Aug 2013 07:36:58 -0700 (PDT)
> From: "Scott Robinson [guest]" <guest at bioconductor.org>
> To: bioconductor at r-project.org, scott.robinson at glasgow.ac.uk
> Subject: [BioC] Can I input ordinal variables into a model in Limma?
>
>
> I am working on some microarray data where the samples come from
> patients with different severities of disease state - something like
> "mild", "moderate", "severe".
>
> I suppose this is an 'ordinal' variable, but only know how to input
> categorical and continuous variables into the model and searching the
> Limma manual for the word 'ordinal' doesn't get me anywhere.
>
> Is it possible to work ordinal variables into my model? If not and I
> still want to use Limma is it best to treat it as categorical or
> continuous? Or is there an alternative package I could use which has
> this functionality?
>
> Many thanks in advance,
>
> Scott
>
> -- output of sessionInfo():
>
>> sessionInfo()
> R version 3.0.1 (2013-05-16)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> locale:
> [1] LC_COLLATE=English_United Kingdom.1252
> [2] LC_CTYPE=English_United Kingdom.1252
> [3] LC_MONETARY=English_United Kingdom.1252
> [4] LC_NUMERIC=C
> [5] LC_TIME=English_United Kingdom.1252
>
> attached base packages:
> [1] parallel stats graphics grDevices utils datasets methods
> [8] base
>
> other attached packages:
> [1] limma_3.16.7 sparcl_1.0.3 lattice_0.20-23
> [4] corrplot_0.71 affyPLM_1.36.0 preprocessCore_1.22.0
> [7] simpleaffy_2.36.1 gcrma_2.32.0 genefilter_1.42.0
> [10] affy_1.38.1 Biobase_2.20.1 BiocGenerics_0.6.0
>
> loaded via a namespace (and not attached):
> [1] affyio_1.28.0 annotate_1.38.0 AnnotationDbi_1.22.6
> [4] BiocInstaller_1.10.3 Biostrings_2.28.0 DBI_0.2-7
> [7] grid_3.0.1 IRanges_1.18.3 RSQLite_0.11.4
> [10] splines_3.0.1 stats4_3.0.1 survival_2.37-4
> [13] XML_3.98-1.1 xtable_1.7-1 zlibbioc_1.6.0
>
______________________________________________________________________
The information in this email is confidential and intend...{{dropped:6}}
More information about the Bioconductor
mailing list