Punto 2, Taller #1 Diseño de Experimentos

Contexto

La tasa de consumo de oxígeno es una medida crucial en la comprensión de los procesos fisiológicos de la gran mayoría de animales, lo cual incluye a un grupo de invertebrados que suele encontrarse más frecuentemente en ambientes marinos: los moluscos. Para poder estudiar la injerencia que puede llegar a tener la concentración/salinidad del agua marina y la morfofisiología de los moluscos en su tasa de consumo de oxígeno, se tomaron 48 datos de esta medida en 2 tipos distintos de molusco (A y B) y bajo 3 porcentajes distintos de concentración del agua (50, 75 y 100%).

require(table1)
require(ggplot2)
require(broom)
require(knitr)
library(dplyr)
require(plotly)
library(car)
require(tidyr)
require(gt)
require(htmltools)
load("~/Archivos R/moluscos.RData")

1. Análisis exploratorio univariado

Previo a la realización de las pruebas estadísticas correspondientes, debe hacerse un proceso de observación y comprensión de cada uno de los indicadores estadísticos principales, tanto de tendencia central, como de dispersión y forma, principalmente. En este estudio se presentan 2 factores de tratamiento (concentración del agua y tipo de molusco) y una sola variable de respuesta. Por tanto, únicamente se realizará el análisis exploratorio de la respuesta, que es la cantidad de consumo de oxígeno, empleando las bibliotecas gt y htmltools para hacer visibles los cálculos de cada indicador:

resumen_m1 <- BD_moluscos %>%
  summarise(
    across(
      cons_o,
      list(
        n = ~sum(!is.na(.x)),
        Media = ~mean(.x, na.rm = TRUE),
        DE = ~sd(.x, na.rm = TRUE),
        Mediana = ~median(.x, na.rm = TRUE),
        Min = ~min(.x, na.rm = TRUE),
        Max = ~max(.x, na.rm = TRUE),
        CV = ~sd(.x, na.rm = TRUE) / mean(.x, na.rm = TRUE) * 100
      )
    )
  ) %>%
  pivot_longer(
    cols = everything(),
    names_to = c("Variable", "Estadistico"),
    names_pattern = "^(.*)_(n|Media|DE|Mediana|Min|Max|CV)$"
  ) %>%
  pivot_wider(
    names_from = Estadistico,
    values_from = value
  ) %>%
  select(Variable, n, Media, DE, Mediana, Min, Max, CV) %>%
  mutate(
    across(where(is.numeric), ~round(.x, 2))
    )
resumen_m1 <- resumen_m1 %>%
  rename(
    `Respuesta` = Variable,
    `n` = n,
    `Media` = Media,
    `DE` = DE,
    `Mediana` = Mediana,
    `Mínimo` = Min,
    `Máximo` = Max,
    `CV (%)` = CV
  )
resumen_m1 <- resumen_m1 %>%
  mutate(
      "Respuesta" = "Consumo oxígeno"
    )
tagList(
  tags$div(
    style = "font-size: 14px; color: #5f6368; margin-bottom: 10px;",
    "Tabla 1.Indicadores estadísticos principales de la medida de consumo de oxígeno"
  ),
  
  resumen_m1 %>%
    gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  ) %>%
    cols_align(
      align = "center",
      columns = everything()
    )
)

Tabla 1.Indicadores estadísticos principales de la medida de consumo de oxígeno

Respuesta	n	Media	DE	Mediana	Mínimo	Máximo	CV (%)
Consumo oxígeno	48	9.3	3.68	9.7	1.8	18.8	39.58

Ahora, además de los datos de los indicadores estadísticos de dispersión y tendencia central, se debe encontrar el nivel de simetría de los datos, para lo cual será necesario realizar un histograma de los datos, empleando la biblioteca ggplot2:

dist_o=ggplot(BD_moluscos, aes(x = cons_o)) +
  geom_histogram(bins = 10, fill = "#4682B4", color = "white") +
  labs(title = "Distribución de las medidas de consumo de oxígeno",
       x = "Consumo de oxígeno", y = "Frecuencia") +theme_minimal()+
  theme(axis.text.x = element_text(angle = 20, hjust = 1),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))
ggplotly(dist_o)

Con esta gráfica, podría afirmarse que los datos son relativamente asímetricos de forma negativa, encontrando la mayor densidad de datos hacia 5,57. Sin embargo, al no ser una diferencia significativa, es posible decir, en cambio, que la distribución de los datos es simétrica.

2. Análisis exploratorio bivariado

Ahora, debido a que este estudio únicamente tiene en cuenta una variable de respuesta, es difícil realizar observaciones acerca de sus datos. Por tanto, se debe hacer un análisis bivariado, en el que primero se comparen uno a uno los tratamientos con la respuesta y, posteriormente, sea posible analizar la relación que existe entre los factores juntos y la respuesta, siguiendo las pautas del diseño factorial.

