[R] Multiple tests on 2 way-ANOVA

Spencer Graves spencer.graves at pdf.com
Sat Jul 15 17:32:42 CEST 2006


<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.

> 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



More information about the R-help mailing list