Entering edit mode
Thanks, Gordon. We will try contrasts.fit. If that does not give us
what
we want, I will consider writing an aov.series routine.
The microarray experiments I see are often factorial, RCB, or split
plot
designs. Particularly with Affy arrays, they are often balanced. I
guess
I am used to thinking in terms of "effects" rather than individual
contrasts, but why would the main effect and interaction idea be
different
for microarray experiments than for other experiments?
e.g. treatment 1-3
genotype 1-3
The interaction is always of interest. However, even in the absence
of
interaction, the genotype or treatment effect is often of interest.
e.g.
the genotype may be various knockouts - we are interested in knowing
what
other genes have been induced or repressed. Why is the test for main
effects for a particular gene less interesting in this context than
e.g.
the change in weight or other univariate measure?
--Naomi
At 12:23 AM 12/23/2004 +1100, Gordon K Smyth wrote:
> > 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
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
