En un diseño de medidas repetidas (RM-ANOVA) se evalúa el efecto de un factor sobre una variable independiente medida en dos o más ocasiones para un mismo individuo o para varios individuos que pertenecen a un mismo grupo (e.g., miembros de una misma familia). Como los múltiples valores en la VD están correlacionados entre ellos ya no es posible usar ANOVA en su forma tradicional, donde los errores de estimación son independientes entre los niveles del factor.

library(readxl)
Base <- read_excel("Bases de datos/Base-AMR.xls", 
    sheet = "Base")
attach(Base)

Explorando base de datos

Los datos corresponden a un experimento en el que se estudia un tratamiento para la depresión. Se ha llevado a cabo un seguimiento de dos grupos de pacientes (1: control / 2: tratamiento) en cinco ocasiones diferentes: dv0: pre-test dv1: post-test después de un mes dv3: 3 meses de seguimiento dv6: 6 meses de seguimiento)

La variable dependiente es la puntuación en depresión.

Tiempo=factor(Tiempo)
View(Base)
Base
## # A tibble: 96 x 3
##    Grupo Tiempo Respuesta
##    <dbl>  <dbl>     <dbl>
##  1     1      0       296
##  2     1      1       175
##  3     1      3       187
##  4     1      6       242
##  5     1      0       376
##  6     1      1       329
##  7     1      3       236
##  8     1      6       126
##  9     1      0       309
## 10     1      1       238
## # ... with 86 more rows
names(Base)
## [1] "Grupo"     "Tiempo"    "Respuesta"
summary(Base)
##      Grupo         Tiempo       Respuesta    
##  Min.   :1.0   Min.   :0.00   Min.   :  6.0  
##  1st Qu.:1.0   1st Qu.:0.75   1st Qu.:111.5  
##  Median :1.5   Median :2.00   Median :174.5  
##  Mean   :1.5   Mean   :2.50   Mean   :194.5  
##  3rd Qu.:2.0   3rd Qu.:3.75   3rd Qu.:291.2  
##  Max.   :2.0   Max.   :6.00   Max.   :447.0
ddply(Base, .(Tiempo), summarize, media = mean(Respuesta), median = median(Respuesta), desviacion=sd(Respuesta), minimo=min(Respuesta), maximo=max(Respuesta))
##   Tiempo    media median desviacion minimo maximo
## 1      0 292.3750  309.5   73.86757    138    447
## 2      1 175.5833  138.5  116.21267     31    402
## 3      3 158.0417  142.5   88.77939     16    334
## 4      6 151.9167  137.0   78.02002      6    358

Grafico de cajas sin ggplot

boxplot(Respuesta~Tiempo, xlab="tiempo", ylab="Respuesta a los tratamientos con el tiempo", main="Respuesta de los pacientes con el paso del tiempo", col="wheat", data=Base)
stripchart(Respuesta~Tiempo,
           method = "jitter",add = T, vertical = T,
           pch=19, col = "black")
medias=tapply(Respuesta, Tiempo, mean)
medias
##        0        1        3        6 
## 292.3750 175.5833 158.0417 151.9167
points(medias,pch=16,col ="Red")###Mostrar las medias
abline(h=mean(Respuesta),lty=2, col ="red")##mostar linea

ANOVA

Modelo1=aov(Respuesta~Tiempo)
summary(Modelo1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Tiempo       3 313918  104639   12.71 5.06e-07 ***
## Residuals   92 757406    8233                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.lm(Modelo1)
## 
## Call:
## aov(formula = Respuesta ~ Tiempo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -154.375  -67.875   -8.979   63.552  226.417 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   292.38      18.52  15.786  < 2e-16 ***
## Tiempo1      -116.79      26.19  -4.459 2.32e-05 ***
## Tiempo3      -134.33      26.19  -5.129 1.61e-06 ***
## Tiempo6      -140.46      26.19  -5.363 6.07e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 90.73 on 92 degrees of freedom
## Multiple R-squared:  0.293,  Adjusted R-squared:   0.27 
## F-statistic: 12.71 on 3 and 92 DF,  p-value: 5.062e-07
TukeyHSD(Modelo1)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Respuesta ~ Tiempo)
## 
## $Tiempo
##           diff        lwr       upr     p adj
## 1-0 -116.79167 -185.32782 -48.25551 0.0001349
## 3-0 -134.33333 -202.86949 -65.79718 0.0000095
## 6-0 -140.45833 -208.99449 -71.92218 0.0000036
## 3-1  -17.54167  -86.07782  50.99449 0.9082451
## 6-1  -23.66667  -92.20282  44.86949 0.8029636
## 6-3   -6.12500  -74.66116  62.41116 0.9954666
plot(TukeyHSD(Modelo1))

