====== ANOVA in R ======
#
# 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 구하기
16*((m.a-mean.tot)^2)
16*((m.b-mean.tot)^2)
16*((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
# 위에서
# 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
# 그러나 정확히 어떤 그룹에서 차이가 나는지는 판단해주지 않음
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)
===== ANOVA in R: Output =====
> #
> # 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
1 28.956987 a
2 24.588298 a
3 29.232040 a
4 15.221730 a
5 18.069720 a
6 26.455426 a
7 28.824505 a
8 27.362726 a
9 11.013442 a
10 28.284603 a
11 31.778215 a
12 28.525725 a
13 29.735641 a
14 33.996331 a
15 22.487061 a
16 31.467550 a
17 26.638215 b
18 25.537955 b
19 30.194374 b
20 22.161849 b
21 23.164369 b
22 17.436099 b
23 34.602104 b
24 12.091427 b
25 36.147829 b
26 21.460025 b
27 30.381254 b
28 24.367170 b
29 17.691419 b
30 26.351577 b
31 11.330276 b
32 24.444055 b
33 10.842903 c
34 20.414861 c
35 11.872108 c
36 29.195858 c
37 26.074787 c
38 18.262033 c
39 18.863621 c
40 25.906536 c
41 4.379624 c
42 2.980549 c
43 20.825446 c
44 22.038722 c
45 24.175933 c
46 25.745030 c
47 23.796968 c
48 18.625020 c
> colnames(comb3)[1] <- "values"
> colnames(comb3)[2] <- "group"
> comb3
values group
1 28.956987 a
2 24.588298 a
3 29.232040 a
4 15.221730 a
5 18.069720 a
6 26.455426 a
7 28.824505 a
8 27.362726 a
9 11.013442 a
10 28.284603 a
11 31.778215 a
12 28.525725 a
13 29.735641 a
14 33.996331 a
15 22.487061 a
16 31.467550 a
17 26.638215 b
18 25.537955 b
19 30.194374 b
20 22.161849 b
21 23.164369 b
22 17.436099 b
23 34.602104 b
24 12.091427 b
25 36.147829 b
26 21.460025 b
27 30.381254 b
28 24.367170 b
29 17.691419 b
30 26.351577 b
31 11.330276 b
32 24.444055 b
33 10.842903 c
34 20.414861 c
35 11.872108 c
36 29.195858 c
37 26.074787 c
38 18.262033 c
39 18.863621 c
40 25.906536 c
41 4.379624 c
42 2.980549 c
43 20.825446 c
44 22.038722 c
45 24.175933 c
46 25.745030 c
47 23.796968 c
48 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
600 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
>
>
> # 위에서
> # 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
>
> # 그러나 정확히 어떤 그룹에서 차이가 나는지는 판단해주지 않음
> 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
95% 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
>
====== Post hoc test ======
[[:post hoc test]]
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)
plot(TukeyHSD(a.res, conf.level=.95), las = 2)
pairwise.t.test(comb3$values, comb3$group, p.adj = "bonf")
===== post hoc test: output =====
> 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
95% 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
{{: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의 불확실 성
* 혹은 획득되는 설명력
SS error
* IV가 고려되었음에도 불구하고
* 끝까지 남는 error
SS total = SS between + SS within
SS between / SS total = IV 가 kicked in 되었을 때 없어지는 uncertainty = IV 의 설명력 = R square value
즉, IV로 uncertainty 가 얼마나 없어질까? 라는 아이디어
이를 살펴보기 위해
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)
===== R square: output =====
> 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
>
>
{{tag> statistics r "analysis of variance" anova}}