Bibliotecas

library(tidyverse)
library(effectsize)
library(readxl)
library(janitor)
library(broom)
library(emmeans)
library(car)

Datos

datos <- read_excel("experimento_ratones.xlsx") 
datos

Tamaño del efecto

Modelos DCA y DBCA

modelo1 <- aov(GST ~ Treatment, data = datos)
modelo2 <- aov(GST ~ Treatment + Block, data = datos)

Anova de dos vías aditivo

  • Hipótesis factor 1:

\[H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4 \\ H_1: \mu_i \neq \mu_j\]

  • Hipótesis factor 2:

\[H_0: \mu_1 = \mu_2 \\ H_1: \mu_1 \neq \mu_1\]

Exploración

datos %>% 
  ggplot(aes(x = Strain, y = GST, color = Treatment, fill = Treatment)) +
  geom_boxplot(alpha = 0.5)

Modelo

modelo3 <- aov(GST ~ Treatment + Strain, data = datos)
modelo3 %>% tidy()

Medias estimadas

  • Tratamiento:
modelo3 %>% 
  emmeans(specs = "Treatment") %>% 
  as.data.frame()
  • Raza:
modelo3 %>% 
  emmeans(specs = "Strain") %>% 
  as.data.frame()

Tamaño del efecto

modelo3 %>% 
  eta_squared()

Residuales

par(mfrow = c(2, 2))
plot(modelo3)

  • Prueba de normalidad:
modelo3 %>% 
  residuals() %>% 
  shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.98395, p-value = 0.9872
  • Prueba de homocedasticidad tratamiento:
datos %>% 
  bartlett.test(data = ., GST ~ Treatment)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  GST by Treatment
## Bartlett's K-squared = 0.00092761, df = 1, p-value = 0.9757
  • Prueba de homocedasticidad raza:
datos %>% 
  bartlett.test(data = ., GST ~ Strain)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  GST by Strain
## Bartlett's K-squared = 1.4085, df = 3, p-value = 0.7035

Anova de dos vías multiplicativo

  • Hipótesis factor 1:

\[H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4 \\ H_1: \mu_i \neq \mu_j\]

  • Hipótesis factor 2:

\[H_0: \mu_1 = \mu_2 \\ H_1: \mu_1 \neq \mu_1\]

  • Hipótesis interacción:

\[H_0: (\alpha\tau)_{11} = \cdots = (\alpha\beta)_{42} \\ H_1: (\alpha\tau)_{ij} \neq (\alpha\tau)_{ij}\]

Exploración

datos %>% 
  group_by(Strain, Treatment) %>% 
  summarise(promedio = mean(GST)) %>% 
  ggplot(aes(x = Strain, y = promedio, color = Treatment)) +
  geom_point() +
  geom_line(aes(group = Treatment))

Modelo

modelo4 <- aov(GST ~ Treatment * Strain, data = datos)
modelo4 %>% tidy()

Medias estimadas

  • Tratamiento:
modelo4 %>% 
  emmeans(specs = "Treatment") %>% 
  as.data.frame()
  • Raza:
modelo4 %>% 
  emmeans(specs = "Strain") %>% 
  as.data.frame()
  • Interacción raza-tratamiento:
modelo4 %>% 
  emmeans(specs = pairwise ~ Treatment | Strain) %>% 
  as.data.frame()

Tamaño del efecto

modelo4 %>% 
  eta_squared()

Residuales

par(mfrow = c(2, 2))
plot(modelo4)

  • Prueba de normalidad:
modelo4 %>% 
  residuals() %>% 
  shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.91331, p-value = 0.1316
  • Prueba de homocedasticidad tratamiento:
datos %>% 
  bartlett.test(data = ., GST ~ Treatment)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  GST by Treatment
## Bartlett's K-squared = 0.00092761, df = 1, p-value = 0.9757
  • Prueba de homocedasticidad raza:
