COMMunication
RESearch.NET

User Tools

Site Tools

You are here: index » anova_note
Trace: • anova_note • types_of_error
anova_note

Differences

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

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
anova_note [2025/12/05 15:02] – [output] hkimscilanova_note [2025/12/09 23:16] (current) – hkimscil
Line 1: Line 1:
 ====== with 2 levels ====== ====== with 2 levels ======
 +
 +t-test를 하는 상황 (2 sample independent t-test)
 <tabbed> <tabbed>
   * :anova note:code01   * :anova note:code01
   * *:anova note:output01   * *:anova note:output01
 </tabbed> </tabbed>
-t-test를 하는 상황 (2 sample independent t-test) 
  
-===== output ===== 
  
 ====== with more than 3 levels ====== ====== with more than 3 levels ======
Line 12: Line 12:
   * :anova note:code02   * :anova note:code02
   * *:anova note:output02   * *:anova note:output02
-</tabed> +</tabbed>
- +
-<code> +
-#  +
-# ANOVA test with 4 levels in IV  +
-# +
-rm(list=ls()) +
-rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) } +
-ss <- function(x) { +
-  sum((x-mean(x))^2) +
-} +
- +
-set.seed(11) +
-n <- 31 +
-na <- nb <- nc <- nd <- n +
-mean.a <- 98 +
-mean.b <- 99 +
-mean.c <- 102 +
-mean.d <- 103 +
- +
-A <- rnorm2(na, mean.a, sqrt(900/(na-1))) +
-B <- rnorm2(nb, mean.b, sqrt(900/(nb-1))) +
-C <- rnorm2(nc, mean.c, sqrt(900/(nc-1))) +
-D <- rnorm2(nd, mean.d, sqrt(900/(nd-1))) +
-ss(A) +
-var(A) +
- +
-# A combined group with group A and B +
-# We call it group total +
-# we can obtain its mean, variance, ss, df, etc. +
-#  +
-comb <- data.frame(A, B, C, D) +
-dat <- stack(comb) +
-head(dat) +
-colnames(dat)[1] <- "values" +
-colnames(dat)[2] <- "group" +
-head(dat) +
- +
-m.tot <- mean(dat$values) +
-m.a <- mean(A) +
-m.b <- mean(B) +
-m.c <- mean(C) +
-m.d <- mean(D) +
- +
-# 그룹 간의 차이에서 나타나는 분산 +
-# 수업시간에 설명을 잘 들을 것 +
-min.x <- min(dat$values) +
-max.x <- max(dat$values) +
-br <- seq(floor(min.x), ceiling(max.x), by = 1)  +
-# Example bin width of 1 +
- +
- +
-hist(A, breaks=br, +
-     xlim = c(min.x-5, max.x+5), col=rgb(1,1,1,0.5), +
-     main = "Histogram of 4 groups") +
-hist(B, breaks=br, add=T, col=rgb(1,1,0,.5)) +
-hist(C, breaks=br, add=T, col=rgb(1,.5,1,.5)) +
-hist(D, breaks=br, add=T, col=rgb(.5,1,1,.5)) +
- +
-abline(v = m.tot, lty=2, lwd=3, col="black")  +
-abline(v = m.a, lty=2, lwd=3, col="blue")  +
-abline(v = m.b, lty=2, lwd=3, col="green")  +
-abline(v = m.c, lty=2, lwd=3, col="red")  +
-abline(v = m.d, lty=2, lwd=3, col="purple") +
- +
-# variance를 ms라고 부르기도 한다 +
-var.tot <- var(dat$values) +
-ms.tot <- var.tot +
- +
-ss.tot <- ss(dat$values) +
-# mean.total 에서 그룹a의 평균까지의 차이를 구한 후 +
-# 이를 제곱하여 그룹 A 멤버의 숫자만큼 더한다 =  +
-# 즉, SS를 구하는 방법.  +
-# 전체평균에서 그룹평균을 뺀 것의 제곱을  +
-# 그룹 구성원 숫자만큼 더하는 것 +
-# 그리고 이들을 다시 모두 더하여  +
-# ss.between에 저장 +
-bet.ta <- (m.tot - m.a)^2 * length(A) +
-bet.tb <- (m.tot - m.b)^2 * length(B) +
-bet.tc <- (m.tot - m.c)^2 * length(C) +
-bet.td <- (m.tot - m.d)^2 * length(D) +
-ss.bet <- bet.ta +  +
-  bet.tb +  +
-  bet.tc +  +
-  bet.td +
- +
-ss.a <- ss(A) +
-ss.b <- ss(B) +
-ss.c <- ss(C) +
-ss.d <- ss(D) +
-ss.wit <- ss.a+ss.b+ss.c+ss.c +
- +
-ss.tot +
-ss.bet +
-ss.wit +
-ss.bet+ss.wit +
- +
-df.tot <- length(dat$values) - 1 +
-df.bet <- nlevels(dat$group) - 1 +
-df.wit <- (length(A)-1) +  +
-  (length(B)-1) +  +
-  (length(C)-1) +  +
-  (length(D)-1) +
-df.tot +
-df.bet +
-df.wit +
- +
-ms.tot <- ss.tot / df.tot +
-ms.bet <- ss.bet / df.bet +
-ms.wit <- ss.wit / df.wit +
- +
-# ms.between은 그룹의 차이때문에 생긴 +
-# 분산으로 IV 혹은 treatment 때문에 생기는 +
-# 차이에 기인하는 분산이고 +
- +
-# ms.within은 각 그룹 내부에서 일어나는 분산이므로 +
-# (variation이므로) 연구자의 관심사와는 상관이 없이 +
-# 나타나는 random한 분산이라고 하면  +
- +
-# t test 때와 마찬가지로  +
-# 그룹의 차이 / 랜덤 차이를 (에러 -> 분산은 에러라고도 했다) +
-# 구해볼 수 있다.  +
- +
-# 즉, 그룹갑분산은 사실 = diff (between groups)  +
-# 그리고 그룹내 분산은 사실 = re +
-# 따라서 우리는 위 둘 간의 비율을 t test와 같이  +
-# 살펴볼 수 있다 +
-# 이것을 f.calculated 이라고 하고 +
- +
-f.cal <- ms.bet / ms.wit +
-f.cal +
- +
-# 컴퓨터 계산이 쉬워지기 전에는 아래처럼 0.5 level  +
-# 에서의 f값을 구한 후 이것과 계산된 f값을 비교해봤었다.  +
-qf(.05, df1 = df.bet, df2 = df.wit, lower.tail = FALSE) +
-f.cal +
-# 위에서 f.calculated > qf값이므로 +
-# f.calculated 값으로 영가설을 부정하고 +
-# 연구가설을 채택하면 판단이 잘못일 확률이  +
-# 0.05보다 작다는 것을 안다.  +
-# 그러나 컴퓨터계산이 용이해지고서는 qf대신에 +
-# pf를 써서 f.cal 값에 해당하는 prob. level을  +
-# 알아본다.  +
- +
-# percentage of f distribution with +
-# df1 and df2 option +
-# 이는 그림의 왼쪽을 나타내므로  +
-# 차이가 점점 커지게 되는 오른쪽을  +
-# 계산하기 위해서는 1-x를 취한다 +
- +
-p.val <- pf(f.cal, df.bet, df.wit, lower.tail=F) +
-p.val +
- +
-f.dat <- aov(dat$values~dat$group, data=dat) +
-summary(f.dat) +
- +
-# graph 로 이해  +
-x <- rf(50000, df1 = df.bet, df2 = df.wit) +
-y.max <- max(df(x, df1=df.bet, df2=df.wit)) +
- +
-hist(x, +
-     breaks = "Scott", +
-     freq = FALSE, +
-     xlim = c(0, f.cal + 1), +
-     ylim = c(0, y.max + .3), +
-     xlab = "", +
-     main = paste("Histogram for a F-distribution with  +
-             df1 = ", df.bet,  +
-          ", df2 = ", df.wit,  +
-          ", F cal = ", round(f.cal,3),  +
-          ", p val = ", round(p.val,5)),  +
-          cex.main = 0.9) +
-curve(df(x, df1 = df.bet, df2 = df.wit),  +
-      from = 0, to = f.cal + 1, n = 5000,  +
-      col = "red", lwd = 2,  +
-      add = T) +
-abline(v=f.cal, col="blue", lwd=2, lty="dotted") +
- +
-f.cal +
-p.val +
-1 - p.val +
- +
-# Now check this +
-ss.tot +
-ss.bet +
-ss.wit +
-ss.tot  +
-ss.bet + ss.wit +
- +
-# 한편 df는  +
-# df.total  30 - 1 +
-df.tot  +
-df.bet +
-df.wit +
-df.tot  +
-df.bet + df.wit +
- +
- +
- +
-################################################## +
-a.res <- aov(values ~ group, data=dat) +
-a.res.sum <- summary(a.res) +
-a.res.sum +
-# 그러나 정확히 어떤 그룹에서 차이가 나는지는 판단해주지 않음  +
-pairwise.t.test(dat$values, dat$group, p.adj = "none") +
-# OR +
-pairwise.t.test(dat$values, dat$group, p.adj = "bonf") +
-pairwise.t.test(dat$values, dat$group, p.adj = "holm") +
- +
-# OR TukeyHSD(anova.output) +
-TukeyHSD(a.res) +
- +
-boxplot(dat$values~dat$group) +
- +
-f.cal +
-p.val +
- +
-boxplot(dat$values~dat$group, main="values by group", +
-        yaxt="n", xlab="value", horizontal=TRUE, +
-        col=terrain.colors(4)) +
-legend("topleft", inset=.05, title="group", +
-       c("A","B","C", "D"), fill=terrain.colors(4), horiz=TRUE) +
-abline(v=mean(dat$values), col="red", lwd=2) +
- +
-         +
-# how much IV explains the DV  +
-# in terms of SS? +
-r.square <- ss.bet / ss.tot +
-eta  <- r.square +
-eta +
-lm.res <- lm(dat$values~dat$group, data = dat) +
-summary(lm.res) +
-summary(a.res) +
- +
-</code> +
- +
-{{:pasted:20250918-210412.png}} +
-{{:pasted:20250918-210425.png}} +
-{{:pasted:20250918-210512.png}} +
- +
-===== output ===== +
-<code> +
-> #  +
-> # ANOVA test with 4 levels in IV  +
-> # +
-> rm(list=ls()) +
-> rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) } +
-> ss <- function(x) { +
-+     sum((x-mean(x))^2) +
-+ } +
->  +
-> set.seed(11) +
-> n <- 31 +
-> na <- nb <- nc <- nd <- n +
-> mean.a <- 98 +
-> mean.b <- 99 +
-> mean.c <- 102 +
-> mean.d <- 103 +
->  +
-> A <- rnorm2(na, mean.a, sqrt(900/(na-1))) +
-> B <- rnorm2(nb, mean.b, sqrt(900/(nb-1))) +
-> C <- rnorm2(nc, mean.c, sqrt(900/(nc-1))) +
-> D <- rnorm2(nd, mean.d, sqrt(900/(nd-1))) +
-> ss(A) +
-[1] 900 +
-> var(A) +
-     [,1] +
-[1,]   30 +
->  +
-> # A combined group with group A and B +
-> # We call it group total +
-> # we can obtain its mean, variance, ss, df, etc. +
-> #  +
-> comb <- data.frame(A, B, C, D) +
-> dat <- stack(comb) +
-> head(dat) +
-     values ind +
-1  96.21237   A +
-2 100.86294   A +
-3  89.24341   A +
-4  90.40224   A +
-5 109.53642   A +
-6  93.62876   A +
-> colnames(dat)[1] <- "values" +
-> colnames(dat)[2] <- "group" +
-> head(dat) +
-     values group +
-1  96.21237     A +
-2 100.86294     A +
-3  89.24341     A +
-4  90.40224     A +
-5 109.53642     A +
-6  93.62876     A +
->  +
-> m.tot <- mean(dat$values) +
-> m.a <- mean(A) +
-> m.b <- mean(B) +
-> m.c <- mean(C) +
-> m.d <- mean(D) +
->  +
-> # 그룹 간의 차이에서 나타나는 분산 +
-> # 수업시간에 설명을 잘 들을 것 +
-> min.x <- min(dat$values) +
-> max.x <- max(dat$values) +
-> br <- seq(floor(min.x), ceiling(max.x), by = 1)  +
-> # Example bin width of 1 +
->  +
->  +
-> hist(A, breaks=br, +
-+      xlim = c(min.x-5, max.x+5), col=rgb(1,1,1,0.5), +
-+      main = "Histogram of 4 groups") +
-> hist(B, breaks=br, add=T, col=rgb(1,1,0,.5)) +
-> hist(C, breaks=br, add=T, col=rgb(1,.5,1,.5)) +
-> hist(D, breaks=br, add=T, col=rgb(.5,1,1,.5)) +
->  +
-> abline(v = m.tot, lty=2, lwd=3, col="black")  +
-> abline(v = m.a, lty=2, lwd=3, col="blue")  +
-> abline(v = m.b, lty=2, lwd=3, col="green")  +
-> abline(v = m.c, lty=2, lwd=3, col="red")  +
-> abline(v = m.d, lty=2, lwd=3, col="purple") +
->  +
-> # variance를 ms라고 부르기도 한다 +
-> var.tot <- var(dat$values) +
-> ms.tot <- var.tot +
->  +
-> ss.tot <- ss(dat$values) +
-> # mean.total 에서 그룹a의 평균까지의 차이를 구한 후 +
-> # 이를 제곱하여 그룹 A 멤버의 숫자만큼 더한다 =  +
-> # 즉, SS를 구하는 방법.  +
-> # 전체평균에서 그룹평균을 뺀 것의 제곱을  +
-> # 그룹 구성원 숫자만큼 더하는 것 +
-> # 그리고 이들을 다시 모두 더하여  +
-> # ss.between에 저장 +
-> bet.ta <- (m.tot - m.a)^2 * length(A) +
-> bet.tb <- (m.tot - m.b)^2 * length(B) +
-> bet.tc <- (m.tot - m.c)^2 * length(C) +
-> bet.td <- (m.tot - m.d)^2 * length(D) +
-> ss.bet <- bet.ta +  +
-+     bet.tb +  +
-+     bet.tc +  +
-+     bet.td +
->  +
-> ss.a <- ss(A) +
-> ss.b <- ss(B) +
-> ss.c <- ss(C) +
-> ss.d <- ss(D) +
-> ss.wit <- ss.a+ss.b+ss.c+ss.c +
->  +
-> ss.tot +
-[1] 4127 +
-> ss.bet +
-[1] 527 +
-> ss.wit +
-[1] 3600 +
-> ss.bet+ss.wit +
-[1] 4127 +
->  +
-> df.tot <- length(dat$values) - 1 +
-> df.bet <- nlevels(dat$group) - 1 +
-> df.wit <- (length(A)-1) +  +
-+     (length(B)-1) +  +
-+     (length(C)-1) +  +
-+     (length(D)-1) +
-> df.tot +
-[1] 123 +
-> df.bet +
-[1] 3 +
-> df.wit +
-[1] 120 +
->  +
-> ms.tot <- ss.tot / df.tot +
-> ms.bet <- ss.bet / df.bet +
-> ms.wit <- ss.wit / df.wit +
->  +
-> # ms.between은 그룹의 차이때문에 생긴 +
-> # 분산으로 IV 혹은 treatment 때문에 생기는 +
-> # 차이에 기인하는 분산이고 +
->  +
-> # ms.within은 각 그룹 내부에서 일어나는 분산이므로 +
-> # (variation이므로) 연구자의 관심사와는 상관이 없이 +
-> # 나타나는 random한 분산이라고 하면  +
->  +
-> # t test 때와 마찬가지로  +
-> # 그룹의 차이 / 랜덤 차이를 (에러 -> 분산은 에러라고도 했다) +
-> # 구해볼 수 있다.  +
->  +
-> # 즉, 그룹갑분산은 사실 = diff (between groups)  +
-> # 그리고 그룹내 분산은 사실 = re +
-> # 따라서 우리는 위 둘 간의 비율을 t test와 같이  +
-> # 살펴볼 수 있다 +
-> # 이것을 f.calculated 이라고 하고 +
->  +
-> f.cal <- ms.bet / ms.wit +
-> f.cal +
-[1] 5.855556 +
->  +
-> # 컴퓨터 계산이 쉬워지기 전에는 아래처럼 0.5 level  +
-> # 에서의 f값을 구한 후 이것과 계산된 f값을 비교해봤었다.  +
-> qf(.05, df1 = df.bet, df2 = df.wit, lower.tail = FALSE) +
-[1] 2.680168 +
-> f.cal +
-[1] 5.855556 +
-> # 위에서 f.calculated > qf값이므로 +
-> # f.calculated 값으로 영가설을 부정하고 +
-> # 연구가설을 채택하면 판단이 잘못일 확률이  +
-> # 0.05보다 작다는 것을 안다.  +
-> # 그러나 컴퓨터계산이 용이해지고서는 qf대신에 +
-> # pf를 써서 f.cal 값에 해당하는 prob. level을  +
-> # 알아본다.  +
->  +
-> # percentage of f distribution with +
-> # df1 and df2 option +
-> # 이는 그림의 왼쪽을 나타내므로  +
-> # 차이가 점점 커지게 되는 오른쪽을  +
-> # 계산하기 위해서는 1-x를 취한다 +
->  +
-> p.val <- pf(f.cal, df.bet, df.wit, lower.tail=F) +
-> p.val +
-[1] 0.0009119191 +
->  +
-> f.dat <- aov(dat$values~dat$group, data=dat) +
-> summary(f.dat) +
-             Df Sum Sq Mean Sq F value   Pr(>F)     +
-dat$group     3    527   175.7   5.856 0.000912 *** +
-Residuals   120   3600    30.0                      +
---- +
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +
->  +
-> # graph 로 이해  +
-> x <- rf(50000, df1 = df.bet, df2 = df.wit) +
-> y.max <- max(df(x, df1=df.bet, df2=df.wit)) +
->  +
-> hist(x, +
-+      breaks = "Scott", +
-+      freq = FALSE, +
-+      xlim = c(0, f.cal + 1), +
-+      ylim = c(0, y.max + .3), +
-+      xlab = "", +
-+      main = paste("Histogram for a F-distribution with  +
-+              df1 = ", df.bet,  +
-+                   ", df2 = ", df.wit,  +
-+                   ", F cal = ", round(f.cal,3),  +
-+                   ", p val = ", round(p.val,5)),  +
-+      cex.main = 0.9) +
-> curve(df(x, df1 = df.bet, df2 = df.wit),  +
-+       from = 0, to = f.cal + 1, n = 5000,  +
-+       col = "red", lwd = 2,  +
-+       add = T) +
-> abline(v=f.cal, col="blue", lwd=2, lty="dotted") +
->  +
-> f.cal +
-[1] 5.855556 +
-> p.val +
-[1] 0.0009119191 +
-> 1 - p.val +
-[1] 0.9990881 +
->  +
-> # Now check this +
-> ss.tot +
-[1] 4127 +
-> ss.bet +
-[1] 527 +
-> ss.wit +
-[1] 3600 +
-> ss.tot  +
-[1] 4127 +
-> ss.bet + ss.wit +
-[1] 4127 +
->  +
-> # 한편 df는  +
-> # df.total  30 - 1 +
-> df.tot  +
-[1] 123 +
-> df.bet +
-[1] 3 +
-> df.wit +
-[1] 120 +
-> df.tot  +
-[1] 123 +
-> df.bet + df.wit +
-[1] 123 +
->  +
->  +
->  +
-> ################################################## +
-> a.res <- aov(values ~ group, data=dat) +
-> a.res.sum <- summary(a.res) +
-> a.res.sum +
-             Df Sum Sq Mean Sq F value   Pr(>F)     +
-group         3    527   175.7   5.856 0.000912 *** +
-Residuals   120   3600    30.0                      +
---- +
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 +
-> # 그러나 정확히 어떤 그룹에서 차이가 나는지는 판단해주지 않음  +
-> pairwise.t.test(dat$values, dat$group, p.adj = "none") +
- +
- Pairwise comparisons using t tests with pooled SD  +
- +
-data:  dat$values and dat$group  +
- +
-  A       B       C       +
-B 0.47366 -       -       +
-C 0.00478 0.03305 -       +
-D 0.00047 0.00478 0.47366 +
- +
-P value adjustment method: none  +
-> # OR +
-> pairwise.t.test(dat$values, dat$group, p.adj = "bonf") +
- +
- Pairwise comparisons using t tests with pooled SD  +
- +
-data:  dat$values and dat$group  +
- +
-  A      B      C      +
-B 1.0000 -      -      +
-C 0.0287 0.1983 -      +
-D 0.0028 0.0287 1.0000 +
- +
-P value adjustment method: bonferroni  +
-> pairwise.t.test(dat$values, dat$group, p.adj = "holm") +
- +
- Pairwise comparisons using t tests with pooled SD  +
- +
-data:  dat$values and dat$group  +
- +
-  A      B      C      +
-B 0.9473 -      -      +
-C 0.0239 0.0991 -      +
-D 0.0028 0.0239 0.9473 +
- +
-P value adjustment method: holm  +
->  +
-> # OR TukeyHSD(anova.output) +
-> TukeyHSD(a.res) +
-  Tukey multiple comparisons of means +
-    95% family-wise confidence level +
- +
-Fit: aov(formula = values ~ group, data = dat) +
- +
-$group +
-    diff        lwr      upr     p adj +
-B-A    1 -2.6246725 4.624673 0.8894474 +
-C-A    4  0.3753275 7.624673 0.0243602 +
-D-A    5  1.3753275 8.624673 0.0026368 +
-C-B    3 -0.6246725 6.624673 0.1415754 +
-D-B    4  0.3753275 7.624673 0.0243602 +
-D-C    1 -2.6246725 4.624673 0.8894474 +
- +
->  +
-> boxplot(dat$values~dat$group) +
->  +
-> f.cal +
-[1] 5.855556 +
-> p.val +
-[1] 0.0009119191 +
->  +
-> boxplot(dat$values~dat$group, main="values by group", +
-+         yaxt="n", xlab="value", horizontal=TRUE, +
-+         col=terrain.colors(4)) +
-> legend("topleft", inset=.05, title="group", +
-+        c("A","B","C", "D"), fill=terrain.colors(4), horiz=TRUE) +
-> abline(v=mean(dat$values), col="red", lwd=2) +
->  +
->  +
-> # how much IV explains the DV  +
-> # in terms of SS? +
-> r.square <- ss.bet / ss.tot +
-> eta  <- r.square +
-> eta +
-[1] 0.1276957 +
-> lm.res <- lm(dat$values~dat$group, data = dat) +
-> summary(lm.res) +
- +
-Call: +
-lm(formula = dat$values ~ dat$group, data = dat) +
- +
-Residuals: +
-     Min       1Q   Median       3Q      Max  +
--12.9462  -3.7708   0.0944   3.1225  13.2340  +
- +
-Coefficients: +
-            Estimate Std. Error t value Pr(>|t|)     +
-(Intercept)  98.0000     0.9837  99.620  < 2e-16 *** +
-dat$groupB    1.0000     1.3912   0.719 0.473664     +
-dat$groupC    4.0000     1.3912   2.875 0.004779 **  +
-dat$groupD    5.0000     1.3912   3.594 0.000474 *** +
---- +
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1+
  
-Residual standard error: 5.477 on 120 degrees of freedom 
-Multiple R-squared:  0.1277, Adjusted R-squared:  0.1059  
-F-statistic: 5.856 on 3 and 120 DF,  p-value: 0.0009119 
  
-> summary(a.res) 
-             Df Sum Sq Mean Sq F value   Pr(>F)     
-group         3    527   175.7   5.856 0.000912 *** 
-Residuals   120   3600    30.0                      
---- 
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
->  
->  
->  
-</code> 
anova_note.1764946924.txt.gz · Last modified: by hkimscil
Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki