Recordatorio de lo que es una vía.

Factores intresujetos.

Factores entresujetos.

Por lo tanto.

\[\text{Tiempo}= \text{Factor - Intrasujeto} \\ \text{FSCA, FSBA, FCCA, FCBA}= \text{Factor - Entrasujeto}\]

Una vía.

Hipótesis.

\[H_0 = \mu_{t1} = \mu_{t2} = \mu_{t3}\]

# Solo hay tiempo y respuesta (Se desea saber qué es lo que sucede en ese cultivo a lo largo del tiempo).
# install.packages("datarium")
# id: Identificación de las parcelas.

library(datarium)

data("selfesteem", package = "datarium")

datos = selfesteem
head(datos, 3)
datos 
boxplot(datos[,-1])

# Construcción de los datos en formato largo (gather), pocas columnas.

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, union
library(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)
boxplot(datos)

Resumen estadísitico.

# CV>20%: Los datos son normales, son buenos para analizar, homogeneidad.

# CV<20%: Heterogeneidad, causante d eproblemas para analizar, se pueden imputar.

datos %>% 
  group_by(tiempo) %>%
  summarise(media = mean(rto),
            desv = sd(rto),
            n = n(),
            cv = 100*desv/media)
boxplot(datos$rto ~ datos$tiempo)

Revisar datos atípicos (Outliers).

library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
datos %>%
  group_by(tiempo) %>%
  identify_outliers(rto)
status=boxplot(datos)

a = status$out
a
## numeric(0)
  • Hay dos datos que no representan confiabilidad, se parecen mucho, pero no son preocupantes, hay que tener cuidado con los exremos, ya que si representan ser datos atípicos.

Análisis inferencial - Probando normalidad.

datos %>%
  group_by(tiempo) %>%
  shapiro_test(rto)
  • Los datos aparentan tener el comportamiento de una distribución normal, puesto a que todos presentan valores >5%, lo que indica que la presencia de datos atípicos es practicamente nula.

Supuesto de esfericidad o Test de Mauchly.

  • Solo se usa para medidas repetidas como la presencia de variabilidad en tiempos sucesivos.
res.aov <- anova_test(data = datos,
                      dv = rto,
                      wid = id,
                      within = tiempo)

# Esfericidad: p-valor > 5%, no se rechaza la H_0, por lo tanto, hay esfericidad.

res.aov$`Mauchly's Test for Sphericity`
get_anova_table(res.aov)

  • Se compara los diferentes tiempos sucesivos entre si, para poder ver su variabilidad (t1 - t2, t1 - t3, t2 - t3).

  • Si es >5% se puede asumir que las varianzas son iguales, por lo que habria que probar otra metodología de análisis.

  • Se rechazó la H_0: el aceite producido en los 3 diferentews tiempos es diferentes, el tiempo de corte tiene efecto.

Comparación entre tiempos aposteriori.

# Se usa el p-valor ajustado (p.adj)

datos %>%
  pairwise_t_test(
    rto ~ tiempo, paired = TRUE,
    p.adjust.method = "bonferroni")
# P-valor ajustado
# P-valor < 5%: Se rechaza H_0, todos son diferentes.

  • Todos los tiempos son diferentes en rendimiento.

  • El mejor es el t3.

  • Los peores son el t1 y t2.

Dos vías.

# Wide format
data("selfesteem2", package = "datarium")

datos2 = selfesteem2
#View(datos2)
datos2 = selfesteem2
datos2$treatment = gl(2,12,24, c('con fert', 'sin fert'))
datos2 = datos2 %>% 
  gather(key='tiempo', value = 'rto',
         t1,t2,t3)
datos2

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.
  • El mejor resultado fue el tiempo 1 con fertilización.

Visualización - análisis descriptivo.

library(ggplot2)

ggplot(datos2)+
  aes(tiempo, rto, fill=treatment)+
  geom_boxplot()

  • Se esperaria que fertilizar fuera significativamente mejor que no fertilizar, sin embargo, vemos que no hay mucha variabilidad entre fertilizar y no fertilizar, esto se puede deber a diversos factores que hayan afectado la efectividad del fertilizante, como lavado por lluvias, vencimiento del principio activo, evaporación por demasiada radiación solar, entre otras. No significa que fertilizar y no fertilizar dé igual, sino que tuvo que pasar algo para que la fertilización no tuviera los resultados afectados.

Prueba de datos atípicos.

