[R] Multiple tests on 2 way-ANOVA
Grathwohl, Dominik, LAUSANNE, NRC-BAS
dominik.grathwohl at rdls.nestle.com
Tue Jul 25 18:26:32 CEST 2006
Spencer Graves <spencer.graves <at> pdf.com> writes:
>
> <comments in line>
>
> Grathwohl, Dominik, LAUSANNE, NRC-BAS wrote:
> > Dear r-helpers,
> >
> > I have a question about multiple testing.
> > Here an example that puzzles me:
> > All matrixes and contrast vectors are presented in treatment contrasts.
> >
> > 1. example:
> > library(multcomp)
> > n<-60; sigma<-20
> > # n = sample size per group
> > # sigma standard deviation of the residuals
> >
> > cov1 <- matrix(c(3/4,-1/2,-1/2,-1/2,1,0,-1/2,0,1), nrow = 3, ncol=3, byrow=TRUE,
> > dimnames = list(c("A", "B", "C"), c("C.1", "C.2", "C.3")))
> > # cov1 = variance covariance matrix of the beta coefficients of a
> > # 2x2 factorial design (see Piantadosi 2005, p. 509)
> >
> > cm1 <- matrix(c(0, 1, 0, 0, 0, 1), nrow = 2, ncol=3, byrow=TRUE,
> > dimnames = list(c("A", "B"), c("C.1", "C.2", "C.3")))
> > # cm1 = contrast matrix for main effects
> >
> > v1 <- csimint(estpar=c(100, 6, 5), df=4*n-3, covm=cov1*sigma^2/n, cmatrix=cm1, conf.level=0.95)
> > summary(v1)
> >
> > The adjusted p-values are almost the Bonferroni p-values.
> > If I understood right: You need not to adjust for multiple testing
> > on main effects in a 2x2 factorial design
> > assuming the absence of interaction.
>
> SG: Where did you get this idea? A p value of 0.05 says that if the
> null hypothesis of no effect is true, a result at least as extreme as
> that observed work actually occur with probability 0.05. Thus, with 2
> independent tests, the probability of getting a result at least that
> extreme in one or both of the tests is 1-(1-0.05)^2 = 0.0975, which is
> almost 2*0.05. Thus, if I were to consider only main effects in a 2x2
> factorial design, this is what I would get from Bonferroni.
>
How I get this idea? There are two viewpoints on multiple tests on factorial designs (lets restrict to a 2x2 factorial and absence of interaction):
1.) A factorial design, you are doing two trials for the price of one. In more detail: If you investigate a treatment A (present or absent) you can conduct a parallel group design with two groups. If you investigate treatment B, you need to conduct a second parallel group design. For the two parallel group designs no adjustment would have been considered. The factorial design is randomized in a way that investigating the two treatment effects can be seen as independent of each other like the two parallel group designs, so no adjustment is necessary.
2.) A experiment with a factorial design is still one experiment thus controlling experiment wise alpha error for two treatments need to be corrected. Since the treatments are randomized a way that they can be seen as independent, the Bonferroni correction is appropriate. And exactly this is doing the csimint function.
3.) My revised viewpoint:
A 2x2 factorial design has four groups:
group 1: non A, non B,
group 2: A, non B,
group 3: non A, B
group 4: A, B
For estimating effect of A as well the effect of B, the "non A, non B" group is involved. So strictly spoken the effects are not completely independent estimated. So csimint is doing a good job.
> > I do not think that there is a bug,
> > I want to understand, why multcomp does adjust for multiple tests
> > having all information about the design of the trial (variance covariance matrix)?
> > Or do I have to introduce somehow more information?
> >
> > 2. example:
> > And I have second question: How do I proper correct for multiple testing
> > if I want to estimate in the presence of interaction the two average main effects.
> > Can some one point me to some literature where I can learn these things?
> > Here the example, 2x2 factorial with interaction, estimation of average main effects:
> >
> > cov2 <- matrix(
> > c(1,-1,-1, 1,
> > -1, 2, 1,-2,
> > -1, 1, 2,-2,
> > 1,-2,-2, 4)
> > , nrow=4, ncol=4, byrow=TRUE)
> > cm2 <- matrix(c(0, 1, 0, 1/2, 0, 0, 1, 1/2), nrow = 2, ncol=4, byrow=TRUE,
> > dimnames = list(c("A", "B"), c("C.1", "C.2", "C.3", "C.4")))
> > v2 <- csimint(estpar=c(100, 6, 5, 2), df=4*n-4, covm=cov2*sigma^2/n, cmatrix=cm2, conf.level=0.95)
> > summary(v2)
>
> SG: The Bonferroni p-value is the observed times the number of rows in
> the contrast matrix. The number of columns is irrelevant to Bonferroni.
>
> SG: I'm not sure, but I believe that the the adjusted p value would
> likely be close to (if not exactly) the rank of the contrast matrix;
> given the rank, it is (I think) independent of the number of rows and
> columns.
>
> SG: These two assertions are consistent with the following example,
> where I increase the number of dimensions by a factor of 4 without
> changing the rank. The Bonferroni p value increased by a factor of 4
> while the adjusted p value did not change, as predicted.
>
> cm2.4 <- rbind(cm2, cm2, cm2, cm2)
> v2.4 <- csimint(estpar=c(100, 6, 5, 2), df=4*n-4,
> covm=cov2*sigma^2/n,
> cmatrix=cm2.4, conf.level=0.95)
> summary(v2.4)
> >
> > I do not believe that this is the most efficient way for doing this,
> > since I made already bad experience with the first example.
>
> SG: I hope this reply converts the "bad experience" to "good". As for
> efficiency, you did very well by including such simple but elegant
> examples. Your post might have more efficiently elicited more and more
> elegant responses sooner with a more carefully chosen Subject, perhaps
> like "Multiple Comparisons Questions". However, the selection of a
> possible better subject might rely on information you didn't have.
>
> Hope this helps.
> Spencer Graves
> >
> > My R.version:
> >
> > platform i386-pc-mingw32
> > arch i386
> > os mingw32
> > system i386, mingw32
> > status
> > major 2
> > minor 2.1
> > year 2005
> > month 12
> > day 20
> > svn rev 36812
> > language R
> >
> > ______________________________________________
> > R-help <at> stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
> ______________________________________________
> R-help <at> stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
>
Dear Spencer,
Thank you very much for this detailed answer to my problem.
I pick up the discussion quite late, because I was in holiday.
For this fast moving times, especially within the R-help group,
a delay for an answer of 10 days is extraordinary.
Hope for your understanding.
Dominik
More information about the R-help
mailing list