
rm(list=ls())
rnorm2 <- function(n,mean,sd) { 
  mean+sd*scale(rnorm(n)) 
}
ss <- function(x) {
  sum((x-mean(x))^2)
}
set.seed(10)
n <- 30
n.o <- n.p <- n
o <- rnorm(n.o, 100, 10)
p <- rnorm(n.p, 104, 10)

# old way
m.o <- mean(o)
m.p <- mean(p)
df.o <- n.o - 1
df.p <- n.p - 1
diff <- m.o - m.p
m.o
m.p
diff
pv <- (ss(o)+ss(p))/(df.o+df.p)
se <- sqrt((pv/n.o) + (pv/n.p))
t.cal <- diff/se
t.cal
pt(t.cal, df.o+df.p) * 2

t.out <- t.test(o, p, var.equal=T)
t.out
t.out$statistic
t.out$p.value

#
comb <- list(o = o, p = p)
o
p
comb
op <- stack(comb)
op
colnames(op)[1] <- "values"
colnames(op)[2] <- "group"
op$group <- factor(op$group)
head(op)
tail(op)
boxplot(op$values~op$group)

# plot(op$values~op$group)
boxplot(op$values~op$group, main="values by group",
        yaxt="n", xlab="value", horizontal=TRUE,
        col=terrain.colors(2))
abline(v=mean(op$values), col="red", lwd=3)
legend("topleft", inset=.05, title="group",
       c("o","p"), fill=terrain.colors(2), horiz=TRUE)

m.tot <- mean(op$values)
m.o <- mean(o)
m.p <- mean(p)
m.tot
m.o 
m.p 

min.x <- min(op$values)
max.x <- max(op$values)
br <- seq(floor(min.x), ceiling(max.x), by = 1) 

hist(o, breaks=br, 
     col=rgb(1,1,1,.5))
abline(v=m.o, col="black", lwd=3)
hist(p, add=T, breaks=br,
     col=rgb(.5,1,1,.5))
abline(v=m.p, col="blue", lwd=3)
abline(v=m.tot, col='red', lwd=3)

hist(o, breaks=br, 
     col=rgb(1,1,1,.5))
hist(p, add=T, breaks=br,
     col=rgb(.5,1,1,.5))
abline(v=m.tot, col='red', lwd=3)

ss.tot <- ss(op$values)
df.tot <- length(op$values)-1
ss.tot/df.tot
var(op$values)
ss.tot
m.tot
m.o
m.p
ss.o <- ss(o)
ss.p <- ss(p)
ss.o
ss.p

hist(o, breaks=br, 
     col=rgb(1,1,1,.5), 
     main = paste0("o mean = ", round(m.o,2), 
                   "\n", "tot mean = ", round(m.tot,2), 
                   ""))
abline(v=m.o, col="black", lwd=3)
abline(v=m.tot, col='red', lwd=3)

hist(p, breaks=br,
     col=rgb(.5,1,1,.5),
     main = paste0("p mean = ", round(m.p,2), 
                   "\n", "tot mean = ", round(m.tot,2), 
                   "")     )
abline(v=m.p, col="blue", lwd=3)
abline(v=m.tot, col='red', lwd=3)

n.o*(m.tot-m.o)^2   # m.tot 에서 o그룹공통 까지의 거리를 제곱해서 모두 더한 값
                    # 위 그림에서 빨간색 선에서 검은색 선까지의 거리를 
                    # 제곱해서 모두 더한 값
n.p*(m.tot-m.p)^2   # m.tot 에서 p그룹공통 까지의 거리를 제곱해서 모두 더한 값
                    # 위 그림에서 빨간색 선에서 파란색 선까지의 
                    # 거리를 제곱해서 모두 더한 값
> ss.bet <- n.o*(m.tot-m.o)^2 +   # 이것은 그룹 (IV, 독립변인) 때문에 
            n.p*(m.tot-m.p)^2     # 생긴 그룹 간 차이이다
> ss.bet            # 따라서 이것을 SS between group이라고 부른다


# ss within group 

hist(o, breaks=br, 
     col=rgb(1,1,1,.5))
abline(v=m.o, col="black", lwd=3)
ss.o <- ss(o)   # o집단의 평균인 검은색 선에서 개인 점수까지의 
                # 거리는 (오차는) 독립변인과 상관없이 랜덤하게 
                # 나타나는 것
ss.o            # o집단의 것을 ss.o라고 부른다

hist(p, breaks=br,
     col=rgb(.5,1,1,.5))
abline(v=m.p, col="blue", lwd=3)
ss.p <- ss(p)   # p집단도 마찬가지이다. 이 집단 내의 sum of square값은 
                # p 집단의 공통특징인 평균에서 개인점수가 랜덤하게 
ss.p            # 나타나는 것이고, 이것을 sum of square p라고 부른다

ss.wit <- ss.o + ss.p # 위 각 집단 내부의 ss값을 합한 값을 ss within group
                      # 이라고 부른다

# 이 둘은 각 그룹의 평균을 중심으로 
# random 하게 나타나는 평균에서의 거리이다 (에러).
# 이것을 sum of square within group 이라고 부른다
# 이것을 우리는 random difference라고 불러도 되겠다
ss.wit

# 반면에 ss between group은 IV의 효과 때문에 나타나는
# 것이다. 우리는 이것을 group difference라고 부른다.
ss.bet
ss.bet + ss.wit
ss.tot

# degrees of freedom
df.tot <- length(op$values)-1
df.bet <- nlevels(op$group) - 1
df.wit <- (n.o-1) + (n.p-1)
df.tot
df.bet
df.wit

# variance of total, between group, and within group
ms.tot <- ss.tot / df.tot
ms.bet <- ss.bet / df.bet
ms.wit <- ss.wit / df.wit

f.cal <- ms.bet / ms.wit
f.cal
p.val <- pf(f.cal, df1=df.bet, df2=df.wit, lower.tail = F)
p.val
summary(aov(op$values~op$group))

# t.cal 값을 확인해보면
t.test(o,p, var.equal = T) # 혹은
t.cal <- diff/se 
p.t.cal <- pt(abs(t.cal), df=df.o+df.p, lower.tail = F)*2
p.t.cal

# 하나는 sd 를 unit으로 계산한 값이고 (t.cal)
# 다른 하나는 variance를 unit으로 계산한 값이지만 
# 원리는 같다. 
t.cal^2
f.cal