a. Relación de la concentración del agua con el consumo de oxígeno

En primera instancia, se realizará un análisis sobre el primer factor de tratamiento, correspondiente a la concentración del agua del mar. Para esto, se tabularán los datos de los indicadores de tendencia central de esta relación a continuación, empleando las bibliotecas Kable y KableExtra.

require(kableExtra)
tabla_cons_o_final <- BD_moluscos %>%
  summarise(
    `50 % - Media` = mean(cons_o[c_agua == 50], na.rm = TRUE),
    `50 % - DE` = sd(cons_o[c_agua == 50], na.rm = TRUE),
    `75 % - Media` = mean(cons_o[c_agua == 75], na.rm = TRUE),
    `75 % - DE` = sd(cons_o[c_agua == 75], na.rm = TRUE),
    `100 % - Media` = mean(cons_o[c_agua == 100], na.rm = TRUE),
    `100 % - DE` = sd(cons_o[c_agua == 100], na.rm = TRUE)
  ) %>%
  mutate(Respuesta = "Consumo de oxígeno") %>%
  select(Respuesta, `50 % - Media`, `50 % - DE`, `75 % - Media`, `75 % - DE`,
    `100 % - Media`, `100 % - DE`) %>%
  mutate(across(where(is.numeric), ~round(.x, 2)))
kable(tabla_cons_o_final,
  col.names = c("Respuesta", "Media", "DE", "Media", "DE", "Media", "DE"
  ), caption = "Tabla 2.Indicadores estadísticos principales entre concentración del agua y consumo de oxígeno",
  align = "ccccccc",
  escape = FALSE
) %>%
  add_header_above(
    c(
      " " = 1,
      "50 %" = 2,
      "75 %" = 2,
      "100 %" = 2
    ),
    align = "c"
  ) %>%
  add_header_above(
    c(" " = 1,
      "Concentración del agua" = 6), align = "c"
  ) %>%
  kable_styling(
    bootstrap_options = "bordered",
    full_width = FALSE,
    position = "center",
    font_size = 14
  ) %>%
  row_spec(0, bold = TRUE, align = "center") %>%
  column_spec(1, bold = TRUE, width = "3.5cm")

Tabla 2.Indicadores estadísticos principales entre concentración del agua y consumo de oxígeno

	Concentración del agua
	50 %		75 %		100 %
Respuesta	Media	DE	Media	DE	Media	DE
Consumo de oxígeno	12.25	3.2	6.99	2.8	8.67	3

Como adicional, se realiza un gráfico boxplot, utilizando ggplot2, para identificar las regiones de esta interacción que están más relacionadas, que se ven por las regiones que se traslapan, además de ver los datos atípicos que tenga la muestra.

dis_int_conc_cons <- ggplot(BD_moluscos, aes(x = as.factor(c_agua), y = cons_o)) +
  geom_boxplot(fill = c("#E76F51", "#2A9D8F", "#4361EE")) +
  labs(x = "Concentración del agua", y = "Consumo de oxígeno")

ggplotly(dis_int_conc_cons) %>%
layout(title = "",
       xaxis = list(title = "Concentración del agua",
                    showgrid = TRUE, gridcolor = "#D9D9D9"),
       yaxis = list(title = "Consumo de oxígeno",
                    showgrid = TRUE, gridcolor = "#D9D9D9"),
       paper_bgcolor = "#FFFFFF",
       plot_bgcolor = "#FFFFFF")

Estos resultados muestran que, en primer lugar, los datos obtenidos de las 3 concentraciones de agua estudiadas no difieren tanto entre sí, es decir, tienen desviaciones estándar muy similares. Además, puede afirmarse que , de los 3 tratamientos posibles, es en el 50% de concentración del agua de mar en donde los resultados varían más con respecto a las otras dos concentraciones estudiadas.

b. Relación del tipo de molusco con el consumo de oxígeno

Para este caso, se realizará el mismo procedimiento del inciso a, pero esta vez analizando la relación que existe entre el tipo de molusco y el consumo de oxígeno:

tabla_mol_final <- BD_moluscos %>%
  summarise(
    `A - Media` = mean(cons_o[molusco == "A"], na.rm = TRUE),
    `A - DE` = sd(cons_o[molusco == "A"], na.rm = TRUE),
    `B - Media` = mean(cons_o[molusco == "B"], na.rm = TRUE),
    `B - DE` = sd(cons_o[molusco == "B"], na.rm = TRUE),
  ) %>%
  mutate(Respuesta = "Consumo de oxígeno") %>%
  select(Respuesta, `A - Media`, `A - DE`, `B - Media`, `B - DE`) %>%
  mutate(across(where(is.numeric), ~round(.x, 2)))
