Factor intrasujetos tiempo Factores entresujetos.

Hipótesis.

\[H_0= \mu_t1= \mu_t2=\mu_t3 \]

library(datarium)
## Warning: package 'datarium' was built under R version 4.2.3
data("selfesteem", package = "datarium")

datos = selfesteem
head(datos, 3)
##   id       t1       t2       t3
## 1  1 4.005027 5.182286 7.107831
## 2  2 2.558124 6.912915 6.308434
## 3  3 3.244241 4.443434 9.778410
boxplot(datos)

library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
data("selfesteem", package = "datarium")

datos = selfesteem
datos = datos  %>%
  gather(key = "tiempo",
         value = "rto",
         t1, t2, t3) %>%
  mutate_at(vars(id, tiempo), as.factor)
datos
## # A tibble: 30 × 3
##    id    tiempo   rto
##    <fct> <fct>  <dbl>
##  1 1     t1      4.01
##  2 2     t1      2.56
##  3 3     t1      3.24
##  4 4     t1      3.42
##  5 5     t1      2.87
##  6 6     t1      2.05
##  7 7     t1      3.53
##  8 8     t1      3.18
##  9 9     t1      3.51
## 10 10    t1      3.04
## # ℹ 20 more rows
boxplot(datos)

datos %>%
  group_by(tiempo) %>%
  summarise(media = mean(rto),
            desv = sd(rto),
            n = n(),
            cv = 100*desv/media)
## # A tibble: 3 × 5
##   tiempo media  desv     n    cv
##   <fct>  <dbl> <dbl> <int> <dbl>
## 1 t1      3.14 0.552    10  17.6
## 2 t2      4.93 0.863    10  17.5
## 3 t3      7.64 1.14     10  15.0
boxplot(datos$rto~datos$tiempo)

library(rstatix)
## Warning: package 'rstatix' was built under R version 4.2.3
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
datos %>%
  group_by(tiempo) %>%
  identify_outliers(rto)
## # A tibble: 2 × 5
##   tiempo id      rto is.outlier is.extreme
##   <fct>  <fct> <dbl> <lgl>      <lgl>     
## 1 t1     6      2.05 TRUE       FALSE     
## 2 t2     2      6.91 TRUE       FALSE

probando normalidad

datos %>%
  group_by(tiempo) %>%
  shapiro_test(rto)
## # A tibble: 3 × 4
##   tiempo variable statistic     p
##   <fct>  <chr>        <dbl> <dbl>
## 1 t1     rto          0.967 0.859
## 2 t2     rto          0.876 0.117
## 3 t3     rto          0.923 0.380

los datos aparentan tener normalidad, puesto que los valores son >5%, lo que indica una presencia de datos atipicos casi nula

supuesto de esfericidad.

res.aov <- anova_test(data = datos,
                      dv = rto,
                      wid = id,
                      within = tiempo)


res.aov$`Mauchly's Test for Sphericity`
##   Effect     W     p p<.05
## 1 tiempo 0.551 0.092
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##   Effect DFn DFd      F        p p<.05   ges
## 1 tiempo   2  18 55.469 2.01e-08     * 0.829

se compara los diferentes tiempos sucesivos entre si, para poder ver su variabilidad

si es >5% se puede asumir que las varianzas son igules.

# Se usa el p-valor ajustado (p.adj)

datos %>%
  pairwise_t_test(
    rto ~ tiempo, paired = TRUE,
    p.adjust.method = "bonferroni")
## # A tibble: 3 × 10
##   .y.   group1 group2    n1    n2 statistic    df           p p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>       <dbl> <dbl> <chr>       
## 1 rto   t1     t2        10    10     -4.97     9 0.000772     2e-3 **          
## 2 rto   t1     t3        10    10    -13.2      9 0.000000334  1e-6 ****        
## 3 rto   t2     t3        10    10     -4.87     9 0.000886     3e-3 **

Medidas repetidas de dos vías.

data("selfesteem2", package = "datarium")

