c:ms:2025:w07.2_factorial_anova
Differences
This shows you the differences between two versions of the page.
| Both sides previous revisionPrevious revisionNext revision | Previous revision | ||
| c:ms:2025:w07.2_factorial_anova [2025/04/20 23:02] – hkimscil | c:ms:2025:w07.2_factorial_anova [2025/05/06 22:25] (current) – hkimscil | ||
|---|---|---|---|
| Line 1: | Line 1: | ||
| + | ====== Two-way ANOVA ====== | ||
| + | Factorial ANOVA | ||
| < | < | ||
| + | rm(list=ls(all=TRUE)) | ||
| + | |||
| ################################################# | ################################################# | ||
| # two-way anova | # two-way anova | ||
| Line 5: | Line 9: | ||
| ################################################# | ################################################# | ||
| - | n.a.group <- 3 # a treatment 숫자 | + | n.a.group <- 3 # a treatment 숫자 |
| - | n.b.group <- 2 # b 그룹 숫자 | + | n.b.group <- 2 # b 그룹 숫자 |
| n.sub <- 30 # 총 샘플 숫자 | n.sub <- 30 # 총 샘플 숫자 | ||
| n.sub/ | n.sub/ | ||
| Line 12: | Line 16: | ||
| # 데이터 생성 | # 데이터 생성 | ||
| set.seed(9) | set.seed(9) | ||
| - | a <- gl(3, 10, n.sub, labels=c(' | + | a <- gl(n.a.group, |
| - | b <- gl(2, 5, n.sub, labels=c(' | + | n.sub/ |
| + | | ||
| + | | ||
| + | b <- gl(n.b.group, | ||
| + | (n.sub/ | ||
| + | | ||
| + | | ||
| a | a | ||
| b | b | ||
| y <- rnorm(30, mean=10) + | y <- rnorm(30, mean=10) + | ||
| - | 3.14 * (a=='a1' & b=='b2') + | + | 3.14 * (a=='drg1' & b=='exerc') + |
| - | 1.43 * (a=='a3' & b=='b2' | + | 1.43 * (a=='drg3' & b=='exerc' |
| y | y | ||
| Line 29: | Line 39: | ||
| table(a,b) | table(a,b) | ||
| tapply(y, list(a,b), mean) # 각 셀에서의 평균 | tapply(y, list(a,b), mean) # 각 셀에서의 평균 | ||
| - | n.within.each <- tapply(y, list(a,b), length) | + | tapply(y, a, mean) |
| + | tapply(y, b, mean) | ||
| + | n.within.each <- tapply(y, list(a,b), length) # the same as table(a, b) | ||
| n.within.each | n.within.each | ||
| df.within.each <- n.within.each - 1 # 각 셀에서의 샘플숫자 | df.within.each <- n.within.each - 1 # 각 셀에서의 샘플숫자 | ||
| Line 63: | Line 75: | ||
| ss.tot <- var.tot * df.tot | ss.tot <- var.tot * df.tot | ||
| ms.tot <- ss.tot/ | ms.tot <- ss.tot/ | ||
| + | |||
| + | mean.tot | ||
| + | var.tot | ||
| + | df.tot | ||
| + | ss.tot | ||
| + | ms.tot | ||
| + | ss.tot/ | ||
| + | |||
| ## between | ## between | ||
| Line 114: | Line 134: | ||
| df.within | df.within | ||
| - | ms.tot <- ss.tot/ | + | ms.tot <- ss.tot / df.tot # we did it above |
| - | ms.bet <- ss.bet | + | ms.bet <- ss.bet |
| ms.a <- ss.a / df.a | ms.a <- ss.a / df.a | ||
| ms.b <- ss.b / df.b | ms.b <- ss.b / df.b | ||
| Line 146: | Line 166: | ||
| # qf(alpha, df.a, df.within, lower.tail = F) | # qf(alpha, df.a, df.within, lower.tail = F) | ||
| # 도 마찬가지 | # 도 마찬가지 | ||
| + | f.a.crit <- qf(ci, df.a, df.within) | ||
| + | f.a > f.a.crit # 만약 f.a가 크다면 | ||
| + | # 혹은 f.a 값에 해당하는 percentage를 구하면 | ||
| pf(f.a, df.a, df.within, lower.tail = F) | pf(f.a, df.a, df.within, lower.tail = F) | ||
| f.b | f.b | ||
| qf(ci, df.b, df.within) | qf(ci, df.b, df.within) | ||
| + | f.b > qf(ci, df.b, df.within) | ||
| pf(f.b, df.b, df.within, lower.tail = F) | pf(f.b, df.b, df.within, lower.tail = F) | ||
| f.ab | f.ab | ||
| qf(ci, df.ab, df.within) | qf(ci, df.ab, df.within) | ||
| + | f.ab > qf(ci, df.ab, df.within) | ||
| pf(f.ab, df.ab, df.within, lower.tail = F) | pf(f.ab, df.ab, df.within, lower.tail = F) | ||
| # aov result | # aov result | ||
| - | summary(aov.dat) | + | aov.dat.all <- aov(y ~ a * b, data = dat) # anova test |
| + | summary(aov.dat.all) # summary of the test output | ||
| + | aov.dat.noint <- aov(y ~ a + b, data = dat) | ||
| + | # interaction이 갖는 SS와 df를 residuals이 갖음을 주목 | ||
| + | summary(aov.dat.noint) | ||
| + | aov.dat.all2 <- aov(y ~ a + b + a:b, data = dat) | ||
| + | summary(aov.dat.all2) | ||
| + | aov.dat.all # what is standard error of residuals | ||
| + | |||
| + | # post-hoc test | ||
| + | TukeyHSD(aov.dat.all) | ||
| + | TukeyHSD(aov.dat.all, | ||
| + | TukeyHSD(aov.dat.all, | ||
| + | TukeyHSD(aov.dat.all, | ||
| + | boxplot(dat$y~dat$a) | ||
| + | boxplot(dat$y~dat$b) | ||
| + | plot(TukeyHSD(aov.dat.all, | ||
| + | las = 2 # rotate x-axis ticks | ||
| + | ) | ||
| + | |||
| + | |||
| + | # | ||
| + | # | ||
| + | # see another e.g. at | ||
| + | # https:// | ||
| + | # | ||
| + | # | ||
| </ | </ | ||
| - | ====== | + | |
| + | ====== | ||
| < | < | ||
| + | > rm(list=ls(all=TRUE)) | ||
| + | > | ||
| > ################################################# | > ################################################# | ||
| > # two-way anova | > # two-way anova | ||
| Line 166: | Line 221: | ||
| > ################################################# | > ################################################# | ||
| > | > | ||
| - | > n.a.group <- 3 # a treatment 숫자 | + | > n.a.group <- 3 # a treatment 숫자 |
| - | > n.b.group <- 2 # b 그룹 숫자 | + | > n.b.group <- 2 # b 그룹 숫자 |
| > n.sub <- 30 # 총 샘플 숫자 | > n.sub <- 30 # 총 샘플 숫자 | ||
| > n.sub/ | > n.sub/ | ||
| Line 174: | Line 229: | ||
| > # 데이터 생성 | > # 데이터 생성 | ||
| > set.seed(9) | > set.seed(9) | ||
| - | > a <- gl(3, 10, n.sub, labels=c(' | + | > a <- gl(n.a.group, |
| - | > b <- gl(2, 5, n.sub, labels=c(' | + | + |
| + | + n.sub, | ||
| + | + labels=c(' | ||
| + | > b <- gl(n.b.group, | ||
| + | + | ||
| + | + n.sub, | ||
| + | + labels=c(' | ||
| > a | > a | ||
| - | | + | |
| - | [19] a2 a2 a3 a3 a3 a3 a3 a3 a3 a3 a3 a3 | + | [22] drg3 drg3 drg3 drg3 drg3 drg3 drg3 drg3 drg3 |
| - | Levels: | + | Levels: |
| > b | > b | ||
| - | | + | |
| - | [19] b2 b2 b1 b1 b1 b1 b1 b2 b2 b2 b2 b2 | + | [19] exerc exerc noex noex noex noex noex exerc exerc exerc exerc exerc |
| - | Levels: | + | Levels: |
| > y <- rnorm(30, mean=10) + | > y <- rnorm(30, mean=10) + | ||
| - | + 3.14 * (a=='a1' & b=='b2') + | + | + 3.14 * (a=='drg1' & b=='exerc') + |
| - | + 1.43 * (a=='a3' & b=='b2' | + | + 1.43 * (a=='drg3' & b=='exerc' |
| > y | > y | ||
| - | | + | |
| - | | + | [11] 11.277571 |
| - | [11] 11.277571 | + | [21] 11.756993 10.182252 |
| - | [16] | + | |
| - | [21] 11.756993 10.182252 | + | |
| - | [26] 14.111990 11.652524 10.723328 11.847213 11.799557 | + | |
| > | > | ||
| > dat <- data.frame(a, | > dat <- data.frame(a, | ||
| > dat | > dat | ||
| - | | + | |
| - | 1 | + | 1 |
| - | 2 | + | 2 |
| - | 3 | + | 3 |
| - | 4 | + | 4 |
| - | 5 | + | 5 |
| - | 6 | + | 6 |
| - | 7 | + | 7 |
| - | 8 | + | 8 |
| - | 9 | + | 9 |
| - | 10 a1 b2 12.777063 | + | 10 drg1 exerc 12.777063 |
| - | 11 a2 b1 11.277571 | + | 11 drg2 noex 11.277571 |
| - | 12 a2 b1 | + | 12 drg2 noex |
| - | 13 a2 b1 10.071054 | + | 13 drg2 noex 10.071054 |
| - | 14 a2 b1 | + | 14 drg2 noex |
| - | 15 a2 b1 11.845257 | + | 15 drg2 noex 11.845257 |
| - | 16 a2 b2 | + | 16 drg2 exerc |
| - | 17 a2 b2 | + | 17 drg2 exerc |
| - | 18 a2 b2 | + | 18 drg2 exerc |
| - | 19 a2 b2 10.887884 | + | 19 drg2 exerc 10.887884 |
| - | 20 a2 b2 | + | 20 drg2 exerc |
| - | 21 a3 b1 11.756993 | + | 21 drg3 noex 11.756993 |
| - | 22 a3 b1 10.182252 | + | 22 drg3 noex 10.182252 |
| - | 23 a3 b1 | + | 23 drg3 noex |
| - | 24 a3 b1 10.926422 | + | 24 drg3 noex 10.926422 |
| - | 25 a3 b1 | + | 25 drg3 noex |
| - | 26 a3 b2 14.111990 | + | 26 drg3 exerc 14.111990 |
| - | 27 a3 b2 11.652524 | + | 27 drg3 exerc 11.652524 |
| - | 28 a3 b2 10.723328 | + | 28 drg3 exerc 10.723328 |
| - | 29 a3 b2 11.847213 | + | 29 drg3 exerc 11.847213 |
| - | 30 a3 b2 11.799557 | + | 30 drg3 exerc 11.799557 |
| > # aov.dat <- aov(y~a*b) # anova test | > # aov.dat <- aov(y~a*b) # anova test | ||
| > # summary(aov.dat) # summary of the test output | > # summary(aov.dat) # summary of the test output | ||
| Line 233: | Line 291: | ||
| > # hand calculation | > # hand calculation | ||
| > table(a,b) | > table(a,b) | ||
| - | | + | |
| - | a b1 b2 | + | a noex exerc |
| - | | + | |
| - | | + | |
| - | | + | |
| > tapply(y, list(a,b), mean) # 각 셀에서의 평균 | > tapply(y, list(a,b), mean) # 각 셀에서의 평균 | ||
| - | | + | |
| - | a1 | + | drg1 |
| - | a2 10.491789 | + | drg2 10.491789 |
| - | a3 10.381089 12.026922 | + | drg3 10.381089 12.026922 |
| > tapply(y, a, mean) | > tapply(y, a, mean) | ||
| - | | + | |
| 11.350981 | 11.350981 | ||
| > tapply(y, b, mean) | > tapply(y, b, mean) | ||
| - | b1 | + | noex exerc |
| 10.18655 11.45709 | 10.18655 11.45709 | ||
| > n.within.each <- tapply(y, list(a,b), length) # the same as table(a, b) | > n.within.each <- tapply(y, list(a,b), length) # the same as table(a, b) | ||
| > n.within.each | > n.within.each | ||
| - | b1 b2 | + | noex exerc |
| - | a1 | + | drg1 |
| - | a2 | + | drg2 |
| - | a3 | + | drg3 |
| > df.within.each <- n.within.each - 1 # 각 셀에서의 샘플숫자 | > df.within.each <- n.within.each - 1 # 각 셀에서의 샘플숫자 | ||
| > df.within.each | > df.within.each | ||
| - | b1 b2 | + | noex exerc |
| - | a1 | + | drg1 |
| - | a2 | + | drg2 |
| - | a3 | + | drg3 |
| > df.within <- sum(df.within.each) # df within | > df.within <- sum(df.within.each) # df within | ||
| > df.within | > df.within | ||
| Line 267: | Line 325: | ||
| > var.within <- tapply(y, list(a,b), var) # var.within | > var.within <- tapply(y, list(a,b), var) # var.within | ||
| > var.within | > var.within | ||
| - | | + | |
| - | a1 0.2628787 0.7362999 | + | drg1 0.2628787 0.7362999 |
| - | a2 1.0308917 1.6504481 | + | drg2 1.0308917 1.6504481 |
| - | a3 0.9510727 1.5677577 | + | drg3 0.9510727 1.5677577 |
| > ss.within.each <- tapply(y, list(a,b), var) * df.within.each | > ss.within.each <- tapply(y, list(a,b), var) * df.within.each | ||
| > ss.within.each | > ss.within.each | ||
| - | b1 b2 | + | noex exerc |
| - | a1 1.051515 2.945200 | + | drg1 1.051515 2.945200 |
| - | a2 4.123567 6.601793 | + | drg2 4.123567 6.601793 |
| - | a3 3.804291 6.271031 | + | drg3 3.804291 6.271031 |
| > ss.within <- sum(ss.within.each) # ss.within | > ss.within <- sum(ss.within.each) # ss.within | ||
| > ss.within | > ss.within | ||
| Line 290: | Line 348: | ||
| > mean.b <- tapply(y, list(b), mean) | > mean.b <- tapply(y, list(b), mean) | ||
| > mean.a | > mean.a | ||
| - | | + | |
| 11.350981 | 11.350981 | ||
| > mean.b | > mean.b | ||
| - | b1 | + | noex exerc |
| 10.18655 11.45709 | 10.18655 11.45709 | ||
| > | > | ||
| Line 299: | Line 357: | ||
| > var.b <- tapply(y, list(b), var) | > var.b <- tapply(y, list(b), var) | ||
| > var.a | > var.a | ||
| - | a1 | + | drg1 |
| 3.521366 1.567182 1.871915 | 3.521366 1.567182 1.871915 | ||
| > var.b | > var.b | ||
| - | | + | |
| 0.7773781 3.7300196 | 0.7773781 3.7300196 | ||
| > | > | ||
| Line 328: | Line 386: | ||
| > mean.each <- tapply(y, list(a,b), mean) | > mean.each <- tapply(y, list(a,b), mean) | ||
| > mean.each | > mean.each | ||
| - | | + | |
| - | a1 | + | drg1 |
| - | a2 10.491789 | + | drg2 10.491789 |
| - | a3 10.381089 12.026922 | + | drg3 10.381089 12.026922 |
| > mean.tot <- mean(y) | > mean.tot <- mean(y) | ||
| > mean.tot | > mean.tot | ||
| Line 339: | Line 397: | ||
| > n.b.each <- tapply(y, list(b), length) | > n.b.each <- tapply(y, list(b), length) | ||
| > n.each | > n.each | ||
| - | b1 b2 | + | noex exerc |
| - | a1 | + | drg1 |
| - | a2 | + | drg2 |
| - | a3 | + | drg3 |
| > n.a.each | > n.a.each | ||
| - | a1 a2 a3 | + | drg1 drg2 drg3 |
| - | 10 10 10 | + | 10 |
| > n.b.each | > n.b.each | ||
| - | b1 b2 | + | noex exerc |
| - | 15 15 | + | |
| > | > | ||
| > | > | ||
| Line 483: | Line 541: | ||
| Residuals | Residuals | ||
| --- | --- | ||
| - | Signif. codes: | + | Signif. codes: |
| - | 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 | + | |
| > aov.dat.noint <- aov(y ~ a + b, data = dat) | > aov.dat.noint <- aov(y ~ a + b, data = dat) | ||
| > # interaction이 갖는 SS와 df를 residuals이 갖음을 주목 | > # interaction이 갖는 SS와 df를 residuals이 갖음을 주목 | ||
| Line 493: | Line 550: | ||
| Residuals | Residuals | ||
| --- | --- | ||
| - | Signif. codes: | + | Signif. codes: |
| - | 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 | + | |
| > aov.dat.all2 <- aov(y ~ a + b + a:b, data = dat) | > aov.dat.all2 <- aov(y ~ a + b + a:b, data = dat) | ||
| > summary(aov.dat.all2) | > summary(aov.dat.all2) | ||
| Line 503: | Line 559: | ||
| Residuals | Residuals | ||
| --- | --- | ||
| - | Signif. codes: | + | Signif. codes: |
| - | 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 | + | |
| > aov.dat.all # what is standard error of residuals | > aov.dat.all # what is standard error of residuals | ||
| Call: | Call: | ||
| Line 517: | Line 572: | ||
| Estimated effects may be unbalanced | Estimated effects may be unbalanced | ||
| > | > | ||
| + | > # post-hoc test | ||
| + | > TukeyHSD(aov.dat.all) | ||
| + | Tukey multiple comparisons of means | ||
| + | 95% family-wise confidence level | ||
| + | |||
| + | Fit: aov(formula = y ~ a * b, data = dat) | ||
| + | |||
| + | $a | ||
| + | diff | ||
| + | drg2-drg1 -1.4405079 -2.575730 -0.3052857 0.0111258 | ||
| + | drg3-drg1 -0.1469757 -1.282198 | ||
| + | drg3-drg2 | ||
| + | |||
| + | $b | ||
| + | | ||
| + | exerc-noex 1.270533 0.5044869 2.03658 0.0022274 | ||
| + | |||
| + | $`a:b` | ||
| + | diff lwr upr p adj | ||
| + | drg2: | ||
| + | drg3: | ||
| + | drg1: | ||
| + | drg2: | ||
| + | drg3: | ||
| + | drg3: | ||
| + | drg1: | ||
| + | drg2: | ||
| + | drg3: | ||
| + | drg1: | ||
| + | drg2: | ||
| + | drg3: | ||
| + | drg2: | ||
| + | drg3: | ||
| + | drg3: | ||
| + | |||
| + | > TukeyHSD(aov.dat.all, | ||
| + | Tukey multiple comparisons of means | ||
| + | 95% family-wise confidence level | ||
| + | |||
| + | Fit: aov(formula = y ~ a * b, data = dat) | ||
| + | |||
| + | $a | ||
| + | diff | ||
| + | drg2-drg1 -1.4405079 -2.575730 -0.3052857 0.0111258 | ||
| + | drg3-drg1 -0.1469757 -1.282198 | ||
| + | drg3-drg2 | ||
| + | |||
| + | > TukeyHSD(aov.dat.all, | ||
| + | Tukey multiple comparisons of means | ||
| + | 95% family-wise confidence level | ||
| + | |||
| + | Fit: aov(formula = y ~ a * b, data = dat) | ||
| + | |||
| + | $b | ||
| + | | ||
| + | exerc-noex 1.270533 0.5044869 2.03658 0.0022274 | ||
| + | |||
| + | > TukeyHSD(aov.dat.all, | ||
| + | Tukey multiple comparisons of means | ||
| + | 95% family-wise confidence level | ||
| + | |||
| + | Fit: aov(formula = y ~ a * b, data = dat) | ||
| + | |||
| + | $`a:b` | ||
| + | diff lwr upr p adj | ||
| + | drg2: | ||
| + | drg3: | ||
| + | drg1: | ||
| + | drg2: | ||
| + | drg3: | ||
| + | drg3: | ||
| + | drg1: | ||
| + | drg2: | ||
| + | drg3: | ||
| + | drg1: | ||
| + | drg2: | ||
| + | drg3: | ||
| + | drg2: | ||
| + | drg3: | ||
| + | drg3: | ||
| + | |||
| + | > boxplot(dat$y~dat$a) | ||
| + | > boxplot(dat$y~dat$b) | ||
| + | > plot(TukeyHSD(aov.dat.all, | ||
| + | + las = 2 # rotate x-axis ticks | ||
| + | + ) | ||
| > | > | ||
| > | > | ||
| + | > # | ||
| + | > # | ||
| + | > # see another e.g. at | ||
| + | > # https:// | ||
| + | > # | ||
| + | > # | ||
| + | > | ||
| + | > | ||
| </ | </ | ||
| + | ====== output interpretation ====== | ||
| + | |||
| + | {{: | ||
| + | {{: | ||
| + | {{: | ||
| + | {{: | ||
| + | {{: | ||
| + | |||
| + | 상호작용효과 포함 해석 | ||
| + | * drug1의 효과에 운동의 역할은 아주 중요한 것으로 밝혀졌다. | ||
| + | * 반면에 drug2, 3에 한해서는 운동의 역할은 제한적이었다. | ||
| + | * 특히 drg3의 경우 운동이 효과를 증대시키기는 하였지만 그것이 통계학적으로 의미있는 것은 아니었다. | ||
| + | * 또한 drg2의 경우 운동이 오히려 효과를 감소하는 것으로 나타났다. 그렇지만 이 또한 통계학적으로 의미가 있는 것은 아니었다. | ||
| + | * 이를 통해서 drg1을 처치받는 환자의 경우에는 운동을 꼭 하도록 하여야 함을 알 수 있고 | ||
| + | * drg2의 경우에는 그 효과가 미미하므로 권장을 하는 정도가 필요함을 알 수 있었다. | ||
| + | 혹은 | ||
| + | * drg1과 운동이 병행 된 경우에는 대부분의 다른 처방보다 월등한 효과가 있는 것으로 밝혀졌다. | ||
| + | * 운동이 병행된 drg1과 drg3의 | ||
| + | * drg2의 운동유무, | ||
| + | |||
| + | {{: | ||
c/ms/2025/w07.2_factorial_anova.1745190156.txt.gz · Last modified: by hkimscil
