design matrix paired data
2
0
Entering edit mode
Guest User ★ 13k
@guest-user-4897
Last seen 10.3 years ago
Dear all, I made a design matrix according to the limma user guide. However, I wanted to check with you, whether it is correct or not. I feel like I made a mistake, but I do not know for sure. I have 3 patients and 2 conditions: df <- data.frame(patient = c("p1","p2","p3", "p1","p2","p3"), treatment = c("control", "control", "control", "treatment", "treatment", "treatment") ) df patient treatment 1 p1 control 2 p2 control 3 p3 control 4 p1 treatment 5 p2 treatment 6 p3 treatment design <- model.matrix(~factor(patient) + factor(treatment, levels = c("control","treatment")), data= df) design (Intercept) factor(patient)p2 factor(patient)p3 1 1 0 0 2 1 1 0 3 1 0 1 4 1 0 0 5 1 1 0 6 1 0 1 factor(treatment, levels = c("control", "treatment"))treatment 1 0 2 0 3 0 4 1 5 1 6 1 attr(,"assign") [1] 0 1 1 2 attr(,"contrasts") attr(,"contrasts")$`factor(patient)` [1] "contr.treatment" attr(,"contrasts")$`factor(treatment, levels = c("control", "treatment"))` [1] "contr.treatment" colnames(design) <- c("intercept", "patient1", "patient2", "treatment") fit <- lmFit(eset, design) fit <- eBayes(fitNoAdj) tt <- topTable(fitNoAdj, coef= "treatment", adjust="fdr") Is this correct? Thank you in advance! Ninni -- output of sessionInfo(): > sessionInfo() R version 3.1.0 (2014-04-10) Platform: x86_64-pc-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] parallel stats graphics grDevices utils datasets methods [8] base other attached packages: [1] annotate_1.42.0 pd.hugene.2.0.st_3.8.1 [3] oligo_1.28.0 Biostrings_2.32.0 [5] XVector_0.4.0 IRanges_1.22.6 [7] oligoClasses_1.26.0 RColorBrewer_1.0-5 [9] xtable_1.7-3 hugene20sttranscriptcluster.db_2.14.0 [11] org.Hs.eg.db_2.14.0 RSQLite_0.11.4 [13] DBI_0.2-7 AnnotationDbi_1.26.0 [15] GenomeInfoDb_1.0.2 Biobase_2.24.0 [17] BiocGenerics_0.10.0 limma_3.20.1 loaded via a namespace (and not attached): [1] affxparser_1.36.0 affyio_1.32.0 BiocInstaller_1.14.2 [4] bit_1.1-12 codetools_0.2-8 ff_2.2-13 [7] foreach_1.4.2 GenomicRanges_1.16.2 iterators_1.0.7 [10] preprocessCore_1.26.0 splines_3.1.0 stats4_3.1.0 [13] XML_3.98-1.1 zlibbioc_1.10.0 -- Sent via the guest posting facility at bioconductor.org.
limma limma • 2.3k views
ADD COMMENT •
0
Entering edit mode
@james-w-macdonald-5106
Last seen 24 minutes ago
United States
Hi Ninni, On 5/8/2014 4:18 AM, Ninni Nahm [guest] wrote: > Dear all, > > I made a design matrix according to the limma user guide. However, I wanted to check with you, whether it is correct or not. I feel like I made a mistake, but I do not know for sure. > > > I have 3 patients and 2 conditions: > > > df <- data.frame(patient = c("p1","p2","p3", "p1","p2","p3"), > treatment = c("control", "control", "control", "treatment", "treatment", "treatment") > ) > > df > patient treatment > 1 p1 control > 2 p2 control > 3 p3 control > 4 p1 treatment > 5 p2 treatment > 6 p3 treatment > > > design <- model.matrix(~factor(patient) + factor(treatment, levels = c("control","treatment")), data= df) > > > > > design > (Intercept) factor(patient)p2 factor(patient)p3 > 1 1 0 0 > 2 1 1 0 > 3 1 0 1 > 4 1 0 0 > 5 1 1 0 > 6 1 0 1 > factor(treatment, levels = c("control", "treatment"))treatment > 1 0 > 2 0 > 3 0 > 4 1 > 5 1 > 6 1 > attr(,"assign") > [1] 0 1 1 2 > attr(,"contrasts") > attr(,"contrasts")$`factor(patient)` > [1] "contr.treatment" > > attr(,"contrasts")$`factor(treatment, levels = c("control", "treatment"))` > [1] "contr.treatment" > > colnames(design) <- c("intercept", "patient1", "patient2", "treatment") > > > fit <- lmFit(eset, design) > fit <- eBayes(fitNoAdj) > tt <- topTable(fitNoAdj, coef= "treatment", adjust="fdr") > > > Is this correct? Yes. Best, Jim > > Thank you in advance! > > Ninni > > > > -- output of sessionInfo(): > >> sessionInfo() > R version 3.1.0 (2014-04-10) > Platform: x86_64-pc-linux-gnu (64-bit) > > locale: > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 > [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 > [7] LC_PAPER=en_US.UTF-8 LC_NAME=C > [9] LC_ADDRESS=C LC_TELEPHONE=C > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C > > attached base packages: > [1] parallel stats graphics grDevices utils datasets methods > [8] base > > other attached packages: > [1] annotate_1.42.0 pd.hugene.2.0.st_3.8.1 > [3] oligo_1.28.0 Biostrings_2.32.0 > [5] XVector_0.4.0 IRanges_1.22.6 > [7] oligoClasses_1.26.0 RColorBrewer_1.0-5 > [9] xtable_1.7-3 hugene20sttranscriptcluster.db_2.14.0 > [11] org.Hs.eg.db_2.14.0 RSQLite_0.11.4 > [13] DBI_0.2-7 AnnotationDbi_1.26.0 > [15] GenomeInfoDb_1.0.2 Biobase_2.24.0 > [17] BiocGenerics_0.10.0 limma_3.20.1 > > loaded via a namespace (and not attached): > [1] affxparser_1.36.0 affyio_1.32.0 BiocInstaller_1.14.2 > [4] bit_1.1-12 codetools_0.2-8 ff_2.2-13 > [7] foreach_1.4.2 GenomicRanges_1.16.2 iterators_1.0.7 > [10] preprocessCore_1.26.0 splines_3.1.0 stats4_3.1.0 > [13] XML_3.98-1.1 zlibbioc_1.10.0 > > -- > Sent via the guest posting facility at bioconductor.org. > > _______________________________________________ > 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 -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099
ADD COMMENT •
0
Entering edit mode
Bernd Klaus ▴ 610
@bernd-klaus-6281
Last seen 6.1 years ago
Germany
Hi Ninni, looks good! However I would suggest removing the "clutter" in front of the design matrix column names. The command tt <- topTable(fitNoAdj, coef= "treatment", adjust="fdr") will probably not work since the last column is called "factor(treatment, levels = c("control", "treatment"))treatment" and not "treatment". Another reason to rename the column names of the design matrix ... Best wishes, Bernd On Thu, 8 May 2014 01:18:36 -0700 (PDT) "Ninni Nahm \[guest\]" <guest at="" bioconductor.org=""> wrote: > Dear all, > > I made a design matrix according to the limma user guide. However, I wanted to check with you, whether it is correct or not. I feel like I made a mistake, but I do not know for sure. > > > I have 3 patients and 2 conditions: > > > df <- data.frame(patient = c("p1","p2","p3", "p1","p2","p3"), > treatment = c("control", "control", "control", "treatment", "treatment", "treatment") > ) > > df > patient treatment > 1 p1 control > 2 p2 control > 3 p3 control > 4 p1 treatment > 5 p2 treatment > 6 p3 treatment > > > design <- model.matrix(~factor(patient) + factor(treatment, levels = c("control","treatment")), data= df) > > > > > design > (Intercept) factor(patient)p2 factor(patient)p3 > 1 1 0 0 > 2 1 1 0 > 3 1 0 1 > 4 1 0 0 > 5 1 1 0 > 6 1 0 1 > factor(treatment, levels = c("control", "treatment"))treatment > 1 0 > 2 0 > 3 0 > 4 1 > 5 1 > 6 1 > attr(,"assign") > [1] 0 1 1 2 > attr(,"contrasts") > attr(,"contrasts")$`factor(patient)` > [1] "contr.treatment" > > attr(,"contrasts")$`factor(treatment, levels = c("control", "treatment"))` > [1] "contr.treatment" > > colnames(design) <- c("intercept", "patient1", "patient2", "treatment") > > > fit <- lmFit(eset, design) > fit <- eBayes(fitNoAdj) > tt <- topTable(fitNoAdj, coef= "treatment", adjust="fdr") > > > Is this correct? > > Thank you in advance! > > Ninni > > > > -- output of sessionInfo(): > > > sessionInfo() > R version 3.1.0 (2014-04-10) > Platform: x86_64-pc-linux-gnu (64-bit) > > locale: > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 > [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 > [7] LC_PAPER=en_US.UTF-8 LC_NAME=C > [9] LC_ADDRESS=C LC_TELEPHONE=C > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C > > attached base packages: > [1] parallel stats graphics grDevices utils datasets methods > [8] base > > other attached packages: > [1] annotate_1.42.0 pd.hugene.2.0.st_3.8.1 > [3] oligo_1.28.0 Biostrings_2.32.0 > [5] XVector_0.4.0 IRanges_1.22.6 > [7] oligoClasses_1.26.0 RColorBrewer_1.0-5 > [9] xtable_1.7-3 hugene20sttranscriptcluster.db_2.14.0 > [11] org.Hs.eg.db_2.14.0 RSQLite_0.11.4 > [13] DBI_0.2-7 AnnotationDbi_1.26.0 > [15] GenomeInfoDb_1.0.2 Biobase_2.24.0 > [17] BiocGenerics_0.10.0 limma_3.20.1 > > loaded via a namespace (and not attached): > [1] affxparser_1.36.0 affyio_1.32.0 BiocInstaller_1.14.2 > [4] bit_1.1-12 codetools_0.2-8 ff_2.2-13 > [7] foreach_1.4.2 GenomicRanges_1.16.2 iterators_1.0.7 > [10] preprocessCore_1.26.0 splines_3.1.0 stats4_3.1.0 > [13] XML_3.98-1.1 zlibbioc_1.10.0 > > -- > Sent via the guest posting facility at bioconductor.org. > > _______________________________________________ > 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
ADD COMMENT •
0
Entering edit mode
Hi Bernd! Thanks again for your help! Maybe you did not see it, but I have this line here: colnames(design) <- c("intercept", "patient1", "patient2", "treatment") to rename the column names, so the tt <- topTable(fitNoAdj, coef= "treatment", adjust="fdr") command should work :) (unless I misunderstood you) Thank you :) Ninni On Thu, May 8, 2014 at 11:03 AM, Bernd Klaus <bernd.klaus@embl.de> wrote: > Hi Ninni, > > looks good! > > However I would suggest removing the > "clutter" in front of the design matrix column > names. > > The command > > tt <- topTable(fitNoAdj, coef= "treatment", adjust="fdr") > > will probably not work since the last column is called > > "factor(treatment, levels = c("control", "treatment"))treatment" > > and not > > "treatment". Another reason to rename the column names of > the design matrix ... > > Best wishes, > > Bernd > > On Thu, 8 May 2014 01:18:36 -0700 (PDT) > "Ninni Nahm \[guest\]" <guest@bioconductor.org> wrote: > > > Dear all, > > > > I made a design matrix according to the limma user guide. However, I > wanted to check with you, whether it is correct or not. I feel like I made > a mistake, but I do not know for sure. > > > > > > I have 3 patients and 2 conditions: > > > > > > df <- data.frame(patient = c("p1","p2","p3", "p1","p2","p3"), > > treatment = c("control", "control", "control", > "treatment", "treatment", "treatment") > > ) > > > > df > > patient treatment > > 1 p1 control > > 2 p2 control > > 3 p3 control > > 4 p1 treatment > > 5 p2 treatment > > 6 p3 treatment > > > > > > design <- model.matrix(~factor(patient) + factor(treatment, levels = > c("control","treatment")), data= df) > > > > > > > > > > design > > (Intercept) factor(patient)p2 factor(patient)p3 > > 1 1 0 0 > > 2 1 1 0 > > 3 1 0 1 > > 4 1 0 0 > > 5 1 1 0 > > 6 1 0 1 > > factor(treatment, levels = c("control", "treatment"))treatment > > 1 0 > > 2 0 > > 3 0 > > 4 1 > > 5 1 > > 6 1 > > attr(,"assign") > > [1] 0 1 1 2 > > attr(,"contrasts") > > attr(,"contrasts")$`factor(patient)` > > [1] "contr.treatment" > > > > attr(,"contrasts")$`factor(treatment, levels = c("control", > "treatment"))` > > [1] "contr.treatment" > > > > colnames(design) <- c("intercept", "patient1", "patient2", "treatment") > > > > > > fit <- lmFit(eset, design) > > fit <- eBayes(fitNoAdj) > > tt <- topTable(fitNoAdj, coef= "treatment", adjust="fdr") > > > > > > Is this correct? > > > > Thank you in advance! > > > > Ninni > > > > > > > > -- output of sessionInfo(): > > > > > sessionInfo() > > R version 3.1.0 (2014-04-10) > > Platform: x86_64-pc-linux-gnu (64-bit) > > > > locale: > > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 > > [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 > > [7] LC_PAPER=en_US.UTF-8 LC_NAME=C > > [9] LC_ADDRESS=C LC_TELEPHONE=C > > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C > > > > attached base packages: > > [1] parallel stats graphics grDevices utils datasets methods > > [8] base > > > > other attached packages: > > [1] annotate_1.42.0 pd.hugene.2.0.st_3.8.1 > > [3] oligo_1.28.0 Biostrings_2.32.0 > > [5] XVector_0.4.0 IRanges_1.22.6 > > [7] oligoClasses_1.26.0 RColorBrewer_1.0-5 > > [9] xtable_1.7-3 > hugene20sttranscriptcluster.db_2.14.0 > > [11] org.Hs.eg.db_2.14.0 RSQLite_0.11.4 > > [13] DBI_0.2-7 AnnotationDbi_1.26.0 > > [15] GenomeInfoDb_1.0.2 Biobase_2.24.0 > > [17] BiocGenerics_0.10.0 limma_3.20.1 > > > > loaded via a namespace (and not attached): > > [1] affxparser_1.36.0 affyio_1.32.0 BiocInstaller_1.14.2 > > [4] bit_1.1-12 codetools_0.2-8 ff_2.2-13 > > [7] foreach_1.4.2 GenomicRanges_1.16.2 iterators_1.0.7 > > [10] preprocessCore_1.26.0 splines_3.1.0 stats4_3.1.0 > > [13] XML_3.98-1.1 zlibbioc_1.10.0 > > > > -- > > Sent via the guest posting facility at bioconductor.org. > > > > _______________________________________________ > > Bioconductor mailing list > > Bioconductor@r-project.org > > https://stat.ethz.ch/mailman/listinfo/bioconductor > > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor > [[alternative HTML version deleted]]
ADD REPLY •
0
Entering edit mode
Hi Ninni, yes I overlooked that line ;-) Bernd On Thu, 8 May 2014 13:38:18 +0200 Ninni Nahm <ninninahm at="" gmail.com=""> wrote: > Hi Bernd! > > Thanks again for your help! > Maybe you did not see it, but I have this line here: > > colnames(design) <- c("intercept", "patient1", "patient2", "treatment") > > to rename the column names, so the > > tt <- topTable(fitNoAdj, coef= "treatment", adjust="fdr") > > command should work :) (unless I misunderstood you) > > Thank you :) > > Ninni > > > > On Thu, May 8, 2014 at 11:03 AM, Bernd Klaus <bernd.klaus at="" embl.de=""> wrote: > > > Hi Ninni, > > > > looks good! > > > > However I would suggest removing the > > "clutter" in front of the design matrix column > > names. > > > > The command > > > > tt <- topTable(fitNoAdj, coef= "treatment", adjust="fdr") > > > > will probably not work since the last column is called > > > > "factor(treatment, levels = c("control", "treatment"))treatment" > > > > and not > > > > "treatment". Another reason to rename the column names of > > the design matrix ... > > > > Best wishes, > > > > Bernd > > > > On Thu, 8 May 2014 01:18:36 -0700 (PDT) > > "Ninni Nahm \[guest\]" <guest at="" bioconductor.org=""> wrote: > > > > > Dear all, > > > > > > I made a design matrix according to the limma user guide. However, I > > wanted to check with you, whether it is correct or not. I feel like I made > > a mistake, but I do not know for sure. > > > > > > > > > I have 3 patients and 2 conditions: > > > > > > > > > df <- data.frame(patient = c("p1","p2","p3", "p1","p2","p3"), > > > treatment = c("control", "control", "control", > > "treatment", "treatment", "treatment") > > > ) > > > > > > df > > > patient treatment > > > 1 p1 control > > > 2 p2 control > > > 3 p3 control > > > 4 p1 treatment > > > 5 p2 treatment > > > 6 p3 treatment > > > > > > > > > design <- model.matrix(~factor(patient) + factor(treatment, levels = > > c("control","treatment")), data= df) > > > > > > > > > > > > > > > design > > > (Intercept) factor(patient)p2 factor(patient)p3 > > > 1 1 0 0 > > > 2 1 1 0 > > > 3 1 0 1 > > > 4 1 0 0 > > > 5 1 1 0 > > > 6 1 0 1 > > > factor(treatment, levels = c("control", "treatment"))treatment > > > 1 0 > > > 2 0 > > > 3 0 > > > 4 1 > > > 5 1 > > > 6 1 > > > attr(,"assign") > > > [1] 0 1 1 2 > > > attr(,"contrasts") > > > attr(,"contrasts")$`factor(patient)` > > > [1] "contr.treatment" > > > > > > attr(,"contrasts")$`factor(treatment, levels = c("control", > > "treatment"))` > > > [1] "contr.treatment" > > > > > > colnames(design) <- c("intercept", "patient1", "patient2", "treatment") > > > > > > > > > fit <- lmFit(eset, design) > > > fit <- eBayes(fitNoAdj) > > > tt <- topTable(fitNoAdj, coef= "treatment", adjust="fdr") > > > > > > > > > Is this correct? > > > > > > Thank you in advance! > > > > > > Ninni > > > > > > > > > > > > -- output of sessionInfo(): > > > > > > > sessionInfo() > > > R version 3.1.0 (2014-04-10) > > > Platform: x86_64-pc-linux-gnu (64-bit) > > > > > > locale: > > > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > > > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 > > > [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 > > > [7] LC_PAPER=en_US.UTF-8 LC_NAME=C > > > [9] LC_ADDRESS=C LC_TELEPHONE=C > > > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C > > > > > > attached base packages: > > > [1] parallel stats graphics grDevices utils datasets methods > > > [8] base > > > > > > other attached packages: > > > [1] annotate_1.42.0 pd.hugene.2.0.st_3.8.1 > > > [3] oligo_1.28.0 Biostrings_2.32.0 > > > [5] XVector_0.4.0 IRanges_1.22.6 > > > [7] oligoClasses_1.26.0 RColorBrewer_1.0-5 > > > [9] xtable_1.7-3 > > hugene20sttranscriptcluster.db_2.14.0 > > > [11] org.Hs.eg.db_2.14.0 RSQLite_0.11.4 > > > [13] DBI_0.2-7 AnnotationDbi_1.26.0 > > > [15] GenomeInfoDb_1.0.2 Biobase_2.24.0 > > > [17] BiocGenerics_0.10.0 limma_3.20.1 > > > > > > loaded via a namespace (and not attached): > > > [1] affxparser_1.36.0 affyio_1.32.0 BiocInstaller_1.14.2 > > > [4] bit_1.1-12 codetools_0.2-8 ff_2.2-13 > > > [7] foreach_1.4.2 GenomicRanges_1.16.2 iterators_1.0.7 > > > [10] preprocessCore_1.26.0 splines_3.1.0 stats4_3.1.0 > > > [13] XML_3.98-1.1 zlibbioc_1.10.0 > > > > > > -- > > > Sent via the guest posting facility at bioconductor.org. > > > > > > _______________________________________________ > > > 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 > >
ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 770 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6