COMMunication
RESearch.NET

This is an old revision of the document!

---

#################################################
# two-way anova
# subject = factor(paste('sub', 1:30, sep=''))
#################################################

n.a.group <- 3 # a treatment 숫자
n.b.group <- 2 # b 그룹 숫자
n.sub <- 30 # 총 샘플 숫자
n.sub/n.a.group

# 데이터 생성
set.seed(9)
a <- gl(3, 10, n.sub, labels=c('a1', 'a2', 'a3'))
b <- gl(2, 5, n.sub, labels=c('b1', 'b2'))
a
b
y <- rnorm(30, mean=10) + 
  3.14 * (a=='a1' & b=='b2') + 
  1.43 * (a=='a3' & b=='b2') 
y

dat <- data.frame(a, b, y)
dat
# aov.dat <- aov(y~a*b) # anova test
# summary(aov.dat) # summary of the test output

# hand calculation
table(a,b)
tapply(y, list(a,b), mean) # 각 셀에서의 평균
n.within.each <- tapply(y, list(a,b), length)
n.within.each
df.within.each <- n.within.each - 1  # 각 셀에서의 샘플숫자
df.within.each
df.within <- sum(df.within.each) # df within
df.within

var.within <- tapply(y, list(a,b), var) # var.within
var.within
ss.within.each <- tapply(y, list(a,b), var) * df.within.each
ss.within.each
ss.within <- sum(ss.within.each) # ss.within
ss.within

ms.within <- ss.within / df.within
ms.within

# interaction.plot(a,b,y)

mean.a <- tapply(y, list(a), mean)
mean.b <- tapply(y, list(b), mean)
mean.a
mean.b

var.a <- tapply(y, list(a), var)
var.b <- tapply(y, list(b), var)
var.a
var.b

mean.tot <- mean(dat$y)
var.tot <- var(dat$y)
df.tot <- n.sub - 1 
ss.tot <- var.tot * df.tot
ms.tot <- ss.tot/df.tot

## between
mean.each <- tapply(y, list(a,b), mean)
mean.each 
mean.tot <- mean(y)
mean.tot
n.each <- tapply(y, list(a,b), length)
n.a.each <- tapply(y, list(a), length)
n.b.each <- tapply(y, list(b), length)
n.each
n.a.each
n.b.each



ss.bet <- sum(n.each*(mean.each-mean.tot)^2)
ss.bet

ss.tot
ss.within
ss.bet
ss.bet + ss.within

ss.a <- sum(n.a.each * ((mean.tot - mean.a)^2))
ss.b <- sum(n.b.each * ((mean.tot - mean.b)^2))
ss.a
ss.b
ss.ab <- ss.bet - (ss.a + ss.b) # ss.ab = ss.interaction
ss.ab

ss.tot
ss.bet
ss.within
ss.a
ss.b
ss.ab

df.tot <- n.sub - 1
df.bet <- (n.a.group * n.b.group) - 1
df.a <- n.a.group - 1
df.b <- n.b.group - 1
df.ab <- df.bet - (df.a + df.b)
df.within <- sum(df.within.each)

df.tot
df.bet
df.a
df.b
df.ab
df.within

ms.tot <- ss.tot/df.tot # we did it above 
ms.bet <- ss.bet - df.bet
ms.a <- ss.a / df.a
ms.b <- ss.b / df.b
ms.ab <- ss.ab / df.ab
ms.within <- ss.within / df.within

ms.tot
ms.bet
ms.within
ms.a
ms.b
ms.ab


f.a <- ms.a / ms.within
f.b <- ms.b / ms.within
f.ab <- ms.ab / ms.within

alpha <- .05
# confidence interval
ci <- 1 - alpha

f.a
# 봐야할 F분포표에서의 F-value
# qt 처럼 qf 사용
# qf(alpha, df.between, df.within, lower.tail=F) 처럼 사용
qf(ci, df.a, df.within) # or
qf(alpha, df.a, df.within, lower.tail = F)
# 혹은 
# qf(alpha, df.a, df.within, lower.tail = F)
# 도 마찬가지
pf(f.a, df.a, df.within, lower.tail = F)

f.b
qf(ci, df.b, df.within)
pf(f.b, df.b, df.within, lower.tail = F)

f.ab
qf(ci, df.ab, df.within)
pf(f.ab, df.ab, df.within, lower.tail = F)

# aov result
summary(aov.dat)