I am analyzing a 2-factor factorial Affy experiment, with 3 d.f. for
each
factor.
I would like to get the F-tests for the main effects and interactions
using
limma.
I have computed all the contrasts, and got the t-tests (both
unadjusted and
eBayes). I do know how to combine these into F-tests "by hand" but I
could
not figure out if there was a simple way to do this using limma.
I had a look at FStat (classifyTestsF). There seems to be a problem,
in
that the matrix tstat is not premultiplied by the contrast matrix when
the
F-statistics are computed. So, if the contrasts are not full-rank, an
error is generated (instead of the F-statistics) because nrow(Q) !=
ncol(tstat).. (See the lines below).
if (fstat.only) {
fstat <- drop((tstat%*% Q)^2 %*% array(1, c(r, 1)))
attr(fstat, "df1") <- r
attr(fstat, "df2") <- df
return(fstat)
}
I figured that before I fiddled with the code, I would check to make
sure
that I didn't miss an existing routine to do this.
Thanks in advance.
Naomi S. Altman 814-865-3791 (voice)
Associate Professor
Bioinformatics Consulting Center
Dept. of Statistics 814-863-7114 (fax)
Penn State University 814-865-1348
(Statistics)
University Park, PA 16802-2111
> Date: Tue, 21 Dec 2004 17:21:19 -0500
> From: Naomi Altman <naomi@stat.psu.edu>
> Subject: [BioC] F-tests for factorial effects - limma
> To: bioconductor@stat.math.ethz.ch
>
> I am analyzing a 2-factor factorial Affy experiment, with 3 d.f. for
each
> factor.
>
> I would like to get the F-tests for the main effects and
interactions using
> limma.
>
> I have computed all the contrasts, and got the t-tests (both
unadjusted and
> eBayes). I do know how to combine these into F-tests "by hand" but
I could
> not figure out if there was a simple way to do this using limma.
limma doesn't have any easy way to deal with main effects and
interactions, at least not with main
effects, interactions are actually simpler. I haven't implemented
this because I've never been
able to figure out what one would do with these things in a microarray
context.
To compute F-tests for main effects and interaction, the easiest way
would probably be to compute
the SS for main effects and interactions by non-limma means, then use
shrinkVar() to adjust the
residual mean squares, i.e., the F-statistic denominators.
If you only want F-tests for interactions, the following code would
work:
X <- model.matrix(~a*b)
fit <- lmFit(eset, X)
p <- ncol(X)
cont.ia <- diag(p)[,attr(X,"assign")==3]
fit.ia <- eBayes(contrasts.fit(fit, cont.ia))
Now fit.ia contains the F-statistic and p-values for the interaction
in fit.ia$F and
fit.ia$F.p.value.
> I had a look at FStat (classifyTestsF). There seems to be a
problem, in
> that the matrix tstat is not premultiplied by the contrast matrix
when the
> F-statistics are computed. So, if the contrasts are not full-rank,
an
> error is generated (instead of the F-statistics) because nrow(Q) !=
> ncol(tstat).. (See the lines below).
No, the code is correct. FStat is quite happy with non full rank
contrasts but the contrast
matrix must be applied using contrasts.fit() before entering FStat().
You should not expect to
see a contrast matrix inside the classifyTestsF() code.
Gordon
> if (fstat.only) {
> fstat <- drop((tstat%*% Q)^2 %*% array(1, c(r, 1)))
> attr(fstat, "df1") <- r
> attr(fstat, "df2") <- df
> return(fstat)
> }
>
> I figured that before I fiddled with the code, I would check to make
sure
> that I didn't miss an existing routine to do this.
>
> Thanks in advance.
>
> Naomi S. Altman 814-865-3791 (voice)
> Associate Professor
> Bioinformatics Consulting Center
> Dept. of Statistics 814-863-7114 (fax)
> Penn State University 814-865-1348
(Statistics)
> University Park, PA 16802-2111