###################################################
###################################################
###################################################

# 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)

# check probability level (pr)
1 - pf(24.75886525, 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).

> # Data preparation
> # Wide format
> # install.packages("datarium")
> # library(datarium)
> 
> data("selfesteem", package = "datarium")
> head(selfesteem, 3)
# A tibble: 3 × 4
     id    t1    t2    t3
  <int> <dbl> <dbl> <dbl>
1     1  4.01  5.18  7.11
2     2  2.56  6.91  6.31
3     3  3.24  4.44  9.78
> 
> # Gather columns t1, t2 and t3 into long format
> # Convert id and time into factor variables
> 
> # for %>% function, dplyr packages needed
> # install.packages("dplyr")
> # library(dplyr)
> 
> # for convert_as_factor function, rstatix needed
> # install.packages("rstatix")
> # library(rstatix)
> 
> selfesteem <- selfesteem %>% 
+   gather(key = "time", value = "score", t1, t2, t3) %>% 
+   convert_as_factor(id, time)
> head(selfesteem, 3)
# A tibble: 3 × 3
  id    time  score
  <fct> <fct> <dbl>
1 1     t1     4.01
2 2     t1     2.56
3 3     t1     3.24
> 
> # get_summary_stats(group_by(selfesteem, time), 
> # score, type = "mean_sd")
> # the above is the same as the below
> selfesteem %>%
+   group_by(time) %>%
+   get_summary_stats(score, type = "mean_sd")
# A tibble: 3 × 5
  time  variable     n  mean    sd
  <fct> <fct>    <dbl> <dbl> <dbl>
1 t1    score       10  3.14 0.552
2 t2    score       10  4.93 0.863
3 t3    score       10  7.64 1.14 
> 
> 
> res.aov <- anova_test(
+   data = selfesteem, 
+   dv = score, 
+   wid = id, 
+   within = time)
> get_anova_table(res.aov)
ANOVA Table (type III tests)

  Effect DFn DFd      F        p p<.05   ges
1   time   2  18 55.469 2.01e-08     * 0.829
> 
> # ges = generalized effect size
> # F
> # (2,18)
> 
> ###################################################
> ###################################################
> ###################################################
> 
> # 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
> 
> # check probability level (pr)
> 1 - pf(24.75886525, 3, 12)
[1] 1.992501e-05
> 
> # 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).
> 
> 

# 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.mov.rev.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)
#view model summary
summary(m.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")
# or
with(movrev, 
     pairwise.t.test(score, movie, 
                     paired = T, 
                     p.adjust.method = "bonf"))