kable(tabla_mol_final,
      col.names = c("Respuesta", "Media", "DE", "Media", "DE"), 
      caption = "Tabla 3.Indicadores estadísticos principales entre el tipo de molusco y el consumo de oxígeno",
      align = "ccccccc",
      escape = FALSE
) %>%
  add_header_above(
    c(
      " " = 1,
      "A" = 2,
      "B" = 2
    ),
    align = "c"
  ) %>%
  add_header_above(
    c(" " = 1,
      "Molusco" = 4), align = "c"
  ) %>%
  kable_styling(
    bootstrap_options = "bordered",
    full_width = FALSE,
    position = "center",
    font_size = 14
  ) %>%
  row_spec(0, bold = TRUE, align = "center") %>%
  column_spec(1, bold = TRUE, width = "3.5cm")

Tabla 3.Indicadores estadísticos principales entre el tipo de molusco y el consumo de oxígeno

	Molusco
	A		B
Respuesta	Media	DE	Media	DE
Consumo de oxígeno	10	3.27	8.61	4

dis_int_mol_cons <- ggplot(BD_moluscos, aes(x = molusco, y = cons_o)) +
  geom_boxplot(fill = c("#E68F71", "#1368EE")) +
  labs(x = "Tipo de molusco", y = "Consumo de oxígeno")

ggplotly(dis_int_mol_cons) %>%
  layout(title = "",
         xaxis = list(title = "Tipo de molusco",
                      showgrid = TRUE, gridcolor = "#D9D9D9"),
         yaxis = list(title = "Consumo de oxígeno",
                      showgrid = TRUE, gridcolor = "#D9D9D9"),
         paper_bgcolor = "#FFFFFF",
         plot_bgcolor = "#FFFFFF")

En este caso, tanto el análisis de indicadores estadísticos como la representación mediante el boxplot indican que los datos de la interacción entre este factor y la respuesta, independientemente del tipo de molusco, indican una relación bastante estrecha entre sí, por lo que este factor no parece ser muy determinante en afectar la cantidad de consumo de oxígeno.

3. Análisis de varianza para la interacción entre variables (ANOVA)

Antes de realizar propiamente el análisis de varianzas, es necesario observar el gráfico que muestra la interacción entre todas las variables del estudio, que puede hacerse gracias a la biblioteca ggplot 2 y se muestra a continuación:

graf_int=ggplot(BD_moluscos, aes(x = c_agua, y = cons_o, colour = molusco))+geom_smooth()+geom_point()+labs(
  title = "Consumo de oxígeno de acuerdo con el tipo de molusco
  y la concentración del agua del mar",
  x = "Concentración del agua (%)",
  y = "Consumo de oxígeno",
  colour = "Tipo de molusco"
) +
  theme(plot.title = element_text(hjust = 0.5))
ggplotly(graf_int)

Ahora bien, una vez hecho todo este análisis exploratorio, se realizó un ANOVA para confirmar cuál de las variables del experimento está más relacionada con el consumo de oxígeno en los especímenes estudiados, obteniendo los siguientes resultados:

anova_mol=aov(cons_o~as.factor(c_agua)+molusco+as.factor(c_agua):molusco, data=BD_moluscos)
tabla=tidy(anova_mol)
tabla$term <- c("Concentración del agua", "Molusco", "Interacción", "Residuals")
tabla <- tabla %>%
  mutate(
    across(c(sumsq, meansq, statistic), ~ round(.x, 1)),
    p.value = ifelse(is.na(p.value), "", format.pval(p.value, digits = 3, eps = 0.001))
  )

kable(tabla,
      col.names = c("Fuente de variación", "Grados de libertad",
                    "Suma de cuadrados", "Cuadrado medio",
                    "Valor F", "p-valor"),
      caption = "Tabla 4. Anova de dos vías mostrando interacción entre  consumo de oxígeno, tipo de molusco y concentración del agua de mar",
      align = "c") %>%
  kable_styling(bootstrap_options = "bordered", full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE)

Tabla 4. Anova de dos vías mostrando interacción entre consumo de oxígeno, tipo de molusco y concentración del agua de mar

Fuente de variación	Grados de libertad	Suma de cuadrados	Cuadrado medio	Valor F	p-valor
Concentración del agua	2	230.8	115.4	13.2	<0.001
Molusco	1	23.2	23.2	2.7	0.111
Interacción	2	15.4	7.7	0.9	0.424
Residuals	42	368.0	8.8	NA

Con la información de la Tabla 4 y la gráfica demostrando la interacción de los factores, pueden reafirmarse diversas ideas previas. Entre estas, se reafirma la poca diferencia existente entre los tratamientos hechos a partir del tipo de molusco, así como su poca incidencia en la cantidad de consumo de oxígeno. Asimismo, puede confirmarse que la variable que más incide en la cantidad de oxígeno consumido por el molusco es la concentración del agua de mar, siendo entonces que entre menor sea esta concentración, más elevada será la cantidad de oxígeno que consuma el animal.

