Bibliotecas

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

Datos

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

Tamaño del efecto

Artículo guía

Anova de una vía (DCA)

Hipótesis

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

Nivel de significancia

Tomaremos un nivel de significancia del 1 % para los ejemplos de este documento.

Resumen descriptivo

datos %>% 
  group_by(Strain) %>% 
  summarise(promedio = mean(GST),
            desviacion = sd(GST),
            mediana = median(GST),
            minimo = min(GST),
            maximo = max(GST))

Exploración

Gráfico 1

datos %>% 
  ggplot(aes(x = Strain, y = GST)) +
  geom_boxplot()

Gráfico 2

datos %>% 
  group_by(Strain) %>% 
  summarise(promedio = mean(GST),
            mediana = median(GST),
            desviacion = sd(GST)) %>% 
  ggplot(aes(x = Strain, y = promedio,
             ymin = promedio - desviacion,
             ymax = promedio + desviacion)) +
  geom_pointrange() +
  geom_point(aes(y = mediana), color = "red")

Modelo

modelo_dca <- aov(GST ~ Strain, data = datos)
modelo_dca %>% tidy()

Suma de cuadrados total (SCT):

422076.50 + 28613.25

## [1] 450689.8

Se calcula de la siguiente manera:

media_total <- mean(datos$GST)

sc_total <- sum((datos$GST - media_total)^2)
sc_total

## [1] 450689.8

Ahora calculamos la suma de cuadrados entre tratamientos:

r_nih <- datos %>% filter(Strain == "NIH")
r_balbc <- datos %>% filter(Strain == "BALB/C")
r_aj <- datos %>% filter(Strain == "A/J")
r_129 <- datos %>% filter(Strain == "129/Ola")

sc_nih <- ((mean(r_nih$GST) - media_total)^2) * 4
sc_balbc <- ((mean(r_balbc$GST) - media_total)^2) * 4
sc_aj <- ((mean(r_aj$GST) - media_total)^2) * 4
sc_129 <- ((mean(r_129$GST) - media_total)^2) * 4

sc_ttos <- sc_nih + sc_balbc + sc_aj + sc_129
sc_ttos

## [1] 28613.25

Finalmente calculamos la suma de cuadrados del residual:

sc_total - sc_ttos

## [1] 422076.5

Medias estimadas

modelo_dca %>% 
  emmeans(specs = "Strain") %>% 
  as.data.frame()

Comparaciones múltiples

t-student

pairwise.t.test(x = datos$GST,
                g = datos$Strain)

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  datos$GST and datos$Strain 
## 
##        129/Ola A/J BALB/C
## A/J    1       -   -     
## BALB/C 1       1   -     
## NIH    1       1   1     
## 
## P value adjustment method: holm

Tukey

modelo_dca %>% 
  TukeyHSD()

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = GST ~ Strain, data = datos)
## 
## $Strain
##                   diff       lwr      upr     p adj
## A/J-129/Ola      84.25 -309.4684 477.9684 0.9185996
## BALB/C-129/Ola  -30.50 -424.2184 363.2184 0.9954768
## NIH-129/Ola      28.75 -364.9684 422.4684 0.9962012
## BALB/C-A/J     -114.75 -508.4684 278.9684 0.8223712
## NIH-A/J         -55.50 -449.2184 338.2184 0.9742353
## NIH-BALB/C       59.25 -334.4684 452.9684 0.9689979

Tamaño del efecto

Eta cuadrado:

modelo_dca %>% 
  eta_squared()

F-Cohen:

modelo_dca %>% 
  cohens_f()

Residuales

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

Prueba de normalidad:

modelo_dca %>% 
  residuals() %>% 
  shapiro.test()

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.96966, p-value = 0.8331

Prueba de homocedasticidad:

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

Diseño en bloques (DBA)

Hipótesis principal:

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

Hipótesis de verificación:

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

Conteo por combinación

datos %>%
  count(Block, Treatment)

Exploración

datos %>% 
  ggplot(aes(x = Treatment, y = GST)) +
  facet_wrap(~Block) +
  geom_boxplot()

Modelo

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

Medias estimadas