datos2 = selfesteem2
datos2
## # A tibble: 24 × 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
## # ℹ 14 more rows
datos2 = selfesteem2
datos2$treatment = gl(2,12,24, c('con fert', 'sin fert'))
datos2 = datos2 %>% 
  gather(key='tiempo', value = 'rto',
         t1,t2,t3)
datos2 %>%
  group_by(treatment, tiempo) %>%
  summarise(media = mean(rto),
            desv = sd(rto),
            n = n(),
            cv = 100*desv/media)
## `summarise()` has grouped output by 'treatment'. You can override using the
## `.groups` argument.
## # A tibble: 6 × 6
## # Groups:   treatment [2]
##   treatment tiempo media  desv     n    cv
##   <fct>     <chr>  <dbl> <dbl> <int> <dbl>
## 1 con fert  t1      88    8.08    12  9.18
## 2 con fert  t2      83.8 10.2     12 12.2 
## 3 con fert  t3      78.7 10.5     12 13.4 
## 4 sin fert  t1      87.6  7.62    12  8.70
## 5 sin fert  t2      87.8  7.42    12  8.45
## 6 sin fert  t3      87.7  8.14    12  9.28

El mejor resultado fue el tiempo 1 con fertilizacion.

visualización.

library(ggplot2)

ggplot(datos2)+
  aes(tiempo, rto, fill=treatment)+
  geom_boxplot()

Pruebas de datos atípicos.

datos2 %>%
  group_by(treatment, tiempo) %>%
  identify_outliers(rto)
## [1] treatment  tiempo     id         rto        is.outlier is.extreme
## <0 rows> (or 0-length row.names)
datos2 %>%
  group_by(treatment, tiempo) %>%
  shapiro_test(rto)
## # A tibble: 6 × 5
##   treatment tiempo variable statistic      p
##   <fct>     <chr>  <chr>        <dbl>  <dbl>
## 1 con fert  t1     rto          0.828 0.0200
## 2 con fert  t2     rto          0.868 0.0618
## 3 con fert  t3     rto          0.887 0.107 
## 4 sin fert  t1     rto          0.919 0.279 
## 5 sin fert  t2     rto          0.923 0.316 
## 6 sin fert  t3     rto          0.886 0.104

Revisión p-valor.

res.aov <- anova_test(
  data = datos2,
  dv = rto,
  wid = id,
  within = c(treatment,
             tiempo)
  )
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           tiempo 1.31 14.37 27.369 5.03e-05     * 0.049
## 3 treatment:tiempo 2.00 22.00 30.424 4.63e-07     * 0.050

Si hay interacción

datos2 %>%
  group_by(treatment, tiempo) %>%
  identify_outliers(rto)
## [1] treatment  tiempo     id         rto        is.outlier is.extreme
## <0 rows> (or 0-length row.names)
datos2 %>%
  group_by(treatment, tiempo) %>%
  shapiro_test(rto)
## # A tibble: 6 × 5
##   treatment tiempo variable statistic      p
##   <fct>     <chr>  <chr>        <dbl>  <dbl>
## 1 con fert  t1     rto          0.828 0.0200
## 2 con fert  t2     rto          0.868 0.0618
## 3 con fert  t3     rto          0.887 0.107 
## 4 sin fert  t1     rto          0.919 0.279 
## 5 sin fert  t2     rto          0.923 0.316 
## 6 sin fert  t3     rto          0.886 0.104

Gráfico de interacción.

interaction.plot(datos2$tiempo,
                 datos2$treatment,
                 datos2$rto)

datos2 %>% 
  group_by(tiempo, treatment) %>% 
  summarise(mean_rto = mean(rto)) %>% 
  ggplot()+
  aes(tiempo, mean_rto,
      color=treatment,
      group=treatment)+
  geom_point(size=5)+
  geom_line(linewidth=3)
## `summarise()` has grouped output by 'tiempo'. You can override using the
## `.groups` argument.
## Warning: Ignoring unknown parameters: linewidth

> En el grafico podemos ver una interacción que es evidente en el tiempo 3