Diseño Medidas Repetidas con Tratamiento

Carlos Jiménez-Gallardo
Estadístico
MSc Infórmatica Educativa
Universidad de La Frontera
carlos.jimenez@ufrontera.cl

Data Scientist
www.innovate.cl cjimenez@innovate.cl

Introducción

Considere la siguiente pregunta: ¿Cuanto tiempo debe transcurrir para que algún tratamiento provoque algún efecto en en el disminución de actividad bacteriana?

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 datos
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

##    tratamiento tiempo replica bacterias id_sujeto
## 1           gc      0       1        94      gc_1
## 2           t1      0       1        98      t1_1
## 3           t2      0       1       116      t2_1
## 4           t3      0       1       101      t3_1
## 5           gc      2       1       105      gc_1
## 6           t1      2       1       127      t1_1
## 7           t2      2       1       111      t2_1
## 8           t3      2       1        95      t3_1
## 9           gc      4       1       101      gc_1
## 10          t1      4       1       116      t1_1
## 11          t2      4       1       124      t2_1
## 12          t3      4       1       120      t3_1

## 'data.frame':    80 obs. of  5 variables:
##  $ tratamiento: Factor w/ 4 levels "gc","t1","t2",..: 1 2 3 4 1 2 3 4 1 2 ...
##  $ tiempo     : num  0 0 0 0 2 2 2 2 4 4 ...
##  $ replica    : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ bacterias  : num  94 98 116 101 105 127 111 95 101 116 ...
##  $ id_sujeto  : chr  "gc_1" "t1_1" "t2_1" "t3_1" ...
##  - attr(*, "out.attrs")=List of 2
##   ..$ dim     : Named int [1:3] 4 4 5
##   .. ..- attr(*, "names")= chr [1:3] "tratamiento" "tiempo" "replica"
##   ..$ dimnames:List of 3
##   .. ..$ tratamiento: chr [1:4] "tratamiento=gc" "tratamiento=t1" "tratamiento=t2" "tratamiento=t3"
##   .. ..$ tiempo     : chr [1:4] "tiempo= 0" "tiempo= 2" "tiempo= 4" "tiempo=14"
##   .. ..$ replica    : chr [1:5] "replica=1" "replica=2" "replica=3" "replica=4" ...

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  5 variables:
##  $ tratamiento: Factor w/ 4 levels "gc","t1","t2",..: 1 2 3 4 1 2 3 4 1 2 ...
##  $ tiempo     : Factor w/ 4 levels "0","2","4","14": 1 1 1 1 2 2 2 2 3 3 ...
##  $ replica    : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ bacterias  : num  94 98 116 101 105 127 111 95 101 116 ...
##  $ id_sujeto  : Factor w/ 20 levels "gc_1","gc_2",..: 1 6 11 16 1 6 11 16 1 6 ...
##  - attr(*, "out.attrs")=List of 2
##   ..$ dim     : Named int [1:3] 4 4 5
##   .. ..- attr(*, "names")= chr [1:3] "tratamiento" "tiempo" "replica"
##   ..$ dimnames:List of 3
##   .. ..$ tratamiento: chr [1:4] "tratamiento=gc" "tratamiento=t1" "tratamiento=t2" "tratamiento=t3"
##   .. ..$ tiempo     : chr [1:4] "tiempo= 0" "tiempo= 2" "tiempo= 4" "tiempo=14"
##   .. ..$ replica    : chr [1:5] "replica=1" "replica=2" "replica=3" "replica=4" ...

identificando posibles outliers

descripYG(datos,bacterias,tratamiento)

## Picking joint bandwidth of 9.42

##   Group  n   Mean Median       SD Kurtosis  Skewness        CV Min Max    P25
## 1    gc 20 109.55  106.5 14.14762 2.538499 0.6113762 0.1291431  89 140 101.00
## 2    t1 20 125.00  117.0 27.94732 2.422913 0.6588415 0.2235786  80 183 106.75
## 3    t2 20 115.85  109.5 16.24897 2.030230 0.7015565 0.1402587  96 146 103.75
## 4    t3 20 120.55  118.0 24.73113 2.255070 0.5766381 0.2051525  85 174 101.00
##      P75   IQR
## 1 117.25 16.25
## 2 137.75 31.00
## 3 127.00 23.25
## 4 140.00 39.00

