COMMunication
RESearch.NET

User Tools

Site Tools

You are here: index » r » repeated_measure_anova
Trace: • anova • schedule • probability • visualization • schedule • social_network_analysis • schedule • schedule
r:repeated_measure_anova

Differences

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

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
r:repeated_measure_anova [2020/06/03 23:49] – removed hkimscilr:repeated_measure_anova [2025/05/12 08:25] (current) – [E.g.2 by hand] hkimscil
Line 1: Line 1:
 +====== E.g. 1 ======
 +{{r:rep.meas.anova.csv}}
 +<code>
 +###################################################
 +###################################################
 +###################################################
  
 +# 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).
 +</code>
 +
 +<code>
 +> # 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
 +> 
 +> 
 +</code>
 +====== E.g. 2 ======
 +{{:r:rep.meas.anova.eg.movie.review.csv}}
 +<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
 +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"))
 +
 +</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     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 
 +> 
 +> 
 +</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     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
 +> 
 +>
 +> 
 +</code>
r/repeated_measure_anova.1591195798.txt.gz · Last modified: 2020/06/03 23:49 by hkimscil
Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki