Differences

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

--- t-test_summing_up [2025/09/18 08:13] – [8] hkimscil
+++ t-test_summing_up [2025/09/18 08:45] (current) – [t-test summing up] hkimscil
@@ Line 1: / Line 1: @@
 ====== t-test summing up ======
 <code>
 rm(list=ls())
@@ Line 210: / Line 209: @@
 dat <- stack(comb)
 head(dat)
+tail(dat)
 m.tot <- mean(s.all)
@@ Line 216: / Line 216: @@
 ss.tot <- ss(s.all)
+bet.s1 <- (m.tot - m.s1)^2 * length(s1)
+bet.s2 <- (m.tot - m.s2)^2 * length(s2)
+ss.bet <- bet.s1 + bet.s2
 ss.s1 <- ss(s1)
 ss.s2 <- ss(s2)
+ss.wit <- ss.s1+ss.s2
-df.tot <- length(s.all)-1
+ss.tot
-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.bet + ss.wit
+df.tot <- length(s.all) - 1
-ss.tot
-library(dplyr)
-# df.bet <- length(unique(dat)) - 1
 df.bet <- nlevels(dat$ind) - 1
-df.wit <- df.s1+df.s2
+df.s1 <- length(s1)-1
+df.s2 <- length(s2)-1
+df.wit <- df.s1 + df.s2
+df.tot
 df.bet
 df.wit
 df.bet+df.wit
-df.tot
+ss.tot/df.tot
+ms.tot <- ss.tot/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)
@@ Line 273: / Line 259: @@
 t.cal.ts
-# this is anova after all.
+# the above is anova after all.
 m1 <- lm(dat$values~dat$ind, data = dat)
@@ Line 287: / Line 273: @@
 sum.m1$fstatistic[1]
 ms.bet/ms.wit
 </code>
 ====== t-test summing up output ======
@@ Line 316: / Line 301: @@
 </WRAP>
-==== 2 ====
+===== 2 =====
 <WRAP group>
@@ Line 357: / Line 342: @@
 </WRAP>
-==== 3 ====
+===== 3 =====
 <WRAP group>
 <WRAP column half>
@@ Line 396: / Line 381: @@
 </WRAP>
-==== 4 ====
+===== 4 =====
 <WRAP group>
 <WRAP column half>
@@ Line 439: / Line 424: @@
 </WRAP>
-==== 5 ====
+===== 5 =====
 <WRAP group>
 <WRAP column half>
@@ Line 505: / Line 490: @@
 </WRAP>
-==== 6 ====
+===== 6 =====
 <WRAP group>
 <WRAP column half>
@@ Line 587: / Line 572: @@
 </WRAP>
-==== 7 ====
+===== 7 =====
 <WRAP group>
 <WRAP column half>
@@ Line 639: / Line 624: @@
-==== 8 ====
+===== 8 t-test and anova =====
 <WRAP group>
 <WRAP column half>
 <code>
