#
# 3 샘플 종류 추출
# A, B, C 학년에 따라서 욕하는 정도가 달라질것이라는 
# 가설
set.seed(201)
rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) }
A <- rnorm2(16, 26, sqrt(600/15))
B <- rnorm2(16, 24, sqrt(750/15))
C <- rnorm2(16, 19, sqrt(900/15))

A <- c(A)
B <- c(B)
C <- c(C)

# 평균구하기
mean(A)
mean(B)
mean(C)

# 3 샘플을 합치기 
# 두번재 컬럼 = group A, B, C 가 되도록
comb3 <- stack(list(a=A, b=B, c=C))
comb3
colnames(comb3)[1] <- "values"
colnames(comb3)[2] <- "group"
comb3

# 전체구성원을 하나로 보고 분산값을 구해본다 
# ms.tot = ss.tot / df.tot
# ss.tot = sum(xi-mean.tot)^2
# df.tot = N - 1
# attach(comb3)
mean.tot <- mean(comb3$values)
ss.tot <- sum((comb3$values-mean.tot)^2)
df.tot <- 48-1
ms.tot <- ss.tot/df.tot
ss.tot
df.tot
ms.tot
var(comb3$values)

# A, B, C 평균
m.a <- mean(A)
m.b <- mean(B)
m.c <- mean(C)
ms.a <- var(A)
ms.b <- var(B)
ms.c <- var(C)
df.a <- 15
df.b <- 15
df.c <- 15
# 사실 우리는 아래 각 그룹의 SS가 
# 600, 750, 900 인것을 알고 있다. 
ss.a <- ms.a*df.a
ss.b <- ms.b*df.b
ss.c <- ms.c*df.c

ss.within <- ss.a + ss.b + ss.c
ss.within <- ss.within

16*((m.a-mean.tot)^2)
16*((m.b-mean.tot)^2)
16*((m.c-mean.tot)^2)
ss.bet <- 16*((m.a-mean.tot)^2) + 16*((m.b-mean.tot)^2) + 16*((m.c-mean.tot)^2)

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

df.tot 
df.within <- df.a + df.b + df.c
df.within
df.bet <- 3-1 
df.bet
df.bet + df.within

ms.bet <- ss.bet/df.bet
ms.within <- ss.within /df.within 
f.calculated <- ms.bet/ms.within
ms.bet
ms.within
f.calculated 
# 이 점수에서의 F critical value =
fp.value <- pf(f.calculated, df1=2, df2=45, lower.tail = F)
fp.value

# f.calculated 와 fp.value로 판단
f.calculated
fp.value


a.res <- aov(values ~ group, data=comb3)
a.res.sum <- summary(a.res)
a.res.sum

# 위에서 
# ssd라는 function을 만들면 

ssd <- function(x) { 
          var(x) * (length(x)-1) }
ss.a1 <- ssd(A)
ss.b2 <- ssd(B)
ss.c1 <- ssd(C)

ss.a == ss.a1
ss.b == ss.b1
ss.c == ss.c1

# 그러나 정확히 어떤 그룹에서 차이가 나는지는 판단해주지 않음 
pairwise.t.test(comb3$values, comb3$group, p.adj = "none")
# OR
pairwise.t.test(comb3$values, comb3$group, p.adj = "bonf")
pairwise.t.test(comb3$values, comb3$group, p.adj = "holm")

# OR TukeyHSD(anova.output)
TukeyHSD(a.res)