COMMunication
RESearch.NET

##
## t-test 이해 확인
## remember t test = diff / rand error
##

set.seed(101)
rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) }

A <- rnorm2(16, 26, sqrt(1160/15))
B <- rnorm2(16, 19, sqrt(1000/15))
A <- c(A)
B <- c(B)

# we know sqrt(1160/15) is A's sdev
# hence, A's var is sqrt(1160/15)^2
# hence, A's SS is sqrt(1160/15)^2 * 15
# this is 1160

# from the above,
# the difference between the A and B means
# remember we try to find
# difference due to the treatment / 
# / random chance of error
diff <- mean(A) - mean(B)

# for se
# we know that the situation refers to
# #2 se two independent samples t-test
# which is sqrt(pooled.var/na + pooled.var/nb)

SSa <- 1160
SSb <- 1000
n.a <- 16
n.b <- 16
df.a <- n.a - 1 
df.b <- n.b - 1

# 확률과통계 수업 학생은 아래에서 
# se값을 구한 결과만 보고 이것이 
# se = sqrt(sigma^2/n)에
# 대응하는 값으로 생각하세요

# we know that we are testing the difference 
# between two independent sample means.
# Hence, we need to use poole variance between 
# the two group. See
# http://commres.net/wiki/t-test#t-test_%EB%B9%84%EA%B5%90
pooled.var <- (SSa + SSb) / (df.a + df.b)
se <- sqrt( (pooled.var/n.a) + (pooled.var/n.b) )

# Remember t test calculation is based on 
# diff / random error 
t.calculated <- diff / se
pooled.var
diff
se
t.calculated

# Now use t.test function for two group
# (independent sample) t-test
# with an assumption that variances of 
# the two gorup are the same.
t.result <- t.test(A, B, var.equal = T)
t.result

# t.result$statistic = t.calculated
# t.result$p.value = probability level of 
# wrong decision with the t calculated value
str(t.result)
t.result$statistic
t.result$p.value

# the above p.value can be obtained with 
# pt function
p.value.calculated <- 2*pt(t.calculated, df = df.a + df.b, lower.tail=FALSE)
p.value.calculated
t.result$p.value

##
## different approach (분산으로 분석)
## 

# A combined group with group A and B
# We call it group total
# we can obtain its mean, variance, ss, df, etc.
# 
A
B
dat <- c(A,B)
dat

mean.total <- mean(dat)
var.total <- var(dat)
# variance를 ms라고 부르기도 한다
# ms = mean square (square 값을 n 으로 나눈 값이 분산이므로)
ms.total <- var.total

df.total <- length(dat) - 1
ss.total <- var.total * df.total
# ss.total 은 아래처럼 구해도 된다
ss.total.check <- sum((dat-mean(dat))^2)

mean.total
var.total
ms.total
df.total 
ss.total
ss.total.check

# Now for each group
mean.a <- mean(A)
mean.b <- mean(B)
mean.a
mean.b

# 그룹 간의 차이에서 나타나는 분산
# 수업시간에 설명을 잘 들을 것

# mean.total 에서 그룹a의 평균까지의 차이를 구한 후
# 이를 제곱하여 그룹 A 멤버의 숫자만큼 곱한다 
# 이를 그룹b의 평균과도 다시 계산하고 
# 이 둘을 더하여
# ss.between에 저장

length(A) * ((mean.total - mean.a)^2)
length(B) * ((mean.total - mean.b)^2)

ss.between <- 
  length(A)*((mean.total - mean.a)^2) + 
  length(B)*((mean.total - mean.b)^2)
ss.between

# df between group은 연구에 사용된 
# 그룹의 숫자에서 1을 뺀 숫자
df.between <- 2 - 1

# 이 그룹 간 차이에 기인하는 분산 값은
ms.between <- ss.between / df.between

# 한편 아래 ss.a 와 ss.b는 각 그룹 내의 
# 분산을 알아보기 위한 방법
ss.a <- var(A) * df.a
ss.b <- var(B) * df.b
ss.within <- ss.a + ss.b

