====== ANOVA in R ======
#
# 3 샘플 종류 추출
# A, B, C 학년에 따라서 욕하는 정도가 달라질것이라는
# 가설
set.seed(201)
ss <- 16
rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) }
A <- rnorm2(ss, 26, sqrt(600/15))
B <- rnorm2(ss, 24, sqrt(750/15))
C <- rnorm2(ss, 19, sqrt(900/15))
A <- c(A)
B <- c(B)
C <- c(C)
A
B
C
# 평균구하기
mean(A)
mean(B)
mean(C)
# just checking
var(A)
var(B)
var(C)
# 각 그룹의 SS값은?
# 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 <- length(comb3$values)-1
ms.tot <- ss.tot / df.tot
ss.tot
df.tot
ms.tot
var(comb3$values)
# 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
# 3그룹이 전체평균에서 얼마나 떨어져서 분포되어 있나 보기 위해서
# 전체평균에서 각 그룹의 평균이 얼마나 떨어져 있는지를 구한 후
# (deviation score 혹은 error score), 이를 제곱한다
# 위의 세 점수를 더하기 전에 각 그룹에 구성원은 16명씩이므로
# (ss) 이를 곱한 후에 모두 더한다. 이것을 SS between group
# 이라 부른다.
ss*((m.a-mean.tot)^2)
ss*((m.b-mean.tot)^2)
ss*((m.c-mean.tot)^2)
# 그리고 우리는 이것이 독립변인이 존재함으로써 밝혀진
# SS값임을 알고 있다
ss.bet <- 16*((m.a-mean.tot)^2) + 16*((m.b-mean.tot)^2) + 16*((m.c-mean.tot)^2)
ss.bet
# 한편 ss.within은 각 그룹 내부에 구성원 간의
# 변화도이므로 (variability) 독립변인이 있음에도
# 불구하고 랜덤하게 나타나는 SS이다
ss.within
# 그리고 ss.total은
# 독립변인 때문에 설명되는 혹은 독립변인이
# 있기 때문에 밝혀진 ss.between 값과
# 랜덤하게 흩어진 그룹구성원 간의 ss.within
# 값으로 구성된다.
ss.tot
ss.bet
ss.within
ss.bet + ss.within
# 이는 df값도 마찬가지이다.
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
# 컴퓨터 계산이 쉬워지기 전에는
# 아래처럼 0.5 level에서의 f값을 구한 후
# 이것과 계산된 f값을 비교해봤었다.
qf(.05, df1 = 2, df2 = 45, lower.tail = FALSE)
f.calculated
# 위에서 f.calculated > qf값이므로
# f.calculated 값으로 영가설을 부정하고
# 연구가설을 채택하면 판단이 잘못일 확률이
# 0.05보다 작다는 것을 안다.
# 그러나 컴퓨터계산이 용이해지고서는 qf대신에
# pf를 써서 f.cal 값에 해당하는 prob. level을
# 알아본다.
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
>
> # 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 <- ss.within
>
> 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.b2 <- 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
sse.ch <- sum(tapply(comb3$values, comb3$group, var)*15)
sse.ch
mse.ch <- sse.ch/45
mse.ch
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
>
>