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)