descripYG(datos,bacterias,tiempo)

## Picking joint bandwidth of 5.07

##   Group  n   Mean Median        SD Kurtosis    Skewness         CV Min Max
## 1     0 20 101.80  103.0  8.036758 4.108857 -0.92744605 0.07894654  80 116
## 2     2 20 107.35  104.5 11.806666 2.668148  0.37177253 0.10998292  85 131
## 3     4 20 113.60  117.0 11.361986 2.428278  0.04051978 0.10001749  93 138
## 4    14 20 148.20  146.0 16.916513 2.506638  0.25538921 0.11414651 117 183
##      P25    P75   IQR
## 1  98.75 107.25  8.50
## 2 101.00 111.50 10.50
## 3 102.75 120.00 17.25
## 4 137.50 158.25 20.75

datos %>%
  group_by(tiempo) %>%
  identify_outliers(bacterias)

## # A tibble: 3 × 7
##   tiempo tratamiento replica bacterias id_sujeto is.outlier is.extreme
##   <fct>  <fct>       <fct>       <dbl> <fct>     <lgl>      <lgl>     
## 1 0      t1          2              80 t1_2      TRUE       FALSE     
## 2 2      t1          5             131 t1_5      TRUE       FALSE     
## 3 2      t3          5              85 t3_5      TRUE       FALSE

Obtención del Modelo

Establecimiento Modelo para analisis de supuesto

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

#Anova
modelo_aov <- aov(bacterias ~ tratamiento * 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)    
## tratamiento  3 2625.0   875.0   14.46 8.09e-05 ***
## Residuals   16  968.2    60.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: id_sujeto:tiempo
##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## tiempo              3  26140    8713  100.40  < 2e-16 ***
## tratamiento:tiempo  9   4007     445    5.13 7.21e-05 ***
## Residuals          48   4166      87                     
## ---
## 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

cat("\n--- Prueba de Normalidad (Shapiro-Wilk por grupo) ---\n")

## 
## --- Prueba de Normalidad (Shapiro-Wilk por grupo) ---

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

## # A tibble: 16 × 4
##    tratamiento tiempo     W p_valor
##    <fct>       <fct>  <dbl>   <dbl>
##  1 gc          0      0.863  0.240 
##  2 gc          2      0.929  0.587 
##  3 gc          4      0.851  0.197 
##  4 gc          14     0.964  0.839 
##  5 t1          0      0.891  0.362 
##  6 t1          2      0.855  0.210 
##  7 t1          4      0.976  0.911 
##  8 t1          14     0.905  0.437 
##  9 t2          0      0.877  0.296 
## 10 t2          2      0.969  0.869 
## 11 t2          4      0.959  0.804 
## 12 t2          14     0.920  0.532 
## 13 t3          0      0.927  0.579 
## 14 t3          2      0.961  0.812 
## 15 t3          4      0.805  0.0885
## 16 t3          14     0.835  0.152

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

datos_ancho <- datos %>%
  pivot_wider(
    id_cols = c(tratamiento, replica, id_sujeto),
    names_from = tiempo,
    values_from = bacterias,
    names_prefix = "t"
  )

# Modelo para prueba de esfericidad
if(require(rstatix, quietly = TRUE)) {
  esfericidad <- datos %>%
    anova_test(dv = bacterias, 
               wid = id_sujeto,
               within = tiempo,
               between = tratamiento)
  print(esfericidad$`Mauchly's Test for Sphericity`)
  
  # Correcciones si se viola esfericidad
  cat("\nCorrecciones (usar si p < 0.05 en Mauchly):\n")
  print(esfericidad$`Sphericity Corrections`)
}

##               Effect     W     p p<.05
## 1             tiempo 0.855 0.805      
## 2 tratamiento:tiempo 0.855 0.805      
## 
## Correcciones (usar si p < 0.05 en Mauchly):
##               Effect   GGe      DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]
## 1             tiempo 0.905 2.71, 43.44 6.21e-19         * 1.108 3.32, 53.18
## 2 tratamiento:tiempo 0.905 8.14, 43.44 1.45e-04         * 1.108 9.97, 53.18
##      p[HF] p[HF]<.05
## 1 1.08e-20         *
## 2 7.21e-05         *

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