datos2 %>%
  group_by(treatment, tiempo) %>%
  identify_outliers(rto)
datos2 %>%
  group_by(treatment, tiempo) %>%
  shapiro_test(rto)

Revisión de p-valor, Varianza.

res.aov <- anova_test(
  data = datos2,
  dv = rto,
  wid = id,
  within = c(treatment,
             tiempo)
  )
get_anova_table(res.aov)
  • Si hay interacción, ya que el p-valor es <5%.
datos2 %>%
  group_by(treatment, tiempo) %>%
  identify_outliers(rto)
datos2 %>%
  group_by(treatment, tiempo) %>%
  shapiro_test(rto)

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.

  • En el tiempo 1 no hay interacción (interacción muy baja, casi nula), y en el tiempo 3 hay una interacción muy amplia (espacio entre una linea y otra).

Tres vías.

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
table(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  12

Descriptivo tabular.

datos3 %>%
  group_by(FQ, FO, tiempo) %>%
  get_summary_stats(rto, type = "mean_sd")

Descriptivo gráfico.

library(ggpubr)
ggboxplot(
  datos3, x = "FQ", y = "rto",
  color = "tiempo", palette = "jco",
  facet.by = "FO", short.panel.labs = FALSE
  )

Identificación de datos atípicos.

datos3 %>% 
  group_by(FQ,FO, tiempo) %>% 
  identify_outliers(rto)

  • Dato extremo puede afectar el rendimiento del experimento, en este caso NO HAY.

  • En caso de hubiera extremos (columna “is.extreme” en “true”), solo hace falta una nota que mencione el dato con el extremo ya que ese dato puede afectar la estabilidad de todos los resultados.

Prueba de normalidad.

datos3
datos3 %>% 
  group_by(FQ,FO,tiempo) %>% 
  shapiro_test(rto)
ggplot(datos3)+
  aes(rto, fill=tiempo)+
  geom_density(alpha=0.5)

  • Todas las densidades parecen tener distribución normal.
# Global
ggplot(datos3)+
  aes(rto)+
  geom_density(alpha=0.5)

  • Se pruebas los diferentes test.
shapiro_test(datos3$rto)
library(nortest)
lillie.test(datos3$rto)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  datos3$rto
## D = 0.10013, p-value = 0.001229
cvm.test(datos3$rto)
## 
##  Cramer-von Mises normality test
## 
## data:  datos3$rto
## W = 0.23031, p-value = 0.002205
sf.test(datos3$rto)
## 
##  Shapiro-Francia normality test
## 
## data:  datos3$rto
## W = 0.96406, p-value = 0.001295
ad.test(datos3$rto)
## 
##  Anderson-Darling normality test
## 
## data:  datos3$rto
## A = 1.432, p-value = 0.001021
pearson.test((datos3$rto))
## 
##  Pearson chi-square normality test
## 
## data:  (datos3$rto)
## P = 28.5, p-value = 0.004672
  • Siempre se dice que cuando el p-valor es <5%, rechaza, y cuando es >5%, no rechaza. Cuando se hicieron las pruebas de normalidad, TODAS dieron <5%, por ende, todas rechazan, todas dicen que no son normales, sin embargo, si se cambia el punto de corte y se establece el 1% en lugar de 5%, en este caso, la prueba simetría cayó en el 3%, como 3% > 1%, entonces no se rechaza, por lo que la conclusión es que los datos son simétricos, aunque no sean normales, es buena señal que tenga cierto grado de simetría, con simetría se puede conformar así no haya normalidad.

Test de simetría.

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 
##                  19

Análisis de varianza - inferencial.

res.aov <- anova_test(
  data = datos3, dv = rto, wid = id,
  within = c(FQ, FO, tiempo)
  )
get_anova_table(res.aov)
  • Hay triple interacción, no hay un tratamiento predominante, significa que los efectos anteriores son “mentirosos”, en este caso, quiere decir que una fertilización es efectiva en un tiempo y luego en otro tiempo ya no lo es, lo cual, no tiene sentido, algunos son buenos, parcialmente sirven en un momento y luego no.

Gráficos de interración.

interaction.plot(datos3$tiempo,
                 datos3$FO,
                 datos3$rto)

interaction.plot(datos3$tiempo,
                 datos3$FQ,
                 datos3$rto)

tbl=tapply(datos3$rto, list(datos3$FQ,
                            datos3$FO,
                            datos3$tiempo),
       mean)
       
addmargins(tbl, 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.68132

Referencia.