DE junio 02
Paula Rodriguez
2023-06-02
MEDIDAS REPETIDAS
https://www.datanovia.com/en/lessons/repeated-measures-anova-in-r/
¿Cuándo me interesa?
- Cuando ingresa un factor intrasujetos (tiempo)
- Factores entresujetos (FSCA, FSBA, FCCA, FCBA) En este momento SON TODOS AQUELLOS MODELOS QUE HEMOS VISTO PERO CON EL TIEMPO (SEGUIMIENTO).
Las medidas repetidas pueden ser de una vía, dos vías, tres vías.
Medidas repetidas una via (Un solo factor: Tiempo)
Medición de aceites (Limonaria)
data("selfesteem", package = "datarium")
datos = selfesteem
head(datos, 3)## id t1 t2 t3
## 1 1 4.005027 5.182286 7.107831
## 2 2 2.558124 6.912915 6.308434
## 3 3 3.244241 4.443434 9.778410Datos: Tiempos equidistantes * 10 parcelas y tres tiempos: - Primera columna -> Número de parcela - Segunda columna -> 30 días después de siembra - Tercera columna -> 60 días después de siembra - Cuarta columna -> 90 días después de siembra
Funcion gather: Convertir en formato largo
¿Cómo quedan los datos al aplicar otra función?
Todos los datos estan ahora en una sola columna, es decir, convertir los datos en formato largo, anteriormente estaban en formato ancho. Función gather reune los datos
library(dplyr)## Warning: package 'dplyr' was built under R version 4.2.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)## Warning: package 'tidyr' was built under R version 4.2.3# install.packages("datarium")
data("selfesteem", package = "datarium")
# Conversión
datos = selfesteem
datos = datos %>%
gather(key = "tiempo",
value = "rto",
t1, t2, t3) %>%
mutate_at(vars(id, tiempo), as.factor)
View(datos)Resumen estadístico de los datos
datos %>%
group_by(tiempo) %>%
summarise(media = mean (rto),
desv = sd(rto),
n = n(),
cv = 100*desv/media)## # A tibble: 3 × 5
## tiempo media desv n cv
## <fct> <dbl> <dbl> <int> <dbl>
## 1 t1 3.14 0.552 10 17.6
## 2 t2 4.93 0.863 10 17.5
## 3 t3 7.64 1.14 10 15.0\[H0 = \mu[t1] = \mu[t2] = \mu[t3]\]
Coeficientes de variación son menores al 20%, demostrando que los datos tienen un comportamiento normal, algo bueno en la agricultura. Si es mayor, puede ser por heterogeneidad, datos atípicos, posiblemente se deba imputar.
Visualización
boxplot(datos$rto ~ datos$tiempo)
#### ¿Tenemos datos atípicos?
library(rstatix)## Warning: package 'rstatix' was built under R version 4.2.3##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filterdatos %>%
group_by(tiempo) %>%
identify_outliers(rto)## # A tibble: 2 × 5
## tiempo id rto is.outlier is.extreme
## <fct> <fct> <dbl> <lgl> <lgl>
## 1 t1 6 2.05 TRUE FALSE
## 2 t2 2 6.91 TRUE FALSEExplicación: En el tiempo 1, parcela, así como en el tiempo 2 parcela 2, el rendimiento es 2.046 y 6.912 respectivamente, pero no son extremos (CV < 20%), es decir, son sospechosos de ser atípicos pero no son preocupantes.
¿Son datos normales? -> Prueba de normalidad
datos %>%
group_by(tiempo) %>%
shapiro_test(rto)## # A tibble: 3 × 4
## tiempo variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 t1 rto 0.967 0.859
## 2 t2 rto 0.876 0.117
## 3 t3 rto 0.923 0.380Normalidad, todos los datos son mayores al 5%, es decir, los datos tienen un comportamiento normal, en caso de que se presente un valor < 5%, se menciona que X dato presenta anormalidad posiblemente asociado a un dato atípico.
Supuesto de esfericidad
Comparar si la variabilidad entre tiempos sucesivos es mas o menos la misma.
res.aov <- anova_test(data = datos,
dv = rto,
wid = id,
within = tiempo)
# Esfericidad -> Prueba Mauchly's
res.aov$`Mauchly's Test for Sphericity`## Effect W p p<.05
## 1 tiempo 0.551 0.092get_anova_table(res.aov)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 tiempo 2 18 55.469 2.01e-08 * 0.829P value es menor al 5% por lo que se rechaza la hipótesis nula, es decir, que el aceite que se está produciendo en los tres tiempos no es el mismo, es decir, que el corte en diferentes tiempos si tiene efecto, el corte 3 parece ser el mejor.
Pruebas de comparación
datos %>%
pairwise_t_test(
rto ~ tiempo, paired = TRUE,
p.adjust.method = "bonferroni")## # A tibble: 3 × 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 rto t1 t2 10 10 -4.97 9 0.000772 2e-3 **
## 2 rto t1 t3 10 10 -13.2 9 0.000000334 1e-6 ****
## 3 rto t2 t3 10 10 -4.87 9 0.000886 3e-3 **P value -> Ajustados (Todos valores pequeños)
Todos los tiempos son diferentes, en este caso el tiempo 3 es la mejor opción.
Medidas repetidas de dos vías (Dos factores: Tiempo y fertilizante (ejemplo))
Seguir el mismo procedimiento
data("selfesteem2", package = "datarium")
datos2 = selfesteem2
datos2$treatment = gl(2,12,24, c('con fert', 'sin fert'))
View(datos2)Convertir en formato largo
datos2 = datos2 %>%
gather(key='tiempo',
value = 'rto',
t1,t2,t3)Resumen estadístico
datos2 %>%
group_by(treatment, tiempo) %>%
summarise(media = mean(rto),
desv = sd(rto),
n = n(),
cv = 100*desv/media)## `summarise()` has grouped output by 'treatment'. You can override using the
## `.groups` argument.## # A tibble: 6 × 6
## # Groups: treatment [2]
## treatment tiempo media desv n cv
## <fct> <chr> <dbl> <dbl> <int> <dbl>
## 1 con fert t1 88 8.08 12 9.18
## 2 con fert t2 83.8 10.2 12 12.2
## 3 con fert t3 78.7 10.5 12 13.4
## 4 sin fert t1 87.6 7.62 12 8.70
## 5 sin fert t2 87.8 7.42 12 8.45
## 6 sin fert t3 87.7 8.14 12 9.28Nuevamente el coeficiente de variación es menor al 20%.
Visualización
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.2.3ggplot(datos2)+
aes(tiempo, rto, fill=treatment)+
geom_boxplot()
Aumento de la fertilización en variabilidad
¿Tenemos datos atípicos?
datos2 %>%
group_by(treatment, tiempo) %>%
identify_outliers(rto)## [1] treatment tiempo id rto is.outlier is.extreme
## <0 rows> (or 0-length row.names)¿Son datos normales? -> Prueba de normalidad
datos2 %>%
group_by(treatment, tiempo) %>%
shapiro_test(rto)## # A tibble: 6 × 5
## treatment tiempo variable statistic p
## <fct> <chr> <chr> <dbl> <dbl>
## 1 con fert t1 rto 0.828 0.0200
## 2 con fert t2 rto 0.868 0.0618
## 3 con fert t3 rto 0.887 0.107
## 4 sin fert t1 rto 0.919 0.279
## 5 sin fert t2 rto 0.923 0.316
## 6 sin fert t3 rto 0.886 0.104ANOVA
res.aov <- anova_test(
data = datos2,
dv = rto,
wid = id,
within = c(treatment,
tiempo)
)
get_anova_table(res.aov)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 treatment 1.00 11.00 15.541 2.00e-03 * 0.059
## 2 tiempo 1.31 14.37 27.369 5.03e-05 * 0.049
## 3 treatment:tiempo 2.00 22.00 30.424 4.63e-07 * 0.050Si hay interacción. p < 5%, no se debe mirar ni el tratamiento ni el tiempo, se observa la interacción. No se aplican comparaciones porque hay interacción.
res.aov$`Mauchly's Test for Sphericity`## Effect W p p<.05
## 1 tiempo 0.469 0.023 *
## 2 treatment:tiempo 0.616 0.089get_anova_table(res.aov)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 treatment 1.00 11.00 15.541 2.00e-03 * 0.059
## 2 tiempo 1.31 14.37 27.369 5.03e-05 * 0.049
## 3 treatment:tiempo 2.00 22.00 30.424 4.63e-07 * 0.050¿Qué se recomienda hacer cuando interacción?
Un gráfico de interacción :)
interaction.plot(datos2$tiempo,
datos2$treatment,
datos2$rto)
datos2 %>%
group_by(tiempo, treatment) %>%
summarise(mean_rto = mean(rto)) %>%
ggplot()+
aes(tiempo, mean_rto,
color=treatment,
group=treatment)+
geom_point(size=5)+
geom_line(linewidth=3)## `summarise()` has grouped output by 'tiempo'. You can override using the
## `.groups` argument.
* Parece ser que el tratamiento sin fertilizante es mejor. * La
interaccón se denota con los cambios evidentes según el tratamiento