Diferencias estadísticas
Análisis estadísticos básicos
Como tenemos varios días de muestreo (1 al 6), primero veamos la estadística básica para el último día de muestreo, el día 6, creemos un subset de datos para solo con los datos de ese día y salvémoslo como un objeto llamado dia6:
> dia6 <- subset(x = df, Dia == 6) # subset con los datos para el dia requerido
Nota. todo lo que está después del signo # no es considerado, son solo comentarios opcionales para describir lo que estamos haciendo.
Obtengamos ahora un resumen estadístico de los datos de Peso para el día 6:
> dia6 %>% group_by(Tratamiento) %>% get_summary_stats(Peso, type = "common") # Stats para el dia ingresado
Veremos que genera datos de estadística básica para los pesos por tratamiento para el día 6:
# A tibble: 3 × 11
Tratamiento variable n min max median iqr mean sd
<chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Control Peso 60 123 285 167 29 171. 27.9
2 T1 Peso 60 105 291 222. 65 217. 43.1
3 T2 Peso 60 4 55 13 11 15.8 11.0
# ℹ 2 more variables: se <dbl>, ci <dbl>
Podemos hacer lo mismo para los datos de Talla:
> dia6 %>% group_by(Tratamiento) %>% get_summary_stats(Talla, type = "common")
# A tibble: 3 × 11
Tratamiento variable n min max median iqr mean sd
<chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Control Talla 60 16 21 18 3 18.4 1.5
2 T1 Talla 60 16 22 20 1 19.4 1.48
3 T2 Talla 60 5 11 7 2 7.2 1.60
# ℹ 2 more variables: se <dbl>, ci <dbl>
Cálculo automático del último día
Para no tener que estar viendo cada vez cual fue el último día de tratamiento, podemos hacer obtenerlo automáticamente y salvarlo como un objeto:
> LastDay <- tail(dataR$Dia, 1)
Asumimos que la tabla está ordenada ascendentemente por el día y la última linea tendrá el valor del último día, por eso usamos tail para cortar el último valor. Así podremos usar este objeto para hacer nuestro subset de datos de ese día y salvarlo en un nuevo objeto llamado Data.Ultimo.Dia:
> Data.Ultimo.Dia <- subset(x = df, Dia == LastDay)
Podemos usar este objeto para todos los análisis posteriores.
Diferencias entre Unidades Experimentales
Veamos primero si tenemos diferencias significativas entre nuestras unidades experimentales en el último día de nuestro experimento; haremos un nuevo objeto llamado dif.kruskal.ExpUnits para el día 6 (dia6, objeto creado ya anteriormente) y calculando diferencias con el test de Kruskal-Wallis (kruskal_test) para pesos por unidad experimental y calcular la significancia:
> dif.kruskal.ExpUnits <- dia6 %>% group_by(Tratamiento) %>% kruskal_test(Peso ~ UExp) %>% adjust_pvalue(method = "bonferroni") %>% add_significance("p.adj")
Si vemos el resultado, podremos ver en una tabla los mismos:
> view(dif.kruskal.ExpUnits)
Para los pesos al día 6 tenemos diferencias significativas entre las unidades experimentales.
Diferencias entre tratamientos
Pesos (mg)
Vamos ahora a calcular diferencias entre el peso en los tratamientos:
> res.kruskal <- dia6 %>% kruskal_test(Peso ~ Tratamiento)
Si hacemos view(res.kruskal) veremos que hay una diferencia altamente significativa (p = 4.61 x 10-30), pero no sabemos entre que tratamientos hay diferencias, calculemos ahora el estadístico de Dunn:
> dia6 %>% dunn_test(Peso ~ Tratamiento, p.adjust.method = "bonferroni")
y obtendremos que hay diferencias altamente significativas (p.adj.) entre todos los tratamientos:
# A tibble: 3 × 9
.y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
1 Peso Control T1 60 60 3.97 7.28e- 5 2.18e- 4 ***
2 Peso Control T2 60 60 -7.48 7.54e-14 2.26e-13 ****
3 Peso T1 T2 60 60 -11.4 2.49e-30 7.47e-30 ****
Metamos este cálculo a un objeto llamado pwcmg:
> pwcmg <- dia6 %>% dunn_test(Peso ~ Tratamiento, p.adjust.method = "bonferroni")
y agregémosle al objeto pwcmg una posición para que en un gráfico aparezcan las diferencias estadísticas:
> pwcmg <- pwcmg %>% add_xy_position(x = "Tratamiento")
Grafiquemos
> ggboxplot(dia6, x = "Tratamiento", y = "Peso") + stat_pvalue_manual(pwcmg, hide.ns = TRUE) + labs(subtitle = get_test_label(res.kruskal, detailed = TRUE),caption = get_pwc_label(pwcmg))
Tallas (mm)
Todo lo anterior fueron cálculos para pesos pero también lo podemos hacer para las tallas:
> pwcmm <- dia6 %>% dunn_test(Talla ~ Tratamiento, p.adjust.method = "bonferroni")
> pwcmm <- pwcmm %>% add_xy_position(x = "Tratamiento")
Grafiquemos
> ggboxplot(dia6, x = "Tratamiento", y = "Talla") + stat_pvalue_manual(pwcmm, hide.ns = TRUE) + labs(subtitle = get_test_label(res.kruskal, detailed = TRUE),caption = get_pwc_label(pwcmm))
Podemos ver que no hay diferencias significativas entre las tallas del Control y el T1 pero entre los demás para el día 6.