[BioC] Anova; can I use lmFit and eBayes from Limma, or should it be aov?

[john herbert]

Dear Bioconductors,
This may be more a question of statistics understanding but it is
related to Limma.

Is it possible to perform an Anova using functions from the Limma
package? If so, how can a post hoc test be applied to the eBayes
"MArrayLM" resulting object?

For me, it is a complex design I have array data, different clients,
time course, stage of disease and treatment. The Limma manual eludes
to Limma being equivalent to Anova in some cases for single channel
arrays (on pages 35 and 41). I am using single channel data.

This code runs;

dd <- data.frame(patient=targets$patient,stage=as.factor(targets$stage),treat=as.factor(targets$treat),time=targets$time)

head(dd)

  patient stage treat time
1       1     A     Y   24
2       1     A     Y   72
3       2     B     N   24
4       2     B     N   72
5       2     B     N    0
6       2     B     Y   72

design <- model.matrix(~patient+stage*treat*time,data=dd)

fit <- lmFit(temp,design)
fit2 <- eBayes(fit)

topTable(fit2)
                       ID X.Intercept.       patient      stageB
stageC       treatY          time stageB.treatY stageC.treatY
stageB.time   stageC.time
A_23_P113716 A_23_P113716     17.57842 -0.0042698998 -0.18798801
-0.18726778 -0.008996453  1.282717e-03    0.38532184    0.18501711
0.0025256233  0.0005692329
A_23_P77779   A_23_P77779     17.07476  0.0015734113 -0.14199010
-0.07687905 -0.132257119  8.288697e-05    0.04068160    0.03375537
0.0024200711  0.0014022044

               treatY.time stageB.treatY.time stageC.treatY.time
AveExpr         F       P.Value     adj.P.Val
A_23_P113716 -1.340200e-03      -0.0072744604      -0.0022265946
17.49773 144347.32 6.481891e-159 2.657575e-154
A_23_P77779  -8.857291e-04       0.0008166484      -0.0004867999
17.01806  98921.73 8.930817e-153 1.830817e-148


How do I interpret the results, if I am performing Anova in Limma
correctly? What are the numbers for X.intercept and stageB etc?

If I take all probes into a variable with
top_data <- topTable(fit2,n=200000000)
and look for ones that are not significant, I don't find any.
which(top_data$adj.P.Val > 0.01)
integer(0)

So it looks like every single probe varies significantly between any
of the conditions, which my instincts say is not correct.

Maybe I should run aov like this instead? But there is still a problem
of post hoc testing as

aof2 <- function(x,y){
  aov(x ~ stage*treat*time+Error(patient), y)
}
anovaresults <- apply(temp, 1, aof2,dd)

lapply(anovaresults,TukeyHSD)
Error in UseMethod("TukeyHSD") :
  no applicable method for 'TukeyHSD' applied to an object of class
"c('aovlist', 'listof')"

If anyone knows how to solve either of these issues, it will be appreciated.

Thanks,

John.