Differences

This shows you the differences between two versions of the page.

--- note.w02 [2025/09/12 10:38] – [9] hkimscil
+++ note.w02 [2025/09/18 10:40] (current) – [T-test sum up] hkimscil
@@ Line 23: / Line 23: @@
 sd(p1)
-p2 <- rnorm2(N.p, m.p+5, sd.p)
+p2 <- rnorm2(N.p, m.p+20, sd.p)
 mean(p2)
 sd(p2)
@@ Line 31: / Line 31: @@
 var(p1)
+hist(p1)
 hist(p1, breaks=50, col = rgb(1, 1, 1, 0.5),
      main = "histogram of p1 and p2",)
@@ Line 41: / Line 41: @@
 hist(p1, breaks=50, col=rgb(0,.5,.5,.5))
 abline(v=mean(p1),lwd=2)
-abline(v=mean(p1)-sd(p1), lwd=2)
+abline(v=m.p1-sd.p1, lwd=2)
 abline(v=mean(p1)+sd(p1), lwd=2)
 abline(v=c(m.p1-2*sd.p1, m.p1+2*sd.p1), lwd=2, col='red')
@@ Line 60: / Line 60: @@
 pnorm(m.p1+3*sd.p1, m.p1, sd.p1) -
   pnorm(m.p1-3*sd.p1, m.p1, sd.p1)
+pnorm(121, 100, 10) - pnorm(85, 100, 10)
 m.p1
@@ Line 69: / Line 71: @@
 pnorm(1)-pnorm(-1)
 pnorm(2)-pnorm(-2)
-pnorm(3)-pnorm(3)
+pnorm(3)-pnorm(-3)
 -pnorm(-2)*2
@@ Line 154: / Line 156: @@
 ################################
-s.size <- 50
+s.size <- 10
 means.temp <- c()
@@ Line 199: / Line 201: @@
 se.z <- sqrt(var(p1)/s.size)
+se.z
 se.z <- c(se.z)
 se.z
@@ Line 244: / Line 247: @@
 c(hi2, hiz2)
 c(hi3, hiz3)
 hist(means, breaks=50,
@@ Line 270: / Line 272: @@
      col = rgb(1, 1, 1, .5))
 abline(v=mean(means), col="black", lwd=3)
-abline(v=m.sample.i.got, col='darkgreen', lwd=3)
+abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3),
+       col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"),
+       lwd=2)
+abline(v=m.sample.i.got, col='red', lwd=3)
 # what is the probablity of getting
@@ Line 301: / Line 306: @@
 sd(means)
-tmp <- mean(means) - (m.s.from.p2 - mean(means))
+m.k <- mean(s.from.p2)
+se.k <- sd(s.from.p2)/sqrt(s.size)
+tmp <- mean(means) - (m.s.from.p2
+                    - mean(means))
 tmp
@@ Line 308: / Line 318: @@
      col = rgb(1, 1, 1, .5))
 abline(v=mean(means), col="black", lwd=3)
+abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3),
+       col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"),
+       lwd=2)
 abline(v=m.s.from.p2, col='blue', lwd=3)
@@ Line 315: / Line 328: @@
 m.s.from.p2
 pnorm(m.s.from.p2, mean(p1), se.z, lower.tail = F)
+pnorm(m.s.from.p2, m.k, se.k, lower.tail = F)
 # then, what is the probabilty of getting
 # greater than m.sample.i.got and
@@ Line 323: / Line 336: @@
 abline(v=tmp, col='red', lwd=3)
 * pnorm(m.s.from.p2, mean(p1), se.z, lower.tail = F)
+* pnorm(m.s.from.p2, m.k, se.k, lower.tail = F)
 se.z
@@ Line 338: / Line 354: @@
 pt(z.cal, 49, lower.tail = F)*2
 t.test(s.from.p2, mu=mean(p1), var.equal = T)
 </code>
@@ Line 753: / Line 772: @@
 <WRAP column half>
 ...........................................................................
+[[:central limit theorem]]
+$$\overline{X} \sim \displaystyle \text{N} \left(\mu, \dfrac{\sigma^{2}}{n} \right)$$
+[[:hypothesis testing]]
+[[:types of error]]
+[[:r:types of error]]
 </WRAP>
 </WRAP>
@@ Line 870: / Line 895: @@
 <WRAP column half>
 ...........................................................................
+{{:pasted:20250915-084136.png}}
 </WRAP>
 </WRAP>
@@ Line 887: / Line 913: @@
 > m.s.from.p2 <- mean(s.from.p2)
 > m.s.from.p2
-[1] 119.0929
+[1] 116.0821
 >
 > se.s
-[1] 3.163228
+[1] 3.166623
 > se.z
 [1] 3.162278
 > sd(means)
-[1] 3.163228
+[1] 3.166623
 >
 > m.k <- mean(s.from.p2)
@@ Line 901: / Line 927: @@
 >
 > tmp <- mean(means) - (m.s.from.p2
-+                     - mean(means))
++                       - mean(means))
 > tmp
-[1] 80.90409
+[1] 83.92356
 >
 > hist(means, breaks=30,
@@ Line 909: / Line 935: @@
 +      col = rgb(1, 1, 1, .5))
 > abline(v=mean(means), col="black", lwd=3)
+> abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3),
++        col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"),
++        lwd=2)
 > abline(v=m.s.from.p2, col='blue', lwd=3)
 >
@@ Line 915: / Line 944: @@
 > # m.sample.i.got?
 > m.s.from.p2
-[1] 119.0929
+[1] 116.0821
 > pnorm(m.s.from.p2, mean(p1), se.z, lower.tail = F)
-[1] 7.816511e-10
+[1] 1.832221e-07
 > pnorm(m.s.from.p2, m.k, se.k, lower.tail = F)
 [1] 0.5
@@ Line 927: / Line 956: @@
 > abline(v=tmp, col='red', lwd=3)
 > 2 * pnorm(m.s.from.p2, mean(p1), se.z, lower.tail = F)
