c:ms:2025:w07_anova_note
This is an old revision of the document!
Table of Contents
ANOVA in R
#
# ANOVA test with 4 levels in IV
#
rm(list=ls())
set.seed(101)
rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) }
n <- 30
na <- n
nb <- n
nc <- n
nd <- n
mean.a <- 26
mean.b <- 25
mean.c <- 23
mean.d <- 22
A <- rnorm2(na, mean.a, sqrt(1160/(na-1)))
B <- rnorm2(nb, mean.b, sqrt(1000/(nb-1)))
C <- rnorm2(nc, mean.c, sqrt(1000/(nc-1)))
D <- rnorm2(nd, mean.d, sqrt(1000/(nd-1)))
# A combined group with group A and B
# We call it group total
# we can obtain its mean, variance, ss, df, etc.
#
A
B
C
D
comb <- data.frame(A, B, C, D)
dat <- stack(comb)
head(dat)
colnames(dat)[1] <- "values"
colnames(dat)[2] <- "group"
head(dat)
mean.total <- mean(dat$values)
var.total <- var(dat$values)
# variance를 ms라고 부르기도 한다
ms.total <- var.total
df.total <- length(dat$values) - 1
ss.total <- var.total * df.total
ss.total.check <- sum(
(dat$values - mean(dat$values))^2
)
mean.total
var.total
ms.total
df.total
ss.total
ss.total.check
# Now for each group
mean.a <- mean(A)
mean.b <- mean(B)
mean.c <- mean(C)
mean.d <- mean(D)
mean.a
mean.b
mean.c
mean.d
# 그룹 간의 차이에서 나타나는 분산
# 수업시간에 설명을 잘 들을 것
hist(dat$values)
abline(v = mean(dat$values), lty=2, lwd=3, col="red")
abline(v = mean.a, lty=2, lwd=3, col="blue")
abline(v = mean.b, lty=2, lwd=3, col="darkgreen")
abline(v = mean.c, lty=2, lwd=3, col="black")
abline(v = mean.d, lty=2, lwd=3, col="orange")
# mean.total 에서 그룹a의 평균까지의 차이를 구한 후
# 이를 제곱하여 그룹 A 멤버의 숫자만큼 더한다 =
# 즉, SS를 구하는 방법.
# 전체평균에서 그룹평균을 뺀 것의 제곱을
# 그룹 구성원 숫자만큼 더하는 것
# 그리고 이들을 다시 모두 더하여
# ss.between에 저장
length(A) * ((mean.total - mean.a)^2)
length(B) * ((mean.total - mean.b)^2)
length(C) * ((mean.total - mean.c)^2)
length(D) * ((mean.total - mean.d)^2)
ss.between <-
length(A) * ((mean.total - mean.a)^2) +
length(B) * ((mean.total - mean.b)^2) +
length(C) * ((mean.total - mean.c)^2) +
length(D) * ((mean.total - mean.d)^2)
ss.between
# df between group은 연구에 사용된
# 그룹의 숫자에서 1을 뺀 숫자
n.groups <- nlevels(dat$group)
df.between <- n.groups - 1
# 이 그룹 간 차이에 기인하는 분산 값은
ms.between <- ss.between / df.between
# 한편 ss.a 와 ss.b는 각 그룹 내의
# 분산을 알아보기 위한 방법
df.a <- length(A) - 1
df.b <- length(B) - 1
df.c <- length(C) - 1
df.d <- length(D) - 1
ss.a <- var(A) * df.a
ss.b <- var(B) * df.b
ss.c <- var(C) * df.c
ss.d <- var(D) * df.d
ss.within <- ss.a + ss.b + ss.c + ss.d
df.within <- df.a + df.b + df.c + df.d
ms.within <- ss.within / df.within
# 여기까지 우리는
# 전체분산
# 그룹간분산
# 그룹내분산 값을
# 구한 것
# ms.between은 그룹의 차이때문에 생긴
# 분산으로 IV 혹은 treatment 때문에 생기는
# 차이에 기인하는 분산이고
# ms.within은 각 그룹 내부에서 일어나는 분산이므로
# (variation이므로) 연구자의 관심사와는 상관이 없이
# 나타나는 random한 분산이라고 하면
# t test 때와 마찬가지로
# 그룹의 차이 / 랜덤 차이를 (에러 -> 분산은 에러라고도 했다)
# 구해볼 수 있다.
# 즉, 그룹갑분산은 사실 = diff (between groups)
# 그리고 그룹내 분산은 사실 = re
# 따라서 우리는 위 둘 간의 비율을 t test와 같이
# 살펴볼 수 있다
# 이것을 f.calculated 이라고 하고
f.calculated <- ms.between / ms.within
# 이 값을 출력해 본다
f.calculated
# 컴퓨터 계산이 쉬워지기 전에는 아래처럼 0.5 level
# 에서의 f값을 구한 후 이것과 계산된 f값을 비교해봤었다.
qf(.05, df1 = df.between, df2 = df.within, lower.tail = FALSE)
f.calculated
# 위에서 f.calculated > qf값이므로
# f.calculated 값으로 영가설을 부정하고
# 연구가설을 채택하면 판단이 잘못일 확률이
# 0.05보다 작다는 것을 안다.
# 그러나 컴퓨터계산이 용이해지고서는 qf대신에
# pf를 써서 f.cal 값에 해당하는 prob. level을
# 알아본다.
# percentage of f distribution with
# df1 and df2 option
# 이는 그림의 왼쪽을 나타내므로
# 차이가 점점 커지게 되는 오른쪽을
# 계산하기 위해서는 1-x를 취한다
f.calculated.pvalue <- 1-pf(f.calculated, df1=df.between, df2=df.within)
f.calculated.pvalue
f.calculated
# graph 로 이해
x <- rf(500000, df1 = df.between, df2 = df.within)
y.max <- max(df(x,df1=df.between, df2=df.within))
hist(x,
breaks = "Scott",
freq = FALSE,
xlim = c(0, f.calculated + 1),
ylim = c(0, y.max + .3),
xlab = "",
main = paste("Histogram for
a F-distribution with
df1 = ", df.between,
"and df2 = ", df.within,
"F calculated value = ", round(f.calculated,3),
"p value = ", f.calculated.pvalue),
cex.main = 0.9
)
curve(df(x, df1 = df.between, df2 = df.within),
from = 0, to = f.calculated + 1, n = 5000,
col = "red", lwd = 2,
add = T)
abline(v=f.calculated, col="blue", lwd=2, lty="dotted")
f.calculated
f.calculated.pvalue
1 - f.calculated.pvalue
# Now check this
ss.total
ss.between
ss.within
ss.total
ss.between + ss.within
# 한편 df는
# df.total 30 - 1
df.total
df.between
df.within
df.total
df.between + df.within
##################################################
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.calculated f.calculated.pvalue # how much IV explains the DV # in terms of SS? r.square <- ss.between / ss.total eta <- r.square eta lm.res <- lm(dat$values~dat$group, data = dat) summary(lm.res) summary(a.res)
ANOVA in R: Output
> #
> # ANOVA test with 4 levels in IV
> #
> rm(list=ls())
> set.seed(101)
> rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) }
> n <- 30
> na <- n
> nb <- n
> nc <- n
> nd <- n
> mean.a <- 26
> mean.b <- 25
> mean.c <- 22
> mean.d <- 20
>
> A <- rnorm2(na, mean.a, sqrt(1160/(na-1)))
> B <- rnorm2(nb, mean.b, sqrt(1000/(nb-1)))
> C <- rnorm2(nc, mean.c, sqrt(1000/(nc-1)))
> D <- rnorm2(nd, mean.d, sqrt(1000/(nd-1)))
>
> # A combined group with group A and B
> # We call it group total
> # we can obtain its mean, variance, ss, df, etc.
> #
> A
[,1]
[1,] 24.15936
[2,] 30.80932
[3,] 21.51824
[4,] 28.24999
[5,] 28.97978
[6,] 35.51391
[7,] 31.31140
[8,] 25.77399
[9,] 33.56897
[10,] 24.93735
[11,] 30.61240
[12,] 20.61063
[13,] 37.43502
[14,] 15.52399
[15,] 24.83574
[16,] 25.16385
[17,] 20.19498
[18,] 27.06992
[19,] 20.43785
[20,] 11.10716
[21,] 25.38778
[22,] 31.99065
[23,] 24.59883
[24,] 15.54592
[25,] 32.26250
[26,] 15.95114
[27,] 30.16291
[28,] 25.72414
[29,] 30.16421
[30,] 30.39809
attr(,"scaled:center")
[1] -0.08287722
attr(,"scaled:scale")
[1] 0.8355107
> B
[,1]
[1,] 30.92777
[2,] 27.40269
[3,] 31.57423
[4,] 13.93715
[5,] 32.61602
[6,] 21.65799
[7,] 26.76631
[8,] 31.07316
[9,] 16.23555
[10,] 28.37195
[11,] 28.56653
[12,] 30.14509
[13,] 12.52765
[14,] 23.17424
[15,] 19.47689
[16,] 28.11125
[17,] 29.06156
[18,] 21.76270
[19,] 24.14405
[20,] 17.31020
[21,] 19.22003
[22,] 24.23347
[23,] 29.11289
[24,] 17.80809
[25,] 30.09268
[26,] 19.78715
[27,] 26.75141
[28,] 32.27229
[29,] 32.52368
[30,] 23.35537
attr(,"scaled:center")
[1] -0.1405676
attr(,"scaled:scale")
[1] 1.025799
> C
[,1]
[1,] 19.93680
[2,] 13.13256
[3,] 17.68193
[4,] 22.13674
[5,] 23.96961
[6,] 23.75823
[7,] 17.40748
[8,] 22.35212
[9,] 21.13146
[10,] 21.02997
[11,] 30.39517
[12,] 31.04547
[13,] 28.28695
[14,] 21.01354
[15,] 10.72282
[16,] 15.34118
[17,] 23.25979
[18,] 13.91989
[19,] 22.28969
[20,] 21.17085
[21,] 32.41776
[22,] 28.04180
[23,] 18.45008
[24,] 18.25799
[25,] 11.25474
[26,] 24.25413
[27,] 21.50399
[28,] 29.43868
[29,] 25.75134
[30,] 30.64723
attr(,"scaled:center")
[1] 0.08931913
attr(,"scaled:scale")
[1] 0.9936574
> D
[,1]
[1,] 27.502532
[2,] 19.249172
[3,] 17.265865
[4,] 15.085873
[5,] 21.168938
[6,] 14.658798
[7,] 16.756660
[8,] 27.742432
[9,] 22.476618
[10,] 14.513907
[11,] 21.084269
[12,] 15.862551
[13,] 32.407055
[14,] 26.575540
[15,] 23.989868
[16,] 18.058006
[17,] 19.989914
[18,] 6.202226
[19,] 16.624780
[20,] 29.690645
[21,] 16.009971
[22,] 19.173433
[23,] 18.504309
[24,] 29.182495
[25,] 24.122752
[26,] 14.773496
[27,] 15.629022
[28,] 14.417493
[29,] 15.869201
[30,] 25.412177
attr(,"scaled:center")
[1] 0.0894333
attr(,"scaled:scale")
[1] 0.9674409
> comb <- data.frame(A, B, C, D)
> dat <- stack(comb)
> head(dat)
values ind
1 24.15936 A
2 30.80932 A
3 21.51824 A
4 28.24999 A
5 28.97978 A
6 35.51391 A
> colnames(dat)[1] <- "values"
> colnames(dat)[2] <- "group"
> head(dat)
values group
1 24.15936 A
2 30.80932 A
3 21.51824 A
4 28.24999 A
5 28.97978 A
6 35.51391 A
>
> mean.total <- mean(dat$values)
> var.total <- var(dat$values)
> # variance를 ms라고 부르기도 한다
> ms.total <- var.total
>
> df.total <- length(dat$values) - 1
> ss.total <- var.total * df.total
> ss.total.check <- sum(
+ (dat$values - mean(dat$values))^2
+ )
>
> mean.total
[1] 23.25
> var.total
[1] 40.69328
> ms.total
[1] 40.69328
> df.total
[1] 119
> ss.total
[1] 4842.5
> ss.total.check
[1] 4842.5
>
> # Now for each group
> mean.a <- mean(A)
> mean.b <- mean(B)
> mean.c <- mean(C)
> mean.d <- mean(D)
> mean.a
[1] 26
> mean.b
[1] 25
> mean.c
[1] 22
> mean.d
[1] 20
>
> # 그룹 간의 차이에서 나타나는 분산
> # 수업시간에 설명을 잘 들을 것
> hist(dat$values)
> abline(v = mean(dat$values), lty=2, lwd=3, col="red")
> abline(v = mean.a, lty=2, lwd=3, col="blue")
> abline(v = mean.b, lty=2, lwd=3, col="darkgreen")
> abline(v = mean.c, lty=2, lwd=3, col="black")
> abline(v = mean.d, lty=2, lwd=3, col="orange")
>
> # mean.total 에서 그룹a의 평균까지의 차이를 구한 후
> # 이를 제곱하여 그룹 A 멤버의 숫자만큼 더한다 =
> # 즉, SS를 구하는 방법.
> # 전체평균에서 그룹평균을 뺀 것의 제곱을
> # 그룹 구성원 숫자만큼 더하는 것
> # 그리고 이들을 다시 모두 더하여
> # ss.between에 저장
>
> length(A) * ((mean.total - mean.a)^2)
[1] 226.875
> length(B) * ((mean.total - mean.b)^2)
[1] 91.875
> length(C) * ((mean.total - mean.c)^2)
[1] 46.875
> length(D) * ((mean.total - mean.d)^2)
[1] 316.875
>
>
> ss.between <-
+ length(A) * ((mean.total - mean.a)^2) +
+ length(B) * ((mean.total - mean.b)^2) +
+ length(C) * ((mean.total - mean.c)^2) +
+ length(D) * ((mean.total - mean.d)^2)
>
> ss.between
[1] 682.5
> # df between group은 연구에 사용된
> # 그룹의 숫자에서 1을 뺀 숫자
> n.groups <- nlevels(dat$group)
> df.between <- n.groups - 1
> # 이 그룹 간 차이에 기인하는 분산 값은
> ms.between <- ss.between / df.between
>
> # 한편 ss.a 와 ss.b는 각 그룹 내의
> # 분산을 알아보기 위한 방법
> df.a <- length(A) - 1
> df.b <- length(B) - 1
> df.c <- length(C) - 1
> df.d <- length(D) - 1
> ss.a <- var(A) * df.a
> ss.b <- var(B) * df.b
> ss.c <- var(C) * df.c
> ss.d <- var(D) * df.d
> ss.within <- ss.a + ss.b + ss.c + ss.d
> df.within <- df.a + df.b + df.c + df.d
> ms.within <- ss.within / df.within
>
> # 여기까지 우리는
> # 전체분산
> # 그룹간분산
> # 그룹내분산 값을
> # 구한 것
>
> # ms.between은 그룹의 차이때문에 생긴
> # 분산으로 IV 혹은 treatment 때문에 생기는
> # 차이에 기인하는 분산이고
>
> # ms.within은 각 그룹 내부에서 일어나는 분산이므로
> # (variation이므로) 연구자의 관심사와는 상관이 없이
> # 나타나는 random한 분산이라고 하면
>
> # t test 때와 마찬가지로
> # 그룹의 차이 / 랜덤 차이를 (에러 -> 분산은 에러라고도 했다)
> # 구해볼 수 있다.
>
> # 즉, 그룹갑분산은 사실 = diff (between groups)
> # 그리고 그룹내 분산은 사실 = re
> # 따라서 우리는 위 둘 간의 비율을 t test와 같이
> # 살펴볼 수 있다
>
>
> # 이것을 f.calculated 이라고 하고
> f.calculated <- ms.between / ms.within
> # 이 값을 출력해 본다
> f.calculated
[,1]
[1,] 6.34375
>
> # 컴퓨터 계산이 쉬워지기 전에는 아래처럼 0.5 level
> # 에서의 f값을 구한 후 이것과 계산된 f값을 비교해봤었다.
> qf(.05, df1 = df.between, df2 = df.within, lower.tail = FALSE)
[1] 2.682809
> f.calculated
[,1]
[1,] 6.34375
> # 위에서 f.calculated > qf값이므로
> # f.calculated 값으로 영가설을 부정하고
> # 연구가설을 채택하면 판단이 잘못일 확률이
> # 0.05보다 작다는 것을 안다.
> # 그러나 컴퓨터계산이 용이해지고서는 qf대신에
> # pf를 써서 f.cal 값에 해당하는 prob. level을
> # 알아본다.
>
> # percentage of f distribution with
> # df1 and df2 option
> # 이는 그림의 왼쪽을 나타내므로
> # 차이가 점점 커지게 되는 오른쪽을
> # 계산하기 위해서는 1-x를 취한다
> f.calculated.pvalue <- 1-pf(f.calculated, df1=df.between, df2=df.within)
> f.calculated.pvalue
[,1]
[1,] 0.0005078937
> f.calculated
[,1]
[1,] 6.34375
>
> # graph 로 이해
> x <- rf(500000, df1 = df.between, df2 = df.within)
> y.max <- max(df(x,df1=df.between, df2=df.within))
>
> hist(x,
+ breaks = "Scott",
+ freq = FALSE,
+ xlim = c(0, f.calculated + 1),
+ ylim = c(0, y.max + .3),
+ xlab = "",
+ main = paste("Histogram for
+ a F-distribution with
+ df1 = ", df.between,
+ "and df2 = ", df.within,
+ "F calculated value = ", round(f.calculated,3),
+ "p value = ", f.calculated.pvalue),
+ cex.main = 0.9
+ )
> curve(df(x, df1 = df.between, df2 = df.within),
+ from = 0, to = f.calculated + 1, n = 5000,
+ col = "red", lwd = 2,
+ add = T)
> abline(v=f.calculated, col="blue", lwd=2, lty="dotted")
>
> f.calculated
[,1]
[1,] 6.34375
> f.calculated.pvalue
[,1]
[1,] 0.0005078937
> 1 - f.calculated.pvalue
[,1]
[1,] 0.9994921
>
>
>
> # Now check this
> ss.total
[1] 4842.5
> ss.between
[1] 682.5
> ss.within
[,1]
[1,] 4160
> ss.total
[1] 4842.5
> ss.between + ss.within
[,1]
[1,] 4842.5
>
> # 한편 df는
> # df.total 30 - 1
> df.total
[1] 119
> df.between
[1] 3
> df.within
[1] 116
> df.total
[1] 119
> df.between + df.within
[1] 119
>
>
> ##################################################
> 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 682 227.50 6.344 0.000508 ***
Residuals 116 4160 35.86
---
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.51908 - -
C 0.01092 0.05478 -
D 0.00017 0.00159 0.19842
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.0655 0.3287 -
D 0.0010 0.0096 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.519 - -
C 0.044 0.164 -
D 0.001 0.008 0.397
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 -5.030484 3.0304845 0.9164944
C-A -4 -8.030484 0.0304845 0.0525532
D-A -6 -10.030484 -1.9695155 0.0009862
C-B -3 -7.030484 1.0304845 0.2171272
D-B -5 -9.030484 -0.9695155 0.0085400
D-C -2 -6.030484 2.0304845 0.5689321
>
> boxplot(dat$values~dat$group)
> f.calculated
[,1]
[1,] 6.34375
> f.calculated.pvalue
[,1]
[1,] 0.0005078937
>
> # how much IV explains the DV
> # in terms of SS?
> r.square <- ss.between / ss.total
> eta <- r.square
> eta
[1] 0.1409396
> 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
-14.8928 -4.1826 -0.3859 4.4517 12.4071
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.000 1.093 23.780 < 2e-16 ***
dat$groupB -1.000 1.546 -0.647 0.519080
dat$groupC -4.000 1.546 -2.587 0.010919 *
dat$groupD -6.000 1.546 -3.880 0.000174 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.988 on 116 degrees of freedom
Multiple R-squared: 0.1409, Adjusted R-squared: 0.1187
F-statistic: 6.344 on 3 and 116 DF, p-value: 0.0005079
> summary(a.res)
Df Sum Sq Mean Sq F value Pr(>F)
group 3 682 227.50 6.344 0.000508 ***
Residuals 116 4160 35.86
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>
Post hoc test
mean.a mean.b mean.c mean.d d.ab <- mean.a - mean.b d.ac <- mean.a - mean.c d.ad <- mean.a - mean.d d.bc <- mean.b - mean.c d.bd <- mean.b - mean.d d.cd <- mean.c - mean.d d.ab d.ac d.ad d.bc d.bd d.cd # mse (ms within) from the a.res.sum output # a.res.sum == summary(aov(values ~ group, data=dat)) a.res.sum # mse = 50 mse <- 35.86 # 혹은 fansy way from dat data.frame # df.a 는 각 그룹의 df (모든 그룹의 df값이 같으므로 df.a를 사용) tapply(dat$values, dat$group, var) # 각 그룹의 분산 tapply(dat$values, dat$group, var) * df.a # 각 그룹의 SS # 이 분산값에 df값을 곱하면 각 그룹의 SS값 # 이 값들을 sum 하면 총 ss.within 값 sse.ch <- sum(tapply(dat$values, dat$group, var) * df.a) sse.ch mse.ch <- sse.ch / (df.a*4) # 이 값에 df.a * 4 하면 ms.within값 mse.ch mse <- 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.ad <- d.ad / se t.bc <- d.bc / se t.bd <- d.bd / se t.cd <- d.cd / se t.ab t.ac t.ad t.bc t.bd t.cd # 이제 위의 점수를 .05 레벨에서 비교할 점수를 찾아 비교한다 # qtukey 펑션을 이용한다 t.crit <- qtukey( .95, nmeans = 3, df = 45) t.crit # 혹은 R이 보통 제시한 거꾸로 방법으로 보면 ptukey(t.ab, nmeans=4, df=df.within, lower.tail = F) ptukey(t.ac, nmeans=4, df=df.within, lower.tail = F) ptukey(t.ad, nmeans=4, df=df.within, lower.tail = F) ptukey(t.bc, nmeans=4, df=df.within, lower.tail = F) ptukey(t.bd, nmeans=4, df=df.within, lower.tail = F) ptukey(t.cd, nmeans=4, df=df.within, lower.tail = F) TukeyHSD(a.res, conf.level=.95)
plot(TukeyHSD(a.res, conf.level=.95), las = 2)
pairwise.t.test(dat$values, dat$group, p.adj = "bonf")
post hoc test: output
> mean.a
[1] 26
> mean.b
[1] 25
> mean.c
[1] 22
> mean.d
[1] 20
>
> d.ab <- mean.a - mean.b
> d.ac <- mean.a - mean.c
> d.ad <- mean.a - mean.d
> d.bc <- mean.b - mean.c
> d.bd <- mean.b - mean.d
> d.cd <- mean.c - mean.d
>
> d.ab
[1] 1
> d.ac
[1] 4
> d.ad
[1] 6
> d.bc
[1] 3
> d.bd
[1] 5
> d.cd
[1] 2
>
> # mse (ms within) from the a.res.sum output
> # a.res.sum == summary(aov(values ~ group, data=dat))
> a.res.sum
Df Sum Sq Mean Sq F value Pr(>F)
group 3 682 227.50 6.344 0.000508 ***
Residuals 116 4160 35.86
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> # mse = 50
> mse <- 35.86
> # 혹은 fansy way from dat data.frame
> # df.a 는 각 그룹의 df (모든 그룹의 df값이 같으므로 df.a를 사용)
> tapply(dat$values, dat$group, var) # 각 그룹의 분산
A B C D
40.00000 34.48276 31.03448 37.93103
> tapply(dat$values, dat$group, var) * df.a # 각 그룹의 SS
A B C D
1160 1000 900 1100
> # 이 분산값에 df값을 곱하면 각 그룹의 SS값
> # 이 값들을 sum 하면 총 ss.within 값
> sse.ch <- sum(tapply(dat$values, dat$group, var) * df.a)
> sse.ch
[1] 4160
> mse.ch <- sse.ch / (df.a*4) # 이 값에 df.a * 4 하면 ms.within값
> mse.ch
[1] 35.86207
> mse <- 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.ad <- d.ad / se
> t.bc <- d.bc / se
> t.bd <- d.bd / se
> t.cd <- d.cd / se
>
> t.ab
[1] 0.9146248
> t.ac
[1] 3.658499
> t.ad
[1] 5.487749
> t.bc
[1] 2.743874
> t.bd
[1] 4.573124
> t.cd
[1] 1.82925
>
> # 이제 위의 점수를 .05 레벨에서 비교할 점수를 찾아 비교한다
> # qtukey 펑션을 이용한다
> t.crit <- qtukey( .95, nmeans = 3, df = 45)
> t.crit
[1] 3.427507
>
> # 혹은 R이 보통 제시한 거꾸로 방법으로 보면
> ptukey(t.ab, nmeans=4, df=df.within, lower.tail = F)
[1] 0.9164944
> ptukey(t.ac, nmeans=4, df=df.within, lower.tail = F)
[1] 0.05255324
> ptukey(t.ad, nmeans=4, df=df.within, lower.tail = F)
[1] 0.0009861641
> ptukey(t.bc, nmeans=4, df=df.within, lower.tail = F)
[1] 0.2171272
> ptukey(t.bd, nmeans=4, df=df.within, lower.tail = F)
[1] 0.008539966
> ptukey(t.cd, nmeans=4, df=df.within, lower.tail = F)
[1] 0.5689321
>
> TukeyHSD(a.res, conf.level=.95)
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 -5.030484 3.0304845 0.9164944
C-A -4 -8.030484 0.0304845 0.0525532
D-A -6 -10.030484 -1.9695155 0.0009862
C-B -3 -7.030484 1.0304845 0.2171272
D-B -5 -9.030484 -0.9695155 0.0085400
D-C -2 -6.030484 2.0304845 0.5689321
> plot(TukeyHSD(a.res, conf.level=.95), las = 2)
> 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.0655 0.3287 -
D 0.0010 0.0096 1.0000
P value adjustment method: bonferroni
>
>
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.total ss.between r.sq <- ss.between / ss.total r.sq # then . . . . lm.res <- lm(values ~ group, data = dat) summary(lm.res) anova(lm.res)
R square: output
> ss.total
[1] 4842.5
> ss.between
[1] 682.5
> r.sq <- ss.between / ss.total
> r.sq
[1] 0.1409396
>
> # then . . . .
>
> lm.res <- lm(values ~ group, data = dat)
> summary(lm.res)
Call:
lm(formula = values ~ group, data = dat)
Residuals:
Min 1Q Median 3Q Max
-14.8928 -4.3342 -0.3732 4.4517 13.0126
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.000 1.093 23.780 < 2e-16 ***
groupB -1.000 1.546 -0.647 0.519080
groupC -4.000 1.546 -2.587 0.010919 *
groupD -6.000 1.546 -3.880 0.000174 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.988 on 116 degrees of freedom
Multiple R-squared: 0.1409, Adjusted R-squared: 0.1187
F-statistic: 6.344 on 3 and 116 DF, p-value: 0.0005079
> anova(lm.res)
Analysis of Variance Table
Response: values
Df Sum Sq Mean Sq F value Pr(>F)
group 3 682.5 227.500 6.3437 0.0005079 ***
Residuals 116 4160.0 35.862
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
c/ms/2025/w07_anova_note.1744267832.txt.gz · Last modified: 2025/04/10 15:50 by hkimscil