Greenhouse–Geisser (GG) = 0.905 y Huynh–Feldt (HFe) = 1.108

Los epsilons sólo se aplican si la esfericidad está violada (p < .05). En este caso, NO lo está, entonces estos valores son informativos y no modifican el análisis.

Analisis del Modelo

summary(modelo_aov)

## 
## Error: id_sujeto
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## tratamiento  3 2625.0   875.0   14.46 8.09e-05 ***
## Residuals   16  968.2    60.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: id_sujeto:tiempo
##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## tiempo              3  26140    8713  100.40  < 2e-16 ***
## tratamiento:tiempo  9   4007     445    5.13 7.21e-05 ***
## Residuals          48   4166      87                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

se observa que hay diferencia en el factor tiempo, así como del factor tratamiento a traves del tiempo, entre 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
resumen <- datos %>%
  group_by(tratamiento, tiempo) %>%
  summarise(
    n = n(),
    media = mean(bacterias),
    sd = sd(bacterias),
    se = sd(bacterias)/sqrt(n()),
    .groups = "drop"
  )

print(resumen)

## # A tibble: 16 × 6
##    tratamiento tiempo     n media    sd    se
##    <fct>       <fct>  <int> <dbl> <dbl> <dbl>
##  1 gc          0          5 101    8.97  4.01
##  2 gc          2          5 105    7.65  3.42
##  3 gc          4          5 103    9.14  4.09
##  4 gc          14         5 129.   8.64  3.87
##  5 t1          0          5  97.8 10.8   4.85
##  6 t1          2          5 120.  10.6   4.75
##  7 t1          4          5 115.   8.35  3.73
##  8 t1          14         5 167.  10.3   4.59
##  9 t2          0          5 108.   5.13  2.29
## 10 t2          2          5 103    5.43  2.43
## 11 t2          4          5 112.  10.3   4.62
## 12 t2          14         5 141.   4.22  1.89
## 13 t3          0          5 101.   4.27  1.91
## 14 t3          2          5 102.  14.0   6.26
## 15 t3          4          5 124    8.49  3.79
## 16 t3          14         5 156.  10.7   4.79

# Estadísticas generales
cat("\n=== ESTADÍSTICAS DESCRIPTIVAS GENERALES ===\n")

## 
## === ESTADÍSTICAS DESCRIPTIVAS GENERALES ===

datos %>%
  group_by(tratamiento) %>%
  summarise(
    n_obs = n(),
    media_general = mean(bacterias),
    sd_general = sd(bacterias),
    .groups = "drop"
  )

## # A tibble: 4 × 4
##   tratamiento n_obs media_general sd_general
##   <fct>       <int>         <dbl>      <dbl>
## 1 gc             20          110.       14.1
## 2 t1             20          125        27.9
## 3 t2             20          116.       16.2
## 4 t3             20          121.       24.7

# VISUALIZACIÓN EXPLORATORIA
# ----------------------------------------------------------------------------
# Gráfico de perfil de medias
ggplot(resumen, aes(x = tiempo, y = media, color = tratamiento, group = tratamiento)) +
  geom_line(size = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = media - se, ymax = media + se), width = 0.2) +
  labs(title = "Perfiles de Medias por Tratamiento",
       x = "Tiempo (días)", 
       y = "Bacterias (media ± SE)",
       color = "Tratamiento") +
  theme_bw() +
  theme(legend.position = "bottom")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Gráfico individual por réplica
ggplot(datos, aes(x = tiempo, y = bacterias, color = tratamiento)) +
  geom_line(aes(group = id_sujeto), alpha = 0.5) +
  geom_point(alpha = 0.7) +
  facet_wrap(~ tratamiento) +
  labs(title = "Trayectorias Individuales por Tratamiento",
       x = "Tiempo (días)", 
       y = "Bacterias") +
  theme_bw()

# Boxplot por tratamiento y tiempo
ggplot(datos, aes(x = tiempo, y = bacterias, fill = tratamiento)) +
  geom_boxplot() +
  labs(title = "Distribución de Bacterias por Tratamiento y Tiempo",
       x = "Tiempo (días)", 
       y = "Bacterias") +
  theme_bw()

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.

Proceso PostHoc

El proceso debe hacerse por cada Factor,

Respecto de los tratamientos

emm_ord <- emmeans(modelo_aov, ~ tratamiento)  # medias marginales por nivel de 'tratamiento'