df.a <- length(A)-1
df.b <- length(B)-1
df.within <- df.a + df.b

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

# 여기까지 우리는 
# 전체분산
# 그룹간분산
# 그룹내분산 값을 
# 구한 것

# ms.between은 그룹의 차이때문에 생긴
# 분산으로 IV 혹은 treatment 때문에 생기는
# 차이에 기인하는 분산이다. 

# ms.within은 각 그룹 내부에서 일어나는 분산이므로
# (variation이므로) 연구자의 관심사와는 상관이 없이
# 나타나는 random한 차이를 알려주는 분산이라고 하면 

# t test 때와 마찬가지로 
# 그룹의 차이 / 랜덤 차이를 (에러 -> 분산은 에러라고도 했다)
# 구해볼 수 있다. 

# 즉, 그룹갑분산은 사실 = diff (between groups) 
# 그리고 그룹내 분산은 사실 = re
# 따라서 우리는 위 둘 간의 비율을 t test와 같이 
# 살펴볼 수 있다

# 이것을 f.calculated 이라고 하고
f.calculated <- ms.between / ms.within

# 이 값을 출력해 본다
f.calculated

# 이 계산은 차이와 랜덤에러의 비율이
# df에 따라서 얼마나 되어야 그 차이가
# 충분히 큰 것인지를 판단하기 위해서
# 쓰인다. 여기서 df에는 두 가지 종류가
# 있다. df.between 그리고 df.within
# 이 정보를 이용하여 f.calculated 값에 대한 
# probability를 구할 수 있다 (t-test때 처럼)

# percentage of f distribution with
# df1 and df2 option
# 이는 그림의 왼쪽을 나타내므로 
# 차이가 점점 커지게 되는 오른쪽을 
# 계산하기 위해서는 1-x를 취한다
f.calculated.pvalue <- 1-pf(f.calculated, df1=df.between, df2=df.within)
f.calculated.pvalue
# 한편,  t test를 했었을 때 (A, B 그룹을 가지고 independent 
# samples t-test를) 아웃 풋은 
t.result 

# 그리고 f 계산에서의 p value는 t test에서의 p.value와 같다
f.calculated.pvalue
t.result$p.value

# 또한
# 여기엣 t 값은 t.result$statistic 으로 프린트아웃할 수 있다
# 이 값이 2.33333 이었다
t.result$statistic
# 혹은 우리가 계산한 값이었던 
# t.calculated
t.calculated

# 그런데 위의 값을 제곱한 값이 바로 f.calculated 값
f.calculated
t.calculated
t.calculated^2

# 혹은 f.calculated 값을 제곱근한 값이 t.calculated
f.calculated
sqrt(f.calculated) 
t.calculated

# 또한 우리는 계산한 ms.total의 분자 부분인 ss.total은 
# 독립변인이 영향을 주는 ss.between 값과
# 독립변인이 있음에도 불구하고 여전히 남는 ss.within
# 으로 이루어짐을 아래처럼 확인한다.
ss.total
ss.between
ss.within
ss.total 
ss.between + ss.within

# 한편 df는 
# df.total  30 - 1
df.total 
df.between
df.within
df.total 
df.between + df.within

# 한 편
# 아래는 위의 분산을 쪼개서 분석하는 펑션인 
# aov를 이용해서 통계분석을 한 결과이다.

A 
B
comb <- stack(list(a=A, b=B))
comb
colnames(comb)[1] <- "values"
colnames(comb)[2] <- "group"
comb

a.res <- aov(values ~ group, data=comb)
a.res.sum <- summary(a.res)
a.res.sum
# 위에서 F value는 5.444 
# 그리고 전체적인 아웃풋을 보면 
# Df group 과 Df Residuals
# Sum Sq group 과 Residuals
# Mean Sq (MS) group 과 MS residuals

# MS.group = 
ms.between

# MS.within = 
ms.within

# F value
ms.between / ms.within 
f.calculated

# 아래는 기존의 아웃풋에서 확인하는 것
str(a.res.sum)
a.res.sum[[1]][1,4]
sqrt(a.res.sum[[1]][1,4])
t.result$statistic