Guide to Univariate & Multivariate Approaches to Mixed-Design 2-Way ANOVA in R
(1 between-subject factor, 1 within-subjects factor [repeated measure])

Jason R. Finley
July 19th, 2009

UNIVARIATE APPROACH TO MIXED-DESIGN 2-Way ANOVA
(1 between-subject factor, 1 within-subjects factor)

Factor A is between-subjects
	(Factor A was presentation method.  it had 2 levels: auditory and visual.  so a=2)
Factor B is within-subjects (repeated measure)
	(Factor B was how long the subject had been going on the task or something like that.  So B1=1hour, B2=2hours, B3=3hours, B4=4hours.  so b=4)
The dependent measure was like response time or something.  Doesn't really matter.

hw9data1a.txt: 
subject	A	B	Y
s1	1	1	3
s1	1	2	4
s1	1	3	7
s1	1	4	7
s2	1	1	6
s2	1	2	5
s2	1	3	8
s2	1	4	8
s3	1	1	3
s3	1	2	4
s3	1	3	7
s3	1	4	9
s4	1	1	3
s4	1	2	3
s4	1	3	6
s4	1	4	8
s5	2	1	1
s5	2	2	2
s5	2	3	5
s5	2	4	10
s6	2	1	2
s6	2	2	3
s6	2	3	6
s6	2	4	10
s7	2	1	2
s7	2	2	4
s7	2	3	5
s7	2	4	9
s8	2	1	2
s8	2	2	3
s8	2	3	6
s8	2	4	11

Load the file into a data structure & assign some variable names to make things easier (taking care to "factorize" A and B):
> filepath <- c("/Users/uniace/Desktop/hw9data1a.txt")
> data91a <- read.table(filepath, header=TRUE)
> A=factor(data91a$A)
> B=factor(data91a$B)
> subject=data91a$subject
> Y=data91a$Y
The following is an alternative to naming out all the variables by hand:
> attach(data91a)

NOTE: YOU MUST FACTORIZE ALL VARIABLES (factors) EXCEPT THE DEPENDENT MEASURE!!!  This includes SUBJECT factor!  The only reason it wasn't factorized above is because the subject numbers used letters, so R automatically treated that as a factor.  Don't risk anything, just factorize all factors!!!

Make overall ANOVA:
> aov91a <- aov(Y~(B*A) + Error(subject/B))
Note:
	•	For a 2-Way Design with 1 between-subjects factor and 1 within-subjects factor (A=between-Ss, B=within-Ss), use:
	⁃	aov(Y~(B*A) + Error(subject/B))
	•	For a 2-Way Design with 2 between-subjects factors, we would use:
	⁃	aov(Y~B*A)
	•	For a 2-Way Design with 2 within-subjects factors, we would use:
	⁃	aov(Y~(B*A) + Error(subject/(B*A))
	•	The order (e.g., B*A vs A*B) shouldn't matter, if you have a balanced design (equal sample size in all between-subjects conditions).  I think.
	•	For a 1-Way Design within-subjects, we would use:
	⁃	aov(Y~B + Error(subject/B))
	•	For a 1-Way Design between-subjects:
	⁃	aov(Y~A)

Now see results of ANOVA!:
> summary(aov91a)

Error: subject
          Df Sum Sq Mean Sq F value Pr(>F)
A          1 3.1250  3.1250       2 0.2070
Residuals  6 9.3750  1.5625               

Error: subject:B
          Df  Sum Sq Mean Sq F value    Pr(>F)    
B          3 194.500  64.833 127.890 2.516e-12 ***
B:A        3  19.375   6.458  12.740 0.0001051 ***
Residuals 18   9.125   0.507                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

We can cut & paste to make one pretty ANOVA table (I did the math myself to make the TOTAL line):
ANOVA TABLE:
          Df  	Sum Sq 	Mean Sq 	F value	Pr(>F) 
A          1 	3.1250  	3.1250    2 		0.2070   
B          3 	194.500  	64.833 	127.890 	2.516e-12 ***
B:A        3  	19.375   	6.458  	12.740 	0.0001051 ***
subject    6 	9.3750  	1.5625   
subject:B 18  	9.125   	0.507  
TOTAL  	31	235.5
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

ADJUSTMENTS FOR VIOLATIONS OF SPHERICITY:
I could not figure out how to get the Greenhouse-Geisser or Huynh-Feldt Epsilon values for 2-way mixed design using R.  Angie found a webpage that says how to do it for 1-way within-subjects design:
http://www.psych.upenn.edu/~baron/rpsych/rpsych.html#SECTION000711000000000000000

These are the relevant notes I have from Stats lecture with David Budescu:
Multiple-Group Sphericity Assumption(s):
1) sphericity holds in each group
2) covariance matrix same in all groups
...
so the approach to Epsilon adjustment:  ..... calculated based on pooled matrix ... & the same adjustment is applied for test of B and A*B

