Introducción

Considere la siguiente pregunta: ¿qué diferencias se genera en la actividad microbiana entre los tiempos de medición ?

Generación del diseño

El diseño en Medidas Repetidas considera que existe un factor dominado por el tiempo. Observandosé el cambio a través de este tiempo, de la misma unidad experimental.

el modelo a considera es

\(Y_{ij}= \mu + \tau_i * \alpha_j +s_i+ \epsilon_{ij}\)

Ejemplo

Librerías

library(tidyverse)    # Manipulación de 📈
library(car)          # ANOVA tipo III y pruebas
library(emmeans)      # Comparaciones post-hoc
library(ggplot2)      # Visualización
library(rstatix)      # Herramientas estadísticas adicionales
library(multcomp)
library(Analitica)

Datos Simulados para ejemplificación

##   id_sujeto tiempo bacterias
## 1         1      0        94
## 2         2      0        98
## 3         3      0       116
## 4         4      0       101
## 5         5      0       101
## 6         6      0       117
## 'data.frame':    80 obs. of  3 variables:
##  $ id_sujeto: Factor w/ 20 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ tiempo   : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ bacterias: num  94 98 116 101 101 117 105 87 93 96 ...
##  - attr(*, "out.attrs")=List of 2
##   ..$ dim     : Named int [1:2] 20 4
##   .. ..- attr(*, "names")= chr [1:2] "id_sujeto" "tiempo"
##   ..$ dimnames:List of 2
##   .. ..$ id_sujeto: chr [1:20] "id_sujeto=1" "id_sujeto=2" "id_sujeto=3" "id_sujeto=4" ...
##   .. ..$ tiempo   : chr [1:4] "tiempo= 0" "tiempo= 2" "tiempo= 4" "tiempo=14"

El proceso de análisis de medidas repetidas con R, el requisito es que las variables necesarias sean declaradas como factor.

Revisamos y corregimos o modificamos lo que necesitamos, toda variable que deba ser factor se corrige, a pesar de ser “CHAR”

# Ajustar tipos de variables
datos <- datos %>%
  mutate(
    tiempo = factor(tiempo),
    id_sujeto = factor(id_sujeto)
  )

str(datos)
## 'data.frame':    80 obs. of  3 variables:
##  $ id_sujeto: Factor w/ 20 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ tiempo   : Factor w/ 4 levels "0","2","4","14": 1 1 1 1 1 1 1 1 1 1 ...
##  $ bacterias: num  94 98 116 101 101 117 105 87 93 96 ...
##  - attr(*, "out.attrs")=List of 2
##   ..$ dim     : Named int [1:2] 20 4
##   .. ..- attr(*, "names")= chr [1:2] "id_sujeto" "tiempo"
##   ..$ dimnames:List of 2
##   .. ..$ id_sujeto: chr [1:20] "id_sujeto=1" "id_sujeto=2" "id_sujeto=3" "id_sujeto=4" ...
##   .. ..$ tiempo   : chr [1:4] "tiempo= 0" "tiempo= 2" "tiempo= 4" "tiempo=14"

identificando posibles outliers

descripYG(datos,bacterias,tiempo)
## Picking joint bandwidth of 4.03

##   Group  n   Mean Median       SD Kurtosis    Skewness         CV Min Max
## 1     0 20 101.40  101.0 9.848323 2.714991 -0.09444535 0.09712350  80 118
## 2     2 20 107.55  106.5 8.319128 2.027274 -0.20389646 0.07735126  91 121
## 3     4 20 117.15  116.0 9.598657 2.652036  0.35030963 0.08193476 101 138
## 4    14 20 154.85  154.0 9.740177 3.349782  0.06931822 0.06290072 133 177
##      P25    P75   IQR
## 1  94.75 105.50 10.75
## 2 101.75 115.25 13.50
## 3 111.75 122.50 10.75
## 4 149.00 160.00 11.00
datos %>%
  group_by(tiempo) %>%
  identify_outliers(bacterias)
## # A tibble: 1 × 5
##   tiempo id_sujeto bacterias is.outlier is.extreme
##   <fct>  <fct>         <dbl> <lgl>      <lgl>     
## 1 14     10              177 TRUE       FALSE

la pregunta que debe hacerse es si corrigo mis datos ya sea eliminado o corrigiendo, en ambos casos debe realizar verificación de datos.

Obtención del Modelo

Establecimiento Modelo para análisis de supuesto

