# z-test 정리
# 
rnorm2 <- function(n,mean,sd) { mean+sd*scale(rnorm(n)) } 
set.seed(101)
# 평균이 80, 표준편차가 10인 모집단을 가정
# 아래 sa 샘플은 모집단의 구성원에게 머리가 
# 좋아지는 약을 투여한 후에 얻은 수학시험
# 점수
# 약의 효과가 있었을까?
p.mean <- 80
p.sigma <- 10 
n <- 100
sa <- rnorm2(n, 83, 10)

# 가설? 
# 영가설?
# confidence interval with 95% certainty?

# 83 != 80
# If 83 = 80 ?
# se for population (when sample size is 100)
se <- sqrt(p.sigma^2/n)
se
se2 <- se*2
# se를 위로 두 개 아래로 두 개 쓴 범위는
p.mean + c(-se2, +se2)
# 영가설이 맞다면 백중 95의 확률로 
# 샘플의 평균이 78-82 사이에서 
# 나와야 함 (왜냐하면 약을 먹은 샘플
# 100명도 모집단의 평범한 사람들일테니까)

# 그런데 샘플의 평균이 83이므로 영가설을 
# 부정함. 이 샘플과 모집단의 평균은 
# 다름. 즉, 샘플의 평균이 모집단에 나오기
# 힘든 평균이므로 이 샘플은 다른 
# 모집단에서 나온 샘플이라고 생각하게 됨
# 따라서 연구가설을 채택하여 
# 두 집단 간에 차이가 있다는 것을 받아들임


# 다른 방법. se를 단위로 하는 표준점수는?
sa.zs <- (mean(sa)-p.mean)/se
sa.zs
# 즉 z-score는 sa.zs 
# 이렇게 따지면 z-score는 +-2 범위 밖이면 
# 영가설을 부정하고 연구가설을 채택할 수 
# 있음
# 우리는 이 점수를 가지고 pnorm 을 구해
# 볼 수 있음
pnorm(sa.zs, lower.tail = F)
# 양방향 테스트일 경우
pnorm(sa.zs, lower.tail = F) * 2
# 즉, 만약에 영가설이 맞다면 
# 샘플 평균 83점 혹은 이에 해당하는 낮은 점수가 
# 나올 확률은 [1] 0.002699796 라는 것. 
# 이것은 5/100보다 작으므로 영가설을 부정하고 
# 연구가설을 채택

# 옛날 방식
qnorm(0.05/2)
abs(qnorm(0.05/2))
sa.zs

# statistical test 
# 두 집단 간의 평균 차이 / random error
# 로 나누어 주는 것
# se = random error 


# t-test 
# sample의 숫자가 작을 때 30미만일 
# 때 구하는 z score 는 부 정확 
# 즉, 이 z score로 구하는 probability도
# 부정확하여 보정이 필요함
# 이 보정되는 probability는 샘플의 크기에
# 따라서 달라짐

# 따라서 우리는 normal distribution에서
# probability를 구하지 않고 (pnorm을 사
# 용하지 않고), 따로 계산하여 구한 
# distribution을 사용함. 이것이 t distribution
# pnorm 처럼 pt를 사용함

set.seed(102)
n <- 16
sb <- rnorm2(n, 86, 10)

se <- sqrt(p.sigma^2/n)
se
se2 <- se*2
p.mean + c(-se2, +se2)

sb.tscore <- (mean(sb)-p.mean)/se
sb.tscore
pt(sb.tscore, df=n-1, lower.tail = F) * 2

# 옛날방식
qt(0.05/2, df=n-1)
qt(0.05/2, df=n-1, lower.tail = F)
sb.tscore

# 샘플의 숫자가 클 때에는?
set.seed(103)
n <- 2500
sc <- rnorm2(n, 80.75, 10)

se <- sqrt(p.sigma^2/n)
se
se2 <- se*2
p.mean + c(-se2, +se2)

sc.tscore <- (mean(sc)-p.mean)/se
sc.tscore
pt(sc.tscore, df = n-1, lower.tail = F) * 2
pnorm(sc.tscore, lower.tail = F) * 2