-[1] 1.563302e-09
+[1] 3.664442e-07
 >
 > 2 * pnorm(m.s.from.p2, m.k, se.k, lower.tail = F)
@@ Line 936: / Line 965: @@
 [1] 3.162278
 > sd(s.from.p2)/sqrt(s.size)
-[1] 3.35771
+[1] 3.209216
 > se.z.adjusted <- sqrt(var(s.from.p2)/s.size)
 > se.z.adjusted
-[1] 3.35771
+[1] 3.209216
 > 2 * pnorm(m.s.from.p2, mean(p1), se.z.adjusted,
 +           lower.tail = F)
-[1] 1.298387e-08
+[1] 5.408382e-07
 >
 > z.cal <- (m.s.from.p2 - mean(p1))/se.z.adjusted
 > z.cal
-[1] 5.686277
+[1] 5.011228
 > pnorm(z.cal, lower.tail = F)*2
-[1] 1.298387e-08
+[1] 5.408382e-07
 >
 >
 > pt(z.cal, 49, lower.tail = F)*2
-[1] 7.095934e-07
+[1] 7.443756e-06
 > t.test(s.from.p2, mu=mean(p1), var.equal = T)
@@ Line 958: / Line 987: @@
 data:  s.from.p2
-t = 5.6863, df = 9, p-value = 0.0002995
+t = 5.0112, df = 9, p-value = 0.0007277
 alternative hypothesis: true mean is not equal to 100
 percent confidence interval:
-.4972 126.6885
+.8224 123.3419
 sample estimates:
 mean of x
-.0929
+.0821
+>
->
 >
 </code>
@@ Line 973: / Line 1001: @@
 ...........................................................................
 {{:pasted:20250912-193249.png}}
+{{:pasted:20250915-083229.png}}
 </WRAP>
 </WRAP>
