En el Anova de medidas repetidas se analizan los resultados obtenidos en un diseño experimental en donde se ha manipulado una única variable independiente (un único factor) con 2 o más niveles pero de forma intra-sujeto. Esto viene a significar, que todos los individuos (o unidades de observación), han pasado por todos los niveles del factor. A este tipo de diseños también se les conoce como diseños de medidas repetidas en el sentido de que a cada sujeto se le repite la medición de la variable dependiente de respuesta en diversas condiciones, tantas como niveles tenga el factor manipulado. También se les conoce como diseños de medidas dependientes debido a que las puntuaciones de un mismo sujeto muestran dependencia estadística entre ellas, están relacionadas.

Se trata de comprobar, por ejemplo, desde un punto de vista de contrastes paramétricos, si las medias de una medición de una variable dependiente continua, son las mismas en distintos momentos del tiempo (3 o más), por ejemplo si un tratamiento de combinación de dieta y ejercicio (las 2 cosas como un todo) tiene la correspondiente contraprestación en pérdida de peso en 4 instantes del tiempo, esto es, se trata de 4 muestras relacionadas, dependientes o pareadas, tomadas sobre los mismos 10 individuos. Es algo parecido a la prueba de la T de Student para muestras relacionadas en un antes y un después, pero esta vez con 3 o más instantes de tiempo, en lugar de 2.

Importar datos

#ejercicio
#https://estamatica.net/anova-de-medidas-repetidas/
library(openxlsx)
url3<-"https://docs.google.com/spreadsheets/d/e/2PACX-1vTSnSZS9JeD0AjKOwD3pNECLaQODO7T0u4E8xKi6_pts5QOAY5R_GtXDiqEl5Ta83Q-e8s0zk0Z7G1G/pub?output=xlsx"
df = read.xlsx(url3,sheet="18_rep")
head(df,7)
##   individuo peso etapa
## 1         1   80     1
## 2         2   77     1
## 3         3   79     1
## 4         4   56     1
## 5         5   68     1
## 6         6   65     1
## 7         7   94     1

Exploracion de datos

boxplot(df$peso~df$etapa, col="white", 
        border =c("#5e35b1","#2196f3","#ffc400","#00bcd4"))

#ver la estructura de los datos
str(df)
## 'data.frame':    40 obs. of  3 variables:
##  $ individuo: num  1 2 3 4 5 6 7 8 9 10 ...
##  $ peso     : num  80 77 79 56 68 65 94 77 76 81 ...
##  $ etapa    : num  1 1 1 1 1 1 1 1 1 1 ...
#transformar a factor
df$individuo<-factor(df$individuo)
df$etapa<-factor(df$etapa)
#ver nuevamente la estructura
str(df)
## 'data.frame':    40 obs. of  3 variables:
##  $ individuo: Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ peso     : num  80 77 79 56 68 65 94 77 76 81 ...
##  $ etapa    : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...

Test de Normalidad

Si el tamaño muestral es lo suficientemente grande (n>30 o n>50, dependiendo del autor), se puede asumir normalidad, independientemente de cómo se comporten las variables en la población, pudiendo proceder de un modo paramético con el Test ANOVA de medidas repetidas con SPSS. En igualdad de condiciones, un test paramétrico (ANOVA) siempre va a ser mucho más potente que un test no paramétrico (FRIEDMAN), es decir, vamos a tener mayor probabilidad de encontrar diferencias estadísticamente significativas, que es lo que se trata, disminuyendo además la probabilidad de que esas diferencias encontradas, sean debidas al azar. 

shapiro.test(df$peso[df$etapa==1])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$peso[df$etapa == 1]
## W = 0.94165, p-value = 0.5715
shapiro.test(df$peso[df$etapa==2])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$peso[df$etapa == 2]
## W = 0.96822, p-value = 0.8739
shapiro.test(df$peso[df$etapa==3])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$peso[df$etapa == 3]
## W = 0.97154, p-value = 0.9048
shapiro.test(df$peso[df$etapa==4])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$peso[df$etapa == 4]
## W = 0.98787, p-value = 0.9935

Desarrollo del modelo

Las variables, cumplen el supuesto

Vamos a realizar el análisis de medidas repetidas ANOVA paramétrico por medio del paquete ez (install.packages(“ez”)). Para ello hay que indicar nuestros datos (data=), nuestra variable dependiente (dv=), nuestros individuos (wid=) y el tiempo (within=):

library(ez)
#Aplicamos un medidas repetidas:
ezANOVA(data=df, dv=peso, wid=individuo, within=etapa)
## $ANOVA
##   Effect DFn DFd        F           p p<.05       ges
## 2  etapa   3  27 41.48314 3.01836e-10     * 0.1266346
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W           p p<.05
## 2  etapa 0.1160743 0.005649013     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2  etapa 0.4582623 1.034209e-05         * 0.5135537 3.525929e-06         *