# 따라서 샘플숫자가 클 때에는 z-test 
# 대신 t-test를 해도 됨. 따라서, 
# 모든 경우 t-test를 하게 됨

# 알파와 베타
# type 1 error 와 type 2 error

set.seed(102)
n <- 16
sb1 <- rnorm2(n, 86, 10)

se <- sqrt(p.sigma^2/n)
se
se2 <- se*2
p.mean + c(-se2, +se2)

sb1.tscore <- (mean(sb1)-p.mean)/se
sb1.tscore
pt(sb1.tscore, df=n-1, lower.tail = F) * 2

# 위에서 범할 수 있는 오류는?

# 그렇다면 아래는 
set.seed(102)
n <- 16
sb2 <- rnorm2(n, 84, 10)

se <- sqrt(p.sigma^2/n)
se
se2 <- se*2
p.mean + c(-se2, +se2)

sb2.tscore <- (mean(sb2)-p.mean)/se
sb2.tscore
pt(sb2.tscore, df=n-1, lower.tail = F) * 2

# see http://commres.net/wiki/types_of_error

######################################
## 두 그룹 간의 차이 비교 테스트
######################################
# see http://commres.net/wiki/t-test

set.seed(109)
ssz <- 100
ts.t1 <- rnorm2(ssz, 50, 5)
ts.t2 <- rnorm2(ssz, 52, 5)
ts.diff <- ts.t2-ts.t1
mean(ts.diff)
sd(ts.diff)
var(ts.diff)
mean(ts.t2) - mean(ts.t1)
m.diff <- mean(ts.diff)
se <- sd(ts.diff)/sqrt(ssz)
t.cal <- m.diff/se
t.cal
pt(t.cal, df=ssz-1, lower.tail = F) * 2
#
t.test(ts.t1, ts.t2, paired = T)

##################################
set.seed(109)
na <- 100
nb <- 100
df.a <- na-1
df.b <- nb-1
ts.a <- rnorm2(na, 103, 10)
ts.b <- rnorm2(nb, 100, 10)

m.diff <- mean(ts.a)-mean(ts.b)
m.diff
ss.a <- var(ts.a) * df.a
ss.b <- var(ts.b) * df.b
ss.a
ss.b

v.pool <- (ss.a + ss.b)/(df.a + df.b)
se <- sqrt((v.pool/na)+(v.pool/nb))

t.cal <- m.diff/se
t.cal
pt(t.cal, df=df.a + df.b, lower.tail = F) * 2
# 
t.test(ts.a, ts.b)

#########################################
set.seed(104)
n <- 30

tsk <- rnorm2(n, 106, 9)
mu <- 100
m.tsk <- mean(tsk)
sd.tsk <- sd(tsk)
m.tsk
sd.tsk

m.diff <- m.tsk - mu
se <- 9/sqrt(n)
m.diff
se

t.cal <- m.diff/se
t.cal
pt(t.cal, df=n-1, lower.tail = F) * 2
#
t.test(tsk, mu=100)

