Modelling interaction using limma package
1
0
Entering edit mode
Fangxin Hong ▴ 810
@fangxin-hong-912
Last seen 10.2 years ago
Does anyone know how to model interaction (two-way ANOVA) using limma package (lmFit)? Also, it seems that eBayes can be applied to other model fit object, like 'lm.series', 'gls.series' and 'rlm.series' beside 'lmFit'. Does anyone use this before? Thanks a lot! I am quite new to bioconductor... Fangxin -- Fangxin Hong, Ph.D. Bioinformatics Specialist Plant Biology Laboratory The Salk Institute 10010 N. Torrey Pines Rd. La Jolla, CA 92037 E-mail: fhong@salk.edu
• 1.1k views
ADD COMMENT
0
Entering edit mode
Ramon Diaz ★ 1.1k
@ramon-diaz-159
Last seen 10.2 years ago
Dear Fangxin, On Tuesday 14 September 2004 19:41, Fangxin Hong wrote: > Does anyone know how to model interaction (two-way ANOVA) using limma > package (lmFit)? You need to set up the design matrix and decide what contrasts you are interested in. At the bottom there is a somewhat long example (you'll have to supply your own data; this is from a script I have); I have added a few checks to make sure I was testing what I wanted to test. Please note that, for some reason I don't recall now, I made some minor modification to classifyTestsF because it was more convenient for my purposes. [Disclaimer: my example is probably overly complicated]. > Also, it seems that eBayes can be applied to other model fit object, like > 'lm.series', 'gls.series' and 'rlm.series' beside 'lmFit'. Does anyone > use this before? Yes. I've used it with lm and rlm, and maybe also with gls. No major differences. But of course, you might want to make sure that using gls or rlm (or lm) makes sense in your case; for instance, the MASS book (by Venables & Ripley) I think contains some cautionary points about using rlm and related methods with too small sample sizes. Hope this helps. ****************************** The example; I can email you privately this as an attachement. ### treat is a factor, and age:centered is continuous. options(contrasts = c("contr.sum", "contr.poly")) d.trt.BY.age <- model.matrix(~ -1 + treat + treat:age.centered) colnames(d.trt.BY.age) <- c("Colon", "Mama", "Normal", "Colon.by.age", "Mama.by.age", "Normal.by.age") contrasts.d.trt.BY.age <- makeContrasts(Colon - Normal, Mama - Normal, Mama - Colon, 0.5*Colon + 0.5* Mama - Normal, ## testing "cancer" vs. normal Colon.by.age - Normal.by.age, Mama.by.age - Normal.by.age, Mama.by.age - Colon.by.age, levels = d.trt.BY.age) lm.trt.BY.age <- lm.series(d.clean, d.trt.BY.age) lm.trt.BY.age.eb <- ebayes(lm.trt.BY.age) lm.trt.BY.age.contrasts <- contrasts.fit(lm.trt.BY.age, contrasts.d.trt.BY.age) lm.trt.BY.age.contrasts.eb <- ebayes(lm.trt.BY.age.contrasts) d.trt.BY.age2 <- model.matrix(~ -1 + treat + age.centered + treat*age.centered) colnames(d.trt.BY.age2) <- c("Colon", "Mama", "Normal", "age", "Colon.by.age", "Mama.by.age") contrasts.d.trt.BY.age2 <- makeContrasts(Colon.by.age, Mama.by.age, levels = d.trt.BY.age2) #### Get the interaction F-statistic tmp <- classifyTestsF(lm.trt.BY.age2.contrasts.eb$t, contrasts = contrasts.d.trt.BY.age2, design = d.trt.BY.age2, df = lm.trt.BY.age2$df + lm.trt.BY.age2.contrasts.eb$df, fstat.only = TRUE) lm.interaction.Fstat <- cbind(F.stat =tmp$Fstatistic, df1 = rep(2, length(tmp$Fstatistic)), df2 = tmp$df2) lm.interaction.Fstat <- cbind(lm.interaction.Fstat, p.value = pf(lm.interaction.Fstat[,1], lm.interaction.Fstat[,2], lm.interaction.Fstat[,3], lower.tail = FALSE)) lm.interaction.Fstat <- cbind(lm.interaction.Fstat, fdr.p.value = p.adjust(lm.interaction.Fstat[, 4], method = "fdr")) summary(lm.interaction.Fstat) ############# Checks ### check appropriate contrasts testing for testing interaction via calssifyTestsF l0 <- lm(d.clean[2, ] ~ -1 + treat + age.centered) l1 <- lm(d.clean[2, ] ~ -1 + treat + age.centered + age.centered*treat) l2 <- lm(d.clean[2, ] ~ -1 + treat + treat:age.centered) anova(l0) anova(l1) anova(l2) anova(l0, l1) ## this and the following are the same anova(l0, l2) anova(l1, l2) ## nothing, as expected ## just for checking!!! Note I modify ts, etc. lm.trt.BY.age2.contrasts.eb$t[2,] <- summary(l1)$coefficients[5:6, 3] classifyTestsF(lm.trt.BY.age2.contrasts.eb$t[2, ], contrasts = contrasts.d.trt.BY.age2[2, ], design = d.trt.BY.age2, df = lm.trt.BY.age2$df[2, ] + lm.trt.BY.age2.contrasts.eb$df, fstat.only = TRUE) ### OK, works as it should ### changes to classifyTestsF ## I change the return value when fstat.only = TRUE classifyTestsF <- function (object, cor.matrix = NULL, design = NULL, contrasts = diag(ncol(design)), df = Inf, p.value = 0.01, fstat.only = FALSE) { if (is.list(object)) { if (is.null(object$t)) stop("tstat cannot be extracted from object") if (missing(design) && !is.null(object$design)) design <- object$design if (missing(contrasts) && !is.null(object$contrasts)) contrasts <- object$contrasts if (missing(df) && !is.null(object$df.prior) && !is.null(object$df.residual)) df <- object$df.prior + object$df.residual tstat <- as.matrix(object$t) } else { tstat <- as.matrix(object) } ngenes <- nrow(tstat) ntests <- ncol(tstat) if (ntests == 1) { qT <- qt(p.value/2, df, lower.tail = FALSE) return(sign(tstat) * (abs(tstat) > qT)) } if (!is.null(cor.matrix) && !is.null(design)) stop("Cannot specify both cor.matrix and design") if (!is.null(design)) { design <- as.matrix(design) R <- chol(crossprod(design)) cor.matrix <- crossprod(backsolve(R, contrasts, transpose = TRUE)) d <- sqrt(diag(cor.matrix)) cor.matrix <- cor.matrix/(d %*% t(d)) } if (is.null(cor.matrix)) { r <- ntests Q <- diag(r)/sqrt(r) } else { E <- eigen(cor.matrix, symmetric = TRUE) r <- sum(E$values/E$values[1] > 1e-08) Q <- matvec(E$vectors[, 1:r], 1/sqrt(E$values[1:r]))/sqrt(r) } if (fstat.only) { fstat <- list() fstat$Fstatistic <- drop((tstat %*% Q)^2 %*% array(1, c(r, 1))) fstat$df1 <- r fstat$df2 <- df return(fstat) } qF <- qf(p.value, r, df, lower.tail = FALSE) if (length(qF) == 1) qF <- rep(qF, ngenes) result <- matrix(0, ngenes, ntests, dimnames = dimnames(tstat)) if (is.null(colnames(tstat)) && !is.null(colnames(contrasts))) colnames(result) <- colnames(contrasts) for (i in 1:ngenes) { x <- tstat[i, ] if (anyis.na(x))) result[i, ] <- NA else if (crossprod(crossprod(Q, x)) > qF[i]) { ord <- order(abs(x), decreasing = TRUE) result[i, ord[1]] <- sign(x[ord[1]]) for (j in 2:ntests) { bigger <- ord[1:(j - 1)] x[bigger] <- sign(x[bigger]) * abs(x[ord[j]]) if (crossprod(crossprod(Q, x)) > qF[i]) result[i, ord[j]] <- sign(x[ord[j]]) else break } } } new("TestResults", result) } -- Ram?n D?az-Uriarte Bioinformatics Unit Centro Nacional de Investigaciones Oncol?gicas (CNIO) (Spanish National Cancer Center) Melchor Fern?ndez Almagro, 3 28029 Madrid (Spain) Fax: +-34-91-224-6972 Phone: +-34-91-224-6900 http://ligarto.org/rdiaz PGP KeyID: 0xE89B3462 (http://ligarto.org/rdiaz/0xE89B3462.asc)
ADD COMMENT
0
Entering edit mode
Dear Ramon; Thank you very much for the detailed answer of my question, it works for me!!! One more question, for one factor, if different individual (gene) has different number of levels (thus different design and constrast matrix). Theoretically, this is still doable using the proposed linear model with Emprical Bayes method according to the paper, but in practice, do you or anybdy on this list know how to do it? what I am thinking is to divide the whole data set into several groups, all individuals in one group has same number of levels, then apply limma package (lmFit and eBayes) as usual. Although this changes the original assumption that variance from all individuals are from the same distribution with the one that variance from individuals in each group are from the same distribution, it should still be valid analysis as each group will still have large amount of genes. Does anyone have better idea beside writing my own code of doing eBayes. Thanks. Fangxin > Dear Fangxin, > > On Tuesday 14 September 2004 19:41, Fangxin Hong wrote: >> Does anyone know how to model interaction (two-way ANOVA) using limma >> package (lmFit)? > > You need to set up the design matrix and decide what contrasts you are > interested in. At the bottom there is a somewhat long example (you'll have > to > supply your own data; this is from a script I have); I have added a few > checks to make sure I was testing what I wanted to test. Please note that, > for some reason I don't recall now, I made some minor modification to > classifyTestsF because it was more convenient for my purposes. > [Disclaimer: > my example is probably overly complicated]. > > >> Also, it seems that eBayes can be applied to other model fit object, >> like >> 'lm.series', 'gls.series' and 'rlm.series' beside 'lmFit'. Does anyone >> use this before? > > > Yes. I've used it with lm and rlm, and maybe also with gls. No major > differences. But of course, you might want to make sure that using gls or > rlm > (or lm) makes sense in your case; for instance, the MASS book (by Venables > & > Ripley) I think contains some cautionary points about using rlm and > related > methods with too small sample sizes. > > > Hope this helps. > > > ****************************** > > The example; I can email you privately this as an attachement. > > ### treat is a factor, and age:centered is continuous. > > > options(contrasts = c("contr.sum", "contr.poly")) > > d.trt.BY.age <- model.matrix(~ -1 + treat + treat:age.centered) > colnames(d.trt.BY.age) <- c("Colon", "Mama", "Normal", > "Colon.by.age", "Mama.by.age", "Normal.by.age") > contrasts.d.trt.BY.age <- > makeContrasts(Colon - Normal, > Mama - Normal, > Mama - Colon, > 0.5*Colon + 0.5* Mama - Normal, ## testing "cancer" vs. > normal > Colon.by.age - Normal.by.age, > Mama.by.age - Normal.by.age, > Mama.by.age - Colon.by.age, > levels = d.trt.BY.age) > > > lm.trt.BY.age <- lm.series(d.clean, d.trt.BY.age) > lm.trt.BY.age.eb <- ebayes(lm.trt.BY.age) > lm.trt.BY.age.contrasts <- contrasts.fit(lm.trt.BY.age, > contrasts.d.trt.BY.age) > lm.trt.BY.age.contrasts.eb <- ebayes(lm.trt.BY.age.contrasts) > > > > > d.trt.BY.age2 <- model.matrix(~ -1 + treat + age.centered + > treat*age.centered) > colnames(d.trt.BY.age2) <- c("Colon", "Mama", "Normal", > "age", "Colon.by.age", "Mama.by.age") > contrasts.d.trt.BY.age2 <- makeContrasts(Colon.by.age, > Mama.by.age, > levels = d.trt.BY.age2) > > > > > #### Get the interaction F-statistic > tmp <- > classifyTestsF(lm.trt.BY.age2.contrasts.eb$t, > contrasts = contrasts.d.trt.BY.age2, > design = d.trt.BY.age2, > df = lm.trt.BY.age2$df + > lm.trt.BY.age2.contrasts.eb$df, > fstat.only = TRUE) > lm.interaction.Fstat <- cbind(F.stat =tmp$Fstatistic, > df1 = rep(2, length(tmp$Fstatistic)), > df2 = tmp$df2) > lm.interaction.Fstat <- cbind(lm.interaction.Fstat, > p.value = pf(lm.interaction.Fstat[,1], > lm.interaction.Fstat[,2], > lm.interaction.Fstat[,3], > lower.tail = FALSE)) > lm.interaction.Fstat <- cbind(lm.interaction.Fstat, > fdr.p.value > p.adjust(lm.interaction.Fstat[, 4], > method = "fdr")) > summary(lm.interaction.Fstat) > > > > ############# Checks > > ### check appropriate contrasts testing for testing interaction via > calssifyTestsF > l0 <- lm(d.clean[2, ] ~ -1 + treat + age.centered) > l1 <- lm(d.clean[2, ] ~ -1 + treat + age.centered + age.centered*treat) > l2 <- lm(d.clean[2, ] ~ -1 + treat + treat:age.centered) > > anova(l0) > anova(l1) > anova(l2) > > anova(l0, l1) ## this and the following are the same > anova(l0, l2) > anova(l1, l2) ## nothing, as expected > > ## just for checking!!! Note I modify ts, etc. > lm.trt.BY.age2.contrasts.eb$t[2,] <- summary(l1)$coefficients[5:6, 3] > > classifyTestsF(lm.trt.BY.age2.contrasts.eb$t[2, ], > contrasts = contrasts.d.trt.BY.age2[2, ], > design = d.trt.BY.age2, > df = lm.trt.BY.age2$df[2, ] + > lm.trt.BY.age2.contrasts.eb$df, > fstat.only = TRUE) > ### OK, works as it should > > > > > ### changes to classifyTestsF > > ## I change the return value when fstat.only = TRUE > classifyTestsF <- > function (object, cor.matrix = NULL, design = NULL, > contrasts = diag(ncol(design)), > df = Inf, p.value = 0.01, fstat.only = FALSE) > { > if (is.list(object)) { > if (is.null(object$t)) > stop("tstat cannot be extracted from object") > if (missing(design) && !is.null(object$design)) > design <- object$design > if (missing(contrasts) && !is.null(object$contrasts)) > contrasts <- object$contrasts > if (missing(df) && !is.null(object$df.prior) && > !is.null(object$df.residual)) > df <- object$df.prior + object$df.residual > tstat <- as.matrix(object$t) > } > else { > tstat <- as.matrix(object) > } > ngenes <- nrow(tstat) > ntests <- ncol(tstat) > if (ntests == 1) { > qT <- qt(p.value/2, df, lower.tail = FALSE) > return(sign(tstat) * (abs(tstat) > qT)) > } > if (!is.null(cor.matrix) && !is.null(design)) > stop("Cannot specify both cor.matrix and design") > if (!is.null(design)) { > design <- as.matrix(design) > R <- chol(crossprod(design)) > cor.matrix <- crossprod(backsolve(R, contrasts, transpose = TRUE)) > d <- sqrt(diag(cor.matrix)) > cor.matrix <- cor.matrix/(d %*% t(d)) > } > if (is.null(cor.matrix)) { > r <- ntests > Q <- diag(r)/sqrt(r) > } > else { > E <- eigen(cor.matrix, symmetric = TRUE) > r <- sum(E$values/E$values[1] > 1e-08) > Q <- matvec(E$vectors[, 1:r], 1/sqrt(E$values[1:r]))/sqrt(r) > } > if (fstat.only) { > fstat <- list() > fstat$Fstatistic <- drop((tstat %*% Q)^2 %*% array(1, c(r, 1))) > fstat$df1 <- r > fstat$df2 <- df > return(fstat) > } > qF <- qf(p.value, r, df, lower.tail = FALSE) > if (length(qF) == 1) > qF <- rep(qF, ngenes) > result <- matrix(0, ngenes, ntests, dimnames = dimnames(tstat)) > if (is.null(colnames(tstat)) && !is.null(colnames(contrasts))) > colnames(result) <- colnames(contrasts) > for (i in 1:ngenes) { > x <- tstat[i, ] > if (anyis.na(x))) > result[i, ] <- NA > else if (crossprod(crossprod(Q, x)) > qF[i]) { > ord <- order(abs(x), decreasing = TRUE) > result[i, ord[1]] <- sign(x[ord[1]]) > for (j in 2:ntests) { > bigger <- ord[1:(j - 1)] > x[bigger] <- sign(x[bigger]) * abs(x[ord[j]]) > if (crossprod(crossprod(Q, x)) > qF[i]) > result[i, ord[j]] <- sign(x[ord[j]]) > else break > } > } > } > new("TestResults", result) > } > > > > > > > > -- > Ram?n D?az-Uriarte > Bioinformatics Unit > Centro Nacional de Investigaciones Oncol?gicas (CNIO) > (Spanish National Cancer Center) > Melchor Fern?ndez Almagro, 3 > 28029 Madrid (Spain) > Fax: +-34-91-224-6972 > Phone: +-34-91-224-6900 > > http://ligarto.org/rdiaz > PGP KeyID: 0xE89B3462 > (http://ligarto.org/rdiaz/0xE89B3462.asc) > > > -- Fangxin Hong, Ph.D. Bioinformatics Specialist Plant Biology Laboratory The Salk Institute 10010 N. Torrey Pines Rd. La Jolla, CA 92037 E-mail: fhong@salk.edu
ADD REPLY
0
Entering edit mode
Dear Fangxin, On Tuesday 21 September 2004 18:38, Fangxin Hong wrote: > Dear Ramon; > Thank you very much for the detailed answer of my question, it works > for me!!! I am glad it was helpful. > > One more question, for one factor, if different individual (gene) has > different number of levels (thus different design and constrast matrix). > Theoretically, this is still doable using the proposed linear model with > Emprical Bayes method according to the paper, but in practice, do you or > anybdy on this list know how to do it? > > what I am thinking is to divide the whole data set into several groups, > all individuals in one group has same number of levels, then apply limma > package (lmFit and eBayes) as usual. Although this changes the original > assumption that variance from all individuals are from the same > distribution with the one that variance from individuals in each group are > from the same distribution, it should still be valid analysis as each > group will still have large amount of genes. Does anyone have better idea > beside writing my own code of doing eBayes. I do not know how I'd approach the issue, and you obviously have thought about it more than I. However, I am having a hard time seeing how you get different levels for different genes for the same array. This factor must be "gene specific", but then, does it make sense to try a "single unified analysis" across all genes? It would even seem that the null you are testing is not the same across genes. For instance, if your factor has levels 1, 2, 3, and in gene A you have levels 1 and 2, and in gene B you have levels 2 and 3, then: for A, H_o: \mu_1 = \mu_2, but for B, H_o: \mu_2 = \mu_3. This seems very strange to me... Best, R. > > Thanks. > > Fangxin > > > Dear Fangxin, > > > > On Tuesday 14 September 2004 19:41, Fangxin Hong wrote: > >> Does anyone know how to model interaction (two-way ANOVA) using limma > >> package (lmFit)? > > > > You need to set up the design matrix and decide what contrasts you are > > interested in. At the bottom there is a somewhat long example (you'll > > have to > > supply your own data; this is from a script I have); I have added a few > > checks to make sure I was testing what I wanted to test. Please note > > that, for some reason I don't recall now, I made some minor modification > > to classifyTestsF because it was more convenient for my purposes. > > [Disclaimer: > > my example is probably overly complicated]. > > > >> Also, it seems that eBayes can be applied to other model fit object, > >> like > >> 'lm.series', 'gls.series' and 'rlm.series' beside 'lmFit'. Does anyone > >> use this before? > > > > Yes. I've used it with lm and rlm, and maybe also with gls. No major > > differences. But of course, you might want to make sure that using gls or > > rlm > > (or lm) makes sense in your case; for instance, the MASS book (by > > Venables & > > Ripley) I think contains some cautionary points about using rlm and > > related > > methods with too small sample sizes. > > > > > > Hope this helps. > > > > > > ****************************** > > > > The example; I can email you privately this as an attachement. > > > > ### treat is a factor, and age:centered is continuous. > > > > > > options(contrasts = c("contr.sum", "contr.poly")) > > > > d.trt.BY.age <- model.matrix(~ -1 + treat + treat:age.centered) > > colnames(d.trt.BY.age) <- c("Colon", "Mama", "Normal", > > "Colon.by.age", "Mama.by.age", "Normal.by.age") > > contrasts.d.trt.BY.age <- > > makeContrasts(Colon - Normal, > > Mama - Normal, > > Mama - Colon, > > 0.5*Colon + 0.5* Mama - Normal, ## testing "cancer" vs. > > normal > > Colon.by.age - Normal.by.age, > > Mama.by.age - Normal.by.age, > > Mama.by.age - Colon.by.age, > > levels = d.trt.BY.age) > > > > > > lm.trt.BY.age <- lm.series(d.clean, d.trt.BY.age) > > lm.trt.BY.age.eb <- ebayes(lm.trt.BY.age) > > lm.trt.BY.age.contrasts <- contrasts.fit(lm.trt.BY.age, > > contrasts.d.trt.BY.age) > > lm.trt.BY.age.contrasts.eb <- ebayes(lm.trt.BY.age.contrasts) > > > > > > > > > > d.trt.BY.age2 <- model.matrix(~ -1 + treat + age.centered + > > treat*age.centered) > > colnames(d.trt.BY.age2) <- c("Colon", "Mama", "Normal", > > "age", "Colon.by.age", "Mama.by.age") > > contrasts.d.trt.BY.age2 <- makeContrasts(Colon.by.age, > > Mama.by.age, > > levels = d.trt.BY.age2) > > > > > > > > > > #### Get the interaction F-statistic > > tmp <- > > classifyTestsF(lm.trt.BY.age2.contrasts.eb$t, > > contrasts = contrasts.d.trt.BY.age2, > > design = d.trt.BY.age2, > > df = lm.trt.BY.age2$df + > > lm.trt.BY.age2.contrasts.eb$df, > > fstat.only = TRUE) > > lm.interaction.Fstat <- cbind(F.stat =tmp$Fstatistic, > > df1 = rep(2, length(tmp$Fstatistic)), > > df2 = tmp$df2) > > lm.interaction.Fstat <- cbind(lm.interaction.Fstat, > > p.value = pf(lm.interaction.Fstat[,1], > > lm.interaction.Fstat[,2], > > lm.interaction.Fstat[,3], > > lower.tail = FALSE)) > > lm.interaction.Fstat <- cbind(lm.interaction.Fstat, > > fdr.p.value > > p.adjust(lm.interaction.Fstat[, 4], > > method = "fdr")) > > summary(lm.interaction.Fstat) > > > > > > > > ############# Checks > > > > ### check appropriate contrasts testing for testing interaction via > > calssifyTestsF > > l0 <- lm(d.clean[2, ] ~ -1 + treat + age.centered) > > l1 <- lm(d.clean[2, ] ~ -1 + treat + age.centered + age.centered*treat) > > l2 <- lm(d.clean[2, ] ~ -1 + treat + treat:age.centered) > > > > anova(l0) > > anova(l1) > > anova(l2) > > > > anova(l0, l1) ## this and the following are the same > > anova(l0, l2) > > anova(l1, l2) ## nothing, as expected > > > > ## just for checking!!! Note I modify ts, etc. > > lm.trt.BY.age2.contrasts.eb$t[2,] <- summary(l1)$coefficients[5:6, 3] > > > > classifyTestsF(lm.trt.BY.age2.contrasts.eb$t[2, ], > > contrasts = contrasts.d.trt.BY.age2[2, ], > > design = d.trt.BY.age2, > > df = lm.trt.BY.age2$df[2, ] + > > lm.trt.BY.age2.contrasts.eb$df, > > fstat.only = TRUE) > > ### OK, works as it should > > > > > > > > > > ### changes to classifyTestsF > > > > ## I change the return value when fstat.only = TRUE > > classifyTestsF <- > > function (object, cor.matrix = NULL, design = NULL, > > contrasts = diag(ncol(design)), > > df = Inf, p.value = 0.01, fstat.only = FALSE) > > { > > if (is.list(object)) { > > if (is.null(object$t)) > > stop("tstat cannot be extracted from object") > > if (missing(design) && !is.null(object$design)) > > design <- object$design > > if (missing(contrasts) && !is.null(object$contrasts)) > > contrasts <- object$contrasts > > if (missing(df) && !is.null(object$df.prior) && > > !is.null(object$df.residual)) > > df <- object$df.prior + object$df.residual > > tstat <- as.matrix(object$t) > > } > > else { > > tstat <- as.matrix(object) > > } > > ngenes <- nrow(tstat) > > ntests <- ncol(tstat) > > if (ntests == 1) { > > qT <- qt(p.value/2, df, lower.tail = FALSE) > > return(sign(tstat) * (abs(tstat) > qT)) > > } > > if (!is.null(cor.matrix) && !is.null(design)) > > stop("Cannot specify both cor.matrix and design") > > if (!is.null(design)) { > > design <- as.matrix(design) > > R <- chol(crossprod(design)) > > cor.matrix <- crossprod(backsolve(R, contrasts, transpose = > > TRUE)) d <- sqrt(diag(cor.matrix)) > > cor.matrix <- cor.matrix/(d %*% t(d)) > > } > > if (is.null(cor.matrix)) { > > r <- ntests > > Q <- diag(r)/sqrt(r) > > } > > else { > > E <- eigen(cor.matrix, symmetric = TRUE) > > r <- sum(E$values/E$values[1] > 1e-08) > > Q <- matvec(E$vectors[, 1:r], 1/sqrt(E$values[1:r]))/sqrt(r) > > } > > if (fstat.only) { > > fstat <- list() > > fstat$Fstatistic <- drop((tstat %*% Q)^2 %*% array(1, c(r, 1))) > > fstat$df1 <- r > > fstat$df2 <- df > > return(fstat) > > } > > qF <- qf(p.value, r, df, lower.tail = FALSE) > > if (length(qF) == 1) > > qF <- rep(qF, ngenes) > > result <- matrix(0, ngenes, ntests, dimnames = dimnames(tstat)) > > if (is.null(colnames(tstat)) && !is.null(colnames(contrasts))) > > colnames(result) <- colnames(contrasts) > > for (i in 1:ngenes) { > > x <- tstat[i, ] > > if (anyis.na(x))) > > result[i, ] <- NA > > else if (crossprod(crossprod(Q, x)) > qF[i]) { > > ord <- order(abs(x), decreasing = TRUE) > > result[i, ord[1]] <- sign(x[ord[1]]) > > for (j in 2:ntests) { > > bigger <- ord[1:(j - 1)] > > x[bigger] <- sign(x[bigger]) * abs(x[ord[j]]) > > if (crossprod(crossprod(Q, x)) > qF[i]) > > result[i, ord[j]] <- sign(x[ord[j]]) > > else break > > } > > } > > } > > new("TestResults", result) > > } > > > > > > > > > > > > > > > > -- > > Ram?n D?az-Uriarte > > Bioinformatics Unit > > Centro Nacional de Investigaciones Oncol?gicas (CNIO) > > (Spanish National Cancer Center) > > Melchor Fern?ndez Almagro, 3 > > 28029 Madrid (Spain) > > Fax: +-34-91-224-6972 > > Phone: +-34-91-224-6900 > > > > http://ligarto.org/rdiaz > > PGP KeyID: 0xE89B3462 > > (http://ligarto.org/rdiaz/0xE89B3462.asc) -- Ram?n D?az-Uriarte Bioinformatics Unit Centro Nacional de Investigaciones Oncol?gicas (CNIO) (Spanish National Cancer Center) Melchor Fern?ndez Almagro, 3 28029 Madrid (Spain) Fax: +-34-91-224-6972 Phone: +-34-91-224-6900 http://ligarto.org/rdiaz PGP KeyID: 0xE89B3462 (http://ligarto.org/rdiaz/0xE89B3462.asc)
ADD REPLY

Login before adding your answer.

Traffic: 610 users visited in the last hour
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