Search
Question: design matrix paired data
0
4.0 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.
modified 4.0 years ago by James W. MacDonald46k • written 4.0 years ago by Guest User12k
0
4.0 years ago by
United States
James W. MacDonald46k 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
0
4.0 years ago by
Bernd Klaus480
Germany
Bernd Klaus480 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
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]]
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 > >