> # z-test
> # 
> rnorm2 <- function(n,mean,sd) { mean+sd*scale(rnorm(n)) } 
> set.seed(101)
> # 평균이 80, 표준편차가 10인 모집단을 가정
> # 아래 sa 샘플은 모집단의 구성원에게 머리가 
> # 좋아지는 약을 투여한 후에 얻은 수학시험
> # 점수
> # 약의 효과가 있었을까?
> p.mean <- 80
> p.sigma <- 10 
> n <- 100
> sa <- rnorm2(n, 83, 10)
> 
> # 가설? 
> # 영가설?
> # confidence interval with 95% certainty?
> 
> # 83 != 80
> # If 83 = 80 ?
> # se for population (when sample size is 100)
> se <- sqrt(p.sigma^2/n)
> se
[1] 1
> se2 <- se*2
> # se를 위로 두 개 아래로 두 개 쓴 범위는
> p.mean + c(-se2, +se2)
[1] 78 82
> # 영가설이 맞다면 백중 95의 확률로 
> # 샘플의 평균이 78-82 사이에서 
> # 나와야 함 (왜냐하면 약을 먹은 샘플
> # 100명도 모집단의 평범한 사람들일테니까)
> 
> # 그런데 샘플의 평균이 83이므로 영가설을 
> # 부정함. 이 샘플과 모집단의 평균은 
> # 다름. 즉, 샘플의 평균이 모집단에 나오기
> # 힘든 평균이므로 이 샘플은 다른 
> # 모집단에서 나온 샘플이라고 생각하게 됨
> # 따라서 연구가설을 채택하여 
> # 두 집단 간에 차이가 있다는 것을 받아들임
> 
> 
> # 다른 방법. se를 단위로 하는 표준점수는?
> sa.zs <- (mean(sa)-p.mean)/se
> sa.zs
[1] 3
> # 즉 z-score는 sa.zs 
> # 이렇게 따지면 z-score는 +-2 범위 밖이면 
> # 영가설을 부정하고 연구가설을 채택할 수 
> # 있음
> # 우리는 이 점수를 가지고 pnorm 을 구해
> # 볼 수 있음
> pnorm(sa.zs, lower.tail = F)
[1] 0.001349898
> # 양방향 테스트일 경우
> pnorm(sa.zs, lower.tail = F) * 2
[1] 0.002699796
> # 즉, 만약에 영가설이 맞다면 
> # 샘플 평균 83점이 나올 확률은 
> # [1] 0.002699796 라는 것. 이것은 5/100보다
> # 작으므로 영가설을 부정하고 연구가설을 채택
> 
> # 옛날 방식
> qnorm(0.05/2)
[1] -1.959964
> abs(qnorm(0.05/2))
[1] 1.959964
> sa.zs
[1] 3
> 
> # statistical test 
> # 두 집단 간의 평균 차이 / random error
> # 로 나누어 주는 것
> # se = random error 
> 
> 
> # t-test 
> # sample의 숫자가 작을 때 30미만일 
> # 때 구하는 z score 는 부 정확 
> # 즉, 이 z score로 구하는 probability도
> # 부정확하여 보정이 필요함
> # 이 보정되는 probability는 샘플의 크기에
> # 따라서 달라짐
> 
> # 따라서 우리는 normal distribution에서
> # probability를 구하지 않고 (pnorm을 사
> # 용하지 않고), 따로 계산하여 구한 
> # distribution을 사용함. 이것이 t distribution
> # pnorm 처럼 pt를 사용함
> 
> set.seed(102)
> n <- 16
> sb <- rnorm2(n, 86, 10)
> 
> se <- sqrt(p.sigma^2/n)
> se
[1] 2.5
> se2 <- se*2
> p.mean + c(-se2, +se2)
[1] 75 85
> 
> sb.tscore <- (mean(sb)-p.mean)/se
> sb.tscore
[1] 2.4
> pt(sb.tscore, df=n-1, lower.tail = F) * 2
[1] 0.02982493
> 
> # 옛날방식
> qt(0.05/2, df=n-1)
[1] -2.13145
> qt(0.05/2, df=n-1, lower.tail = F)
[1] 2.13145
> sb.tscore
[1] 2.4
> 
> # 샘플의 숫자가 클 때에는?
> set.seed(103)
> n <- 2500
> sc <- rnorm2(n, 80.75, 10)
> 
> se <- sqrt(p.sigma^2/n)
> se
[1] 0.2
> se2 <- se*2
> p.mean + c(-se2, +se2)
[1] 79.6 80.4
> 
> sc.tscore <- (mean(sc)-p.mean)/se
> sc.tscore
[1] 3.75
> pt(sc.tscore, df = n-1, lower.tail = F) * 2
[1] 0.0001808498
> pnorm(sc.tscore, lower.tail = F) * 2
[1] 0.0001768346
> 
> 
> # 따라서 샘플숫자가 클 때에는 z-test 
> # 대신 t-test를 해도 됨. 따라서, 
> # 모든 경우 t-test를 하게 됨
> 
> # 알파와 베타
> # type 1 error 와 type 2 error
> 
> set.seed(102)
> n <- 16
> sb1 <- rnorm2(n, 86, 10)
> 
> se <- sqrt(p.sigma^2/n)
> se
[1] 2.5
> se2 <- se*2
> p.mean + c(-se2, +se2)
[1] 75 85
> 
> sb1.tscore <- (mean(sb1)-p.mean)/se
> sb1.tscore
[1] 2.4
> pt(sb1.tscore, df=n-1, lower.tail = F) * 2
[1] 0.02982493
> 
> # 위에서 범할 수 있는 오류는?
> 
> # 그렇다면 아래는 
> set.seed(102)
> n <- 16
> sb2 <- rnorm2(n, 84, 10)
> 
> se <- sqrt(p.sigma^2/n)
> se
[1] 2.5
> se2 <- se*2
> p.mean + c(-se2, +se2)
[1] 75 85
> 
> sb2.tscore <- (mean(sb2)-p.mean)/se
> sb2.tscore
[1] 1.6
> pt(sb2.tscore, df=n-1, lower.tail = F) * 2
[1] 0.130445
> 
> # see http://commres.net/wiki/types_of_error
> 
> ######################################
> ## 두 그룹 간의 차이 비교 테스트
> ######################################
> # see http://commres.net/wiki/t-test
> 
> set.seed(109)
> ssz <- 100
> ts.t1 <- rnorm2(ssz, 50, 5)
> ts.t2 <- rnorm2(ssz, 52, 5)
> ts.diff <- ts.t2-ts.t1
> mean(ts.diff)
[1] 2
> sd(ts.diff)
[1] 7.010937
> var(ts.diff)
         [,1]