modelo_aov, considera el efecto del tratamiento a través del tiempo

#Anova
modelo_aov <- aov(
  bacterias ~ tiempo + Error(id_sujeto / tiempo),
  data = datos
)

summary(modelo_aov) # revisar si esta correctamente realizado el modelo, fijarse en los grados de libertad
## 
## Error: id_sujeto
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 19   1723    90.7               
## 
## Error: id_sujeto:tiempo
##           Df Sum Sq Mean Sq F value Pr(>F)    
## tiempo     3  34468   11489   131.3 <2e-16 ***
## Residuals 57   4988      88                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los Df deben equivaler al numero de elementos del factor menos 1 ( factor tiempo son 4 por tanto Df =3)

Análisis de Supuestos

NORMALIDAD

normalidad <- datos %>%
  group_by(tiempo) %>%
  summarise(
    W = shapiro.test(bacterias)$statistic,
    p_valor = shapiro.test(bacterias)$p.value,
    .groups = "drop"
  )

print(normalidad)
## # A tibble: 4 × 3
##   tiempo     W p_valor
##   <fct>  <dbl>   <dbl>
## 1 0      0.969   0.742
## 2 2      0.959   0.517
## 3 4      0.973   0.823
## 4 14     0.982   0.954

para cada uno de los grupos analizados el p-value resulta ser mayor al alfa (5%), por ende se puede continuar con el ANOVA (en caso contrario aplicar test FRIEDMAN)

Esfericidad

es uno de los supuestos más importantes del ANOVA de medidas repetidas. Su objetivo es verificar si la variabilidad de las diferencias entre todos los pares de mediciones repetidas es similar.

Lo que se busca es que la varianza entre cada transición sea similar.

Cuando la esfericidad no se cumple, el ANOVA tiende a producir valores F demasiado grandes, lo que aumenta la probabilidad de encontrar diferencias significativas cuando en realidad no las hay (incremento del error Tipo I).

esfericidad <- datos %>%
  anova_test(
    dv = bacterias,
    wid = id_sujeto,
    within = tiempo
  )
print(esfericidad$`Mauchly's Test for Sphericity`)
##   Effect     W     p p<.05
## 1 tiempo 0.855 0.734

La hipótesis nula de Mauchly establece que se cumple la esfericidad. Al revisar vemos que p = 0.734 > 0.05, estableciendo que no se rechaza la hipótesis nula. Por tanto, el supuesto de esfericidad se cumple con DF(3,57)

Si se rechazará habría que hacer la correción

print(esfericidad$`Sphericity Corrections`)
##   Effect   GGe      DF[GG]    p[GG] p[GG]<.05   HFe   DF[HF]    p[HF] p[HF]<.05
## 1 tiempo 0.917 2.75, 52.27 1.23e-23         * 1.088 3.26, 62 1.44e-25         *

Esto significa que puedes analizar tu ANOVA de medidas repetidas sin aplicar correcciones (los grados de libertad originales son válidos).

Greenhouse–Geisser (GGe) = 0.917 y Huynh–Feldt (HFe) = 1.088 en caso contrario los nuevos DF (GGe3, GGe57), igualmente para la correccion Huynh-Feldt DF(HFe3.HFe57)

Analisis del Modelo

summary(modelo_aov)
## 
## Error: id_sujeto
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 19   1723    90.7               
## 
## Error: id_sujeto:tiempo
##           Df Sum Sq Mean Sq F value Pr(>F)    
## tiempo     3  34468   11489   131.3 <2e-16 ***
## Residuals 57   4988      88                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

se observa que hay diferencia en el factor tiempo, con un p-valor de 0. Luego realizamos post hoc

Análisis Post Hoc

Una recomendación para saber si se realiza bien los procedimientos es mirar que eventualmente debería encontrar como resultados. Para ello aplique un análisis exploratorio como el siguiente

### Revisar comportamiento esperable