modelo2=aov(Respuesta~Tiempo+Grupo)
summary(modelo2)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Tiempo       3 313918  104639   18.07 2.78e-09 ***
## Grupo        1 230496  230496   39.81 9.98e-09 ***
## Residuals   91 526910    5790                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.lm(modelo2)
## 
## Call:
## aov(formula = Respuesta ~ Tiempo + Grupo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -191.375  -47.677    5.188   46.969  177.417 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   439.38      28.00  15.691  < 2e-16 ***
## Tiempo1      -116.79      21.97  -5.317 7.49e-07 ***
## Tiempo3      -134.33      21.97  -6.115 2.37e-08 ***
## Tiempo6      -140.46      21.97  -6.394 6.81e-09 ***
## Grupo         -98.00      15.53  -6.309 9.98e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 76.09 on 91 degrees of freedom
## Multiple R-squared:  0.5082, Adjusted R-squared:  0.4866 
## F-statistic: 23.51 on 4 and 91 DF,  p-value: 2.291e-13
TukeyHSD(modelo2)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: Grupo
## Warning in TukeyHSD.aov(modelo2): 'which' specified some non-factors which will
## be dropped
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Respuesta ~ Tiempo + Grupo)
## 
## $Tiempo
##           diff        lwr       upr     p adj
## 1-0 -116.79167 -174.28067 -59.30266 0.0000044
## 3-0 -134.33333 -191.82234 -76.84433 0.0000001
## 6-0 -140.45833 -197.94734 -82.96933 0.0000000
## 3-1  -17.54167  -75.03067  39.94734 0.8548972
## 6-1  -23.66667  -81.15567  33.82234 0.7041351
## 6-3   -6.12500  -63.61401  51.36401 0.9923814
Repeated.Measures.ANOVA <- aov(Respuesta ~ Tiempo + Grupo + Tiempo*Grupo)
summary(Repeated.Measures.ANOVA)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## Tiempo        3 313918  104639  19.713 7.31e-10 ***
## Grupo         1 230496  230496  43.423 3.13e-09 ***
## Tiempo:Grupo  3  59792   19931   3.755   0.0137 *  
## Residuals    88 467118    5308                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ahora podemos analizar las dos tablas siguientes. En primer lugar, las pruebas sobre los efectos inter-sujetos que muestran el efecto de la variable grupo en el conjunto de datos sin tener en cuenta las repeticiones (o medidas). Vemos que el grupo tiene un impacto significativo en la puntuación en depresión. A continuación, las pruebas sobre los efectos intra-sujetos muestran el impacto del tiempo (de las diferentes medidas) sobre la variable dependiente. Puede ser útil examinar los términos de interacción entre la repetición y los factores explicativos. Vemos que el factor de repetición tiene un impacto significativo en la puntuación en depresión; la interacción tiene también un impacto significativo.

Este estudio ha mostrado que tanto el tiempo como el tratamiento tienen un impacto significativo en la puntuación en depresión.

Supuestos del Modelo

Normalidad

qqnorm(modelo2$residuals) 
qqline(modelo2$residuals,col = rainbow(4))

par(mfrow = c(1,2))
plot(modelo2, which = 1:4)

shapiro.test(modelo2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo2$residuals
## W = 0.99478, p-value = 0.9738
library(nortest)

lillie.test(residuals(modelo2))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuals(modelo2)
## D = 0.057551, p-value = 0.6049

Varianza constante test de Breush Pagan

library(lmtest)
## Warning: package 'lmtest' was built under R version 4.0.5
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(modelo2)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo2
## BP = 9.9276, df = 4, p-value = 0.04167

De la salida anterior se observa que el valor-P es menor que el nivel de significancia usual de 5%, por lo tanto, hay evidencias para decir que no se cumple la homocedasticidad

Independencia

#independencia
library(lmtest)
dwtest(modelo2,Grupo)
## 
##  Durbin-Watson test
## 
## data:  modelo2
## DW = 1.1633, p-value = 1.164e-05
## alternative hypothesis: true autocorrelation is greater than 0
dwtest(modelo2,Tiempo)
## 
##  Durbin-Watson test
## 
## data:  modelo2
## DW = 1.9712, p-value = 0.2921
## alternative hypothesis: true autocorrelation is greater than 0
plot(modelo2$residuals)

dwtest(modelo2)
## 
##  Durbin-Watson test
## 
## data:  modelo2
## DW = 1.1633, p-value = 1.164e-05
## alternative hypothesis: true autocorrelation is greater than 0

En el ANOVA clásico inter-sujetos hemos visto que debe cumplirse el supuesto de que todas las puntuaciones de diferentes condiciones sean independientes. Cuando tenemos medidas repetidas no podemos cumplir con este supuesto ya que las mediciones están relacionadas porque provienen de un mismo sujeto. Por lo tanto, debemos “cambiar” este supuesto por otro válido.

Test for sphericity – Mauchly’s test

neo <- matrix(Respuesta, nrow = 25, ncol = 4)
## Warning in matrix(Respuesta, nrow = 25, ncol = 4): la longitud de los datos [96]
## no es un submúltiplo o múltiplo del número de filas [25] en la matriz
neo
##       [,1] [,2] [,3] [,4]
##  [1,]  296  364  225    6
##  [2,]  175  270  134  292
##  [3,]  187  358  317  139
##  [4,]  242  447   31  104
##  [5,]  376  402   85  184
##  [6,]  329  294  120  275
##  [7,]  236  266  362   94
##  [8,]  126  220  104  135
##  [9,]  309   70  144  137
## [10,]  238   95  114  150
## [11,]  150  137  338   48
## [12,]  173  375  132   20
## [13,]  222  335   91   85
## [14,]   60  334   77  319
## [15,]   82  129  263   68
## [16,]  135  310   94   67
## [17,]  150  300  141   12
## [18,]  271  253  142  300
## [19,]  250  170  138  138
## [20,]  266  310   38  114
## [21,]  316  245   16  174
## [22,]  291  200   95  296
## [23,]  238  170  329  175
## [24,]  194  282   62  187
## [25,]  321  186   62  242
mauchly.test(lm(neo ~ 1), X = ~ 1)
## 
##  Mauchly's test of sphericity
##  Contrasts orthogonal to
##  ~1
## 
## 
## data:  SSD matrix from lm(formula = neo ~ 1)
## W = 0.81581, p-value = 0.4636

La salida devuelve un valor p superior a 0,05, lo que significa que se acepta la hipótesis nula que establece que no hay diferencia en las varianzas para todas las comparaciones de grupos por pares.