Differences

This shows you the differences between two versions of the page.

--- r:repeated_measure_anova [2024/05/08 08:22] – [E.g. 2] hkimscil
+++ r:repeated_measure_anova [2025/05/12 08:25] (current) – [E.g.2 by hand] hkimscil
@@ Line 30: / Line 30: @@
 #view model summary
 summary(m.aov)
+# Error(patient/drug) = patient error embedded in drug
+# the same thing
+# the latter should be preferred
+m2.aov <- aov(response ~ drug
+              + Error(patient/drug),
+              data = df)
+summary(m2.aov)
+# check this
+m3.aov <- aov(response ~ drug, data = df)
+summary(m3.aov)
 # check probability level (pr)
 - pf(24.75886525, 3, 12)
+# or
+library(broom)
+fit <- tidy(m2.aov)
+fit
+fit$statistic[2]
+- pf(fit$statistic[2], 3, 12)
+# A one-way repeated measures ANOVA was conducted
+# on five individuals to examine the effect that
+# four different drugs had on response time.
+# Results showed that the type of drug used lead
+# to statistically significant differences in
+# response time (F(3, 12) = 24.76, p < 0.001).
 # A one-way repeated measures ANOVA was conducted
@@ Line 44: / Line 69: @@
 <code>
-> ###################################################
-> ###################################################
-> ###################################################
->
 > # data
 > df <- data.frame(patient=rep(1:5, each=4),
@@ Line 103: / Line 124: @@
 Residuals 12  112.8     9.4
 ---
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+>
+> # Error(patient/drug) = patient error embedded in drug
+> # the same thing
+> # the latter should be preferred
+> m2.aov <- aov(response ~ drug
++               + Error(patient/drug),
++               data = df)
+> summary(m2.aov)
+Error: patient
+          Df Sum Sq Mean Sq F value Pr(>F)
+Residuals  4  680.8   170.2
+Error: patient:drug
+          Df Sum Sq Mean Sq F value   Pr(>F)
+drug       3  698.2   232.7   24.76 1.99e-05 ***
+Residuals 12  112.8     9.4
+---
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+>
+> # check this
+> m3.aov <- aov(response ~ drug, data = df)
+> summary(m3.aov)
+            Df Sum Sq Mean Sq F value Pr(>F)
+drug         3  698.2   232.7   4.692 0.0155 *
+Residuals   16  793.6    49.6
+---
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 >
 > # check probability level (pr)
 > 1 - pf(24.75886525, 3, 12)
 [1] 1.992501e-05
->
+> # or
-> # A one-way repeated measures ANOVA was conducted
+> library(broom)
-> # on five individuals to examine the effect that
+> fit <- tidy(m2.aov)
-> # four different drugs had on response time.
+> fit
->
+# A tibble: 3 × 7
-> # Results showed that the type of drug used lead
+  stratum      term     df sumsq meansq statistic  p.value
-> # to statistically significant differences in
+  <chr>        <chr> <dbl> <dbl>  <dbl>     <dbl>    <dbl>
-> # response time (F(3, 12) = 24.76, p < 0.001).
+patient      Resi…     4  681.  170.       NA   NA
+patient:drug drug      3  698.  233.       24.8  1.99e-5
+patient:drug Resi…    12  113.    9.4      NA   NA
+> fit$statistic[2]
+[1] 24.75887
+> 1 - pf(fit$statistic[2], 3, 12)
+[1] 1.992501e-05
 >
 >
@@ Line 126: / Line 184: @@
                  movie=rep(1:3, times=5),
                  score=c(88, 84, 92,
 , 78, 90,
 , 94, 95,
 , 83, 88,
 , 90, 99))
 #view data
@@ Line 143: / Line 201: @@
              + Error(reviewer),
              data = movrev)
+m2.aov <- aov(score ~ movie
+              + Error(reviewer/movie),
+              data = movrev)
 #view model summary
 summary(m.aov)
+summary(m2.aov)
 # pairwise.t.test(movrev$score, movrev$movie, paired = T, p.adjust.method = "bonf")
 attach(movrev)
 pairwise.t.test(score, movie, paired = T, p.adjust.method = "bonf")