-====== T-test sum up ======
-<code>
-rm(list=ls())
-rnorm2 <- function(n,mean,sd){
-  mean+sd*scale(rnorm(n))
-}
-se <- function(sample) {
-  sd(sample)/sqrt(length(sample))
-}
-ss <- function(x) {
-  sum((x-mean(x))^2)
-}
-N.p <- 1000000
-m.p <- 100
-sd.p <- 10
-p1 <- rnorm2(N.p, m.p, sd.p)
-mean(p1)
-sd(p1)
-p2 <- rnorm2(N.p, m.p+10, sd.p)
-mean(p2)
-sd(p2)
-hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5),
-     main = "histogram of p1 and p2",)
-abline(v=mean(p1), col="black", lwd=3)
-hist(p2, breaks=50, add=TRUE, col = rgb(0, 0, 1, 0.5))
-abline(v=mean(p2), col="red", lwd=3)
-s.size <- 1000
-s2 <- sample(p2, s.size)
-mean(s2)
-sd(s2)
-se.z <- sqrt(var(p1)/s.size)
-se.z <- c(se.z)
-se.z.range <- c(-2*se.z,2*se.z)
-se.z.range
-mean(p1)+se.z.range
-mean(s2)
-z.cal <- (mean(s2) - mean(p1)) / se.z
-z.cal
-pnorm(z.cal, lower.tail = F) * 2
-z.cal
-# principles . . .
-# distribution of sample means
-iter <- 100000
-means <- c()
-for (i in 1:iter) {
-  m.of.s <- mean(sample(p1, s.size))
-  means <- append(means, m.of.s)
-}
-hist(means,
-     xlim = c(mean(means)-3*sd(means), mean(s2)+5),
-     col = rgb(1, 0, 0, 0.5))
-abline(v=mean(p1), col="black", lwd=3)
-abline(v=mean(s2),
-       col="blue", lwd=3)
-lo1 <- mean(p1)-se.z
-hi1 <- mean(p1)+se.z
-lo2 <- mean(p1)-2*se.z
-hi2 <- mean(p1)+2*se.z
-abline(v=c(lo1, hi1, lo2, hi2),
-       col=c("green","green", "brown", "brown"),
-       lwd=2)
-se.z
-c(lo2, hi2)
-pnorm(z.cal, lower.tail = F) * 2
-# Note that we use sqrt(var(s2)/s.size)
-# as our se value instread of
-# sqrt(var(p1)/s.size)
-# This is a common practice for R
-# In fact, some z.test (made by someone)
-# function uses the former rather than
-# latter.
-sqrt(var(p1)/s.size)
-se.z
-sqrt(var(s2)/s.size)
-se(s2)
-t.cal.os <- (mean(s2) - mean(p1))/ se(s2)
-z.cal <- (mean(s2) - mean(p1))/ se.z
-t.cal.os
-z.cal
-df.s2 <- length(s2)-1
-df.s2
-p.t.os <- pt(abs(t.cal.os), df.s2, lower.tail = F) * 2
-p.t.os
-t.out <- t.test(s2, mu=mean(p1))
-library(BSDA)
-z.out <- z.test(s2, p1, sigma.x = sd(s2), sigma.y = sd(p1))
-z.out$statistic # se.z 대신에 se(s2) 값으로 구한 z 값
-z.cal # se.z으로 (sqrt(var(p1)/s.size)값) 구한 z 값
-t.out$statistic # se(s2)를 분모로 하여 구한 t 값
-t.cal.os # se(s2)를 이용하여 손으로 구한 t 값
-# But, after all, we use t.test method regardless of
-# variation
-hist(means,
-     xlim = c(mean(means)-3*sd(means), mean(s2)+5),
-     col = rgb(1, 0, 0, 0.5))
-abline(v=mean(p1), col="black", lwd=3)
-abline(v=mean(s2),
-       col="blue", lwd=3)
-lo1 <- mean(p1)-se.z
-hi1 <- mean(p1)+se.z
-lo2 <- mean(p1)-2*se.z
-hi2 <- mean(p1)+2*se.z
-abline(v=c(lo1, hi1, lo2, hi2),
-       col=c("green","green", "brown", "brown"),
-       lwd=2)
-# difference between black and blue line
-# divided by
-# se(s2) (= random difference)
-# t.value
-mean(s2)
-mean(p1)
-diff <- mean(s2)-mean(p1)
-diff
-se(s2)
-diff/se(s2)
-t.cal.os
-########################
-# 2 sample t-test
-########################
-# 가정. 아래에서 추출하는 두
-# 샘플의 모집단의 파라미터를
-# 모른다.
-s1 <- sample(p1, s.size)
-s2 <- sample(p2, s.size)
-mean(s1)
-mean(s2)
-ss(s1)
-ss(s2)
-df.s1 <- length(s1)-1
-df.s2 <- length(s2)-1
-df.s1
-df.s2
-pooled.variance <- (ss(s1)+ss(s2))/(df.s1+df.s2)
-pooled.variance
-se.ts <- sqrt((pooled.variance/length(s1))+(pooled.variance/length(s2)))
-se.ts
-t.cal.ts <- (mean(s1)-mean(s2))/se.ts
-t.cal.ts
-p.val.ts <- pt(abs(t.cal.ts), df=df.s1+df.s2, lower.tail = F) * 2
-p.val.ts
-t.test(s1, s2, var.equal = T)
-se(s1)
-se(s2)
-mean(s1)+c(-se(s1)*2, se(s1)*2)
-mean(s2)+c(-se(s2)*2, se(s2)*2)
-mean(p1)
-mean(p2)
-hist(s1, breaks=50,
-     col = rgb(1, 0, 0, 0.5))
-hist(s2, breaks=50, add=T, col=rgb(0,0,1,1))
-abline(v=mean(s1), col="green", lwd=3)
-# hist(s2, breaks=50, add=TRUE, col = rgb(0, 0, 1, 0.5))
-abline(v=mean(s2), col="lightblue", lwd=3)
-diff <- mean(s1)-mean(s2)
-se.ts
-diff/se.ts
-####
-# repeated measure t-test
-# we can use the above case
-# pop paramter unknown
-# two consecutive measurement
-# for the same sample
-t1 <- s1
-t2 <- s2
-mean(t1)
-mean(t2)
-diff.s <- t1 - t2
-diff.s
-t.cal.rm <- mean(diff.s)/se(diff.s)
-t.cal.rm
-p.val.rm <- pt(abs(t.cal.rm), length(s1)-1, lower.tail = F) * 2
-p.val.rm
-t.test(s1, s2, paired = T)
-# create multiple histogram
-s.all <- c(s1,s2)
-mean(s.all)
-hist(s1, col='grey', breaks=50, xlim=c(50, 150))
-hist(s2, col='darkgreen', breaks=50, add=TRUE)
-abline(v=c(mean(s.all)),
-       col=c("red"), lwd=3)
-abline(v=c(mean(s1), mean(s2)),
-       col=c("black", "green"), lwd=3)
-comb <- data.frame(s1,s2)
-dat <- stack(comb)
-head(dat)
-m.tot <- mean(s.all)
-m.s1 <- mean(s1)
-m.s2 <- mean(s2)
-ss.tot <- ss(s.all)
-ss.s1 <- ss(s1)
-ss.s2 <- ss(s2)
-df.tot <- length(s.all)-1
-df.s1 <- length(s1)-1
-df.s2 <- length(s2)-1
-ms.tot <- var(s.all)
-ms.tot
-ss.tot/df.tot
-var(s1)
-ss.s1 / df.s1
-var(s2)
-ss.s2 / df.s2
-ss.b.s1 <- length(s1) * ((m.tot - m.s1)^2)
-ss.b.s2 <- length(s2) * ((m.tot - m.s1)^2)
-ss.bet <- ss.b.s1+ss.b.s2
-ss.bet
-ss.wit <- ss.s1 + ss.s2
-ss.wit
-ss.bet + ss.wit
-ss.tot
-library(dplyr)
-# df.bet <- length(unique(dat)) - 1
-df.bet <- nlevels(dat$ind) - 1
-df.wit <- df.s1+df.s2
-df.bet
-df.wit
-df.bet+df.wit
-df.tot
-ms.bet <- ss.bet / df.bet
-ms.wit <- ss.wit / df.wit
-ms.bet
-ms.wit
-f.cal <- ms.bet / ms.wit
-f.cal
-pf(f.cal, df.bet, df.wit, lower.tail = F)
-f.test <-  aov(dat$values~ dat$ind, data = dat)
-summary(f.test)
-sqrt(f.cal)
-t.cal.ts
-# this is anova after all.
-</code>
-====== output ======
-<code>
-> rm(list=ls())
-> rnorm2 <- function(n,mean,sd){
-+   mean+sd*scale(rnorm(n))
-+ }
-> se <- function(sample) {
-+   sd(sample)/sqrt(length(sample))
-+ }
-> ss <- function(x) {
-+   sum((x-mean(x))^2)
-+ }
->
-> N.p <- 1000000
-> m.p <- 100
-> sd.p <- 10
->
->
-> p1 <- rnorm2(N.p, m.p, sd.p)
-> mean(p1)
-[1] 100
-> sd(p1)
-[1] 10
->
-> p2 <- rnorm2(N.p, m.p+10, sd.p)
-> mean(p2)
-[1] 110
-> sd(p2)
-[1] 10
->
-> hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5),
-+      main = "histogram of p1 and p2",)
-> abline(v=mean(p1), col="black", lwd=3)
-> hist(p2, breaks=50, add=TRUE, col = rgb(0, 0, 1, 0.5))
-> abline(v=mean(p2), col="red", lwd=3)
->
-> s.size <- 1000
-> s2 <- sample(p2, s.size)
-> mean(s2)
-[1] 110.4892
-> sd(s2)
-[1] 10.36614
->
-> se.z <- sqrt(var(p1)/s.size)
-> se.z <- c(se.z)
-> se.z.range <- c(-2*se.z,2*se.z)
-> se.z.range
-[1] -0.6324555  0.6324555
->
-> mean(p1)+se.z.range
-[1]  99.36754 100.63246
-> mean(s2)
-[1] 110.4892
->
-> z.cal <- (mean(s2) - mean(p1)) / se.z
-> z.cal
-[1] 33.16976
-> pnorm(z.cal, lower.tail = F) * 2
-[1] 2.93954e-241
->
-> z.cal
-[1] 33.16976
->
-> # principles . . .
-> # distribution of sample means
-> iter <- 100000
-> means <- c()
-> for (i in 1:iter) {
-+   m.of.s <- mean(sample(p1, s.size))
-+   means <- append(means, m.of.s)
-+ }
->
-> hist(means,
-+      xlim = c(mean(means)-3*sd(means), mean(s2)+5),
-+      col = rgb(1, 0, 0, 0.5))
-> abline(v=mean(p1), col="black", lwd=3)
-> abline(v=mean(s2),
-+        col="blue", lwd=3)
-> lo1 <- mean(p1)-se.z
-> hi1 <- mean(p1)+se.z
-> lo2 <- mean(p1)-2*se.z
-> hi2 <- mean(p1)+2*se.z
-> abline(v=c(lo1, hi1, lo2, hi2),
-+        col=c("green","green", "brown", "brown"),
-+        lwd=2)
-> se.z
-[1] 0.3162278
-> c(lo2, hi2)
-[1]  99.36754 100.63246
-> pnorm(z.cal, lower.tail = F) * 2
-[1] 2.93954e-241
->
->
-> # Note that we use sqrt(var(s2)/s.size)
-> # as our se value instread of
-> # sqrt(var(p1)/s.size)
-> # This is a common practice for R
-> # In fact, some z.test (made by someone)
-> # function uses the former rather than
-> # latter.
->
-> sqrt(var(p1)/s.size)
-          [,1]
-[1,] 0.3162278
-> se.z
-[1] 0.3162278
->
-> sqrt(var(s2)/s.size)
-[1] 0.3278061
-> se(s2)
-[1] 0.3278061
->
-> t.cal.os <- (mean(s2) - mean(p1))/ se(s2)
-> z.cal <- (mean(s2) - mean(p1))/ se.z
-> t.cal.os
-[1] 31.99818
-> z.cal
-[1] 33.16976
->
-> df.s2 <- length(s2)-1
-> df.s2
-[1] 999
-> p.t.os <- pt(abs(t.cal.os), df.s2, lower.tail = F) * 2
-> p.t.os
-[1] 3.162139e-155
-> t.out <- t.test(s2, mu=mean(p1))
->
-> library(BSDA)
-필요한 패키지를 로딩중입니다: lattice
-다음의 패키지를 부착합니다: ‘BSDA’
-The following object is masked from ‘package:datasets’:
-    Orange
-경고메시지(들):
-패키지 ‘BSDA’는 R 버전 4.3.3에서 작성되었습니다
-> z.out <- z.test(s2, p1, sigma.x = sd(s2), sigma.y = sd(p1))
->
-> z.out$statistic # se.z 대신에 se(s2) 값으로 구한 z 값
-      z
-.9833
-> z.cal # se.z으로 (sqrt(var(p1)/s.size)값) 구한 z 값
-[1] 33.16976
->
-> t.out$statistic # se(s2)를 분모로 하여 구한 t 값
-       t
-.99818
-> t.cal.os # se(s2)를 이용하여 손으로 구한 t 값
-[1] 31.99818
->
-> # But, after all, we use t.test method regardless of
-> # variation
->
-> hist(means,
-+      xlim = c(mean(means)-3*sd(means), mean(s2)+5),
-+      col = rgb(1, 0, 0, 0.5))
-> abline(v=mean(p1), col="black", lwd=3)
-> abline(v=mean(s2),
-+        col="blue", lwd=3)
-> lo1 <- mean(p1)-se.z
-> hi1 <- mean(p1)+se.z
-> lo2 <- mean(p1)-2*se.z
-> hi2 <- mean(p1)+2*se.z
-> abline(v=c(lo1, hi1, lo2, hi2),
-+        col=c("green","green", "brown", "brown"),
-+        lwd=2)
->
-> # difference between black and blue line
-> # divided by
-> # se(s2) (= random difference)
-> # t.value
-> mean(s2)
-[1] 110.4892
-> mean(p1)
-[1] 100
-> diff <- mean(s2)-mean(p1)
-> diff
-[1] 10.4892
-> se(s2)
-[1] 0.3278061
-> diff/se(s2)
-[1] 31.99818
->
-> t.cal.os
-[1] 31.99818
->
-> ########################
-> # 2 sample t-test
-> ########################
-> # 가정. 아래에서 추출하는 두
-> # 샘플의 모집단의 파라미터를
-> # 모른다.
-> s1 <- sample(p1, s.size)
-> s2 <- sample(p2, s.size)
->
-> mean(s1)
-[1] 100.7138
-> mean(s2)
-[1] 109.8804
-> ss(s1)
-[1] 108288.9
-> ss(s2)
-[1] 100751
-> df.s1 <- length(s1)-1
-> df.s2 <- length(s2)-1
-> df.s1
-[1] 999
-> df.s2
-[1] 999
->
-> pooled.variance <- (ss(s1)+ss(s2))/(df.s1+df.s2)
-> pooled.variance
-[1] 104.6246
-> se.ts <- sqrt((pooled.variance/length(s1))+(pooled.variance/length(s2)))
-> se.ts
-[1] 0.4574376
-> t.cal.ts <- (mean(s1)-mean(s2))/se.ts
-> t.cal.ts
-[1] -20.03892
-> p.val.ts <- pt(abs(t.cal.ts), df=df.s1+df.s2, lower.tail = F) * 2
-> p.val.ts
-[1] 1.522061e-81
->
-> t.test(s1, s2, var.equal = T)
-	Two Sample t-test
-data:  s1 and s2
-t = -20.039, df = 1998, p-value < 2.2e-16
-alternative hypothesis: true difference in means is not equal to 0
-percent confidence interval:
- -10.06366  -8.26945
-sample estimates:
-mean of x mean of y
-.7138  109.8804
->
-> se(s1)
-[1] 0.3292374
-> se(s2)
-[1] 0.3175719
->
-> mean(s1)+c(-se(s1)*2, se(s1)*2)
-[1] 100.0553 101.3723
-> mean(s2)+c(-se(s2)*2, se(s2)*2)
-[1] 109.2452 110.5155
->
-> mean(p1)
-[1] 100
-> mean(p2)
-[1] 110
->
-> hist(s1, breaks=50,
-+      col = rgb(1, 0, 0, 0.5))
-> hist(s2, breaks=50, add=T, col=rgb(0,0,1,1))
-> abline(v=mean(s1), col="green", lwd=3)
-> # hist(s2, breaks=50, add=TRUE, col = rgb(0, 0, 1, 0.5))
-> abline(v=mean(s2), col="lightblue", lwd=3)
->
-> diff <- mean(s1)-mean(s2)
-> se.ts
-[1] 0.4574376
-> diff/se.ts
-[1] -20.03892
->
-> ####
-> # repeated measure t-test
-> # we can use the above case
-> # pop paramter unknown
-> # two consecutive measurement
-> # for the same sample
->
-> t1 <- s1
-> t2 <- s2
-> mean(t1)
-[1] 100.7138
-> mean(t2)
-[1] 109.8804
-> diff.s <- t1 - t2
-> diff.s
-   [1] -11.22445947  -3.01532631  -3.47460077 -14.77286374  -9.74734969   2.33544374  -1.53761723
-   [8]   5.16258644 -40.75608052  -4.94498515  13.05099703  -0.24847628  -6.26581498  -8.14116058
-  [15]   5.49152742   6.15799638  -4.65589524  -7.93104859 -15.28488198  -5.99182327 -15.35591992
-  [22]  -1.00704696  18.27285355 -12.09936359  -4.72393289  -4.98301049 -20.95452153   1.62363928
-  [29] -10.58244069 -21.50608157 -53.61230898  -3.85198479 -42.61929736  -6.80266370 -22.92704580
-  [36]   3.01745740 -19.37131113 -27.82551902 -10.05485425 -25.12701225 -12.93162558  -7.55706006
-  [43] -19.16855657  -4.81878797   1.64397602 -28.64658004  16.36241227   8.73170802 -11.56090742
-  [50]   3.21799642 -39.37233043  -9.96051946 -19.11232333 -34.53077051  -4.85780005  -9.52501389
-  [57]  -8.28743679 -38.33995044 -50.60884456  -3.43450084  -0.85381393 -13.30971467  10.13049966
-  [64]   8.65616917 -29.75453733 -25.40674843 -24.98197786 -12.92901371  15.80168803  -8.67599446
-  [71] -20.50728324 -19.37275012 -23.27866089 -11.74962053 -34.35317908 -26.10460199   8.59009957
-  [78] -24.79252799   3.09475727 -19.13505970 -11.72561867 -33.79775614  -6.00167910 -25.03263480
-  [85] -23.66447959  -8.54416282  -6.89905337 -10.45234583 -34.67182776 -37.20205500   6.60270378
-  [92] -21.22842221 -11.68774036  -8.71535203 -15.55746542   2.88009050  -5.51543509   1.21606420
-  [99]  19.63435733 -14.40578016 -11.24172079  10.21723258 -23.41564885  27.60247565  -8.28684078
- [106]  17.72472594  -0.29977586 -14.84142327 -20.25713391 -37.99518419 -27.68545647 -19.78976153
- [113] -10.23092427   0.40875267 -17.36077213 -24.53979674 -24.18070810 -24.13321556 -17.36615616
- [120] -31.50478963  -2.47101725 -22.14003910 -33.63875270 -19.91485505  -3.64251563 -30.62407901
- [127]  -7.91406849  -6.25389594   4.40651820 -19.08031290 -14.01489366 -31.30542657 -27.05443597
- [134]  -6.60642443   1.29892762 -11.73908399   1.96878666  -7.53666283 -42.78007247  -9.04715952
- [141]  14.05326537  -6.85724091 -16.02305236 -24.38581030  -9.91759245   4.75488243  -2.83181250
- [148]  -5.38371023 -19.62202451  -1.55000290 -14.86541899   2.95279749   3.82076747  -3.38868029
- [155]  -0.65101074  18.42861149 -20.55548459   3.02438240  21.23539589   3.32810800  -9.66847192
- [162]   2.48687983  -5.40571073  -1.53265391 -12.93542011  19.87564176 -10.69228781  12.80629134
- [169]   6.51132022 -16.96586244  10.88690774  -0.37382590 -14.69255590 -26.66941722  -3.67854181
- [176]   6.29443209 -10.77182585 -25.65187453   0.09825251   9.03908176  10.25414786  -0.45340800
- [183] -38.63379525 -16.58800530 -17.62672134   5.43887886  -4.95532070  -6.46372777  -3.14562140
- [190] -22.25276559   2.05941880 -44.33676979  -3.34190343 -19.04858920  -8.21990394 -25.53625536
- [197] -17.21120659  28.96404607   6.25767994  -6.17593254 -34.33503461   6.65350479  11.42897662
- [204]  -0.83715976 -28.46397824 -40.67262397  -6.01225907 -11.22598108 -10.92756008 -15.45671946
- [211]  -4.57060131 -27.19860432   5.68618678 -27.70257611 -27.77374648  -5.93312400  -9.37871992
- [218] -24.41623403  16.94244832 -30.46760860  -4.91996788   2.89031604   5.27074167 -23.91666746
- [225] -13.27091592  -7.99640540  15.26148582 -26.01138488 -28.57927092 -17.29274303 -17.07704891
- [232] -11.52528966  -5.50387909  -0.66159232 -11.50347650 -19.90680762 -11.09595230 -11.02710712
- [239] -13.34969773 -37.98584006  -0.95289265 -15.00431567  -1.22592809   7.40922588  11.45790664
- [246] -24.07983488   4.54606079  -0.54863357   6.56528626  -4.04491250 -17.13525622 -22.85976576
- [253] -12.30101864   7.01445235  -6.77058075 -10.69023765 -14.21289974   0.68743488 -22.13964282
- [260]  -4.93155960  -5.32992121  -9.24699990 -34.21542676 -28.10074867 -14.64483350 -21.94636738
- [267]   0.57190289 -32.90838279 -23.39341251  -8.52122572   7.61461839 -12.60688433  -8.17161329
- [274]  -7.73981345  -4.86979671  14.97509924 -17.25493992 -48.85010339  10.16448581  -4.34608694
- [281]  -4.73924884  19.07076764  -5.23571728 -28.73076387  -1.01530521  -1.14387890  -6.96197277
- [288] -10.10160879 -11.59352210   5.83294359 -10.03967522  -9.59761019 -15.88867160  -7.13643475
- [295]  -7.29391701 -31.93027109 -15.56526408  13.22678162 -11.18996097 -25.00719650 -12.64524490
- [302] -35.18133234  13.41791211 -10.87228845 -31.95732546 -18.32496165 -17.76804212 -14.54640557
- [309]  21.57859526   8.22556833 -11.51111550  -7.19669828 -28.98417620  -8.16801696  -3.70794736
- [316]  -6.13751897   8.57292816   0.83289680 -11.25989874  11.25387495 -21.86409747   8.11413189
- [323] -31.81200194 -12.67414744  -7.24837917  -9.76383734  -9.95068555 -17.43835347 -21.88852882
- [330]  -7.57724957 -39.97109963 -18.22405067  -3.35304894  -2.01422861  12.65855868  19.68084692
- [337] -17.60471709   0.17064467  -2.41707219  -8.99719317  -6.91454803   6.10676201  -6.24933210
- [344] -10.52672545  -8.28580388 -21.74741007 -20.10217608 -27.73408273 -15.76758951 -16.05537057
- [351]   7.28801425  14.00094463 -32.21492728 -26.60734855   5.34256687  -2.73678515 -11.37626522
- [358] -23.18199896 -18.58488442 -26.09162223  -8.11379506 -16.24314767  -7.63202196  -7.70819353
- [365]  15.27344158  -6.89298171 -14.69008686 -10.61871285 -15.97716535   2.04244654 -20.72377059
- [372]  -3.33710996  -8.70719328   1.11210610  17.29240572   3.03254095  -0.17471362  -0.03964601
- [379] -37.25886118   1.80635404 -18.13062850 -12.08353080   3.17736630  11.60663212  14.40230421
- [386] -10.17057116 -11.45984282 -11.29130871 -15.60613249   3.94447987  21.39005539  38.66131508
- [393] -21.77681110  12.35832123  20.97222557 -30.89221977  26.68094540  13.28471640   5.74956285
- [400]  -7.60656480  17.47154758   2.14422376 -25.38609321  -7.75483605  12.58051687  -0.12179751
- [407] -17.52668230   1.92572455  -0.40532799 -11.06357792 -11.11874742 -10.35835894 -18.01080311
- [414] -14.30821541 -10.88427284   3.35560021 -13.46512417   6.04047559 -14.68666762   1.68951313
- [421]   4.31329242   7.50742041 -26.79115307 -32.01096580  28.21960400   3.60371811  14.05530287
- [428] -23.35919604   4.65486577 -42.05164227   5.90447867 -22.06140361 -36.42899364 -37.85288612
- [435]  -7.63484440  -7.99917501   5.73410720 -17.24124704  25.18142781 -22.27250215   3.31690400
- [442] -19.75974381  13.54264541 -18.86771849 -37.61540899   1.88862667 -10.85162959 -15.90454635
- [449]  -2.21038773  -1.64057323  -5.05984647 -11.33623396  -3.23993046  -7.07393112 -12.97757803
- [456]   2.81189108 -14.48709633 -26.49136729 -27.92148265  -3.32407602 -32.81165304  12.60241883
- [463] -19.77773032  -7.21862902 -30.05644966 -34.95205067  11.36207603  -1.88732702  -0.76534355
- [470]  -2.96997060 -26.06871359 -35.94190403  15.37335724 -29.72925340 -16.82830601 -32.14796579
- [477]  -7.49787668 -18.90903047  -4.95224428 -15.73171387  13.33385011   6.12192173 -12.11076372
- [484] -18.62856464 -25.58573774 -11.07888550 -16.77332405  -0.02979250   0.28045921   6.30825053
- [491]   2.66879530 -12.26308670  -3.34380103   0.57984817  -6.07862084  -0.47454107   0.18277721
- [498]  -1.49326079 -10.71424083 -18.55815468 -12.12523481 -10.32886876   1.21618496 -26.03417587
- [505]   4.21652961  -7.19317598  -6.49276222 -24.15078462 -16.62200664  18.65907236 -14.55468004
- [512]  12.03275099  13.63985493   1.94955005   4.24574372  -1.87352631  -9.90476053   4.92904187
- [519] -27.58434557  -3.22307117 -17.00433045 -24.92647628  -9.86853681  -3.12201969  -3.23781294
- [526] -44.00039083 -11.40638340 -10.42772214  11.56731373  34.18935566   0.22493781 -21.75192741
- [533] -27.11548294 -26.69534077  11.86361692 -18.87587506 -19.76960989  -6.72305933 -16.37844743
- [540]  28.42329924  -9.05042270  -5.94334753   0.80277394 -10.71835822  -7.39049563  -2.78087519
- [547] -18.87000360  -2.56690042  -2.31189557 -15.13438755   2.47712830  -7.13675100  -0.16789790
- [554] -17.69850124 -11.92887098  -7.41497303 -12.51272098   8.45361690 -30.63879338 -22.19654933
- [561]   0.37289814 -27.30901665 -42.04696582 -12.04972054 -20.74642541  -7.29054494 -11.52182026
- [568]   8.38372199 -18.57905859   7.62453900 -15.52892307  -1.83154260 -16.24286276  -5.74980635
- [575]   5.80843070   5.89082648 -12.38211486 -30.94945977 -34.27141760 -16.68862227 -12.74228148
- [582] -12.87197389 -27.14630137 -14.53291112  -8.24059195  -5.93523740   3.67320111 -14.73114416
- [589]   7.12333661 -20.54493654 -16.23259278  -7.23953837  13.66853283 -13.35719133  -3.19469244
- [596]  -8.63453073  -4.38937208  -4.25683530  -3.01229288  -9.46716386 -14.49889042  -4.80320886
- [603]  -5.93786510  -5.58582436 -38.52292640  -7.31871320 -10.10261534   3.55750183  -0.27859803
- [610]  13.47634673 -29.68104705  -4.28539812   3.67080445  -5.84424964  13.18400307  -7.89086592
- [617]  -1.33920025  -0.95308433 -29.91743460 -18.66903033 -20.15406028  -0.89348972 -12.41231804
- [624]   8.08680667 -34.26161156 -33.98947051  -7.80666121 -10.20912967 -28.95896039  -5.89738802
- [631]   2.32928389 -38.93521867 -15.56868160  -6.70412659  -2.57471692   2.46451660 -12.81558476
- [638]  19.76309298   2.12194484  24.84734206 -13.55620029 -17.87794609 -19.54477100 -19.73182713
- [645] -42.02480771 -13.42077241   9.19843679 -15.49829521 -35.16885995 -17.77559877  -2.28384147
- [652] -17.68303368   6.17812431   0.66249967   5.00542555  -2.26766815 -11.73064973  -6.75727090
- [659] -19.48957914 -28.88885682 -10.59583907 -32.20537872 -11.95661953 -23.77887711 -23.51551361
- [666] -17.32443690 -10.86310515  -8.51040349 -27.57610667 -26.85875525  -3.70329222 -20.50863173
- [673] -12.41671968  -5.15745402 -19.62849573   4.85082387 -29.88557391 -20.24914582 -27.73847105
- [680]  -0.12605203  -9.20839167 -51.56244236   3.17764075  13.96786787 -12.74398310   6.86534410
- [687] -21.75000477  14.10169236 -24.69667641 -20.26619614 -11.67028168  -5.14496708 -21.84000650
- [694] -34.30010114  -4.30214907 -15.81158253  -2.54412477   0.17601622 -22.22290730  -6.51460318
- [701] -19.46561809   9.62212347  -5.62354822  16.70312068  -7.16879691  -5.77420998  -2.36157455
- [708] -27.91638644 -12.61381331   7.45329002  -1.78749631  10.24888993  -2.76665687  -6.47189694
- [715] -20.55376627   1.03372077  -6.75380336 -21.29024889 -16.12342903 -36.81337018   1.75482644
- [722] -14.38944775  -4.27006397 -10.21581755   1.97016866   1.50969462  32.31451580   3.32233756
- [729]  -7.85868267   8.83356066  -3.54004596 -17.21481071 -29.58350979 -22.72248706 -27.86169027
- [736] -20.25705972  -1.67627671 -23.02237081 -20.13752529   6.07661361  -0.84839297  19.31619624
- [743] -13.32818441  -1.51206927 -16.05469364 -18.19320869  21.19248327 -26.85398142   2.82896396
- [750]   1.90853566 -10.76451371  -0.16368097  -5.02703204 -23.15483742  -5.12822113  -4.84245502
- [757] -19.16011286 -13.91801221 -11.66649472  15.04676653 -20.54651422  32.67150526 -11.37626079
- [764] -14.75337241  -0.46891630 -12.09854921 -10.75658868 -18.03655441  10.95871312  27.68695247
- [771]  -1.22676076 -15.78897397 -13.68374038  16.98138996 -15.57048042 -17.53983895 -32.33929466
- [778] -34.01869977  -2.46227514  -6.56500408 -17.04103052 -17.24440339  -6.68381805  -8.43674456
- [785]  -3.88372407 -29.28134174 -40.86613320 -18.64259952  -2.78880196  11.13938536  14.40875929
- [792]  -6.52858972 -38.92485161 -23.57441915 -25.59877817  19.51852353 -20.06650023  -4.32313864
- [799]   4.62056706  -5.18094527  -0.76429848  -2.98392902 -18.50300318 -29.63778693 -23.63235389
- [806]   5.68017122  14.33220812 -13.31317005 -10.81641168   5.22445430 -35.46250519   1.56327962
- [813] -14.60867384  10.29147319 -13.28593538   6.63825469 -14.52606348  17.45903705 -38.05095694
- [820] -13.93270232 -20.21993468  -1.96308711   1.80444271  12.16855255   9.52956342 -14.88194384
- [827] -29.04193544 -24.04102844 -14.01878071  -2.03269506   2.67151865  -5.20017572 -41.14943705
- [834]   8.96661691  28.12018815  -2.37196235 -13.46669223 -15.34687871 -19.92157033  16.10283716
- [841]   5.71060454  -1.87210810  17.82634786   6.46299261 -23.56325888  -9.31538158  13.08900119
- [848]   8.81863004 -16.87823373   2.57469446 -19.81326240  -1.08297141 -15.99656489 -12.78570251
- [855]  -8.53943328 -37.51286174  -5.80934175  29.94051347 -12.29916397  -0.84174744  -5.89053659
- [862] -30.93593593  -6.24638974   6.71567898   0.33777483  -6.43007412  -6.10032287 -18.90969351
- [869]   6.22885535   2.29565188 -25.72416278  -4.48305502  -5.77922453 -13.55585021 -23.84825362
- [876] -14.65449874 -25.51320775 -35.73124575   2.27482359 -23.67720440   7.67981459 -30.63388731
- [883]  -2.12532769  -6.06248123 -14.67967251   2.92695069 -32.55242308 -19.80182640 -12.70340060
- [890]  -0.36473422 -24.33804299 -33.53505272   6.38777505  -7.65940679 -45.22813407  -5.91512961
- [897]   4.82722129 -32.55034135  -5.64002137  -1.85377692 -24.82298659 -11.91899896   5.90226019
- [904]  -6.67799556 -18.39929702 -14.70709248  20.16465530  -2.37785503 -19.40544013 -43.64259489
- [911]  -4.80310727 -20.67267587  -6.12960286  14.76051916 -24.94995895 -21.55367734   3.51347606
- [918] -21.82098554   4.68892318  22.32281743   3.01554647 -14.22391287 -28.44488042   0.32000549
- [925] -26.29548705 -28.39677088  -6.06084948  -2.71491964   1.69227810  -8.71016310  -6.16547536
- [932] -11.67413566 -11.59680714   9.40984647  -9.93669428 -17.84745893   2.04601218 -19.80104095
- [939]  -0.49341925   6.14760676 -22.21183010  13.50485022  -5.22057307 -17.82539558 -24.46518962
- [946]  13.48666595 -11.80484500 -20.85173165 -31.08852302  -4.75860232 -36.88054918   2.98370161
- [953]  -0.37405972  17.58168372   0.48972051   9.79846880 -45.38152568 -25.01098967  -0.10839639
- [960]   0.14270290  10.28628008  21.54101027 -17.49673999 -12.24470665 -11.99999059 -24.05651849
- [967] -18.49661342 -12.30226946  -6.43942483 -23.07187196  -1.29013458 -15.53084359   0.86159642
- [974] -11.26368153  -4.91450784  -9.73711163 -14.70304307  -4.39913543 -21.76712550  -7.49898974
- [981] -25.17421015 -10.35273557  -9.64669400 -14.19622872 -13.91668603 -24.13258717 -15.08519499
- [988]   1.35746984  10.40157841  -2.47480562 -35.30199191 -25.52554695  -2.31850569 -10.24616931
- [995] -22.27223290  -4.57167529  -7.75456863  -2.13869306 -17.30982789 -24.04030147
-> t.cal.rm <- mean(diff.s)/se(diff.s)
-> t.cal.rm
-[1] -19.82846
-> p.val.rm <- pt(abs(t.cal.rm), length(s1)-1, lower.tail = F) * 2
-> p.val.rm
-[1] 4.836389e-74
-> t.test(s1, s2, paired = T)
-	Paired t-test
-data:  s1 and s2
-t = -19.828, df = 999, p-value < 2.2e-16
-alternative hypothesis: true mean difference is not equal to 0
-percent confidence interval:
- -10.073732  -8.259379
-sample estimates:
-mean difference
-      -9.166555
->
-> # create multiple histogram
-> s.all <- c(s1,s2)
-> mean(s.all)
-[1] 105.2971
-> hist(s1, col='grey', breaks=50, xlim=c(50, 150))
-> hist(s2, col='darkgreen', breaks=50, add=TRUE)
-> abline(v=c(mean(s.all)),
-+        col=c("red"), lwd=3)
-> abline(v=c(mean(s1), mean(s2)),
-+        col=c("black", "green"), lwd=3)
->
-> comb <- data.frame(s1,s2)
-> dat <- stack(comb)
-> head(dat)
-     values ind
-  93.17788  s1
-103.00254  s1
-104.53388  s1
-  88.59698  s1
-105.67789  s1
-112.72657  s1
->
-> m.tot <- mean(s.all)
-> m.s1 <- mean(s1)
-> m.s2 <- mean(s2)
->
-> ss.tot <- ss(s.all)
-> ss.s1 <- ss(s1)
-> ss.s2 <- ss(s2)
->
-> df.tot <- length(s.all)-1
-> df.s1 <- length(s1)-1
-> df.s2 <- length(s2)-1
->
-> ms.tot <- var(s.all)
-> ms.tot
-[1] 125.5892
-> ss.tot/df.tot
-[1] 125.5892
->
-> var(s1)
-[1] 108.3973
-> ss.s1 / df.s1
-[1] 108.3973
->
-> var(s2)
-[1] 100.8519
-> ss.s2 / df.s2
-[1] 100.8519
->
-> ss.b.s1 <- length(s1) * ((m.tot - m.s1)^2)
-> ss.b.s2 <- length(s2) * ((m.tot - m.s1)^2)
-> ss.bet <- ss.b.s1+ss.b.s2
-> ss.bet
-[1] 42012.87
->
-> ss.wit <- ss.s1 + ss.s2
-> ss.wit
-[1] 209039.9
->
-> ss.bet + ss.wit
-[1] 251052.8
-> ss.tot
-[1] 251052.8
->
-> library(dplyr)
-다음의 패키지를 부착합니다: ‘dplyr’
-The following objects are masked from ‘package:stats’:
-    filter, lag
-The following objects are masked from ‘package:base’:
-    intersect, setdiff, setequal, union
-> # df.bet <- length(unique(dat)) - 1
-> df.bet <- nlevels(dat$ind) - 1
-> df.wit <- df.s1+df.s2
-> df.bet
-[1] 1
-> df.wit
-[1] 1998
-> df.bet+df.wit
-[1] 1999
-> df.tot
-[1] 1999
->
-> ms.bet <- ss.bet / df.bet
-> ms.wit <- ss.wit / df.wit
-> ms.bet
-[1] 42012.87
-> ms.wit
-[1] 104.6246
->
-> f.cal <- ms.bet / ms.wit
-> f.cal
-[1] 401.5582
-> pf(f.cal, df.bet, df.wit, lower.tail = F)
-[1] 1.522061e-81
->
->
-> f.test <-  aov(dat$values~ dat$ind, data = dat)
-> summary(f.test)
-              Df Sum Sq Mean Sq F value Pr(>F)
-dat$ind        1  42013   42013   401.6 <2e-16 ***
-Residuals   1998 209040     105
----
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
->
-> sqrt(f.cal)
-[1] 20.03892
-> t.cal.ts
-[1] -20.03892
->
-> # this is anova after all.
->
-</code>