COMMunication
RESearch.NET

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== ANOVA in R ======
이 문서는 http://commres.net/wiki/r/anova 의 내용을 변형해서 
다시 만든 문서입니다. http://commres.net/wiki/r/anova 은 
세 그룹 간의 차이를 보는 것이고 이 문서는 4그룹 간의 차이를 보는
문서입니다. 
<code>
# 
# 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(900/(na-1)))
B <- rnorm2(nb, mean.b, sqrt(900/(nb-1)))
C <- rnorm2(nc, mean.c, sqrt(900/(nc-1)))
D <- rnorm2(nd, mean.d, sqrt(900/(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

</code>
<code>
# 그러나 정확히 어떤 그룹에서 차이가 나는지는 판단해주지 않음 
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)
</code>

===== ANOVA in R: Output =====
<code>
> # 
> # 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(900/(na-1)))
> B <- rnorm2(nb, mean.b, sqrt(900/(nb-1)))
> C <- rnorm2(nc, mean.c, sqrt(900/(nc-1)))
> D <- rnorm2(nd, mean.d, sqrt(900/(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.37871
 [2,] 30.23619
 [3,] 22.05233
 [4,] 27.98186
 [5,] 28.62468
 [6,] 34.38014
 [7,] 30.67844
 [8,] 25.80092
 [9,] 32.66698
[10,] 25.06399
[11,] 30.06274
[12,] 21.25288
[13,] 36.07231
[14,] 16.77241
[15,] 24.97448
[16,] 25.26349
[17,] 20.88676
[18,] 26.94242
[19,] 21.10069
[20,] 12.88194
[21,] 25.46073
[22,] 31.27674
[23,] 24.76580
[24,] 16.79173
[25,] 31.51620
[26,] 17.14866
[27,] 29.66682
[28,] 25.75701
[29,] 29.66796
[30,] 29.87397
attr(,"scaled:center")
[1] -0.08287722
attr(,"scaled:scale")
[1] 0.8355107
> B
          [,1]
 [1,] 30.62357
 [2,] 27.27939
 [3,] 31.23686
 [4,] 14.50486
 [5,] 32.22519
 [6,] 21.82949
 [7,] 26.67567
 [8,] 30.76150
 [9,] 16.68531
[10,] 28.19891
[11,] 28.38350
[12,] 29.88106
[13,] 13.16769
[14,] 23.26793
[15,] 19.76032
[16,] 27.95159
[17,] 28.85313
[18,] 21.92882
[19,] 24.18798
[20,] 17.70481
[21,] 19.51664
[22,] 24.27280
[23,] 28.90183
[24,] 18.17715
[25,] 29.83134
[26,] 20.05465
[27,] 26.66153
[28,] 31.89910
[29,] 32.13759
[30,] 23.43977
attr(,"scaled:center")
[1] -0.1405676
attr(,"scaled:scale")
[1] 1.025799
> C
          [,1]
 [1,] 21.04268
 [2,] 14.58761
 [3,] 18.90352
 [4,] 23.12972
 [5,] 24.86854
 [6,] 24.66800
 [7,] 18.64315
 [8,] 23.33405
 [9,] 22.17603
[10,] 22.07975
[11,] 30.96436
[12,] 31.58129
[13,] 28.96433
[14,] 22.06416
[15,] 12.30153
[16,] 16.68289
[17,] 24.19514
[18,] 15.33454
[19,] 23.27483
[20,] 22.21340
[21,] 32.88316
[22,] 28.73176
[23,] 19.63225
[24,] 19.45001
[25,] 12.80615
[26,] 25.13846
[27,] 22.52944
[28,] 30.05695
[29,] 26.55883
[30,] 31.20348
attr(,"scaled:center")
[1] 0.08931913
attr(,"scaled:scale")
[1] 0.9936574
> D
           [,1]
 [1,] 29.117527
 [2,] 21.287702
 [3,] 19.406172
 [4,] 17.338050
 [5,] 23.108952
 [6,] 16.932891
 [7,] 18.923097
 [8,] 29.345116
 [9,] 24.349527
[10,] 16.795435
[11,] 23.028628
[12,] 18.074871
[13,] 33.770366
[14,] 28.238105
[15,] 25.785121
[16,] 20.157663
[17,] 21.990432
[18,]  8.910282
[19,] 18.797986
[20,] 31.193353
[21,] 18.214726
[22,] 21.215850
[23,] 20.581063
[24,] 30.711279
[25,] 25.911186
[26,] 17.041703
[27,] 17.853327
[28,] 16.703969
[29,] 18.081180
[30,] 27.134442
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.37871   A
2 30.23619   A
3 22.05233   A
4 27.98186   A
5 28.62468   A
6 34.38014   A
> colnames(dat)[1] <- "values"
> colnames(dat)[2] <- "group"
> head(dat)
    values group
1 24.37871     A
2 30.23619     A
3 22.05233     A
4 27.98186     A
5 28.62468     A
6 34.38014     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] 24
> var.total
[1] 32.77311
> ms.total
[1] 32.77311
> df.total 
[1] 119
> ss.total
[1] 3900
> ss.total.check
[1] 3900
> 
> # 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] 23
> mean.d
[1] 22
> 
> # 그룹 간의 차이에서 나타나는 분산
> # 수업시간에 설명을 잘 들을 것
> 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] 120
> length(B) * ((mean.total - mean.b)^2)
[1] 30
> length(C) * ((mean.total - mean.c)^2)
[1] 30
> length(D) * ((mean.total - mean.d)^2)
[1] 120
> 
> 
> 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] 300
> # 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,] 3.222222
> 
> # 컴퓨터 계산이 쉬워지기 전에는 아래처럼 0.5 level 
> # 에서의 f값을 구한 후 이것과 계산된 f값을 비교해봤었다. 
> qf(.05, df1 = df.between, df2 = df.within, lower.tail = FALSE)
[1] 2.682809
> f.calculated
         [,1]
[1,] 3.222222
> # 위에서 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.02527283
> f.calculated
         [,1]
[1,] 3.222222
> 
> # 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,] 3.222222
> f.calculated.pvalue
           [,1]
[1,] 0.02527283
> 1 - f.calculated.pvalue
          [,1]
[1,] 0.9747272
> 
> 
> 
> # Now check this
> ss.total
[1] 3900
> ss.between
[1] 300
> ss.within
     [,1]
[1,] 3600
> ss.total 
[1] 3900
> ss.between + ss.within
     [,1]
[1,] 3900
> 
> # 한편 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    300  100.00   3.222 0.0253 *
Residuals   116   3600   31.03                 
---
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.4883 -      -     
C 0.0392 0.1671 -     
D 0.0063 0.0392 0.4883

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.000 -     -    
C 0.235 1.000 -    
D 0.038 0.235 1.000

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.977 -     -    
C 0.196 0.501 -    
D 0.038 0.196 0.977

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 -4.749401  2.7494006 0.8987322
C-A   -3 -6.749401  0.7494006 0.1638896
D-A   -4 -7.749401 -0.2505994 0.0316929
C-B   -2 -5.749401  1.7494006 0.5078699
D-B   -3 -6.749401  0.7494006 0.1638896
D-C   -1 -4.749401  2.7494006 0.8987322

> 
> boxplot(dat$values~dat$group)
> f.calculated
         [,1]
[1,] 3.222222
> f.calculated.pvalue
           [,1]
[1,] 0.02527283
> 
> 
> # how much IV explains the DV 
> # in terms of SS?
> r.square <- ss.between / ss.total
> eta  <- r.square
> eta
[1] 0.07692308
> 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 
-13.1181  -3.9308  -0.3568   3.9491  11.7704 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   26.000      1.017  25.563  < 2e-16 ***
dat$groupB    -1.000      1.438  -0.695  0.48831    
dat$groupC    -3.000      1.438  -2.086  0.03920 *  
dat$groupD    -4.000      1.438  -2.781  0.00633 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.571 on 116 degrees of freedom
Multiple R-squared:  0.07692,	Adjusted R-squared:  0.05305 
F-statistic: 3.222 on 3 and 116 DF,  p-value: 0.02527

> summary(a.res)
             Df Sum Sq Mean Sq F value Pr(>F)  
group         3    300  100.00   3.222 0.0253 *
Residuals   116   3600   31.03                 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> 
></code>

{{:c:ms:2025:pasted:20250414-084011.png}}
{{:c:ms:2025:pasted:20250414-084020.png}}
{{:c:ms:2025:pasted:20250414-084030.png}}
====== Post hoc test ======
[[:post hoc test]]
<code>
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 = 4, df = 30 * 4)
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)
</code>

