Clase 18 y 20
Maria Jose Duarte
2023-06-12
13, 16 quiz, trabajo del diseño de hoy
DISEÑO EN MEDIDAS REPETIDAS
https://www.datanovia.com/en/lessons/repeated-measures-anova-in-r/
El tiempo juega un papel en el modelo, reconocer qué pasa cuando se tienen mediciones en el tiempo, hacer evaluaciones temporales.
Factor intrasujetos = tiempo
Factores entrasujetos = otros factores (al azar, el bloques, simples)
Antes de comenzar, paquetes y funciones importantes:
rstatix: para analisis estadístico.
datarium: contiene los datos requeridos para este caso.
anova_test: dv = variable respuesta, within = factor intrasujetos (tiempo), get_anova_table()
Tres casos de este diseño
MEDIDAS REPETIDAS DE UNA VÍA
El único factor involucrado es el tiempo, se recomienda que sean equidistantes.
Habrá una sola respuesta
Ejemplo: variable respuesta Aceite de limonaria.
#Preparación de los datos
library(dplyr)##
## 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)
# install.packages("datarium")
data("selfesteem", package = "datarium")
datos = selfesteem
datos = datos %>%
gather(key = "tiempo",
value = "rto",
t1, t2, t3) %>%
mutate_at(vars(id, tiempo), as.factor)
View(datos)
#La función gather junta los datos de rendimiento en una sola columna, se convierte formato ancho en largo
#id y tiempo se debe indicar que son factores con as.factor.
# %>% se utiliza para ordenar las funciones en lugar de agruparlas, cual va primero. \[H_0: \mu_{t1} = \mu_{t2} = \mu_{t3}\]
Resumen estadístico de los datos
#vamos a calcular media, número de los datos, desviación
#n: número de datos, cv: coeficientes de variación
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.0Todos los coeficientes de variación <20%, lo que nos dice que se puede confiar en los datos.
Visualización por diagrama de caja
boxplot(datos$rto~ datos$tiempo)
Detección de datos atípicos
library(rstatix)##
## 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 FALSENo se detectan EXTREMOS, no son preocupantemente atípicos de los 30 datos que tenemos 2 son sospechosos.
(is.extreme) da falso.
Normalidad de datos
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.380Observamos pvalores, todos son mayores al 5%, los datos de cada tiempo son normales.
Si no se tiene comportamiento normal, solo anotamos que no lo son
Analisis de varianza
res.aov <- anova_test(data = datos,
dv =rto,
wid = id,
within = tiempo)
get_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.829ges: tamaño del efecto
Solo observo el pvalor, rechazo la H0
Corte 3 mayor rendimiento de aceite
Supuesto de esfericidad
Test de Mauchly
Igualdad de varianzas, tenemos tiempos que son muy distantes t1 -t2. E ste compara si la variabilidad de los datos en tiempos SUCEVOS es igual
res.aov$`Mauchly's Test for Sphericity`## Effect W p p<.05
## 1 tiempo 0.551 0.092Puedo asumir variables iguales
La varianza en tiempos sucesivos no son iguales, se debe buscar otra metodología
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 **Observo pvalores ajustados, observo que todos son menores al 5% Todos los tiempos son diferentes, ya que este observa la relación entre cada tiempo.
MEDIDAS REPETIDAS DE DOS VÍAS
Como tenemos dos factores, se debe analizar la interacción de los factores.
Se añade un factor al tiempo, por ejemplo: variedad, lugar, procedencia de la semilla, etc
Para este ejemplo: fertilizante.
data("selfesteem2", package = "datarium")
datos2 = selfesteem2
datos2$treatment = gl(2,12,24, c('con fert', 'sin fert'))
View(datos2)
#Un control es sin fert y otro con fert#Vamos a convertirlos de formato ancho a largo
datos2 = datos2 %>%
gather(key='tiempo', value = 'rto',
t1,t2,t3)
#De nuevo usamos la función gatherResumen estadístico de los datos
#vamos a calcular media, número de los datos, desviación
#n: número de datos, cv: coeficientes de variación
datos2 %>%
group_by(treatment, tiempo) %>%
summarise(media = mean(rto),
mediana = median(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 × 7
## # Groups: treatment [2]
## treatment tiempo media mediana desv n cv
## <fct> <chr> <dbl> <dbl> <dbl> <int> <dbl>
## 1 con fert t1 88 92 8.08 12 9.18
## 2 con fert t2 83.8 88 10.2 12 12.2
## 3 con fert t3 78.7 81 10.5 12 13.4
## 4 sin fert t1 87.6 90 7.62 12 8.70
## 5 sin fert t2 87.8 90 7.42 12 8.45
## 6 sin fert t3 87.7 89.5 8.14 12 9.28Todos los cv son menores al 20% , el mejor es media = 88, con fert en t1.
Visualización por diagrama de caja
library(ggplot2)
ggplot(datos2)+
aes(tiempo, rto, fill=treatment)+
geom_boxplot()
Sin fertilizantes son mejores para t2 y t3 que con fertilizantes.
Detección de 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)No hay atipicos, 0
Supuesto de normalidad
En la pag se hace de esta forma, normalmente se extraen los residuales pero esta función no permite extrear los residuales.
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.104Analisis de varianza
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.050Primero observamos que SI hay interacción por el pvalor.
Procedimiento con interacción
Vamos a hacer un grafico de interacción.
interaction.plot(datos2$tiempo,
datos2$treatment,
datos2$rto)
#Otra forma.
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.
Lo que causa la interacción es el con fert, porque si presenta cambios a lo largo del tiempo mientras que el sin fert no-
Mejor tratamiento es t1 con fert.
Prueba de estericidad
res.aov$`Mauchly's Test for Sphericity`## Effect W p p<.05
## 1 tiempo 0.469 0.023 *
## 2 treatment:tiempo 0.616 0.089MEDIDAS REPETIDAS DE TRES VÍAS
FO: fertilizante organico
FQ: fertilizante químico
library(rstatix)
library(datarium)
set.seed(123)
data("weightloss", package = "datarium")
datos3 = weightloss %>%
rename(FQ = diet,
FO = exercises) %>%
pivot_longer(c(t1, t2, t3),
names_to = "tiempo",
values_to = "rto")
datos3## # A tibble: 144 × 5
## id FQ FO tiempo rto
## <fct> <fct> <fct> <chr> <dbl>
## 1 1 no no t1 10.4
## 2 1 no no t2 13.2
## 3 1 no no t3 11.6
## 4 2 no no t1 11.6
## 5 2 no no t2 10.7
## 6 2 no no t3 13.2
## 7 3 no no t1 11.4
## 8 3 no no t2 11.1
## 9 3 no no t3 11.4
## 10 4 no no t1 11.1
## # ℹ 134 more rowstable(datos3$FQ, datos3$FO, datos3$tiempo) ## , , = t1
##
##
## no yes
## no 12 12
## yes 12 12
##
## , , = t2
##
##
## no yes
## no 12 12
## yes 12 12
##
## , , = t3
##
##
## no yes
## no 12 12
## yes 12 12Resumen estadístico de los datos
datos3 %>%
group_by(FQ, FO, tiempo) %>%
get_summary_stats(rto, type = "mean_sd")## # A tibble: 12 × 7
## FQ FO tiempo variable n mean sd
## <fct> <fct> <chr> <fct> <dbl> <dbl> <dbl>
## 1 no no t1 rto 12 10.9 0.868
## 2 no no t2 rto 12 11.6 1.30
## 3 no no t3 rto 12 11.4 0.935
## 4 no yes t1 rto 12 10.8 1.27
## 5 no yes t2 rto 12 13.4 1.01
## 6 no yes t3 rto 12 16.8 1.53
## 7 yes no t1 rto 12 11.7 0.938
## 8 yes no t2 rto 12 12.4 1.42
## 9 yes no t3 rto 12 13.8 1.43
## 10 yes yes t1 rto 12 11.4 1.09
## 11 yes yes t2 rto 12 13.2 1.22
## 12 yes yes t3 rto 12 14.7 0.625#Para obtene un resumen estadístico de la variable respuesta Observamos que la media mayor es de 14.7 a partir de utilizar FQ + FO. en el t3.
El valor menor es 10.8 solo utilizar FO en el t1.
Visualización
library(ggpubr)
bxp <- ggboxplot(
datos3, x = "FQ", y = "rto",
color = "tiempo", palette = "jco",
facet.by = "FO",
short.panel.labs = FALSE
)
bxp
La caja mas alta es t3 FO si FQ no, los tratamientos con rto mas altos son el t3 (grises).
El uso de fertilizante organico parece ser mejor.
Detección de datos atipicos
datos3 %>%
group_by(FQ, FO, tiempo) %>%
identify_outliers(rto)## # A tibble: 5 × 7
## FQ FO tiempo id rto is.outlier is.extreme
## <fct> <fct> <chr> <fct> <dbl> <lgl> <lgl>
## 1 no no t3 2 13.2 TRUE FALSE
## 2 yes no t1 1 10.2 TRUE FALSE
## 3 yes no t1 3 13.2 TRUE FALSE
## 4 yes no t1 4 10.2 TRUE FALSE
## 5 yes no t2 10 15.3 TRUE FALSENos interesa la columna que dice IS.EXTREME. Concluimos que no hay ningún extremo.
En el caso de que se obtengan añadimos una nota: este datos puede afectar la estabilidad de los resultados
Supuesto de normalidad
Por Shapiro test, sin embargo vamos a usar otro método por grafico para observar de forma general la distribución de los datos.
Con las otras pruebas tampoco obtenemos normalidad, usualmente la normalidad se realiza con los residuales, la función que utuliza la pag web no permite extraerlos
datos3 %>%
group_by(FQ, FO, tiempo) %>%
shapiro_test(rto)## # A tibble: 12 × 6
## FQ FO tiempo variable statistic p
## <fct> <fct> <chr> <chr> <dbl> <dbl>
## 1 no no t1 rto 0.917 0.264
## 2 no no t2 rto 0.957 0.743
## 3 no no t3 rto 0.965 0.851
## 4 no yes t1 rto 0.922 0.306
## 5 no yes t2 rto 0.912 0.229
## 6 no yes t3 rto 0.953 0.674
## 7 yes no t1 rto 0.942 0.528
## 8 yes no t2 rto 0.982 0.989
## 9 yes no t3 rto 0.931 0.395
## 10 yes yes t1 rto 0.914 0.241
## 11 yes yes t2 rto 0.947 0.598
## 12 yes yes t3 rto 0.937 0.464#Global
ggplot(datos3)+
aes(rto)+
geom_density(alpha=0.5)
#Separados por tiempo
ggplot(datos3)+
aes(rto, fill = tiempo)+
geom_density(alpha=0.5)
#Vamos a hacer solo 1 prueba de normalidad
shapiro.test(datos3$rto)##
## Shapiro-Wilk normality test
##
## data: datos3$rto
## W = 0.96377, p-value = 0.0007474#Otras formas de hacer pruebas de normalidad
library(nortest)
lillie.test(datos3$rto)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: datos3$rto
## D = 0.10013, p-value = 0.001229cvm.test(datos3$rto)##
## Cramer-von Mises normality test
##
## data: datos3$rto
## W = 0.23031, p-value = 0.002205sf.test(datos3$rto)##
## Shapiro-Francia normality test
##
## data: datos3$rto
## W = 0.96406, p-value = 0.001295#Prueba de simetria, vamos a asegurar simetria ya que no tenemos normalidad.
library(lawstat)
symmetry.test(datos3$rto)##
## m-out-of-n bootstrap symmetry test by Miao, Gel, and Gastwirth (2006)
##
## data: datos3$rto
## Test statistic = 2.2399, p-value = 0.032
## alternative hypothesis: the distribution is asymmetric.
## sample estimates:
## bootstrap optimal m
## 19Observamos tanto para el grafico global como separados por tiempo distribución relativamente “normal” de los datos, un poco asimetrica a la derecha.
Para la prueba de simetria obtenemos pvalor del 3%, que en este caso lo comparamos con el 1%, por lo que asumimos simetría. 3% mayor que 1%.
Analisis de varianza
res.aov <- anova_test(
data = datos3, dv = rto, wid = id,
within = c(FQ, FO, tiempo)
)
get_anova_table(res.aov)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 FQ 1.00 11.00 6.021 3.20e-02 * 0.028
## 2 FO 1.00 11.00 58.928 9.65e-06 * 0.284
## 3 tiempo 2.00 22.00 110.942 3.22e-12 * 0.541
## 4 FQ:FO 1.00 11.00 75.356 2.98e-06 * 0.157
## 5 FQ:tiempo 1.38 15.17 0.603 5.01e-01 0.013
## 6 FO:tiempo 2.00 22.00 20.826 8.41e-06 * 0.274
## 7 FQ:FO:tiempo 2.00 22.00 14.246 1.07e-04 * 0.147#siendo dv = variable dependiente. Miramos el pvalor, empezamos por la triple interaccion: tiempo-factores.
Si tenemos interacción no se puede interpretar de ahí hacia arriba, por tanto hacemos un grafico de interacción.
En presencia de interacción no hacemos comparaciones de media.
interaction.plot(datos3$tiempo,
datos3$FQ,
datos3$rto)
interaction.plot(datos3$tiempo,
datos3$FO,
datos3$rto)
Tabla de medias
tabla = tapply(datos3$rto, list(datos3$FQ,
datos3$FO,
datos3$tiempo),
mean)
addmargins(tabla, FUN = mean)## Margins computed over dimensions
## in the following order:
## 1:
## 2:
## 3:## , , t1
##
## no yes mean
## no 10.90917 10.79417 10.85167
## yes 11.74250 11.39333 11.56792
## mean 11.32583 11.09375 11.20979
##
## , , t2
##
## no yes mean
## no 11.56583 13.42083 12.49333
## yes 12.41583 13.22500 12.82042
## mean 11.99083 13.32292 12.65688
##
## , , t3
##
## no yes mean
## no 11.45000 16.81750 14.13375
## yes 13.78667 14.65500 14.22083
## mean 12.61833 15.73625 14.17729
##
## , , mean
##
## no yes mean
## no 11.30833 13.67750 12.49292
## yes 12.64833 13.09111 12.86972
## mean 11.97833 13.38431 12.68132Unir la fila y la columna de los margenes no se obtiene lo mismo que los valores de la tabla.
Al cruzar vertical y horizontalmente el valor de media mas alta es diferente, esto es señal de interacción, por tanto los tiempos 2 y 3 son los que están causando la interacción.
En el caso de que se presenten otras interacciones
Si no se presenta triple interacción pero si, por ejemplo, de otros dos factores, F1 y T, por tanto F2 SI se puede interpretar. Por tanto hacemos grafico de interacción para concluir sobre F1. Para F2 tukey o prueba de comparación de medias.
Primer parcial: 35%