# t-test 이해 확인
pre <- c(3,0,6,7,4,3,2,8,4)
post <- c(5,2,5,7,10,9,7,11,8)
mean(pre)
mean(post)
sd(pre)
sd(post)

diff.prepost <- pre-post
mean.diff <- mean(diff.prepost)
mean.diff 
sd.diff <- sd(diff.prepost)
sd.diff

#
# remember t test = diff / rand error
#
se <- sd.diff/sqrt(length(diff.prepost))
se                          
t.value <- mean.diff/se
t.value
df.value <- length(diff.prepost)-1
df.value

# pt function is for getting percentage of 
# t score with df value
pt(t.value, df=df.value)
2*pt(t.value, df=df.value)

qt(.975, df=df.value)
qt(.025, df=df.value)
test.output <- t.test(pre,post, paired=T)
test.output
# for the reference
# test.output$
str(test.output)



# from the quiz questions
# stu should understand the logic of the ttest 

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 <- 26 - 19

# 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

pooled.var <- (SSa+SSb)/(df.a+df.b)
se <- sqrt(pooled.var/n.a + pooled.var/n.b)
t.calculated <- diff/se
pooled.var
se
t.calculated

t.result <- t.test(A, B, var.equal = T)
t.result
t.result$statistic
t.result$p.value

p.value <- 2*pt(-t.result$statistic, df=df.a+df.b)
p.value 
str(t.result)
t.result$p.value

##

## different approach
# 

A
B
dat <- c(A,B)
dat

var.total <- var(dat)
df.total <- length(dat)-1
ss.total <- var.total*df.total
ss.total.check <- sum((dat-mean(dat))^2)
ss.total
ss.total.check
mean.total <- mean(dat)
mean.total

mean.a <- mean(A)
mean.b <- mean(B)

# mean.total 에서 그룹a의 평균까지의 차이를 구한 후
# 이를 제곱하여 A의 숫자만큼 더한다 = 
# 즉, SS를 구하는 방법. 
# 전체평균에서 그룹평균을 뺀 것의 제곱을 
# 그룹 구성원 숫자만큼 더하는 것


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

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

# Now check this
ss.total
ss.between
ss.within
ss.total == ss.between + ss.within

# 한편 df는 
# df.total  30 - 1
df.between <- 2-1  # 그룹숫자 - 1
df.a <- length(A)-1 # a 구성원 - 1
df.b <- length(B)-1 # b 구성원 - 1
df.within <- df.a + df.b

df.total 
df.between
df.within
df.total == df.between + df.within

# 분산을 구하는 방법은 SS/df 이므로
# 분산을 ms 라고 표기하면 우리는 
# ms.total, ms.between, ms.within을 구할 수 있다

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

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

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

# t test 때와 마찬가지로 
# 그룹의 차이 / 랜덤 차이를 (에러 -> 분산은 에러라고도 했다)
# 구해볼 수 있다. 이것을 f.calculated 이라고 하고
# 이를 프린트아웃 한다

f.calculated <- ms.between / ms.within
f.calculated

# 한편,  t test를 했었을 때 (A, B 그룹을 가지고 independent 
# samples t-test를) 아웃 풋은 
t.result 

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

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

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


# 한 편

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

> ##
> 
> ## different approach
> # 
> 
> A
 [1] 20.994218 31.148068 16.961481 27.240217 28.354539 38.331534
 [7] 31.914700 23.459605 35.361796 22.182136 30.847396 15.575648
[13] 41.264878  7.808831 22.026979 22.527973
> B
 [1] 12.941146 21.270062 13.235378  1.931364 19.232163 27.231465
 [7] 18.276359  7.308871 27.560815  7.799787 25.017185 19.639663
[13] 25.018756 25.302096 28.941002 23.293888
> dat <- c(A,B)
> dat
 [1] 20.994218 31.148068 16.961481 27.240217 28.354539 38.331534
 [7] 31.914700 23.459605 35.361796 22.182136 30.847396 15.575648
