###################################################
###################################################
###################################################
# data
df <- data.frame(patient=rep(1:5, each=4),
drug=rep(1:4, times=5),
response=c(30, 28, 16, 34,
14, 18, 10, 22,
24, 20, 18, 30,
38, 34, 20, 44,
26, 28, 14, 30))
#view data
df
write.csv(df, file="rep.meas.anova.csv")
#fit repeated measures ANOVA model
df$drug <- factor(df$drug)
df$patient <- factor(df$patient)
# Error(patient) = patient error should be isolated
m.aov <- aov(response ~ drug
+ Error(patient),
data = df)
#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)
1 - pf(24.75886525, 3, 12)
# or
library(broom)
fit <- tidy(m2.aov)
fit
fit$statistic[2]
1 - 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
# 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).
> # data
> df <- data.frame(patient=rep(1:5, each=4),
+ drug=rep(1:4, times=5),
+ response=c(30, 28, 16, 34,
+ 14, 18, 10, 22,
+ 24, 20, 18, 30,
+ 38, 34, 20, 44,
+ 26, 28, 14, 30))
>
> #view data
> df
patient drug response
1 1 1 30
2 1 2 28
3 1 3 16
4 1 4 34
5 2 1 14
6 2 2 18
7 2 3 10
8 2 4 22
9 3 1 24
10 3 2 20
11 3 3 18
12 3 4 30
13 4 1 38
14 4 2 34
15 4 3 20
16 4 4 44
17 5 1 26
18 5 2 28
19 5 3 14
20 5 4 30
> write.csv(df, file="rep.meas.anova.csv")
>
>
> #fit repeated measures ANOVA model
> df$drug <- factor(df$drug)
> df$patient <- factor(df$patient)
>
> # Error(patient) = patient error should be isolated
> m.aov <- aov(response ~ drug
+ + Error(patient),
+ data = df)
> #view model summary
> summary(m.aov)
Error: patient
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 4 680.8 170.2
Error: Within
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 ‘***’ 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 ‘***’ 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 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> # check probability level (pr)
> 1 - pf(24.75886525, 3, 12)
[1] 1.992501e-05
> # or
> library(broom)
> fit <- tidy(m2.aov)
> fit
# A tibble: 3 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 patient Resi… 4 681. 170. NA NA
2 patient:drug drug 3 698. 233. 24.8 1.99e-5
3 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
>
>
rep.meas.anova.eg.movie.review.csv
# 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
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,
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,
76, 78, 90,
78, 94, 95,
80, 83, 88,
82, 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,
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,
+ 76, 78, 90,
+ 78, 94, 95,
+ 80, 83, 88,
+ 82, 90, 99))
>
> #view data
> movrev
reviewer movie score
1 1 1 88
2 1 2 84
3 1 3 92
4 2 1 76
5 2 2 78
6 2 3 90
7 3 1 78
8 3 2 94
9 3 3 95
10 4 1 80
11 4 2 83
12 4 3 88
13 5 1 82
14 5 2 90
15 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 ‘***’ 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 ‘***’ 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
1 2
2 0.628 -
3 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
1 2
2 0.628 -
3 0.029 0.060
P value adjustment method: bonferroni
>
>
#### # 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)
> ####
> # by hand
> ####
> movrev
reviewer movie score
1 1 1 88
2 1 2 84
3 1 3 92
4 2 1 76
5 2 2 78
6 2 3 90
7 3 1 78
8 3 2 94
9 3 3 95
10 4 1 80
11 4 2 83
12 4 3 88
13 5 1 82
14 5 2 90
15 5 3 99
> tapply(movrev$score, list(reviewer, movie), mean) # 각 셀의 평균값
1 2 3
1 88 84 92
2 76 78 90
3 78 94 95
4 80 83 88
5 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 ‘***’ 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 ‘***’ 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 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>