pdf("kuehl_06_06.pdf")

A <- rep(1:7,each=7)
B <- rep(1:7,times=7)

pcthear16 <-c(2.1,6.8,8.4,1.4,14.6,7.9,4.8,
            1.7,8.1,8.4,1.4,12.0,3.7,4.5,
            14.4,14.8,27.0,30.9,36.5,36.4,31.4,
            57.4,62.4,37.4,63.3,65.5,65.6,59.8,
            66.2,81.7,53.3,80.7,79.7,80.8,82.4,
            75.2,94.0,74.3,87.9,93.3,87.8,80.5,
            4.1,10.2,10.7,5.5,18.1,11.4,6.1)


audio <- data.frame(A,B,pcthear16)

attach(audio)

meanA <- tapply(pcthear16,A,mean)
meanB <- tapply(pcthear16,B,mean)
meanAB <- mean(pcthear16)
A_levels <- max(A)
B_levels <- max(B)

SSE <- 0
SSA <- 0
SSB <- 0
P <- 0

#for (j in 1:B_levels) {
for (i in 1:A_levels) {
for (j in 1:B_levels) {
P <- P + (meanA[[i]] - meanAB)*(meanB[[j]]-meanAB)*pcthear16[(i-1)*B_levels+j]
SSE <- SSE + (pcthear16[(i-1)*B_levels+j]-meanA[[i]]-meanB[[j]]+meanAB)^2
SSA <- SSA + (meanA[[i]]-meanAB)^2
SSB <- SSB + (meanB[[j]]-meanAB)^2
}
}

P
SSE
SSA
SSB

SSNA <- P^2/((SSA/B_levels)*(SSB/A_levels))
SSNA

SSRES <- SSE-SSNA
SSRES

MSA <- SSA/(A_levels-1)
MSB <- SSB/(B_levels-1)
MSE <- SSE/((A_levels-1)*(B_levels-1))
MSNA <- SSNA/1
MSRES <- SSRES/((A_levels-1)*(B_levels-1)-1)

F.NA <- MSNA/MSRES
P_F.NA <- 1-pf(F.NA,1,(A_levels-1)*(B_levels-1)-1)

F.NA
P_F.NA



dev.off()