R help - Jul 2006 - Multiple tests on 2 way-ANOVA

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

I do not believe that this is the most efficient way for doing this, 
since I made already bad experience with the first example.

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

<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

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