So i got the 2 covariance matrices for Factor B (one at each level of A), but then i didnt know what to do from there, like how to pool it into one matrix.  I found some little info online on how to pool matrices, but didn't really understand what/how to do:
http://tolstoy.newcastle.edu.au/R/e2/help/06/09/1214.html
If anyone figures it out, please let me know!

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

MULTIVARIATE APPROACH TO MIXED-DESIGN 2-Way ANOVA
(1 between-subject factor, 1 within-subjects factor)

Factor A is between-subjects
	(Factor A was presentation method.  it had 2 levels: auditory and visual.  so a=2)
Factor B is within-subjects (repeated measure)
	(Factor B was how long the subject had been going on the task or something like that.  So B1=1hour, B2=2hours, B3=3hours, B4=4hours.  so b=4)
The dependent measure was like response time or something.  Doesn't really matter.

hw9data1b.txt: 
subject	A	B1	B2	B3	B4
s1	auditory	3	4	7	7
s2	auditory	6	5	8	8
s3	auditory	3	4	7	9
s4	auditory	3	3	6	8
s5	visual	1	2	5	10
s6	visual	2	3	6	10
s7	visual	2	4	5	9
s8	visual	2	3	6	11

Load the file into a data structure:
> filepath <- c("/Users/uniace/Desktop/hw9data1b.txt")
> data91b <- read.table(filepath, header=TRUE)

Test for effect of Factor A:
This is done separately from the tests for effect of Factor B and test for interactions.  It is done in the traditional (univariate) way, since Factor A is between-subjects so there's no need to get all fancy and multivariate.  That said, it's important that we test Factor A in the context of the full model (which includes Factor B and Subject factor... I think).  This means that we specify the error term just like we did with the univariate overall ANOVA:
> aov91a <- aov(Y~A + Error(subject/B))
> summary(aov91a)
          Df Sum Sq Mean Sq F value Pr(>F)
A          1 3.1250  3.1250       2 0.2070
Residuals  6 9.3750  1.5625 

Test for effect of FACTOR B (within-Ss) & INTERACTION A:B
The following R commands and table give the results of a MANOVA using the Wilks test, including the test result for B and for the interaction A:B.  We must specify a set of b-1 linearly independent contrasts that we'll test to see if they're all simultaneously equal to zero (for testing effect of Factor B) and test to see if they're all simultaneously equal to each other across levels of A (for testing interaction of A & B). For illustration we'll first use the set of polynomial contrasts. 
> contr.poly(4)
             .L   .Q         .C
[1,] -0.6708204  0.5 -0.2236068
[2,] -0.2236068 -0.5  0.6708204
[3,]  0.2236068 -0.5 -0.6708204
[4,]  0.6708204  0.5  0.2236068
NOTES:
	•	.L = LINEAR contrast
	•	.Q = QUADRATIC contrast
	•	.C = CUBIC contrast
	•	Note how it's turned on its side from how we would normally write it out!

> man<-manova(cbind(data91b$B1,data91b$B2,data91b$B3,data91b$B4) %*% contr.poly(4) ~ A , data=data91b)
Note:
	•	If this was just a 1-Way within-subjects design, we would replace "~ A" with "~ 1".  I'm not sure what we'd do if we had a 2-Way design with both factors within-subjects
> summary(man, test="Wilks", intercept=TRUE)
            Df  Wilks approx F num Df den Df   Pr(>F)   
(Intercept)  1  0.027   47.192      3      4 0.001403 **
A            1  0.144    7.906      3      4 0.037115 * 
Residuals    6                                          
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

NOTES:
	•	The "(Intercept)" line amounts to the multivariate test for a sig. effect of Factor B (within-Ss). 
	•	The "A" line is for the interaction A:B.
	•	%*% means matrix multiplication
	•	contr.poly(4)    gives you the set of 3 polynomial contrasts (orthogonal)
	•	you MUST include intercept=TRUE  when calling summary.  otherwise it won't show the intercept line, and that's what you're really interested in!
	•	all the tests you can use are: "Pillai", "Wilks", "Hotelling-Lawley", "Roy