datos %>% 
  bartlett.test(data = ., GST ~ Strain)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  GST by Strain
## Bartlett's K-squared = 1.4085, df = 3, p-value = 0.7035
  • Prueba de homocedasticidad interacción:
datos %>% 
  leveneTest(data = ., GST ~ Treatment * Strain)

Dos modelos más

Modelo aditivo con bloque

  • Hipótesis factor 1:

\[H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4 \\ H_1: \mu_i \neq \mu_j\]

  • Hipótesis factor 2:

\[H_0: \mu_1 = \mu_2 \\ H_1: \mu_1 \neq \mu_1\]

  • Hipótesis de verificación:

\[ H_0: \beta_1 = \beta_2 \\ H_A: \beta_1 \neq \beta_2 \]

modelo5 <- aov(GST ~ Treatment + Strain + Block, data = datos)
modelo5 %>% tidy()

Modelo multiplicativo con bloque

  • Hipótesis factor 1:

\[H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4 \\ H_1: \mu_i \neq \mu_j\]

  • Hipótesis factor 2:

\[H_0: \mu_1 = \mu_2 \\ H_1: \mu_1 \neq \mu_1\]

  • Hipótesis interacción:

\[H_0: (\alpha\tau)_{11} = \cdots = (\alpha\beta)_{42} \\ H_1: (\alpha\tau)_{ij} \neq (\alpha\tau)_{ij}\]

  • Hipótesis de verificación:

\[ H_0: \beta_1 = \beta_2 \\ H_A: \beta_1 \neq \beta_2 \]

modelo6 <- aov(GST ~ Treatment * Strain + Block, data = datos)
modelo6 %>% tidy()

Comparación de modelos

AIC

AIC(modelo1, modelo2, modelo3, modelo4, modelo5, modelo6)

BIC

BIC(modelo1, modelo2, modelo3, modelo4, modelo5, modelo6)

Mejor modelo

Medias estimadas

  • Tratamiento:
modelo6 %>% 
  emmeans(specs = "Treatment") %>% 
  as.data.frame()
  • Raza:
modelo6 %>% 
  emmeans(specs = "Strain") %>% 
  as.data.frame()
  • Interacción raza-tratamiento:
modelo6 %>% 
  emmeans(specs = pairwise ~ Treatment | Strain) %>% 
  as.data.frame()

Comparaciones múltiples

  • Tabla:
TukeyHSD(modelo6, which = "Treatment:Strain") %>% 
  tidy()
  • Gráfico:
TukeyHSD(modelo6, which = "Treatment:Strain") %>% 
  tidy() %>% 
  ggplot(aes(x = fct_reorder(contrast, estimate), y = estimate,
             ymin = conf.low, ymax = conf.high)) +
  geom_point() +
  geom_errorbar() +
  geom_hline(yintercept = 0, lty = 2, color = "red") +
  coord_flip()

Tamaño del efecto

modelo6 %>% 
  eta_squared()

Residuales

par(mfrow = c(2, 2))
plot(modelo6)

  • Prueba de normalidad:
modelo6 %>% 
  residuals() %>% 
  shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.92236, p-value = 0.1841
  • Prueba de homocedasticidad tratamiento:
datos %>% 
  bartlett.test(data = ., GST ~ Treatment)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  GST by Treatment
## Bartlett's K-squared = 0.00092761, df = 1, p-value = 0.9757
  • Prueba de homocedasticidad raza:
datos %>% 
  bartlett.test(data = ., GST ~ Strain)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  GST by Strain
## Bartlett's K-squared = 1.4085, df = 3, p-value = 0.7035
  • Prueba de homocedasticidad bloque:
datos %>% 
  bartlett.test(data = ., GST ~ Block)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  GST by Block
## Bartlett's K-squared = 0.14923, df = 1, p-value = 0.6993
  • Prueba de homocedasticidad interacción:
datos %>% 
  leveneTest(data = ., GST ~ Treatment * Strain)