[1,] 49.15323
> mean(ts.t2) - mean(ts.t1)
[1] 2
> m.diff <- mean(ts.diff)
> se <- sd(ts.diff)/sqrt(ssz)
> t.cal <- m.diff/se
> t.cal
[1] 2.852686
> pt(t.cal, df=ssz-1, lower.tail = F) * 2
[1] 0.005278446
> #
> t.test(ts.t1, ts.t2, paired = T)

	Paired t-test

data:  ts.t1 and ts.t2
t = -2.8527, df = 99, p-value = 0.005278
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -3.3911219 -0.6088781
sample estimates:
mean difference 
             -2 

> 
> ##################################
> set.seed(109)
> na <- 100
> nb <- 100
> df.a <- na-1
> df.b <- nb-1
> ts.a <- rnorm2(na, 103, 10)
> ts.b <- rnorm2(nb, 100, 10)
> 
> m.diff <- mean(ts.a)-mean(ts.b)
> m.diff
[1] 3
> ss.a <- var(ts.a) * df.a
> ss.b <- var(ts.b) * df.b
> ss.a
     [,1]
[1,] 9900
> ss.b
     [,1]
[1,] 9900
> 
> v.pool <- (ss.a + ss.b)/(df.a + df.b)
> se <- sqrt((v.pool/na)+(v.pool/nb))
> 
> t.cal <- m.diff/se
> t.cal
        [,1]
[1,] 2.12132
> pt(t.cal, df=df.a + df.b, lower.tail = F) * 2
           [,1]
[1,] 0.03513866
> # 
> t.test(ts.a, ts.b)

	Welch Two Sample t-test

data:  ts.a and ts.b
t = 2.1213, df = 198, p-value = 0.03514
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2111461 5.7888539
sample estimates:
mean of x mean of y 
      103       100 

> 
> #########################################
> set.seed(104)
> n <- 30
> 
> tsk <- rnorm2(n, 106, 9)
> mu <- 100
> m.tsk <- mean(tsk)
> sd.tsk <- sd(tsk)
> m.tsk
[1] 106
> sd.tsk
[1] 9
> 
> m.diff <- m.tsk - mu
> se <- 9/sqrt(n)
> m.diff
[1] 6
> se
[1] 1.643168
> 
> t.cal <- m.diff/se
> t.cal
[1] 3.651484
> pt(t.cal, df=n-1, lower.tail = F) * 2
[1] 0.001021295
> #
> t.test(tsk, mu=100)

	One Sample t-test

data:  tsk
t = 3.6515, df = 29, p-value = 0.001021
alternative hypothesis: true mean is not equal to 100
95 percent confidence interval:
 102.6393 109.3607
sample estimates:
mean of x 
      106 

> 
> 
>