Anova de 2 via Medidas Repetidas

library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.5     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.0.2     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggpubr)
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
data("selfesteem2",package="datarium")
selfesteem2
## # A tibble: 24 x 5
##    id    treatment    t1    t2    t3
##    <fct> <fct>     <dbl> <dbl> <dbl>
##  1 1     ctr          83    77    69
##  2 2     ctr          97    95    88
##  3 3     ctr          93    92    89
##  4 4     ctr          92    92    89
##  5 5     ctr          77    73    68
##  6 6     ctr          72    65    63
##  7 7     ctr          92    89    79
##  8 8     ctr          92    87    81
##  9 9     ctr          95    91    84
## 10 10    ctr          92    84    81
## # ... with 14 more rows
selfesteem2<- selfesteem2 %>% gather(key="time",value="score",t1,t2,t3) %>% convert_as_factor(id,time)
selfesteem2
## # A tibble: 72 x 4
##    id    treatment time  score
##    <fct> <fct>     <fct> <dbl>
##  1 1     ctr       t1       83
##  2 2     ctr       t1       97
##  3 3     ctr       t1       93
##  4 4     ctr       t1       92
##  5 5     ctr       t1       77
##  6 6     ctr       t1       72
##  7 7     ctr       t1       92
##  8 8     ctr       t1       92
##  9 9     ctr       t1       95
## 10 10    ctr       t1       92
## # ... with 62 more rows

Resumen Estadístico

selfesteem2 %>% group_by(treatment,time) %>% get_summary_stats(score,type="mean_sd")
## # A tibble: 6 x 6
##   treatment time  variable     n  mean    sd
##   <fct>     <fct> <chr>    <dbl> <dbl> <dbl>
## 1 ctr       t1    score       12  88    8.08
## 2 ctr       t2    score       12  83.8 10.2 
## 3 ctr       t3    score       12  78.7 10.5 
## 4 Diet      t1    score       12  87.6  7.62
## 5 Diet      t2    score       12  87.8  7.42
## 6 Diet      t3    score       12  87.7  8.14

Visualización de los datos

bxp<-ggboxplot(selfesteem2,x="time",y="score",color="treatment",palette="jco")
bxp

## Supuestos del Modelo

Outliers

selfesteem2 %>% group_by(treatment,time) %>% identify_outliers(score)
## [1] treatment  time       id         score      is.outlier is.extreme
## <0 rows> (or 0-length row.names)

Prueba de Normalidad

selfesteem2 %>% group_by(treatment,time) %>% shapiro_test(score)
## # A tibble: 6 x 5
##   treatment time  variable statistic      p
##   <fct>     <fct> <chr>        <dbl>  <dbl>
## 1 ctr       t1    score        0.828 0.0200
## 2 ctr       t2    score        0.868 0.0618
## 3 ctr       t3    score        0.887 0.107 
## 4 Diet      t1    score        0.919 0.279 
## 5 Diet      t2    score        0.923 0.316 
## 6 Diet      t3    score        0.886 0.104
ggqqplot(selfesteem2,"score",ggtheme=theme_bw()) + facet_grid(time~treatment,labeller="label_both")

Cálculo del ANOVA

res.aov<-anova_test(data=selfesteem2,dv=score,wid=id,within=c(treatment,time))
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##           Effect  DFn   DFd      F        p p<.05   ges
## 1      treatment 1.00 11.00 15.541 2.00e-03     * 0.059
## 2           time 1.31 14.37 27.369 5.03e-05     * 0.049
## 3 treatment:time 2.00 22.00 30.424 4.63e-07     * 0.050

La intteracción entre el factor treatment y time es estadísticamente significativa F(2,22)=30.4, p<0.0001

Pruebas Post-Hoc

  1. Si la interacción es significativa
  1. Efectos principales simples: correr un anova de una via con la primera variable
  2. Comparación pareada simple: correr multiples comparaciones pareadas para determinar cuales grupos son diferentes.
  1. Si la interacción no es significativa

Se determina cuales de los efectos principales son significativos. Ver resultados del Anova

one.way<- selfesteem2 %>% group_by(time) %>% anova_test(dv=score,wid=id,within=treatment) %>% get_anova_table() %>% adjust_pvalue(method="bonferroni")
one.way
## # A tibble: 3 x 9
##   time  Effect      DFn   DFd      F       p `p<.05`      ges   p.adj
## * <fct> <chr>     <dbl> <dbl>  <dbl>   <dbl> <chr>      <dbl>   <dbl>
## 1 t1    treatment     1    11  0.376 0.552   ""      0.000767 1      
## 2 t2    treatment     1    11  9.03  0.012   "*"     0.052    0.036  
## 3 t3    treatment     1    11 30.9   0.00017 "*"     0.199    0.00051
pwc<- selfesteem2 %>% 
      group_by(time) %>% 
      pairwise_t_test(
     score~treatment,paired=TRUE,
     p.adjust.method="bonferroni"
)
pwc
## # A tibble: 3 x 11
##   time  .y.   group1 group2    n1    n2 statistic    df       p   p.adj
## * <fct> <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>   <dbl>   <dbl>
## 1 t1    score ctr    Diet      12    12     0.613    11 0.552   0.552  
## 2 t2    score ctr    Diet      12    12    -3.00     11 0.012   0.012  
## 3 t3    score ctr    Diet      12    12    -5.56     11 0.00017 0.00017
## # ... with 1 more variable: p.adj.signif <chr>
pwc<-pwc %>% add_xy_position(x="time")
bxp+
  stat_pvalue_manual(pwc,tip.length=0,hide.ns=TRUE)+
  labs(
     subtitle=get_test_label(res.aov,detailed=TRUE),
     caption=get_pwc_label(pwc)
  )