[BioC] Can I input ordinal variables into a model in Limma?
Gordon K Smyth
smyth at wehi.EDU.AU
Fri Mar 14 01:36:48 CET 2014
Dear Scott,
The linear term is interpreted as linear trend. It is on a standardized
scale, so that the log2-fold-change between the severe and mild levels is
sqrt(0.5) times the reported coefficient (logFC) for the linear term.
Best wishes
Gordon
> Date: Wed, 12 Mar 2014 15:04:11 +0000
> From: Scott Robinson <Scott.Robinson at glasgow.ac.uk>
> To: Gordon K Smyth <smyth at wehi.EDU.AU>
> Cc: Bioconductor mailing list <bioconductor at r-project.org>
> Subject: Re: [BioC] Can I input ordinal variables into a model in
> Limma?
>
> 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:4}}
More information about the Bioconductor
mailing list