[13] 41.264878  7.808831 22.026979 22.527973 12.941146 21.270062
[19] 13.235378  1.931364 19.232163 27.231465 18.276359  7.308871
[25] 27.560815  7.799787 25.017185 19.639663 25.018756 25.302096
[31] 28.941002 23.293888
> 
> var.total <- var(dat)
> df.total <- length(dat)-1
> ss.total <- var.total*df.total
> ss.total.check <- sum((dat-mean(dat))^2)
> ss.total
[1] 2552
> ss.total.check
[1] 2552
> mean.total <- mean(dat)
> mean.total
[1] 22.5
> 
> mean.a <- mean(A)
> mean.b <- mean(B)
> 
> # mean.total 에서 그룹a의 평균까지의 차이를 구한 후
> # 이를 제곱하여 A의 숫자만큼 더한다 = 
> # 즉, SS를 구하는 방법. 
> # 전체평균에서 그룹평균을 뺀 것의 제곱을 
> # 그룹 구성원 숫자만큼 더하는 것
> 
> 
> length(A)*((mean.total - mean.a)^2)
[1] 196
> length(B)*((mean.total - mean.b)^2)
[1] 196
> ss.between <- 
+   length(A)*((mean.total - mean.a)^2) + 
+   length(B)*((mean.total - mean.b)^2)
> 
> ss.between
[1] 392
> 
> # 한편 ss.a 와 ss.b는 각 그룹 내의 
> # 분산을 알아보기 위한 방법
> ss.a <- var(A) * df.a
> ss.b <- var(B) * df.b
> ss.within <- ss.a + ss.b
> 
> # Now check this
> ss.total
[1] 2552
> ss.between
[1] 392
> ss.within
[1] 2160
> ss.total == ss.between + ss.within
[1] FALSE
> 
> # 한편 df는 
> # df.total  30 - 1
> df.between <- 2-1  # 그룹숫자 - 1
> df.a <- length(A)-1 # a 구성원 - 1
> df.b <- length(B)-1 # b 구성원 - 1
> df.within <- df.a + df.b
> 
> df.total 
[1] 31
> df.between
[1] 1
> df.within
[1] 30
> df.total == df.between + df.within
[1] TRUE
> 
> # 분산을 구하는 방법은 SS/df 이므로
> # 분산을 ms 라고 표기하면 우리는 
> # ms.total, ms.between, ms.within을 구할 수 있다
> 
> ms.total <- ss.total / df.total
> ms.between <- ss.between / df.between
> ms.within <- ss.within / df.within
> 
> # 위에서 ms.between은 그룹의 차이때문에 생긴
> # 분산으로 IV 혹은 treatment 때문에 생기는
> # 차이에 기인하는 분산이고
> 
> # ms.within은 각 그룹 내부에서 일어나는 분산이므로
> # (variation이므로) 연구자의 관심사와는 상관이 없이
> # 나타나는 random한 분산이라고 하면 
> 
> # t test 때와 마찬가지로 
> # 그룹의 차이 / 랜덤 차이를 (에러 -> 분산은 에러라고도 했다)
> # 구해볼 수 있다. 이것을 f.calculated 이라고 하고
> # 이를 프린트아웃 한다
> 
> f.calculated <- ms.between / ms.within
> f.calculated
[1] 5.444444
> 
> # 한편,  t test를 했었을 때 (A, B 그룹을 가지고 independent 
> # samples t-test를) 아웃 풋은 
> t.result 

	Two Sample t-test

data:  A and B
t = 2.3333, df = 30, p-value = 0.02652
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  0.8731826 13.1268174
sample estimates:
mean of x mean of y 
       26        19 

> 
> # 여기엣 t 값은 t.result$statistic 으로 프린트아웃할 수 있다
> # 이 값이 2.33333 이었다
> t.result$statistic
       t 