+detach(movrev)
+# or
+with(movrev,
+     pairwise.t.test(score, movie,
+                     paired = T,
+                     p.adjust.method = "bonf"))
+# the second
+movrev <- data.frame(reviewer=rep(1:5, each=3),
+                 movie=rep(1:3, times=5),
+                 score=c(88, 84, 92,
+, 78, 90,
+, 94, 95,
+, 83, 88,
+, 90, 99))
+#view data
+movrev
+write.csv(movrev, file="rep.meas.anova.eg.movie.review.csv")
+movrev$movie <- factor(movrev$movie)
+movrev$reviewer <- factor(movrev$reviewer)
+# Error(reviewer) = reviewer error should be isolated
+# The above is the same as Error(reviewer/movie)
+m.aov <- aov(score ~ movie
+             + Error(reviewer),
+             data = movrev)
+m2.aov <- aov(score ~ movie
+              + Error(reviewer/movie),
+              data = movrev)
+#view model summary
+summary(m.aov)
+summary(m2.aov)
+# pairwise.t.test(movrev$score, movrev$movie, paired = T, p.adjust.method = "bonf")
+attach(movrev)
+pairwise.t.test(score, movie, paired = T, p.adjust.method = "bonf")
+detach(movrev)
 # or
 with(movrev,
@@ Line 157: / Line 259: @@
 </code>
+<code>
+> # the second
+> movrev <- data.frame(reviewer=rep(1:5, each=3),
++                  movie=rep(1:3, times=5),
++                  score=c(88, 84, 92,
++                             76, 78, 90,
++                             78, 94, 95,
++                             80, 83, 88,
++                             82, 90, 99))
+>
+> #view data
+> movrev
+   reviewer movie score
+         1     1    88
+         1     2    84
+         1     3    92
+         2     1    76
+         2     2    78
+         2     3    90
+         3     1    78
+         3     2    94
+         3     3    95
+        4     1    80
+        4     2    83
+        4     3    88
+        5     1    82
+        5     2    90
+        5     3    99
+> write.csv(movrev, file="rep.meas.anova.eg.movie.review.csv")
+>
+> movrev$movie <- factor(movrev$movie)
+> movrev$reviewer <- factor(movrev$reviewer)
+>
+> # Error(reviewer) = reviewer error should be isolated
+> # The above is the same as Error(reviewer/movie)
+> m.aov <- aov(score ~ movie
++              + Error(reviewer),
++              data = movrev)
+> m2.aov <- aov(score ~ movie
++               + Error(reviewer/movie),
++               data = movrev)
+>
+> #view model summary
+> summary(m.aov)
+Error: reviewer
+          Df Sum Sq Mean Sq F value Pr(>F)
+Residuals  4  173.7   43.43
+Error: Within
+          Df Sum Sq Mean Sq F value  Pr(>F)
+movie      2  363.3  181.67   10.19 0.00632 **
+Residuals  8  142.7   17.83
+---
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+> summary(m2.aov)
+Error: reviewer
+          Df Sum Sq Mean Sq F value Pr(>F)
+Residuals  4  173.7   43.43
+Error: reviewer:movie
+          Df Sum Sq Mean Sq F value  Pr(>F)
+movie      2  363.3  181.67   10.19 0.00632 **
+Residuals  8  142.7   17.83
+---
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+>
+> # pairwise.t.test(movrev$score, movrev$movie, paired = T, p.adjust.method = "bonf")
+> attach(movrev)
+The following objects are masked from movrev (pos = 12):
+    movie, reviewer, score
+The following objects are masked from movrev (pos = 14):
+    movie, reviewer, score
+> pairwise.t.test(score, movie, paired = T, p.adjust.method = "bonf")
+	Pairwise comparisons using paired t tests
+data:  score and movie
+     2
+0.628 -
+0.029 0.060
+P value adjustment method: bonferroni
+> detach(movrev)
+> # or
+> with(movrev,
++      pairwise.t.test(score, movie,
++                      paired = T,
++                      p.adjust.method = "bonf"))
+	Pairwise comparisons using paired t tests
+data:  score and movie
+     2
+0.628 -
+0.029 0.060
+P value adjustment method: bonferroni
+>
+>
+</code>
+====== E.g.2 by hand ======
+<code>
+####
+# by hand
+####
+movrev
+tapply(movrev$score, list(reviewer, movie), mean) # 각 셀의 평균값
+m.tot <- mean(movrev$score)
+m.by.movie <- tapply(movrev$score, list(movie), mean)
+m.by.er <- tapply(movrev$score, list(reviewer), mean)
+ss.bet <- sum(5*(m.tot-m.by.movie)^2)
+df.bet <- 3-1 # 영화 가짓수 - 1
+ms.bet <- ss.bet / df.bet
+var.by.movie <- tapply(movrev$score, list(movie), var)
+ss.with <- sum(var.by.movie*4)
+df.with <- 4 + 4 + 4
+ms.with <- ss.with / df.with
+ms.bet
+ms.with
+ss.sub <- sum(3 * (m.tot-m.by.er)^2) # 3 treatments (movies)
+df.sub <- 5 - 1 # 5 persons - 1
+ms.sub <- ss.sub/df.sub
+ms.sub
+ss.res <- ss.with - ss.sub
+ss.res
+df.res <- 4 * 2 # (N-1)*(k-1)
+ms.res <- ss.res/df.res
+ms.res
+f.cal <- ms.bet/ms.res
+f.cal
+pf(f.cal, 2, 8, lower.tail = F)
+ss.bet
+ss.with
+ss.sub
+ss.res
+df.bet
+df.with
+df.sub
+df.res
+ms.bet
+ms.with
+ms.sub
+ms.res
+# check
+ss.with == ss.sub + ss.res
+df.with == df.sub + df.res
+summary(m.aov)
+#
+# no embedded
+# understand the meaning of embedded
+m0.aov <- aov(score ~ movie, data = movrev)
+summary(m0.aov)
+summary(m.aov)
+</code>
+<code>
+> ####
+> # by hand
+> ####
+> movrev
+   reviewer movie score
+         1     1    88
+         1     2    84
+         1     3    92
+         2     1    76
+         2     2    78
+         2     3    90
+         3     1    78
+         3     2    94
+         3     3    95
+        4     1    80
+        4     2    83
+        4     3    88
+        5     1    82
+        5     2    90
+        5     3    99
+> tapply(movrev$score, list(reviewer, movie), mean) # 각 셀의 평균값
+  2  3
+88 84 92
+76 78 90
+78 94 95
+80 83 88
+82 90 99
+> m.tot <- mean(movrev$score)
+>
+>
+> m.by.movie <- tapply(movrev$score, list(movie), mean)
+> m.by.er <- tapply(movrev$score, list(reviewer), mean)
+>
+> ss.bet <- sum(5*(m.tot-m.by.movie)^2)
+> df.bet <- 3-1 # 영화 가짓수 - 1
+> ms.bet <- ss.bet / df.bet
+>
+> var.by.movie <- tapply(movrev$score, list(movie), var)
+> ss.with <- sum(var.by.movie*4)
+> df.with <- 4 + 4 + 4
+> ms.with <- ss.with / df.with
+> ms.bet
+[1] 181.6667
+> ms.with
+[1] 26.36667
+>
+> ss.sub <- sum(3 * (m.tot-m.by.er)^2) # 3 treatments (movies)
+> df.sub <- 5 - 1 # 5 persons - 1
+> ms.sub <- ss.sub/df.sub
+> ms.sub
+[1] 43.43333
+>
+> ss.res <- ss.with - ss.sub
+> ss.res
+[1] 142.6667
+> df.res <- 4 * 2 # (N-1)*(k-1)
+> ms.res <- ss.res/df.res
+> ms.res
+[1] 17.83333
+>
+> f.cal <- ms.bet/ms.res
+> f.cal
+[1] 10.18692
+> pf(f.cal, 2, 8, lower.tail = F)
+[1] 0.006319577
+>
+> ss.bet
+[1] 363.3333
+> ss.with
+[1] 316.4
+> ss.sub
+[1] 173.7333
+> ss.res
+[1] 142.6667
+> df.bet
+[1] 2
+> df.with
+[1] 12
+> df.sub
+[1] 4
+> df.res
+[1] 8
+> ms.bet
+[1] 181.6667
+> ms.with
+[1] 26.36667
+> ms.sub
+[1] 43.43333
+> ms.res
+[1] 17.83333
+>
+> # check
+> ss.with == ss.sub + ss.res
+[1] TRUE
+> df.with == df.sub + df.res
+[1] TRUE
+> summary(m.aov)
+Error: reviewer
+          Df Sum Sq Mean Sq F value Pr(>F)
+Residuals  4  173.7   43.43
+Error: Within
+          Df Sum Sq Mean Sq F value  Pr(>F)
+movie      2  363.3  181.67   10.19 0.00632 **
+Residuals  8  142.7   17.83
+---
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+> #
+> # no embedded
+> # understand the meaning of embedded
+> m0.aov <- aov(score ~ movie, data = movrev)
+> summary(m0.aov)
+            Df Sum Sq Mean Sq F value Pr(>F)
+movie        2  363.3  181.67    6.89 0.0102 *
+Residuals   12  316.4   26.37
+---
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+> summary(m.aov)
+Error: reviewer
+          Df Sum Sq Mean Sq F value Pr(>F)
+Residuals  4  173.7   43.43
+Error: Within
+          Df Sum Sq Mean Sq F value  Pr(>F)
+movie      2  363.3  181.67   10.19 0.00632 **
+Residuals  8  142.7   17.83
+---
+Signif. codes:
+‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+>
+>
+>
+</code>