# Resumen por tratamiento y tiempo
emm_tiempo <- emmeans(modelo_aov, ~ tiempo)
## Note: re-fitting model with sum-to-zero contrasts
pairs(emm_tiempo, adjust = "holm")
##  contrast           estimate   SE df t.ratio p.value
##  tiempo0 - tiempo2     -6.15 2.96 57  -2.079  0.0421
##  tiempo0 - tiempo4    -15.75 2.96 57  -5.324 <0.0001
##  tiempo0 - tiempo14   -53.45 2.96 57 -18.069 <0.0001
##  tiempo2 - tiempo4     -9.60 2.96 57  -3.245  0.0039
##  tiempo2 - tiempo14   -47.30 2.96 57 -15.990 <0.0001
##  tiempo4 - tiempo14   -37.70 2.96 57 -12.745 <0.0001
## 
## P value adjustment: holm method for 6 tests
pairs(emm_tiempo, adjust = "hochberg")
##  contrast           estimate   SE df t.ratio p.value
##  tiempo0 - tiempo2     -6.15 2.96 57  -2.079  0.0421
##  tiempo0 - tiempo4    -15.75 2.96 57  -5.324 <0.0001
##  tiempo0 - tiempo14   -53.45 2.96 57 -18.069 <0.0001
##  tiempo2 - tiempo4     -9.60 2.96 57  -3.245  0.0039
##  tiempo2 - tiempo14   -47.30 2.96 57 -15.990 <0.0001
##  tiempo4 - tiempo14   -37.70 2.96 57 -12.745 <0.0001
## 
## P value adjustment: hochberg method for 6 tests

Contrastes contra el tiempo

Los contrastes son comparaciones planificadas entre medias, estos permiten responder preguntas específicas planteadas antes del análisis.

En un ANOVA de medidas repetidas, primero se responde la pregunta general:

¿qué efecto se genera entre los tiempos de medición?

Si el ANOVA es significativo, los contrastes permiten responder preguntas particulares, por ejemplo:

¿El tiempo inicial difiere del tiempo 2? ¿El tiempo inicial difiere del tiempo 4? ¿El tiempo inicial difiere del tiempo 14?

contrast(
  emm_tiempo,
  method = list(
    "t0 vs t2"  = c(1, -1, 0, 0),
    "t0 vs t4"  = c(1, 0, -1, 0),
    "t0 vs t14" = c(1, 0, 0, -1)
  ),
  adjust = "holm"
)
##  contrast  estimate   SE df t.ratio p.value
##  t0 vs t2     -6.15 2.96 57  -2.079  0.0421
##  t0 vs t4    -15.75 2.96 57  -5.324 <0.0001
##  t0 vs t14   -53.45 2.96 57 -18.069 <0.0001
## 
## P value adjustment: holm method for 3 tests
plot(pairs(emm_tiempo, adjust = "holm"))

cld(emm_tiempo, Letters = letters, adjust = "holm")
##  tiempo emmean  SE df lower.CL upper.CL .group
##  0         101 2.1 76       96      107  a    
##  2         108 2.1 76      102      113   b   
##  4         117 2.1 76      112      123    c  
##  14        155 2.1 76      149      160     d 
## 
## Warning: EMMs are biased unless design is perfectly balanced 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 4 estimates 
## P value adjustment: holm method for 6 tests 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

Alternativa no paramétrica: Friedman

# ------------------------------------------------------------
# 
# Usar si no se cumplen los supuestos
# ------------------------------------------------------------

friedman.test(
  bacterias ~ tiempo | id_sujeto,
  data = datos
)
## 
##  Friedman rank sum test
## 
## data:  bacterias and tiempo and id_sujeto
## Friedman chi-squared = 45.091, df = 3, p-value = 8.85e-10
datos %>%
  pairwise_wilcox_test(
    bacterias ~ tiempo,
    paired = TRUE,
    p.adjust.method = "holm"
  )
## # A tibble: 6 × 9
##   .y.       group1 group2    n1    n2 statistic          p    p.adj p.adj.signif
## * <chr>     <chr>  <chr>  <int> <int>     <dbl>      <dbl>    <dbl> <chr>       
## 1 bacterias 0      2         20    20      45.5 0.0286      2.86e-2 *           
## 2 bacterias 0      4         20    20      15.5 0.000275    8.24e-4 ***         
## 3 bacterias 0      14        20    20       0   0.00000191  1.14e-5 ****        
## 4 bacterias 2      4         20    20      27   0.00223     4.46e-3 **          
## 5 bacterias 2      14        20    20       0   0.00000191  1.14e-5 ****        
## 6 bacterias 4      14        20    20       0   0.00000191  1.14e-5 ****

Se observa que el tratamiento 2 presenta una misma tendencia que el grupo control y que los tratamiento 1 y 3 presentaría similitud en sus resultados. Así mismo se observa que el último periodo de medición podría presentar una diferencia significativa con respecto al estado inicial al menos, en algunos tratamientos.