2.333333 
> # 혹은 우리가 계산한 값이었던 
> # t.calculated
> t.calculated
[1] 2.333333
> 
> # 그런데 위의 값을 제곱한 값이 바로 f.calculated 값
> f.calculated
[1] 5.444444
> t.calculated^2
[1] 5.444444
> 
> # 혹은 f.calculated 값을 제곱근한 값이 t.calculated
> sqrt(f.calculated) 
[1] 2.333333
> t.calculated
[1] 2.333333
> 
> 
> # 한 편
> 
> A 
 [1] 20.994218 31.148068 16.961481 27.240217 28.354539 38.331534
 [7] 31.914700 23.459605 35.361796 22.182136 30.847396 15.575648
[13] 41.264878  7.808831 22.026979 22.527973
> B
 [1] 12.941146 21.270062 13.235378  1.931364 19.232163 27.231465
 [7] 18.276359  7.308871 27.560815  7.799787 25.017185 19.639663
[13] 25.018756 25.302096 28.941002 23.293888
> comb <- stack(list(a=A, b=B))
> comb
      values ind
1  20.994218   a
2  31.148068   a
3  16.961481   a
4  27.240217   a
5  28.354539   a
6  38.331534   a
7  31.914700   a
8  23.459605   a
9  35.361796   a
10 22.182136   a
11 30.847396   a
12 15.575648   a
13 41.264878   a
14  7.808831   a
15 22.026979   a
16 22.527973   a
17 12.941146   b
18 21.270062   b
19 13.235378   b
20  1.931364   b
21 19.232163   b
22 27.231465   b
23 18.276359   b
24  7.308871   b
25 27.560815   b
26  7.799787   b
27 25.017185   b
28 19.639663   b
29 25.018756   b
30 25.302096   b
31 28.941002   b
32 23.293888   b
> colnames(comb)[1] <- "values"
> colnames(comb)[2] <- "group"
> comb
      values group
1  20.994218     a
2  31.148068     a
3  16.961481     a
4  27.240217     a
5  28.354539     a
6  38.331534     a
7  31.914700     a
8  23.459605     a
9  35.361796     a
10 22.182136     a
11 30.847396     a
12 15.575648     a
13 41.264878     a
14  7.808831     a
15 22.026979     a
16 22.527973     a
17 12.941146     b
18 21.270062     b
19 13.235378     b
20  1.931364     b
21 19.232163     b
22 27.231465     b
23 18.276359     b
24  7.308871     b
25 27.560815     b
26  7.799787     b
27 25.017185     b
28 19.639663     b
29 25.018756     b
30 25.302096     b
31 28.941002     b
32 23.293888     b
> 
> a.res <- aov(values ~ group, data=comb)
> a.res.sum <- summary(a.res)
> a.res.sum
            Df Sum Sq Mean Sq F value Pr(>F)  
group        1    392     392   5.444 0.0265 *
Residuals   30   2160      72                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> # 위에서 F value는 5.444 
> # 그리고 전체적인 아웃풋을 보면 
> # Df group 과 Df Residuals
> # Sum Sq group 과 Residuals
> # Mean Sq (MS) group 과 MS residuals
> 
> # MS.group = 
> ms.between
[1] 392
> 
> # MS.within = 
> ms.within
[1] 72
> 
> # F value
> ms.between / ms.within 
[1] 5.444444
> f.calculated
[1] 5.444444
> 
> # 아래는 기존의 아웃풋에서 확인하는 것
> str(a.res.sum)
List of 1
 $ :Classes ‘anova’ and 'data.frame':	2 obs. of  5 variables:
  ..$ Df     : num [1:2] 1 30
  ..$ Sum Sq : num [1:2] 392 2160
  ..$ Mean Sq: num [1:2] 392 72
  ..$ F value: num [1:2] 5.44 NA
  ..$ Pr(>F) : num [1:2] 0.0265 NA
 - attr(*, "class")= chr [1:2] "summary.aov" "listof"
> a.res.sum[[1]][1,4]
[1] 5.444444
> sqrt(a.res.sum[[1]][1,4])
[1] 2.333333
> t.result$statistic
       t 
2.333333 
>