### Howell Table 13.2 ###
## A x B anova  Eysenck example
# Read in data and factor() Age and Condition
# This example has equal sample sizes, so is balanced

###  This solves the problem a couple of ways, not all of which agree. The
###  point is to show that the type = 3 solution for Anova gives a wrong answer here.
###  Well, you get the right answer if you set sum contrasts in a particular way.
###  The comments were written as I went along, and I left them that way.

par(mfrow=c(2,2))
Eysenck <- read.table("http://www.uvm.edu/~dhowell/methods7/DataFiles/Tab13-2.dat",
 header = TRUE)
Eysenck$subj <- factor(1:100)
Eysenck$Condition <- factor(Eysenck$Condition, levels = 1:5, labels = c("Counting",
     "Rhyming", "Adjective", "Imagery","Intention"))
Eysenck$Age <- factor(Eysenck$Age, levels = 1:2, labels = c("Old","Young"))
attach(Eysenck)
#Compute the anova and then print out the result
result <- anova(aov(Recall~Condition*Age, data = Eysenck))
print(result)
plot(Recall~Condition)
interaction.plot(x.factor = Condition, trace.factor = Age, response = Recall, fun = mean, legend = TRUE)

#better Yet!  Maybe
# This should allow me to get Type III SS. But when I use "type = "III" " I get a
# strange answer even though the design is orthogonal. If I use "type = II" with
# this orthogonal designs I get the right answer. Odd because this should give
# legitimate type III SS whether the design is orthogonal or not. And drop1 also
# gives the wrong answer--though it agrees with type III.
library(car)
options(contrasts = c("contr.sum","contr.poly"))
resultsCar <- lm(Recall~Age*Condition, data = Eysenck )
type2 <- Anova(resultsCar, type = "II")
type3 <- Anova(resultsCar, type = 3)     #type3 and the next are wrong if you type contrasts(Conditions) = contr.sum.
                                         #Perhaps Fox uses drop1
drop1(resultsCar, scope = ~Age*Condition, test = "F")
                                          # I GOT IT. Typing contrasts(Conditions) <- contr.sum
                                          # does not work, but typing
                                          # options(contrasts = c('contr.sum', 'contr.poly')) does work

################################################################################
# Now do it via regression with dummy variables
A <- rep(c(1,-1), each = 50)
C1 <- rep(c(1,0,0,0,-1,1,0,0,0,-1), each = 10)
C2 <- rep(c(0,1,0,0,-1,0,1,0,0,-1), each = 10)
C3 <- rep(c(0,0,1,0,-1,0,0,1,0,-1), each = 10)
C4 <- rep(c(0,0,0,1,-1,0,0,0,1,-1), each = 10)
AC1 <- A*C1
AC2 <- A*C2
AC3 <- A*C3
AC4 <- A*C4
data <- cbind(Recall, A, C1, C2, C3, C4, AC1, AC2, AC3, AC4)
SStotal <- sum(summary(ModelFull)[[1]][[2]][1:10])
ModelFull <- aov(Recall~A+C1+C2+C3+C4+AC1+AC2+AC3+AC4)
SSresid.full <-  sum(summary(ModelFull)[[1]][[2]][10])
ModelA <- aov(Recall ~ C1+C2+C3+C4+AC1+AC2+AC3+AC4)
SSA <- summary(ModelA)[[1]][[2]][9]-SSresid.full
ModelC <- aov(Recall~A+AC1+AC2+AC3+AC4)
SSC <- summary(ModelB)[[1]][[2]][6]-SSresid.full
ModelAC <- aov(Recall~A+C1+C2+C3+C4)
SSAC <- summary(ModelAC)[[1]][[2]][6]-SSresid.full
SumSq <- list(SSA, SSC, SSAC, SSresid.full,SStotal)
dfA <- 1; dfC <- 4; dfAC <- 4; dfResid <- 90; dfTotal <- 99
options(digits = 7)
MSA <- SSA/dfA
MSC <- SSC/dfC
MSAC <- SSAC/dfAC
MSresid <- SSresid.full/dfResid
FA <-MSA/MSresid
FC <- MSC/MSresid
FAC <- MSAC/MSresid
cat("         df   SS       MS         F           p  \n"," A      ",   dfA, SSA, '  ',MSA,'  ',FA,'  ',1- pf(FA,dfA,dfResid),' \n',
                                          " C      ",  dfC, SSC,' ', MSC,' ',FC,'  ',1- pf(FC,dfC,dfResid), '\n',
                                          " AC     ", dfAC, SSAC,'   ', MSAC,'  ',FAC,'  ',1- pf(FAC,dfAC,dfResid), '\n',
                                          "Resid  ",  dfResid, SSresid.full, '   ', MSresid,'\n'  )
cat('\n \n')

###############################################################################
# Now for effect size measures
library(MBESS)

C.means <- tapply(Recall, Condition, mean)
A.means <- tapply(Recall, Age, mean)
AC.means <- tapply(Recall, list(Age, Condition),mean)
a <- length(levels(Age))
c <- length(levels(Condition))
n <- length(Recall)/(a*c)
sigma2A <- (a-1)*(MSA-MSresid)/(n*a*c)
sigma2C <- (c-1)*(MSC-MSresid)/(n*a*c)
sigma2AC <- (a-1)*(c-1)*(MSAC-MSresid)/(n*a*c)
sigma2total <- sigma2A + sigma2C + sigma2AC + MSresid
omega2A <- sigma2A/sigma2total
omega2C <- sigma2C/sigma2total
omega2AC <- sigma2AC/sigma2total

cat("Omega squared as the effect size measure \n")
cat("Effect     Omega Sq \n")
cat(" A        ", omega2A, '\n',"C        ",omega2C,'\n',"AC       ",omega2AC, '\n \n')

# Now for partial effects
partialA <- sigma2A/(sigma2A + MSresid)
partialC <- sigma2C/(sigma2C + MSresid)
partialAC <- sigma2AC/(sigma2AC + MSresid)

cat("Partial omega squared as the effect size measure \n")
cat("Effect    Partial Omega Sq \n")
cat(" A        ", partialA, '\n',"C        ",partialC,'\n',"AC       ",partialAC, '\n \n')

################################################################################

# Now to get Simple Effects

# I am going to test each effect with its own error term, though you could
# pool the error terms if you wish.


old.dat <- subset(Eysenck, Eysenck$Age == "Old")
young.dat <- subset(Eysenck, Eysenck$Age == "Young" )

results1 <-anova(aov(Recall ~ Condition, data = old.dat) )
cat("\n\n\t\tResults for older subjects\n")
print(results1)
plot(Recall~Condition, data = old.dat, main = "Older Subjects")

results2 <-anova(aov(Recall ~ Condition, data = young.dat) )
cat("\n\n\t\tResults for younger subjects\n")
print(results2)
plot(Recall~Condition, data = young.dat, main = "Younger Subjects")


################################################################################