Differences

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

--- r:anova [2019/09/18 07:59] – hkimscil
+++ r:anova [2024/04/17 08:30] (current) – [ANOVA in R: Output] hkimscil
@@ Line 1: / Line 1: @@
-[[:r:Oneway ANOVA]]
+====== ANOVA in R ======
-[[:r:Twoway ANOVA]]
+<code>
+#
+# 3 샘플 종류 추출
+# A, B, C 학년에 따라서 욕하는 정도가 달라질것이라는
+# 가설
+set.seed(201)
+rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) }
+A <- rnorm2(16, 26, sqrt(600/15))
+B <- rnorm2(16, 24, sqrt(750/15))
+C <- rnorm2(16, 19, sqrt(900/15))
+A <- c(A)
+B <- c(B)
+C <- c(C)
+# 평균구하기
+mean(A)
+mean(B)
+mean(C)
+# 3 샘플을 합치기
+# 두번재 컬럼 = group A, B, C 가 되도록
+comb3 <- stack(list(a=A, b=B, c=C))
+comb3
+colnames(comb3)[1] <- "values"
+colnames(comb3)[2] <- "group"
+comb3
+# 전체구성원을 하나로 보고 분산값을 구해본다
+# ms.tot = ss.tot / df.tot
+# ss.tot = sum(xi-mean.tot)^2
+# df.tot = N - 1
+# attach(comb3)
+mean.tot <- mean(comb3$values)
+ss.tot <- sum((comb3$values-mean.tot)^2)
+df.tot <- 48-1
+ms.tot <- ss.tot/df.tot
+ss.tot
+df.tot
+ms.tot
+var(comb3$values)
+# within part 구하기
+# 간단한 방법
+tapply(comb3$values, comb3$group, var)*15
+sse <- sum(tapply(comb3$values, comb3$group, var)*15)
+sse
+# 이해한 개념대로 얼른 구하기
+# A, B, C 평균
+m.a <- mean(A)
+m.b <- mean(B)
+m.c <- mean(C)
+ms.a <- var(A)
+ms.b <- var(B)
+ms.c <- var(C)
+df.a <- 15
+df.b <- 15
+df.c <- 15
+# 사실 우리는 아래 각 그룹의 SS가
+# 600, 750, 900 인것을 알고 있다.
+ss.a <- ms.a*df.a
+ss.b <- ms.b*df.b
+ss.c <- ms.c*df.c
+ss.within <- ss.a + ss.b + ss.c
+ss.within == sse # check both are the same
+# between part 구하기
+*((m.a-mean.tot)^2)
+*((m.b-mean.tot)^2)
+*((m.c-mean.tot)^2)
+ss.bet <- 16*((m.a-mean.tot)^2) + 16*((m.b-mean.tot)^2) + 16*((m.c-mean.tot)^2)
+ss.tot
+ss.bet
+ss.within
+ss.bet + ss.within
+df.tot
+df.within <- df.a + df.b + df.c
+df.within
+df.bet <- 3-1
+df.bet
+df.bet + df.within
+ms.bet <- ss.bet/df.bet
+ms.within <- ss.within /df.within
+f.calculated <- ms.bet/ms.within
+ms.bet
+ms.within
+f.calculated
+# 이 점수에서의 F critical value =
+fp.value <- pf(f.calculated, df1=2, df2=45, lower.tail = F)
+fp.value
+# f.calculated 와 fp.value로 판단
+f.calculated
+fp.value
+a.res <- aov(values ~ group, data=comb3)
+a.res.sum <- summary(a.res)
+a.res.sum
+</code>
+<code>
+# 위에서
+# ssd라는 function을 만들면
+ssd <- function(x) {
+          var(x) * (length(x)-1) }
+ss.a1 <- ssd(A)
+ss.b2 <- ssd(B)
+ss.c1 <- ssd(C)
+ss.a == ss.a1
+ss.b == ss.b1
+ss.c == ss.c1
+</code>
+<code>
+# 그러나 정확히 어떤 그룹에서 차이가 나는지는 판단해주지 않음
+pairwise.t.test(comb3$values, comb3$group, p.adj = "none")
+# OR
+pairwise.t.test(comb3$values, comb3$group, p.adj = "bonf")
+pairwise.t.test(comb3$values, comb3$group, p.adj = "holm")
+# OR TukeyHSD(anova.output)
+TukeyHSD(a.res)
+</code>
+===== ANOVA in R: Output =====
+<code>
+> #
+> # 3 샘플 종류 추출
+> # A, B, C 학년에 따라서 욕하는 정도가 달라질것이라는
+> # 가설
+> set.seed(201)
+> rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) }
+> A <- rnorm2(16, 26, sqrt(600/15))
+> B <- rnorm2(16, 24, sqrt(750/15))
+> C <- rnorm2(16, 19, sqrt(900/15))
+>
+> A <- c(A)
+> B <- c(B)
+> C <- c(C)
+>
+> # 평균구하기
+> mean(A)
+[1] 26
+> mean(B)
+[1] 24
+> mean(C)
+[1] 19
+>
+> # 3 샘플을 합치기
+> # 두번재 컬럼 = group A, B, C 가 되도록
+> comb3 <- stack(list(a=A, b=B, c=C))
+> comb3
+      values ind
+  28.956987   a
+  24.588298   a
+  29.232040   a
+  15.221730   a
+  18.069720   a
+  26.455426   a
+  28.824505   a
+  27.362726   a
+  11.013442   a
+28.284603   a
+31.778215   a
+28.525725   a
+29.735641   a
+33.996331   a
+22.487061   a
+31.467550   a
+26.638215   b
+25.537955   b
+30.194374   b
+22.161849   b
+23.164369   b
+17.436099   b
+34.602104   b
+12.091427   b
+36.147829   b
+21.460025   b
+30.381254   b
+24.367170   b
+17.691419   b
+26.351577   b
+11.330276   b
+24.444055   b
+10.842903   c
+20.414861   c
+11.872108   c
+29.195858   c
+26.074787   c
+18.262033   c
+18.863621   c
+25.906536   c
+  4.379624   c
+  2.980549   c
+20.825446   c
+22.038722   c
+24.175933   c
+25.745030   c
+23.796968   c
+18.625020   c
+> colnames(comb3)[1] <- "values"
+> colnames(comb3)[2] <- "group"
+> comb3
+      values group
+  28.956987     a
+  24.588298     a
+  29.232040     a
+  15.221730     a
+  18.069720     a
+  26.455426     a
+  28.824505     a
+  27.362726     a
+  11.013442     a
+28.284603     a
+31.778215     a
+28.525725     a
+29.735641     a
+33.996331     a
+22.487061     a
+31.467550     a
+26.638215     b
+25.537955     b
+30.194374     b
+22.161849     b
+23.164369     b
+17.436099     b
+34.602104     b
+12.091427     b
+36.147829     b
+21.460025     b
+30.381254     b
+24.367170     b
+17.691419     b
+26.351577     b
+11.330276     b
+24.444055     b
+10.842903     c
+20.414861     c
+11.872108     c
+29.195858     c
+26.074787     c
+18.262033     c
+18.863621     c
+25.906536     c
+  4.379624     c
+  2.980549     c
+20.825446     c
+22.038722     c
+24.175933     c
+25.745030     c
+23.796968     c
+18.625020     c
+>
+> # 전체구성원을 하나로 보고 분산값을 구해본다
+> # ms.tot = ss.tot / df.tot
+> # ss.tot = sum(xi-mean.tot)^2
+> # df.tot = N - 1
+> # attach(comb3)
+> mean.tot <- mean(comb3$values)
+> ss.tot <- sum((comb3$values-mean.tot)^2)
+> df.tot <- 48-1
+> ms.tot <- ss.tot/df.tot
+> ss.tot
+[1] 2666
+> df.tot
+[1] 47
+> ms.tot
+[1] 56.7234
+> var(comb3$values)
+[1] 56.7234
+>
+> # within part 구하기
+> # 간단한 방법
+> tapply(comb3$values, comb3$group, var)*15
+  a   b   c
+750 900
+> sse <- sum(tapply(comb3$values, comb3$group, var)*15)
+> sse
+[1] 2250
+>
+> # 이해한 개념대로 얼른 구하기
+> # A, B, C 평균
+> m.a <- mean(A)
+> m.b <- mean(B)
+> m.c <- mean(C)
+> ms.a <- var(A)
+> ms.b <- var(B)
+> ms.c <- var(C)
+> df.a <- 15
+> df.b <- 15
+> df.c <- 15
+> # 사실 우리는 아래 각 그룹의 SS가
+> # 600, 750, 900 인것을 알고 있다.
+> ss.a <- ms.a*df.a
+> ss.b <- ms.b*df.b
+> ss.c <- ms.c*df.c
+>
+> ss.within <- ss.a + ss.b + ss.c
+> ss.within == sse # check both are the same
+[1] TRUE
+>
+>
+> # between part 구하기
+> 16*((m.a-mean.tot)^2)
+[1] 144
+> 16*((m.b-mean.tot)^2)
+[1] 16
+> 16*((m.c-mean.tot)^2)
+[1] 256
+> ss.bet <- 16*((m.a-mean.tot)^2) + 16*((m.b-mean.tot)^2) + 16*((m.c-mean.tot)^2)
+>
+> ss.tot
+[1] 2666
+> ss.bet
+[1] 416
+> ss.within
+[1] 2250
+> ss.bet + ss.within
+[1] 2666
+>
+> df.tot
+[1] 47
+> df.within <- df.a + df.b + df.c
+> df.within
+[1] 45
+> df.bet <- 3-1
+> df.bet
+[1] 2
+> df.bet + df.within
+[1] 47
+>
+> ms.bet <- ss.bet/df.bet
+> ms.within <- ss.within /df.within
+> f.calculated <- ms.bet/ms.within
+> ms.bet
+[1] 208
+> ms.within
+[1] 50
+> f.calculated
+[1] 4.16
+> # 이 점수에서의 F critical value =
+> fp.value <- pf(f.calculated, df1=2, df2=45, lower.tail = F)
+> fp.value
+[1] 0.02199143
+>
+> # f.calculated 와 fp.value로 판단
+> f.calculated
+[1] 4.16
+> fp.value
+[1] 0.02199143
+>
+>
+> a.res <- aov(values ~ group, data=comb3)
+> a.res.sum <- summary(a.res)
+> a.res.sum
+            Df Sum Sq Mean Sq F value Pr(>F)
+group        2    416     208    4.16  0.022 *
+Residuals   45   2250      50
+---
+Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+>
+>
+</code>
+<code>
+> # 위에서
+> # ssd라는 function을 만들면
+>
+> ssd <- function(x) {
++     var(x) * (length(x)-1) }
+> ss.a1 <- ssd(A)
+> ss.b1 <- ssd(B)
+> ss.c1 <- ssd(C)
+>
+> ss.a == ss.a1
+[1] TRUE
+> ss.b == ss.b1
+[1] TRUE
+> ss.c == ss.c1
+[1] TRUE
+>
+</code>
+<code>
+> # 그러나 정확히 어떤 그룹에서 차이가 나는지는 판단해주지 않음
+> pairwise.t.test(comb3$values, comb3$group, p.adj = "none")
+	Pairwise comparisons using t tests with pooled SD
+data:  comb3$values and comb3$group
+  a      b
+b 0.4279 -
+c 0.0075 0.0516
+P value adjustment method: none
+> # OR
+> pairwise.t.test(comb3$values, comb3$group, p.adj = "bonf")
+	Pairwise comparisons using t tests with pooled SD
+data:  comb3$values and comb3$group
+  a     b
+b 1.000 -
+c 0.023 0.155
+P value adjustment method: bonferroni
+> pairwise.t.test(comb3$values, comb3$group, p.adj = "holm")
+	Pairwise comparisons using t tests with pooled SD
+data:  comb3$values and comb3$group
+  a     b
+b 0.428 -
+c 0.023 0.103
+P value adjustment method: holm
+>
+> # OR TukeyHSD(anova.output)
+> TukeyHSD(a.res)
+  Tukey multiple comparisons of means
+% family-wise confidence level
+Fit: aov(formula = values ~ group, data = comb3)
+$group
+    diff        lwr       upr     p adj
+b-a   -2  -8.059034  4.059034 0.7049466
+c-a   -7 -13.059034 -0.940966 0.0201250
+c-b   -5 -11.059034  1.059034 0.1238770
+>
+</code>
+====== Post hoc test ======
+[[:post hoc test]]
+<code>
+m.a
+m.b
+m.c
+d.ab <- m.a - m.b
+d.ac <- m.a - m.c
+d.bc <- m.b - m.c
+d.ab
+d.ac
+d.bc
+# mse (ms within) from the a.res.sum output
+# a.res.sum == summary(aov(values ~ group, data=comb3))
+a.res.sum
+# mse = 50
+mse <- 50
+# 혹은 fansy way from comb3 data.frame
+# 15 는 각 그룹의 df
+# 각 그룹에 따라서 values에 대한 variance를 구하여
+# df 값이 15를 곱해서 SS within_group 값을 모두 구하고
+# 이를 합산한다
+sse.ch <- sum(tapply(comb3$values, comb3$group, var)*15)
+sse.ch
+#### 사실 위의 값은 먼저 구한 ss.within 값
+ss.within
+mse.ch <- sse.ch/45
+mse.ch
+# http://commres.net/wiki/post_hoc_tes 을 보면
+# ms.error 를 그 그룹의 샘플 숫자로 나눈다 (length(A) = 16)
+se <- sqrt(mse/length(A))
+# now t scores for two compared groups
+t.ab <- d.ab / se
+t.ac <- d.ac / se
+t.bc <- d.bc / se
+t.ab
+t.ac
+t.bc
+# 이제 위의 점수를 .05 레벨에서 비교할 점수를 찾아 비교한다
+# qtukey 펑션을 이용한다
+t.crit <- qtukey( .95, nmeans = 3, df = 45)
+t.crit
+# t.ac만이 큰 것을 알 수 있다.
+# 따라서 a 와 c 가 서로 다른 그룹
+# 즉, 1학년과 3학년이 서로 다른 그룹
+# 혹은 R이 보통 제시한 거꾸로 방법으로 보면
+ptukey(t.ab, nmeans=3, df=45, lower.tail = F)
+ptukey(t.ac, nmeans=3, df=45, lower.tail = F)
+ptukey(t.bc, nmeans=3, df=45, lower.tail = F)
+TukeyHSD(a.res, conf.level=.95)
+</code>
+<code>
+plot(TukeyHSD(a.res, conf.level=.95), las = 2)
+</code>
+<code>
+pairwise.t.test(comb3$values, comb3$group, p.adj = "bonf")
+</code>
+===== post hoc test: output =====
+<code>
+> m.a
+[1] 26
+> m.b
+[1] 24
+> m.c
+[1] 19
+>
+> d.ab <- m.a - m.b
+> d.ac <- m.a - m.c
+> d.bc <- m.b - m.c
+>
+> d.ab
+[1] 2
+> d.ac
+[1] 7
+> d.bc
+[1] 5
+>
+> # se
+> # from the a.res.sum output
+> a.res.sum
+            Df Sum Sq Mean Sq F value Pr(>F)
+group        2    416     208    4.16  0.022 *
+Residuals   45   2250      50
+---
+Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+> # mse = 50
+> mse <- 50
+>
+> se <- sqrt(mse/length(A))
+>
+> # now t scores for two compared groups
+> t.ab <- d.ab / se
+> t.ac <- d.ac / se
+> t.bc <- d.bc / se
+>
+> t.ab
+[1] 1.131371
+> t.ac
+[1] 3.959798
+> t.bc
+[1] 2.828427
+>
+> # 이제 위의 점수를 .05 레벨에서 비교할 점수를 찾아 비교한다
+> # qtukey 펑션을 이용한다
+> t.crit <- qtukey( .95, nmeans = 3, df = 45)
+> t.crit
+[1] 3.427507
+>
+> # t.ac만이 큰 것을 알 수 있다.
+> # 따라서 a 와 c 가 서로 다른 그룹
+> # 즉, 1학년과 3학년이 서로 다른 그룹
+>
+> # 혹은 R이 보통 제시한 거꾸로 방법으로 보면
+> ptukey(t.ab, nmeans=3, df=45, lower.tail = F)
+[1] 0.7049466
+> ptukey(t.ac, nmeans=3, df=45, lower.tail = F)
+[1] 0.02012498
+> ptukey(t.bc, nmeans=3, df=45, lower.tail = F)
+[1] 0.123877
+>
+> TukeyHSD(a.res, conf.level=.95)
+  Tukey multiple comparisons of means
+% family-wise confidence level
+Fit: aov(formula = values ~ group, data = comb3)
+$group
+    diff        lwr       upr     p adj
+b-a   -2  -8.059034  4.059034 0.7049466
+c-a   -7 -13.059034 -0.940966 0.0201250
+c-b   -5 -11.059034  1.059034 0.1238770
+</code>
+{{:c:ms:2023:pasted:20230418-223608.png}}
+====== R square or Eta square ======
+SS toal
+  * = Y 변인만으로 Y를 예측했을 때의 오차의 제곱
+  * Y 변인만으로 = Y의 평균을 가지고 Y값을 예측한 것
+  * 학습 초기에 에러의 제곱의 합으로 설명된 것
+SS between
+  * X 변인 (independent variable) 정보가 고려 되었을 때
+  * 독립변인이 고려되었을 때 (됨으로써)
+  * 없어지는 SS total의 불확실 성
+  * 혹은 획득되는 <fc #ff0000>설명력</fc>
+SS error
+  * IV가 고려되었음에도 불구하고
+  * 끝까지 남는 error
+SS total = SS between + SS within
+SS between / SS total = IV 가 kicked in 되었을 때 없어지는 uncertainty = IV 의 설명력 = <fc #ff0000>R square value</fc>
+즉, IV로 uncertainty 가 얼마나 없어질까? 라는 아이디어
+이를 살펴보기 위해
+<code>
+ss.tot
+ss.bet
+r.sq <- ss.bet / ss.tot
+r.sq
+# then . . . .
+lm.res <- lm(values ~ group, data = comb3)
+summary(lm.res)
+anova(lm.res)
+</code>
+===== R square: output =====
+<code>
+> ss.tot
+[1] 2666
+> ss.bet
+[1] 416
+> r.sq <- ss.bet / ss.tot
+> r.sq
+[1] 0.156039
+>
+> # then . . . .
+>
+> lm.res <- lm(values ~ group, data = comb3)
+> summary(lm.res)
+Call:
+lm(formula = values ~ group, data = comb3)
+Residuals:
+    Min      1Q  Median      3Q     Max
+-16.020  -2.783   1.476   4.892  12.148
+Coefficients:
+            Estimate Std. Error t value Pr(>|t|)
+(Intercept)   26.000      1.768   14.71   <2e-16 ***
+groupb        -2.000      2.500   -0.80   0.4279
+groupc        -7.000      2.500   -2.80   0.0075 **
+---
+Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+Residual standard error: 7.071 on 45 degrees of freedom
+Multiple R-squared:  0.156,	Adjusted R-squared:  0.1185
+F-statistic:  4.16 on 2 and 45 DF,  p-value: 0.02199
+> anova(lm.res)
+Analysis of Variance Table
+Response: values
+          Df Sum Sq Mean Sq F value  Pr(>F)
+group      2    416     208    4.16 0.02199 *
+Residuals 45   2250      50
+---
+Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+>
+>
+</code>
 {{tag> statistics r "analysis of variance" anova}}