Podemos ver en nuestros resultados que el tiempo es significativo ($ANOVA el p es menor de 0.05), es decir el rendimiento escolar cambia con el tiempo. Pero en este caso nuestros datos violan uno de los requisitos que es la esfericidad (Mauchly’s test es significativo) , debemos fiarnos de la p de la ‘Sphericity corrections’ que nos confirma lo que hemos deducido al principio.

#Aplicamos un medidas repetidas:
ezANOVA(data=df, dv=peso, wid=individuo, within=etapa, detailed = TRUE)
## $ANOVA
##        Effect DFn DFd        SSn      SSd         F           p p<.05       ges
## 1 (Intercept)   1   9 202635.225 2899.525 628.97096 1.22402e-09     * 0.9854413
## 2       etapa   3  27    434.075   94.175  41.48314 3.01836e-10     * 0.1266346
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W           p p<.05
## 2  etapa 0.1160743 0.005649013     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2  etapa 0.4582623 1.034209e-05         * 0.5135537 3.525929e-06         *

Hay diferencias estadísticamente significativas (p-valor <0,05), entre 4 los instantes de medición del peso en Kgs de manera global (p-valores de Prueba Multivariante), como se pretendía demostrar en la investigación de eficacia de los tratamientos de pérdida de peso. La Hipótesis nula de Mauchly es de homogeneidad de varianzas de las diferencias medias, al resultar esta prueba significativa (hay heterogeneidad de varianzas), se tiene que tener en cuenta las pruebas de correciones de esfericidad, como lo es GGe (Greenhouse-Geisser epsilon) la cual sí es significativa.

#Aplicamos un medidas repetidas:
ezANOVA(data=df, dv=peso, wid=individuo, within=etapa, detailed = TRUE, return_aov = TRUE)
## $ANOVA
##        Effect DFn DFd        SSn      SSd         F           p p<.05       ges
## 1 (Intercept)   1   9 202635.225 2899.525 628.97096 1.22402e-09     * 0.9854413
## 2       etapa   3  27    434.075   94.175  41.48314 3.01836e-10     * 0.1266346
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W           p p<.05
## 2  etapa 0.1160743 0.005649013     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2  etapa 0.4582623 1.034209e-05         * 0.5135537 3.525929e-06         *
## 
## $aov
## 
## Call:
## aov(formula = formula(aov_formula), data = data)
## 
## Grand Mean: 71.175
## 
## Stratum 1: individuo
## 
## Terms:
##                 Residuals
## Sum of Squares   2899.525
## Deg. of Freedom         9
## 
## Residual standard error: 17.94908
## 
## Stratum 2: individuo:etapa
## 
## Terms:
##                   etapa Residuals
## Sum of Squares  434.075    94.175
## Deg. of Freedom       3        27
## 
## Residual standard error: 1.867609
## Estimated effects may be unbalanced
#Más medidas
#https://www.rdocumentation.org/packages/ez/versions/4.4-0/topics/ezANOVA
#eta cuadrado parcila
#https://jllopisperez.com/2013/06/26/la-eta-cuadrada-y-la-eta-cuadrada-parcial/

Y con estas pruebas, una vez se ha rechazado la homodedasticidad, se comprueba que efectivamente hay diferencias estadísticamente significativas entre los instantes de las mediciones de los pesos, pues el valor de probabilidad tanto de Greenhouse-Geisser, como de HFe medidas, es sensiblemente menor que 0,05, como queríamos comprobar con la investigación. Además, con el Eta cuadrado de 0.12 (Generalized Eta-Squared measure of effect size see in references below: Bakeman, 2005.)

Suele considerarse que una eta cuadrada en torno a 0,01 es poco efecto, que una eta cuadrada en torno a 0,06 indica un efecto medio y que una eta cuadrada superior a 0,14 es ya un efecto grande.

Para obtener el eta cuadrado parcial = 434.075/(434.075 + 94.175)

El eta parcial al cuadrado se obtiene que el modelo explica el 82,2% del comportamiento de la variable dependiente, o lo que es lo mismo, el 82% de la pérdida de peso es debida a los tratamientos.

Segundo Método

aov(data = df, peso ~ etapa + Error(individuo/etapa))  # method 2
## 
## Call:
## aov(formula = peso ~ etapa + Error(individuo/etapa), data = df)
## 
## Grand Mean: 71.175
## 
## Stratum 1: individuo
## 
## Terms:
##                 Residuals
## Sum of Squares   2899.525
## Deg. of Freedom         9
## 
## Residual standard error: 17.94908
## 
## Stratum 2: individuo:etapa
## 
## Terms:
##                   etapa Residuals
## Sum of Squares  434.075    94.175
## Deg. of Freedom       3        27
## 
## Residual standard error: 1.867609
## Estimated effects may be unbalanced
anova1<-aov(data = df, peso ~ etapa + Error(individuo/etapa))  # method 2
summary(anova1)
## 
## Error: individuo
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9   2900   322.2               
## 
## Error: individuo:etapa
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## etapa      3  434.1  144.69   41.48 3.02e-10 ***
## Residuals 27   94.2    3.49                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tercer Método