Por último, para terminar de confirmar estas observaciones, se realizará un postANOVA, o prueba post-hoc, realizando la prueba de Shapiro-Wilk (con su respectivo histograma y gráfica Q-Q) y de homogeneidad de varianzas, las cuales se presentarán a continuación:

residuales_m <- residuals(anova_mol)
shapiro_residualesm <- shapiro.test(residuales_m)

tabla_shapiro_residm <- data.frame(
  Prueba = "Normalidad de residuales (Anova interacción concentración del agua y tipo de molusco)",
  Estadistico_W = round(unname(shapiro_residualesm$statistic), 3),
  Valor_p = round(shapiro_residualesm$p.value, 4),
  Normalidad = ifelse(shapiro_residualesm$p.value > 0.05, "No se rechaza", "Se rechaza")
)

kable(tabla_shapiro_residm,
      col.names = c("Prueba", "Estadístico W", "Valor p", "Normalidad (α = 0.05)"),
      caption = "Prueba de Shapiro-Wilk sobre los residuales del modelo",
      row.names = FALSE,
      align = "lccc")

Prueba de Shapiro-Wilk sobre los residuales del modelo

Prueba	Estadístico W	Valor p	Normalidad (α = 0.05)
Normalidad de residuales (Anova interacción concentración del agua y tipo de molusco)	0.958	0.0857	No se rechaza

par(mfrow = c(1,2))
hist(residuales_m, main = "Histograma de residuales", xlab = "Residuales", col = "lightblue")
qqnorm(residuales_m); qqline(residuales_m, col = "red")

levene_conc <- leveneTest(cons_o ~ as.factor(c_agua), data = BD_moluscos)

tabla_levenem <- as.data.frame(levene_conc)
rownames(tabla_levenem) <- NULL
tabla_levenem$`Pr(>F)` <- format.pval(tabla_levenem$`Pr(>F)`, digits = 3, eps = 0.001)
tabla_levenem$`F value` <- round(tabla_levenem$`F value`, 3)

kable(tabla_levenem,
      col.names = c("GL", "F", "Valor p"),
      caption = "Prueba de homogeneidad de varianzas de concentración del agua entre hábitats",
      align = "ccc",
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE,
                position = "center",
                font_size = 14)

Prueba de homogeneidad de varianzas de concentración del agua entre hábitats

GL	F	Valor p
2	0.217	0.806
45	NA	NA

Estas pruebas demuestran que no existe variabilidad grande entre los grupos de los factores, ni siquiera en los tratamientos de mayor diversidad como lo son los de la concentración del agua del mar. Sin embargo, esto sólo dice que las varianzas entre estos grupos son muy similares, no habla como tal de la injerencia que tengan sobre la variable de respuesta.

Para poder conocer la respuesta a la última pregunta, es necesario realizar un test LSD y graficar estos resultados, comparando nuevamente a los valores de consumo de oxígeno con los de concentración del agua de mar.

require(agricolae)
lsd_mol <- LSD.test(anova_mol, "as.factor(c_agua)", p.adj = "none", console = FALSE)
plot(lsd_mol,
     main = "Comparación de consumo de oxígeno por concentración del agua (LSD)",
     xlab = "Concentración del agua", ylab = "Consumo de oxígeno")

En este gráfico se observa entonces que sí existe una variación entre los tratamientos del factor de concentración del agua de mar, donde los tratamientos del 75 y 100% de concentración tienen el mismo grupo asignado (b) al tener medias similares; mientras que, el tratamiento de 50% de concentración de agua de mar queda aislado en el grupo a por tener una media mucho mayor a la del resto de tratamientos, tal como el análisis exploratorio y el ANOVA habían previsto.

Conclusiones

Tras las pruebas realizadas y el análisis de todos los datos estudiados junto a sus interacciones, es posible concluir que:

1. No existe una relación real entre el tipo de molusco y un cambio en la cantidad de consumo de oxígeno del animal.
2. Tampoco existen cambios drásticos al estudiar en conjunto, como interacción, el tipo de molusco y la concentración del agua de mar con la cantidad de consumo de oxígeno.
3. Existe una fuerte relación entre una baja concentración del agua de mar y una alta cantidad de consumo de oxígeno, haciendo de este el único factor incidente para la variable de respuesta estudiada. Se esperaría entonces que los moluscos que viven en zonas con concentración de agua marina menores al 50% tiendan a necesitar una mucho mayor cantidad de oxígeno en comparación a otros moluscos que vivan en zonas con concentraciones mayores al 50%, aproximadamente.