En la clase del diseño split-plot se mencionó que los diseños de medidas repetidas pueden se considerados estadísticamente análogos, pues los tratamientos están anidados en las unidades experimentales. La principal diferencia entre ambos diseños radica en que el anidamiento en los diseños split-plot es espacial, mientras que en los diseños de medidas repetidas el anidamiento es temporal. Las mediciones repetidas en las unidades experimentales pueden corresponder a la aplicación secuencial de tratamientos, o a la medición del efecto acumulado del mismo tratamiento. Debido a este componente temporal, los diseños de medidas repetidas tienden a durar más tiempo que otros diseños donde se aplican todos los tratamientos al mismo tiempo pero en unidades experimentales diferentes.
Como podemos ver, en los diseños de medidas repetidas hay una serie de trueques. Por un lado, estos diseños permiten que los experimentos se realicen con menos unidades experimentales, pero requieren de más tiempo. Existen además importantes consideraciones biológicas, pues en necesario que los tratamientos no tendrán un efecto sobre el comportamiento de la unidad experimental que ya recibió otro tratamiento.
Resulta evidente entonces que, como en todos los diseños experimentales, la decisión de hacer mediciones repetidas debe estar basada en en una serie de criterios:
El análisis de los diseños con medidas repetidas es igual al de los diseños split-plot, porque en esencia consiste de un ANOVA anidado. Por lo tanto, la mayor complicación del análisis es la identificación correcta de la estructura de anidación —¿qué tratamientos están contenidos dentro de cada unidad experimental?
Aquí haremos otro análisis anidado de un diseño de tres vías con medidas repetidas. La base analizada está contenida en datarium
bajo el nombre weightloss
.
Weightloss
Comenzaremos, nuevamente, cargando la base de datos, viendo en qué formato está y si es necesario transformarla.
weight <- datarium::weightloss
knitr::kable(head(weight)) # Imprimiendo las primeras seis filas de weightloss
id | diet | exercises | t1 | t2 | t3 |
---|---|---|---|---|---|
1 | no | no | 10.43 | 13.21 | 11.59 |
2 | no | no | 11.59 | 10.66 | 13.21 |
3 | no | no | 11.35 | 11.12 | 11.35 |
4 | no | no | 11.12 | 9.50 | 11.12 |
5 | no | no | 9.50 | 9.73 | 12.28 |
6 | no | no | 9.50 | 12.74 | 10.43 |
Al igual que los otros diseños anidados que vimos (slefesteem2
), está en formato ancho en el factor tiempo, por lo que hay que utilizar reshape2::melt
para darle formato largo. Para que melt
funcione como queremos (mantener una columna para cada factor experimental), necesitamos especificar qué columnas permanecerán así en el argumento id.vars
:
weight.l <- reshape2::melt(weight, id.vars = c("id", "diet", "exercises"))
names(weight.l) <- c("id", "diet", "exercises", "time", "loss")
knitr::kable(head(weight.l)) #Imprimiendo las primeras seis filas de la tabla resultante
id | diet | exercises | time | loss |
---|---|---|---|---|
1 | no | no | t1 | 10.43 |
2 | no | no | t1 | 11.59 |
3 | no | no | t1 | 11.35 |
4 | no | no | t1 | 11.12 |
5 | no | no | t1 | 9.50 |
6 | no | no | t1 | 9.50 |
Al igual que selfesteem2
, weightloss
contiene tres mediciones de la pérdida de peso por cada tratamiento aplicado a los individuos. Los factores diet
y exercises
fueron aplicados simultáneamente, por lo que sólo time
está anidado en ambos, lo cual constituye la principal diferencia con selfesteem2
que sólo incluye dos factores anidados en los individuos. Con base en esta estructura experimental, podemos contemplar la existencia de interacciones entre diet
y exercises
. La interacción con tiempo debe existir porque está anidado en ambos factores.
Antes de continuar con el análisis hagamos una rápida verificación de la distribución de loss
, esta vez, utilizando únicamente la pruebe de shapiro para cada grupo experimental:
library(tidyverse); library(rstatix)
shap.test <- weight.l %>%
group_by(diet, exercises, time) %>%
shapiro_test(loss)
knitr::kable(shap.test, caption = "Resultados de la prueba de Shapiro dentro de cada grupo experimental.", align = "c")
diet | exercises | time | variable | statistic | p |
---|---|---|---|---|---|
no | no | t1 | loss | 0.9173082 | 0.2643971 |
no | no | t2 | loss | 0.9571979 | 0.7432089 |
no | no | t3 | loss | 0.9649090 | 0.8509028 |
no | yes | t1 | loss | 0.9223457 | 0.3059440 |
no | yes | t2 | loss | 0.9123400 | 0.2285804 |
no | yes | t3 | loss | 0.9525447 | 0.6744326 |
yes | no | t1 | loss | 0.9422786 | 0.5281443 |
yes | no | t2 | loss | 0.9816773 | 0.9894564 |
yes | no | t3 | loss | 0.9313724 | 0.3948864 |
yes | yes | t1 | loss | 0.9141565 | 0.2411125 |
yes | yes | t2 | loss | 0.9472878 | 0.5977346 |
yes | yes | t3 | loss | 0.9373215 | 0.4641723 |
Los resultados sugieren que podemos aceptar la hipótesis nula de no diferencias con una distribución normal.
Entonces procedamos con el análisis de weightloss
. Para comenzar, sabemos que puede haber interacción diet:exercises
, por lo que comenzaremos con el modelo más complicado:
library(lme4)
m1 <- lmer(loss ~ (diet/time) * (exercises/time) + (1|id/exercises/time) + (1|id/diet/time), weight.l)
car::Anova(m1, type = 2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: loss
## Chisq Df Pr(>Chisq)
## diet 3.8387 1 0.050084 .
## exercises 57.7072 1 3.042e-14 ***
## diet:time 16.3963 4 0.002531 **
## time:exercises 54.6825 2 1.336e-12 ***
## diet:exercises 31.8969 1 1.626e-08 ***
## diet:time:exercises 29.3907 2 4.149e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: loss ~ (diet/time) * (exercises/time) + (1 | id/exercises/time) +
## (1 | id/diet/time)
## Data: weight.l
##
## REML criterion at convergence: 441.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.24021 -0.60718 0.05678 0.60456 2.35978
##
## Random effects:
## Groups Name Variance Std.Dev.
## time..diet.id. (Intercept) 1.422e-01 3.771e-01
## time..exercises.id. (Intercept) 9.305e-02 3.050e-01
## diet.id (Intercept) 0.000e+00 0.000e+00
## exercises.id (Intercept) 8.895e-10 2.983e-05
## id (Intercept) 7.457e-02 2.731e-01
## id.1 (Intercept) 4.236e-03 6.509e-02
## Residual 1.047e+00 1.023e+00
## Number of obs: 144, groups:
## time:(diet:id), 72; time:(exercises:id), 72; diet:id, 24; exercises:id, 24; id, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.9092 0.3368 32.391
## dietyes 0.8333 0.4452 1.872
## exercisesyes -0.1150 0.4359 -0.264
## dietno:timet2 0.6567 0.4623 1.420
## dietyes:timet2 0.6733 0.4623 1.456
## dietno:timet3 0.5408 0.4623 1.170
## dietyes:timet3 2.0442 0.4623 4.422
## timet2:exercisesyes 1.9700 0.6165 3.196
## timet3:exercisesyes 5.4825 0.6165 8.893
## dietyes:exercisesyes -0.2342 0.5908 -0.396
## dietyes:timet2:exercisesyes -0.8117 0.8355 -0.971
## dietyes:timet3:exercisesyes -4.2650 0.8355 -5.105
##
## Correlation of Fixed Effects:
## (Intr) dietys exrcss dtn:t2 dtys:2 dtn:t3 dtys:3 tmt2:x tmt3:x
## dietyes -0.661
## exercisesys -0.647 0.450
## dietno:tmt2 -0.686 0.482 0.471
## dietys:tmt2 -0.050 -0.482 0.038 0.073
## dietno:tmt3 -0.686 0.482 0.471 0.500 0.036
## dietys:tmt3 -0.050 -0.482 0.038 0.036 0.500 0.073
## tmt2:xrcssy 0.458 -0.318 -0.707 -0.667 -0.054 -0.333 -0.027
## tmt3:xrcssy 0.458 -0.318 -0.707 -0.333 -0.027 -0.667 -0.054 0.500
## dtys:xrcssy 0.439 -0.663 -0.678 -0.319 0.319 -0.319 0.319 0.479 0.479
## dtys:tmt2:x -0.310 0.469 0.479 0.452 -0.452 0.226 -0.226 -0.678 -0.339
## dtys:tmt3:x -0.310 0.469 0.479 0.226 -0.226 0.452 -0.452 -0.339 -0.678
## dtys:x dty:2:
## dietyes
## exercisesys
## dietno:tmt2
## dietys:tmt2
## dietno:tmt3
## dietys:tmt3
## tmt2:xrcssy
## tmt3:xrcssy
## dtys:xrcssy
## dtys:tmt2:x -0.707
## dtys:tmt3:x -0.707 0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Esta especificación del modelo, con tiempo
anidado en diet
y exercises
parece tener todos los términos relevantes para los datos, por ahora sólo falta examinar los residuales. Como ya se habrán dado cuenta, las gráficas de diagnóstico del paquete lmer
sólo muestran la primera gráfica de aov
y lm
, pero ésta es suficiente:
plot(m1)
Con este diagnóstico, podemos entonces rechazar todas las hipótesis nulas. Para ver rápidamente los efectos, podemos ahora sí, ver los grupos experimentales de loss
:
boxplot(loss ~ diet + exercises + time, weight.l)
Lo que sugiere que los efectos de dieta, ejercicio y tiempo, son todos positivos, es decir, que la dieta, ejercicio y tiempo sobre la pérdida de peso dependen de que todos estén presentes para que haya una mayor pérdida de peso.