library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
modelo<-anova_test(data=df,dv=peso,wid=individuo,within=etapa)
get_anova_table(modelo)
## ANOVA Table (type III tests)
## 
##   Effect  DFn   DFd      F        p p<.05   ges
## 1  etapa 1.37 12.37 41.483 1.03e-05     * 0.127
#con eta cuadrada parcial
modelo<-anova_test(data=df,dv=peso,wid=individuo,within=etapa, effect.size = "pes")
get_anova_table(modelo)
## ANOVA Table (type III tests)
## 
##   Effect  DFn   DFd      F        p p<.05   pes
## 1  etapa 1.37 12.37 41.483 1.03e-05     * 0.822
modelo$`Mauchly's Test for Sphericity`
##   Effect     W     p p<.05
## 1  etapa 0.116 0.006     *
modelo$`Sphericity Corrections`
##   Effect   GGe      DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]    p[HF]
## 1  etapa 0.458 1.37, 12.37 1.03e-05         * 0.514 1.54, 13.87 3.53e-06
##   p[HF]<.05
## 1         *
modelo$ANOVA
##   Effect DFn DFd      F        p p<.05   pes
## 1  etapa   3  27 41.483 3.02e-10     * 0.822
pairwise_t_test(data=df, formula=peso~etapa, paired=T,p.adjust.method = "bonferroni")
## # A tibble: 6 x 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 peso  1      2         10    10      3.93     9   3   e-3 2.1 e-2 *           
## 2 peso  1      3         10    10      4.36     9   2   e-3 1.1 e-2 *           
## 3 peso  1      4         10    10      8.02     9   2.17e-5 1.3 e-4 ***         
## 4 peso  2      3         10    10      3.34     9   9   e-3 5.2 e-2 ns          
## 5 peso  2      4         10    10      9.80     9   4.24e-6 2.54e-5 ****        
## 6 peso  3      4         10    10      9.78     9   4.32e-6 2.59e-5 ****

Todos son significativos

Comparación por parejas

 bonf1 <- pairwise.t.test(df$peso,df$etapa, method="bonf")  #pd.adj="bonferroni"
 print(bonf1)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  df$peso and df$etapa 
## 
##   1    2    3   
## 2 1.00 -    -   
## 3 0.88 1.00 -   
## 4 0.21 0.63 1.00
## 
## P value adjustment method: holm

Esta es la forma correcta

#pairwise comparisons
#library(asbio)
 bonf <- pairwise.t.test(df$peso,df$etapa, pd.adj="bonferroni", paired = T)  #pd.adj="sidak"
 print(bonf)
## 
##  Pairwise comparisons using paired t tests 
## 
## data:  df$peso and df$etapa 
## 
##   1       2       3      
## 2 0.0070  -       -      
## 3 0.0054  0.0086  -      
## 4 8.7e-05 2.5e-05 2.5e-05
## 
## P value adjustment method: holm

Todos aquellos tratamientos en los que en la comparativa 2 a 2 presentan un p-valor (Sig,) menor que 0,05, muestran diferencias estadísticamente significativas en media en el cruce de los 2 momentos en el tiempo, esto es, respecto a la pérdida de peso, objeto de la investigación. La prueba Post Hoc de Bonferroni es más conservadora, pero quizá Sidak sea la más utilizada en este tipo de ANOVA.

summary(df[df$etapa==1,])
##    individuo      peso       etapa 
##  1      :1   Min.   :56.00   1:10  
##  2      :1   1st Qu.:70.00   2: 0  
##  3      :1   Median :77.00   3: 0  
##  4      :1   Mean   :75.30   4: 0  
##  5      :1   3rd Qu.:79.75         
##  6      :1   Max.   :94.00         
##  (Other):4

Friedman

Friedman: Equivalente no paramétrico

En el caso de no cumplirse el supuesto de Normalidad en los datos, o de tamaños muestrales menores de 50, se procede desde un punto de vista no paramétrico con el Test de Friedman, el más conveniente para muestras relacionadas/dependientes, donde se viola este punto de partida:

friedman.test(df$peso, df$etapa, df$individuo)
## 
##  Friedman rank sum test
## 
## data:  df$peso, df$etapa and df$individuo
## Friedman chi-squared = 25.68, df = 3, p-value = 1.113e-05

Como la probabilidad asociada al estadístico de Friedman resulta estadísticamente significativa (valor p menor de 0,05), se rechaza el que el comportamiento en Mediana sea el mismo para los 4 grupos, por lo habría que testear donde se muestran tales diferencias con el Test de Wilcoxon en comparativas 2 a 2, equivalente no paramétrico a la prueba T de Student de muestras pareadas/dependientes (un antes y un después), o con la correspondiente prueba de comparaciones múltiples no paramétrica.

pairwise.wilcox.test(df$peso, df$etapa, paired = TRUE, p.adjust.method = "holm")
## 
##  Pairwise comparisons using Wilcoxon signed rank test with continuity correction 
## 
## data:  df$peso and df$etapa 
## 
##   1     2     3    
## 2 0.033 -     -    
## 3 0.033 0.033 -    
## 4 0.033 0.033 0.033
## 
## P value adjustment method: holm