Differences

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

--- note.w02 [2025/09/11 10:36] – hkimscil
+++ note.w02 [2025/09/12 19:39] (current) – [Sampling Distribution and z-test] 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 301: / Line 304: @@
 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 315: / Line 323: @@
 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 331: @@
 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 349: @@
 pt(z.cal, 49, lower.tail = F)*2
 t.test(s.from.p2, mu=mean(p1), var.equal = T)
 </code>
 ====== output ======
+===== 1 =====
+<WRAP group>
+<WRAP column half>
 <code>
 > rm(list=ls())
@@ Line 355: / Line 373: @@
 + }
 >
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 2 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > ################################
 > N.p <- 1000000
@@ Line 386: / Line 415: @@
 >
 >
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 3 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > hist(p1, breaks=50, col=rgb(0,.5,.5,.5))
 > abline(v=mean(p1),lwd=2)
@@ Line 413: / Line 454: @@
 [1] 0.9973002
 >
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 4 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > m.p1
 [1] 100
@@ Line 463: / Line 516: @@
 [1] 0.03593032
 >
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 5 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > z.p1 <- (p1-mean(p1))/sd(p1)
 > mean(z.p1)
@@ Line 493: / Line 558: @@
 >
 > #
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 6 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > hist(p1, breaks=50, col=rgb(.9,.9,.9,.9))
 > abline(v=mean(p1),lwd=2)
@@ Line 553: / Line 630: @@
 >
 >
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 7 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > ################################
 > s.size <- 50
@@ Line 568: / Line 657: @@
 [1]  98.76098  99.90935  99.29643  99.66014 101.93822
 >
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 8 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > iter <- 1000000
 > # means <- c()
@@ Line 600: / Line 701: @@
 +        lwd=2)
 >
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 9 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > # meanwhile . . . .
 > se.s
 [1] 1.414035
->
+# se.s = sd(means)
+# The below is from CLT
+# see http://commres.org/wiki/central limit theorem
+#
 > se.z <- sqrt(var(p1)/s.size)
 > se.z <- c(se.z)
@@ Line 645: / Line 762: @@
 > se.s
 [1] 1.414035
-> se.z
-[1] 1.414214
 >
-> # because CLT
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 10 =====
+<WRAP group>
+<WRAP column half>
+<code>
+> # because of CLT we can use the
+> # below instead of
+> # mean(means)+-se.s
+> #
 > loz1 <- mean(p1)-se.z
 > hiz1 <- mean(p1)+se.z
@@ Line 671: / Line 801: @@
 >
 >
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 11 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > hist(means, breaks=50,
 +      xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)),
@@ Line 694: / Line 836: @@
 [1]  96 104
 >
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 12 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > m.sample.i.got <- mean(means)+ 1.5*sd(means)
 > m.sample.i.got
@@ Line 726: / Line 880: @@
 [1] 102.1193
 >
+</code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+</WRAP>
+</WRAP>
+===== 13 =====
+<WRAP group>
+<WRAP column half>
+<code>
 > ### one more time
 > # this time, with a story
 > mean(p2)
-[1] 105
+[1] 120
 > sd(p2)
 [1] 10
@@ Line 735: / Line 901: @@
 > m.s.from.p2 <- mean(s.from.p2)
 > m.s.from.p2
-[1] 103.1283
+[1] 119.0929
 >
 > se.s
-[1] 1.414035
+[1] 3.163228
 > se.z
-[1] 1.414214
+[1] 3.162278
 > sd(means)
-[1] 1.414035
+[1] 3.163228
 >
-> 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
-[1] 96.86822
+[1] 80.90409
 >
 > hist(means, breaks=30,
@@ Line 758: / Line 929: @@
 > # m.sample.i.got?
 > m.s.from.p2
-[1] 103.1283
+[1] 119.0929
 > pnorm(m.s.from.p2, mean(p1), se.z, lower.tail = F)
-[1] 0.01348266
+[1] 7.816511e-10
->
+> pnorm(m.s.from.p2, m.k, se.k, lower.tail = F)
+[1] 0.5
 > # then, what is the probabilty of getting
 > # greater than m.sample.i.got and
@@ Line 769: / Line 941: @@
 > abline(v=tmp, col='red', lwd=3)
 > 2 * pnorm(m.s.from.p2, mean(p1), se.z, lower.tail = F)
-[1] 0.02696533
+[1] 1.563302e-09
+>
+> 2 * pnorm(m.s.from.p2, m.k, se.k, lower.tail = F)
+[1] 1
+>
 >
 > se.z
-[1] 1.414214
+[1] 3.162278
 > sd(s.from.p2)/sqrt(s.size)
-[1] 1.414296
+[1] 3.35771
 > se.z.adjusted <- sqrt(var(s.from.p2)/s.size)
 > se.z.adjusted
-[1] 1.414296
+[1] 3.35771
 > 2 * pnorm(m.s.from.p2, mean(p1), se.z.adjusted,
 +           lower.tail = F)
-[1] 0.02697421
+[1] 1.298387e-08
 >
 > z.cal <- (m.s.from.p2 - mean(p1))/se.z.adjusted
 > z.cal
-[1] 2.211891
+[1] 5.686277
 > pnorm(z.cal, lower.tail = F)*2
-[1] 0.02697421
+[1] 1.298387e-08
 >
 >
 > pt(z.cal, 49, lower.tail = F)*2
-[1] 0.03166797
+[1] 7.095934e-07
 > t.test(s.from.p2, mu=mean(p1), var.equal = T)
@@ Line 796: / Line 972: @@
 data:  s.from.p2
-t = 2.2119, df = 49, p-value = 0.03167
+t = 5.6863, df = 9, p-value = 0.0002995
 alternative hypothesis: true mean is not equal to 100
 percent confidence interval:
-.2861 105.9704
+.4972 126.6885
 sample estimates:
 mean of x
-.1283
+.0929
+>
 >
 </code>
+</WRAP>
+<WRAP column half>
+...........................................................................
+{{:pasted:20250912-193249.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>