Test the individual contrasts that we specified (& their interactions with Factor A):
> aovman<-aov(man)
> summary(aovman, intercept=TRUE)
 Response .L :
            Df  Sum Sq Mean Sq F value    Pr(>F)    
(Intercept)  1 184.900 184.900 182.617 1.018e-05 ***
A            1  13.225  13.225  13.062   0.01118 *  
Residuals    6   6.075   1.013                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

 Response .Q :
            Df Sum Sq Mean Sq F value   Pr(>F)   
(Intercept)  1 8.0000  8.0000    25.6 0.002311 **
A            1 3.1250  3.1250    10.0 0.019509 * 
Residuals    6 1.8750  0.3125                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

 Response .C :
            Df  Sum Sq Mean Sq F value  Pr(>F)   
(Intercept)  1 1.60000 1.60000  8.1702 0.02886 * 
A            1 3.02500 3.02500 15.4468 0.00771 **
Residuals    6 1.17500 0.19583                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

NOTES:
	•	.L = LINEAR contrast (tests for linear trend)
	•	.Q = QUADRATIC contrast
	•	.C = CUBIC contrast
	•	The "(Intercept)" line is the test for whether the particular contrast is significant.
	•	The "A" line is for the interaction of the contrast with Factor A.
	•	you MUST include intercept=TRUE  when calling summary.  otherwise it won't show the intercept line, and that's what you're really interested in!

Try it with a different set of contrasts:
Let's try using a different set of contrasts for the multivariate test.  Let's try profile contrast set.
> profile <- matrix(c(1,-1,0,0, 0,1,-1,0, 0,0,1,-1), nrow = 3, ncol=4, byrow=TRUE)
> profile
     [,1] [,2] [,3] [,4]
[1,]    1   -1    0    0
[2,]    0    1   -1    0
[3,]    0    0    1   -1 
Note: but we need to transpose this first!
> profile=t(profile)
> profile
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]   -1    1    0
[3,]    0   -1    1
[4,]    0    0   -1

> man2<-manova(cbind(data91b$B1,data91b$B2,data91b$B3,data91b$B4) %*% profile ~ A , data=data91b)
> summary(man2, test="Wilks", intercept=TRUE)
            Df  Wilks approx F num Df den Df   Pr(>F)   
(Intercept)  1  0.027   47.192      3      4 0.001403 **
A            1  0.144    7.906      3      4 0.037115 * 
Residuals    6                                          
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Note how we get the same values for the Wilks test!  This is expected, because the set of contrasts you use doesn't matter as long as they're linearly independent.  But now let's take a look at the individual contrast results with our new profile contrast set:

> aovman2<-aov(man2)
> summary(aovman2, intercept=TRUE)
 Response 1 :
            Df Sum Sq Mean Sq F value  Pr(>F)  
(Intercept)  1 4.5000  4.5000  7.7143 0.03210 *
A            1 2.0000  2.0000  3.4286 0.11353  
Residuals    6 3.5000  0.5833                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

 Response 2 :
            Df Sum Sq Mean Sq F value    Pr(>F)    
(Intercept)  1   60.5    60.5     121 3.355e-05 ***
A            1    0.5     0.5       1    0.3559    
Residuals    6    3.0     0.5                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

 Response 3 :
            Df Sum Sq Mean Sq F value    Pr(>F)    
(Intercept)  1 60.500  60.500    72.6 0.0001432 ***
A            1 24.500  24.500    29.4 0.0016295 ** 
Residuals    6  5.000   0.833                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Notes:
	•	Each of the three contrasts compares 2 levels of B.  There's probably same way to have made these better labeled, but here's what they mean:
	•	Response 1 compares B1 to B2
	•	Response 2 compares B2 to B3
	•	Response 3 compares B3 to B4
	•	The "Intercept" line shows the result of testing the contrast collapsed across all levels of A (I think).
	•	The "A" line shows the result of testing for interaction between the contrast and Factor A.  
	•	The significant interaction in Response 3 means that the difference between B3 and B4 (subjects' performance at 3 hours vs 4 hours) is not the same across the two levels of Factor A (auditory vs visual).