## Note: re-fitting model with sum-to-zero contrasts

## NOTE: Results may be misleading due to involvement in interactions

pairs(emm_ord, adjust = "holm")

##  contrast estimate   SE df t.ratio p.value
##  gc - t1    -15.45 2.46 16  -6.281  0.0001
##  gc - t2     -6.30 2.46 16  -2.561  0.0628
##  gc - t3    -11.00 2.46 16  -4.472  0.0019
##  t1 - t2      9.15 2.46 16   3.720  0.0075
##  t1 - t3      4.45 2.46 16   1.809  0.1483
##  t2 - t3     -4.70 2.46 16  -1.911  0.1483
## 
## Results are averaged over the levels of: tiempo 
## P value adjustment: holm method for 6 tests

Análisis ajuste Hochberg

pairs(emm_ord, adjust = "hochberg")

##  contrast estimate   SE df t.ratio p.value
##  gc - t1    -15.45 2.46 16  -6.281  0.0001
##  gc - t2     -6.30 2.46 16  -2.561  0.0628
##  gc - t3    -11.00 2.46 16  -4.472  0.0019
##  t1 - t2      9.15 2.46 16   3.720  0.0075
##  t1 - t3      4.45 2.46 16   1.809  0.0893
##  t2 - t3     -4.70 2.46 16  -1.911  0.0893
## 
## Results are averaged over the levels of: tiempo 
## P value adjustment: hochberg method for 6 tests

de acuerdo con la respuesta con el ajuste Hochberg, los tratamientos 1 y 3 presentarían una diferencia con el grupo control, no así el tratamiento 2. Entre los tratamientos 1 y 3 tieneden a presentar similar respuesta

El método del contraste, permite analizar la diferencia como ud estime conveniente, generalmente se usa para analizar un grupo base, comunmente llamado CONTROL.

contrast(emm_ord,
         method = list("0 vs 1" = c( 1,-1, 0, 0),
                       "0 vs 2" = c( 1, 0,-1, 0),
                       "0 vs 3" = c( 1, 0, 0,-1)),
         adjust = "holm")

##  contrast estimate   SE df t.ratio p.value
##  0 vs 1      -15.4 2.46 16  -6.281  <.0001
##  0 vs 2       -6.3 2.46 16  -2.561  0.0209
##  0 vs 3      -11.0 2.46 16  -4.472  0.0008
## 
## Results are averaged over the levels of: tiempo 
## P value adjustment: holm method for 3 tests

en una mirada más específica respecto del grupo control, pareciera que al final del proceso (14 dias), los tratmiento provocan una diferencia significativa contra el grupo control.

plot(pairs(emm_ord, adjust = "holm"))

## Warning: `aes_()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`
## ℹ The deprecated feature was likely used in the emmeans package.
##   Please report the issue at <https://github.com/rvlenth/emmeans/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

cld(emm_ord, Letters = letters, adjust = "holm")

##  tratamiento emmean   SE df lower.CL upper.CL .group
##  gc             110 1.74 16      105      114  a    
##  t2             116 1.74 16      111      121  ab   
##  t3             121 1.74 16      116      125   bc  
##  t1             125 1.74 16      120      130    c  
## 
## Results are averaged over the levels of: tiempo 
## 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.

Respecto de los tiempo

emm_ord <- emmeans(modelo_aov, ~ tiempo)  # medias marginales por nivel de 'tiempo'

## Note: re-fitting model with sum-to-zero contrasts

## NOTE: Results may be misleading due to involvement in interactions

pairs(emm_ord, adjust = "holm")

##  contrast           estimate   SE df t.ratio p.value
##  tiempo0 - tiempo2     -5.55 2.95 48  -1.884  0.0781
##  tiempo0 - tiempo4    -11.80 2.95 48  -4.005  0.0006
##  tiempo0 - tiempo14   -46.40 2.95 48 -15.750  <.0001
##  tiempo2 - tiempo4     -6.25 2.95 48  -2.122  0.0781
##  tiempo2 - tiempo14   -40.85 2.95 48 -13.866  <.0001
##  tiempo4 - tiempo14   -34.60 2.95 48 -11.745  <.0001
## 
## Results are averaged over the levels of: tratamiento 
## P value adjustment: holm method for 6 tests