Search
Question: design matrix paired data
0
gravatar for Guest User
3.5 years ago by
Guest User12k
Guest User12k 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.
ADD COMMENTlink modified 3.5 years ago by James W. MacDonald44k • written 3.5 years ago by Guest User12k
0
gravatar for James W. MacDonald
3.5 years ago by
United States
James W. MacDonald44k wrote:
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 COMMENTlink written 3.5 years ago by James W. MacDonald44k
0
gravatar for Bernd Klaus
3.5 years ago by
Bernd Klaus460
Germany
Bernd Klaus460 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 COMMENTlink written 3.5 years ago by Bernd Klaus460
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 REPLYlink written 3.5 years ago by Ninni Nahm50
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 REPLYlink written 3.5 years ago by Bernd Klaus460
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 342 users visited in the last hour