+> # 2 샘플 히스토그램을 이용해서
+> # anova 설명하기
+> s.all <- c(s1,s2)
+> mean(s.all) # 두 집단을 합친
+[1] 104.7682
+> # 종속변인의 평균 (아래 히스토
+> # 그램에서 검정색).
+> hist(s1, col='grey', breaks=30, xlim=c(50, 150),
++      main = "histogram of s1, s2")
+> hist(s2, col='yellow', breaks=30, add=TRUE)
+> abline(v=c(mean(s.all)),
++        col=c("black"), lwd=3)
+> abline(v=c(mean(s1), mean(s2)),
++        col=c("green", "red"), lwd=2)
+>
+</code>
+</WRAP>
+<WRAP column half>
+.....................................
+{{:pasted:20250918-081521.png}}
+</WRAP>
+</WRAP>
-rm(list=ls())
+===== 9 =====
-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 <- 100000
+<WRAP group>
-m.p <- 100
+<WRAP column half>
-sd.p <- 10
+<code>
+> comb <- data.frame(s1,s2)
+> dat <- stack(comb)
+> head(dat)
+     values ind
+106.26706  s1
+101.89180  s1
+  98.64593  s1
+  93.82977  s1
+  80.35249  s1
+100.47300  s1
+> tail(dat)
+       values ind
+111.33398  s2
+115.99752  s2
+109.85412  s2
+126.80207  s2
+  99.13395  s2
+130.01417  s2
+>
+</code>
+</WRAP>
+<WRAP column half>
+.....................................
+</WRAP>
+</WRAP>
+===== 10 =====
+<WRAP group>
+<WRAP column half>
+<code>
+> m.tot <- mean(s.all)
+> m.s1 <- mean(s1)
+> m.s2 <- mean(s2)
+>
+> ss.tot <- ss(s.all)
+> bet.s1 <- (m.tot - m.s1)^2 * length(s1)
+> bet.s2 <- (m.tot - m.s2)^2 * length(s2)
+> ss.bet <- bet.s1 + bet.s2
+> ss.s1 <- ss(s1)
+> ss.s2 <- ss(s2)
+> ss.wit <- ss.s1+ss.s2
+>
+> ss.tot
+[1] 25203.89
+> ss.bet
+[1] 4219.463
+> ss.wit
+[1] 20984.43
+> ss.bet+ss.wit
+[1] 25203.89
+>
+> df.tot <- length(s.all) - 1
+> df.bet <- nlevels(dat$ind) - 1
+> df.s1 <- length(s1)-1
+> df.s2 <- length(s2)-1
+> df.wit <- df.s1 + df.s2
+>
+> df.tot
+[1] 199
+> df.bet
+[1] 1
+> df.wit
+[1] 198
+> df.bet+df.wit
+[1] 199
+>
+</code>
+</WRAP>
+<WRAP column half>
+.............................
+{{:pasted:20250918-081521.png}}
+  * ss.between 은 두개의 집단과 그 집단의 평균으로 이루어진 분산에서의 ss값을 말한다.
+  * 그림에서 녹색선 위에 s1의 원소들이 모여있고
+  * 붉은색 선위에 s2원소들이 모여 있다고치면
+  * 이 분산은 그룹간 분산이 (between group variance) 된다.
+  * 이를 구하기 위해서 ss.bet의 공식을 사용한다.
+  * ss.within의 경우에는 간단히 s1과 s2 집단의 ss값을 구하는 것으로
+  * 스크립트 초반에 언급된 ss값을 이용해서 구한후 더한다.
+  * df.between은 그룹의 숫자에서 1을 뺀 숫자
+  * df.within은 각 그룹내의 구성원에서 1을 뺀 후, 모두 더한 숫자이다
+  * df.total은 전체 구성원에서 1을 빼거나, df.bet + df.wit  를 더한 숫자
+</WRAP>
+</WRAP>
-p1 <- rnorm2(N.p, m.p, sd.p)
+===== 11 =====
-mean(p1)
-sd(p1)
-p2 <- rnorm2(N.p, m.p+10, sd.p)
+<WRAP group>
-mean(p2)
+<WRAP column half>
-sd(p2)
+<code>
+> ss.tot/df.tot
+[1] 126.6527
+> ms.tot <- ss.tot/df.tot
+>
+> ms.bet <- ss.bet / df.bet
+> ms.wit <- ss.wit / df.wit
+>
+> f.cal <- ms.bet / ms.wit
+> f.cal
+[1] 39.81302
+> pf(f.cal, df.bet, df.wit, lower.tail = F)
+[1] 1.792613e-09
+>
+> f.test <-  aov(dat$values~ dat$ind, data = dat)
+> sum.f.test <- summary(f.test)
+> sum.f.test
+             Df Sum Sq Mean Sq F value   Pr(>F)
+dat$ind       1   4219    4219   39.81 1.79e-09 ***
+Residuals   198  20984     106
+---
+Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+>
+>
+> sum.f.test[[1]][1,4]
+[1] 39.81302
+> sum.f.test[[1]][1,5]
+[1] 1.792613e-09
+> sqrt(f.cal)
+[1] 6.309756
+> t.cal.ts
+[1] -6.309756
+>
+> # the above is anova after all.
+>
+</code>
+</WRAP>
+<WRAP column half>
+.............................
+  * 위에서 구한 ss 값들을 가지고 분산값을 구할 수 있다 (ms = mean square = variance)
+  * ms.between은 ss.bet/df.bet 값을 말하는 것으로 그룹이 나누어짐으로써 구할 수 있었던 분산이다. 이 부분은 그룹 때문에 생긴 (treatment때문에 생긴) 분산이다.
+  * ms.within은 ss.wit/df.wit 값을 말하는 것으로 그룹의 구성 때문에 구할 수 있었던 ss.between 를 제외한 나머지 ss값을 df.wit값으로 나누어준 값을 말한다. 이 부분이 random 파트이다.
+  * 이 둘의 비율을 F값이라고 하고, 이는 기본적으로
+  * group difference 를 random difference 로 나눈 값을 말하고 이는
+  * t-test의 difference/random error와 같은 형태를 갖는다
-hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5),
+  * 따라서, 분산을 이용해서 구한 F값은 표준편차를 이용해서 구한 t값의 제곱값을 갖게 된다.
-     main = "histogram of p1 and p2",)
+</WRAP>
-abline(v=mean(p1), col="black", lwd=3)
+</WRAP>
-hist(p2, breaks=50, add=TRUE, col = rgb(0, 0, 1, 0.5))
-abline(v=mean(p2), col="red", lwd=3)
-# but suppose that we do not
+===== 12 =====
-# know the parameter of p2.
-# But, the real effect of the
-# drug is 10 point improvement
-# (mean = 110).
-hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5),
+<WRAP group>
-     main = "histogram of p1",)
+<WRAP column half>
-abline(v=mean(p1), col="black", lwd=3)
+<code>
+> m1 <- lm(dat$values~dat$ind, data = dat)
+> sum.m1 <- summary(m1)
+> sum.m1
-# we take 100 random sample
+Call:
-# (from p2)
+lm(formula = dat$values ~ dat$ind, data = dat)
-s.size <- 100
-s2 <- sample(p2, s.size)
-mean(s2)
-sd(s2)
-# The sample statistics are
+Residuals:
-# our clue to prove the
+     Min       1Q   Median       3Q      Max
-# effectiveness of the drug.
+-31.4177  -5.7058   0.4976   6.2952  29.2586
-# Let's assume that the drug is
+Coefficients:
-# NOT effective. And we also assume
+            Estimate Std. Error t value Pr(>|t|)
-# the distribution of sample means
+(Intercept)  100.175      1.029   97.31  < 2e-16 ***
-# whose sample size is 100 (in this
+dat$inds2      9.186      1.456    6.31 1.79e-09 ***
-# case).
+---
+Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-se.z <- sqrt(var(p1)/s.size)
+Residual standard error: 10.29 on 198 degrees of freedom
-se.z <- c(se.z)
+Multiple R-squared:  0.1674,	Adjusted R-squared:  0.1632
+F-statistic: 39.81 on 1 and 198 DF,  p-value: 1.793e-09
-z.cal <- (mean(s2) - mean(p1)) / se.z
+>
-z.cal
+> sum.m1$r.square
-pnorm(z.cal, lower.tail = F) * 2
+[1] 0.1674131
+> ss.bet # 독립변인인 s1, s2의 그룹
-# So,
+[1] 4219.463
-iter <- 100000
+> # 구분으로 설명된 ss.tot 중의 일부
-means <- rep(NA, iter)
+> ss.tot # y 전체의 uncertainty
-# means <- c()
+[1] 25203.89
-for (i in 1:iter) {
+> ss.bet/ss.tot
-  m.of.s <- mean(sample(p1, s.size))
+[1] 0.1674131
-  # means <- append(means, m.of.s)
+>
-  means[i] <- m.of.s
+> sum.m1$fstatistic[1]
-}
+   value
+.81302
-hist(means,
+> ms.bet/ms.wit
-     xlim = c(mean(means)-3*sd(means),
+[1] 39.81302
-              mean(means)+10*sd(means)),
+>
-     col = rgb(0, 0, 0, 0))
+>
-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)
-se.s <- se(s2)
-t.cal.os <- (mean(s2) - mean(p1))/ se.s
-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
-t.out <- t.test(s2, mu=mean(p1))
-t.out
-t.cal.os
-p.t.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.out <- t.test(s1, s2, var.equal = T)
-t.out
-hist(s1, breaks=30,
-     col = rgb(1, 0, 0, 0.5))
-hist(s2, breaks=30, 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))
-diff
-se.ts
-diff/se.ts
-t.out$statistic
-####
-# 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
-head(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)
-# 2 샘플 히스토그램을 이용해서
-# anova 설명하기
-s.all <- c(s1,s2)
-mean(s.all) # 두 집단을 합친
-# 종속변인의 평균 (아래 히스토
-# 그램에서 검정색).
-hist(s1, col='grey', breaks=30, xlim=c(50, 150),
-     main = "histogram of s1, s2")
-hist(s2, col='yellow', breaks=30, add=TRUE)
-abline(v=c(mean(s.all)),
-       col=c("black"), lwd=3)
-abline(v=c(mean(s1), mean(s2)),
-       col=c("green", "red"), lwd=2)
-comb <- data.frame(s1,s2)
-dat <- stack(comb)
-head(dat)
-tail(dat)
-m.tot <- mean(s.all)
-m.s1 <- mean(s1)
-m.s2 <- mean(s2)
-ss.tot <- ss(s.all)
-bet.s1 <- (m.tot - m.s1)^2 * length(s1)
-bet.s2 <- (m.tot - m.s2)^2 * length(s2)
-ss.bet <- bet.s1 + bet.s2
-ss.s1 <- ss(s1)
-ss.s2 <- ss(s2)
-ss.wit <- ss.s1+ss.s2
-ss.tot
-ss.bet
-ss.wit
-df.tot <- length(s.all) - 1
-df.bet <- nlevels(dat$ind) - 1
-df.s1 <- length(s1)-1
-df.s2 <- length(s2)-1
-df.wit <- df.s1 + df.s2
-df.tot
-df.bet
-df.wit
-ss.tot/df.tot
-ms.tot <- ss.tot/df.tot
-ms.bet <- ss.bet / df.bet
-ms.wit <- ss.wit / df.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)
-sum.f.test <- summary(f.test)
-sum.f.test
-sum.f.test[[1]][1,4]
-sum.f.test[[1]][1,5]
-sqrt(f.cal)
-t.cal.ts
-# the above is anova after all.
-m1 <- lm(dat$values~dat$ind, data = dat)
-sum.m1 <- summary(m1)
-sum.m1
-sum.m1$r.square
-ss.bet # 독립변인인 s1, s2의 그룹
-# 구분으로 설명된 ss.tot 중의 일부
-ss.tot # y 전체의 uncertainty
-ss.bet/ss.tot
-sum.m1$fstatistic[1]
-ms.bet/ms.wit
 </code>
 </WRAP>
 <WRAP column half>
 ................................
+  * ss.tot 은 s1그룹과 s2그룹의 값을 한줄에 배열하여 구하는 종속변인의 (dat$values) ss값을 말한다. 즉, 독립변이 없을 때의 종속변인의 uncertainty를 말한다.
+  * ss.bet은 독립변인이 고려됨으로써 설명된 종속부분의 일부를 말한다.
+  * ss.tot = ss.bet + ss.wit
+  * total uncertainty (a) = treatment effect (b) + random effect
+  * b/a를 R square값이라고 부른다.
 </WRAP>
 </WRAP>