Differences

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

--- c:ma:2023:lecturer_note [2023/09/17 19:48] – created hkimscil
+++ c:ma:2023:lecturer_note [2023/09/17 19:49] (current) – hkimscil
@@ Line 1: / Line 1: @@
 ====== ANOVA notes ======
 <code R>
-temp <- c(1,2)
+# ANOVA
+set.seed(1024)
+n <-30
+s1 <- round(rnorm(n, 5.6, 1),0)
+s2 <- round(rnorm(n, 5.9, 1),0)
+s3 <- round(rnorm(n, 6.5, 1),0)
+s1
+s2
+s3
+sam <- c(s1,s2,s3)
+sam
+group <- c(rep("g1",n), rep("g2",n), rep("g3",n))
+group
+sa <- data.frame(sam,group)
+sa
+attach(sa)
+aov.sa <- aov(sam ~ group, data=sa)
+summary(aov.sa)
+df.tot <- (3*n)-1
+ss.tot <- var(sam)*df.tot
+var.tot <- var(sam)
+df.tot
+ss.tot
+var.tot
+df.bet <- 3 - 1
+df.s1 <- n-1
+df.s2 <- n-1
+df.s3 <- n-1
+df.within <- df.s1+df.s2+df.s3
+df.within
+ss.s1 <- var(s1)*df.s1
+ss.s2 <- var(s2)*df.s2
+ss.s3 <- var(s3)*df.s3
+ss.s1
+ss.s2
+ss.s3
+ss.within <- ss.s1+ss.s2+ss.s3
+ss.within
+mean.grand <- mean(sam)
+mean.s1 <- mean(s1)
+mean.s2 <- mean(s2)
+mean.s3 <- mean(s3)
+mean.grand
+mean.s1
+mean.s2
+mean.s3
+diff.g1 <- n * (mean.grand - mean.s1)^2
+diff.g2 <- n * (mean.grand - mean.s2)^2
+diff.g3 <- n * (mean.grand - mean.s3)^2
+diff.g1
+diff.g2
+diff.g3
+ss.bet <-  diff.g1 + diff.g2 + diff.g3
+ss.bet
+# sumup
+ss.tot
+ss.bet
+ss.within
+ss.bet+ss.within
+df.tot
+df.bet
+df.within
+df.bet + df.within
+ms.bet <- ss.bet / df.bet
+ms.wit <- ss.within / df.within
+ms.bet
+ms.wit
+fvalue <- ms.bet/ms.wit
+fvalue
+f.res <- (aov(sa$sam~sa$group))
+summary(f.res)
+##########################
+# http://commres.net/wiki/r/oneway_anova
+#
+##########################
+x1 <- c(50.5, 52.1, 51.9, 52.4, 50.6, 51.4, 51.2, 52.2, 51.5, 50.8)
+x2 <- c(47.5, 47.7, 46.6, 47.1, 47.2, 47.8, 45.2, 47.4, 45.0, 47.9)
+x3 <- c(46.0, 47.1, 45.6, 47.1, 47.2, 46.4, 45.9, 47.1, 44.9, 46.2)
+xc <- data.frame(x1,x2,x3)
+xs <- stack(xc)
+xs
+colnames(xs) <- c("score", "temp")
+xs
+levels(xs$temp) <- c("low", "mid", "hi")
+xs
+str(xs) # 데이터 구조 체크
+attach(xs)
+aov.xs <- aov(score ~ temp, data = xs)
+summary(aov.xs)
+mean.tot <- mean(score)
+mean.tot
+df.tot <- length(score)-1
+ss.tot <- var(score) * df.tot
+ss.tot
+df.tot
+var.tot <- var(score)
+var.tot2 <- ss.tot/df.tot
+var.tot == var.tot2
+mean.each <- tapply(score, list(temp), mean)
+mean.each
+var.each <- tapply(score, list(temp), var)
+var.each
+df.within <- 3*(10-1)
+ss.each <- var.each * 9
+ss.within <- sum(ss.each)
+ss.within
+ss.tot - ss.within
+ss.bet <- sum(10*(mean.tot - mean.each)^2)
+ss.bet
+df.bet <- 3-1
+ss.tot
+ss.bet
+ss.within
+ss.bet + ss.within
+df.tot
+df.bet
+df.within
+df.bet + df.within
+ms.bet <- ss.bet / df.bet
+ms.within <- ss.within / df.within
+f.value <- ms.bet / ms.within
+ss.bet
+ss.within
+df.bet
+df.within
+ms.bet
+ms.within
+f.value
+alpha <- .05
+qf(alpha, df.bet, df.within, lower.tail = F)
+pf(f.value, df.bet, df.within, lower.tail = F)
+aov.xs <- aov(score ~ temp, data = xs)
+summary(aov.xs)
+#################################################
+# two-way anova
+# subject = factor(paste('sub', 1:30, sep=''))
+#################################################
+n.a.group <- 3 # a treatment 숫자
+n.b.group <- 2 # b 그룹 숫자
+n.sub <- 30 # 총 샘플 숫자
+# 데이터 생성
+set.seed(9)
+a <- gl(3, 10, n.sub, labels=c('a1', 'a2', 'a3'))
+b <- gl(2, 5, n.sub, labels=c('b1', 'b2'))
+a
+b
+y <- rnorm(30, mean=10) + 3.14 * (a=='a1' & b=='b2') + 1.43 * (a=='a3' & b=='b2')
+y
+dat <- data.frame(a, b, y)
+aov.dat <- aov(y~a*b) # anova test
+summary(aov.dat) # summary of the test output
+# hand calculation
+table(a,b)
+tapply(y, list(a,b), mean) # 각 셀에서의 평균
+df.within.each <- tapply(y, list(a,b), length) -1  # 각 셀에서의 샘플숫자
+n.within.each <- df.within.each + 1
+df.within <- sum(df.within.each) # df within
+var.within <- tapply(y, list(a,b), var) # var.within
+ss.within.each <- tapply(y, list(a,b), var) * df.within.each
+ss.within.each
+ss.within <- sum(ss.within.each) # ss.within
+ss.within
+interaction.plot(a,b,y)
+mean.a <- tapply(y, list(a), mean)
+mean.b <- tapply(y, list(b), mean)
+mean.a
+mean.b
+var.a <- tapply(y, list(a), var)
+var.b <- tapply(y, list(b), var)
+mean.tot <- mean(dat$y)
+var.tot <- var(dat$y)
+df.tot <- n.sub - 1
+ss.tot <- var.tot * df.tot
+## between
+mean.each <- tapply(y, list(a,b), mean)
+mean.each
+mean.tot <- mean(y)
+mean.tot
+n.each <- tapply(y, list(a,b), length)
+n.each
+n.a.each <- tapply(y, list(a), length)
+n.b.each <- tapply(y, list(b), length)
+ss.bet <- sum(n.each*(mean.each-mean.tot)^2)
+ss.bet
+ss.tot
+ss.within
+ss.bet
+ss.bet + ss.within
+ss.a <- sum(n.a.each * ((mean.tot - mean.a)^2))
+ss.b <- sum(n.b.each * ((mean.tot - mean.b)^2))
+ss.a
+ss.b
+ss.ab <- ss.bet - (ss.a + ss.b)
+ss.ab
+ss.tot
+ss.bet
+ss.within
+ss.a
+ss.b
+ss.ab
+df.tot <- n.sub - 1
+df.bet <- (n.a.group * n.b.group) - 1
+df.a <- n.a.group - 1
+df.b <- n.b.group - 1
+df.ab <- df.bet - (df.a + df.b)
+df.within <- sum(df.within.each)
+df.tot
+df.bet
+df.a
+df.b
+df.ab
+df.within
+ms.a <- ss.a / df.a
+ms.b <- ss.b / df.b
+ms.ab <- ss.ab / df.ab
+ms.within <- ss.within / df.within
+ms.a
+ms.b
+ms.ab
+ms.within
+f.a <- ms.a / ms.within
+f.b <- ms.b / ms.within
+f.ab <- ms.ab / ms.within
+alpha <- .05
+f.a
+# 봐야할 F분포표에서의 F-value
+# qt 처럼 qf 사용
+# qf(alpha, df.between, df.within, lower.tail=F) 처럼 사용
+qf(alpha, df.a, df.within, lower.tail = F)
+pf(f.a, df.a, df.within, lower.tail = F)
+f.b
+qf(alpha, df.b, df.within, lower.tail = F)
+pf(f.b, df.b, df.within, lower.tail = F)
+f.ab
+qf(alpha, df.ab, df.within, lower.tail = F)
+pf(f.ab, df.ab, df.within, lower.tail = F)
+# aov result
+summary(aov.dat)
 </code>