<code>
plot(TukeyHSD(a.res, conf.level=.95), las = 2)
</code>

<code>
pairwise.t.test(dat$values, dat$group, p.adj = "bonf")
</code>
===== post hoc test: output =====
<code>
> 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 = 4, df = 30 * 4)
> t.crit
[1] 3.684589
> 
> # 혹은 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 
> 
> 

</code>

{{:c:ms:2025:pasted:20250410-154903.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.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)
</code>

===== R square: output =====
<code>
> ss.total
[1] 3900
> ss.between
[1] 300
> r.sq <- ss.between / ss.total
> r.sq
[1] 0.07692308
> 
> # then . . . . 
> 
> lm.res <- lm(values ~ group, data = dat)
> summary(lm.res)

Call:
lm(formula = values ~ group, data = dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.1181  -3.9308  -0.3568   3.9491  11.7704 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   26.000      1.017  25.563  < 2e-16 ***
groupB        -1.000      1.438  -0.695  0.48831    
groupC        -3.000      1.438  -2.086  0.03920 *  
groupD        -4.000      1.438  -2.781  0.00633 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.571 on 116 degrees of freedom
Multiple R-squared:  0.07692,	Adjusted R-squared:  0.05305 
F-statistic: 3.222 on 3 and 116 DF,  p-value: 0.02527

> anova(lm.res)
Analysis of Variance Table

Response: values
           Df Sum Sq Mean Sq F value  Pr(>F)  
group       3    300 100.000  3.2222 0.02527 *
Residuals 116   3600  31.034                  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 

</code>