modelo_dbca %>% 
  emmeans(specs = "Treatment") %>% 
  as.data.frame()

Comparaciones múltiples

t-student

pairwise.t.test(x = datos$GST,
                g = datos$Treatment)

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  datos$GST and datos$Treatment 
## 
##         Control
## Treated 0.002  
## 
## P value adjustment method: holm

Tukey

modelo_dbca %>% 
  TukeyHSD()

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = GST ~ Treatment + Block, data = datos)
## 
## $Treatment
##                  diff      lwr      upr     p adj
## Treated-Control 238.5 144.2819 332.7181 0.0001077
## 
## $Block
##       diff      lwr      upr    p adj
## B-A 176.25 82.03193 270.4681 0.001399

Tamaño del efecto

Eta cuadrado:

modelo_dbca %>% 
  eta_squared()

F-Cohen:

modelo_dbca %>% 
  cohens_f()

Residuales

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

Prueba de normalidad:

modelo_dbca %>% 
  residuals() %>% 
  shapiro.test()

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.84408, p-value = 0.01116

Prueba de homocedasticidad tratamiento:

datos %>% 
  leveneTest(data = ., GST ~ Treatment)

Prueba de homocedasticidad bloque:

datos %>% 
  leveneTest(data = ., GST ~ Block)

Potencia estadística

DCA

wp.anova(k = 4,
         n = 16,
         f = 0.26,
         alpha = 0.01,
         type = "two.sided")

## Power for One-way ANOVA
## 
##     k  n    f alpha     power
##     4 16 0.26  0.01 0.0476529
## 
## NOTE: n is the total sample size (contrast, two.sided)
## URL: http://psychstat.org/anova

DBCA

wp.anova(k = 6,
         n = 16,
         f = 1.52,
         alpha = 0.01,
         type = "two.sided")

## Power for One-way ANOVA
## 
##     k  n    f alpha     power
##     6 16 1.52  0.01 0.9921112
## 
## NOTE: n is the total sample size (contrast, two.sided)
## URL: http://psychstat.org/anova

Diseño de experimentos

DCA

Diseño completamente al azar:

tratamientos <- c("A", "B", "C", "D")
repeticiones <- 15
dis_dca <- design.crd(
  trt = tratamientos,
  r = repeticiones,
  seed = 2023,
  randomization = TRUE
)

dis_dca$book

¿Este es el número de repeticiones necesaria?

wp.anova(
  k = 4,
  power = 0.8,
  f = 0.25, # tamaño de efecto medio
  alpha = 0.01,
  type = "two.sided"
)

## Power for One-way ANOVA
## 
##     k       n    f alpha power
##     4 190.243 0.25  0.01   0.8
## 
## NOTE: n is the total sample size (contrast, two.sided)
## URL: http://psychstat.org/anova

DBCA

Diseño en bloques completos al azar:

tratamientos <- c("A", "B", "C", "D")
bloques <- 3
dis_dbca <- design.rcbd(
  trt = tratamientos,
  r = bloques,
  seed = 2023,
  randomization = TRUE
)

dis_dbca$book

¿Este es el número de repeticiones necesaria?

wp.anova(k = 7,
         power = 0.8,
         f = 1.25, # tamaño de efecto medio
         alpha = 0.01,
         type = "two.sided")

## Power for One-way ANOVA
## 
##     k        n    f alpha power
##     7 13.58406 1.25  0.01   0.8
## 
## NOTE: n is the total sample size (contrast, two.sided)
## URL: http://psychstat.org/anova

Diseño Experimental

Análisis de experimentos

Edimer David Jaramillo

Bibliotecas

Datos

Tamaño del efecto

Anova de una vía (DCA)

Hipótesis

Nivel de significancia

Resumen descriptivo

Exploración

Gráfico 1

Gráfico 2

Modelo

Medias estimadas

Comparaciones múltiples

t-student

Tukey

Tamaño del efecto

Residuales

Diseño en bloques (DBA)

Conteo por combinación

Exploración

Modelo

Medias estimadas

Comparaciones múltiples

t-student

Tukey

Tamaño del efecto

Residuales

Potencia estadística

DCA

DBCA

Diseño de experimentos

DCA

DBCA