12  Técnicas no paramétricas

12.1 Librerías

library(ggplot2)
library(dplyr)
library(moments)
library(DescTools)
library(patchwork)
theme_set(see::theme_lucid() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()))

12.2 Introducción

Hasta ahora habíamos tratado con pruebas y técnicas que requieren ciertas características, particularmente normalidad y homocedasticidad, para dar resultados confiables. Pero, ¿qué pasa cuando nuestros datos no se ajustan a una distribución normal o nuestras varianzas no son homogéneas? Una opción es recurrir a algún método de transformación de datos, pero eso lo abordaremos en otra ocasión. Hoy hablaremos sobre una familia de técnicas cuyos requerimientos sobre las distribuciones de los datos no son tan estrictos: Las pruebas NO paramétricas o libres de Distribución. Para explicar en qué consisten y cuales son sus diferencias me gustaría hacerlo en forma de una sere de preguntas y respuestas:

  1. ¿Por qué se llaman no paramétricas? Estas técnicas reciben esos nombres debido a que no están basadas en los Parámetros de una Población (\(\mu\) o \(\sigma^2\)) y que en consecuencia no requieren que los conozcamos o estimemos. También podemos decir que son no paramétricas porque no trabajan con los parámetros de ninguna distribución, sino que trabajan con la distribución empírica de los datos; es decir, con la distribución propia de los datos.

  2. ¿Por qué no aplicarlas siempre? Una deducción rápida puede ser que, si estas pruebas no requieren de características especiales de nuestros datos, podemos utilizarlas siempre y dejar de lado las pruebas paramétricas. Esto se debe a que en general, las técnicas paramétricas son más precisas y tienen un mayor poder estadístico; es decir, que son mejores para encontrar un efecto realmente significativo (son más sensibles).

  3. ¿En qué casos aplicarlas?

    1. Cuando sepamos que la distribución de la variable en la población no es normal. Esta parte tiene sus bemoles porque el cubrir este único requisito dependerá de le exhaustividad de nuestro muestreo y si es realmente representativo de lo que sucede a un nivel poblacional. Dejando el tecnicismo de lado, recurriremos a estas técnicas en el mismo escenario en el que aplicaríamos alguna transformación a los datos: cuando tengamos evidencia lo suficientemente grande de que nuestros datos (muestras) no se ajustan a una distribución normal. Como ya habíamos comentado antes debemos sustentar nuestras conclusiones a partir de diversas evidencias y no solo un valor de p. Un análisis gráfico de la distribución de frecuencias (gráfico QQ, KDE o un gráfico de frecuencias), acompañado de una prueba de normalidad (e.g. Shapiro-Wilks) deberían ser lo mínimo realizable para evaluar la normalidad de nuestros datos. Adicionalmente podemos estimar la curtosis (altura de la curva) y el sesgo (simetría de la distribución) y ver qué tan diferentes son de el valor de 1 de una distribución normal.
    2. Otros casos a tener en cuenta son:
      1. Nuestros datos están en escala nominal u ordinal.
      2. El tamaño de muestra es muy pequeño para realizar una prueba paramétrica. ¿Qué tan pequeño? menos de 10-15 datos por muestra. Ën este punto es importante mencionar que con un n<5 no es posible identificar desviaciones de una distribución normal.
      3. El objetivo es realizar inferencias alrededor de la mediana y no de la media (útil con distribuciones muy sesgadas).
  4. ¿Qué pasa si, aún violando los supuestos, empleo una prueba paramétrica? En general, la consecuencia es que se obtienen resultados no confiables. Para ejemplificarlo tomemos al Análisis de la Varianza (ANOVA) paramétrico. Recordemos que en este análisis comparamos la Varianza dentro de los grupos contra la Varianza entre los grupos, por lo que si uno o más grupos presentan desviaciones muy importantes de la normalidad el resultado será un análisis sesgado, causado por las diferencias en la distribución global de frecuencias. Adicionalmente, si no cumplimos con el supuesto de homogeneidad de varianzas, el resultado será que la varianza dentro de uno o más grupos será mayor y, en consecuencia, la razón de las mismas tenderá a ser más pequeña y el resultado será un mayor nivel de significancia (valor de p).

Figura 12.1: Paramétrico vs. no paramétrico

Habiendo tocado esos puntos, exploremos las alternativas no paramétricas a algunas de las pruebas que hemos visto en el curso.

En esta sesión utilizaremos los archivos Datos1 2.csv, Datos_bloques.csv y Abundancias.csv, por lo que primero tendremos que cargarlos:

# Carga de datos
guppys <- read.csv("datos/Datos1 2.csv",
               header = F,
               skip = 1,
               sep = ",") |> select(1:5)
colnames(guppys) <- c("Dieta", "Periodo", "Rep", "LT", "PT")
head(guppys)
bloques <- read.csv("datos/Datos_bloques.csv")
head(bloques)
abundancias <- read.csv("datos/Abundancias.csv")
abundancias

Comencemos entonces a explorar algunas de las distintas herramientas que tenemos a nuestra disposición.

12.3 Comparación de muestras independientes: U de Mann-Whitney

También llamada prueba de Mann-Whitney-Wilcoxon o prueba de suma de rangos de Wilcoxon. Esta prueba es la alternativa a la prueba T de Student para muestras independientes, en donde la comparación se realiza en términos de las medianas (\(\tilde{x}\)) y por tanto la hipótesis nula será \(H_0: \tilde{X}_1 = \tilde{X}_2\).

Antes de pasar al modo de calcular el estadístico de prueba, hablemos sobre las consideraciones que debemos de tener para aplicar esta prueba:

  1. La variable dependiente deberá de estar en escala ordinal o continua.
  2. La variable independiente es una variable categórica con dos grupos.
  3. Las muestras deberán ser independientes; es decir, no deberá de existir una relación entre o dentro de los grupos (medidas repetidas)
  4. Aunque los datos no se ajusten a una distribución normal, deberán de tener la misma forma (ambos con forma de campana y sesgo negativo, o ambos con sesgo positivo, etc.)

El estadístico de prueba es el estadístico U, el cual será el valor mínimo de alguno de los siguientes cálculos, donde R es la suma de rangos. El rango de cada dato se obtiene ordenando todos los datos de mayor a menor, teniendo cuidado de no perder a qué grupo pertenece cada uno de ellos, y numerarlos de 1 hasta el número de datos.

  • \(U_1 = R_1 - \frac{n_1(n_1+1)}{2}\)

  • \(U_2 = R_2 - \frac{n_2(n_2 +1)}{2}\)

Para aplicarla en R utilizaremos la función wilcox.test(formula), pero primero exploremos nuestros datos. Nuestra variable de interés es la variable LT, por lo que realizaremos la evaluación de la normalidad, tanto de manera cuantitativa (prueba de Shapiro-Wilks, sesgo y curtosis), como de manera gráfica (Gráficos de densidad). Generemos una función que nos permita reunir estos datos en una sola tabla de manera rápida:

guppys |> group_by(Dieta) |> rstatix::shapiro_test(LT)
norms <- data.frame(Dieta = NA, k = NA, s = NA)
grupos <- unique(guppys$Dieta)

for (i in 1:length(grupos)) {
  sub_set <- na.omit(guppys$LT[guppys$Dieta == grupos[i]])

  norms[i,] <- c(as.character(grupos[i]),
                 round(kurtosis(sub_set), 2),
                 round(skewness(sub_set), 2)
                 )
}
norms

Y gráficamente:

kde_plots <- ggplot(data = guppys, aes(LT, color = Dieta)) + 
             geom_density(show.legend = F) + 
             facet_wrap(~Dieta, nrow = 2) +
             labs(title = "Longitud Total de Guppys bajo tres dietas",
                  subtitle = "Gráficos de densidad",
                  caption = "Archivo: Datos1 2.csv",
                  x = element_blank(),
                  y = element_blank()) +
             scale_y_continuous(breaks = NULL)
kde_plots

Debido a que todas las distribuciones fueron leptocúrticas (k > 1) y con sesgo negativo (s < 1), todas las pruebas de normalidad tuvieron niveles de significancia pequeños y los gráficos KDE mostraron tendencias hacia la bimodalidad, concluimos que los datos no se ajustan a una distribución normal, ¿no? Sí y no. En efecto la distribución de cada grupo no es normal, pero recordarás que la decisión del tipo de prueba a aplicar no está en función de la distribución de cada grupo, sino de la normalidad de los residuales del modelo ANOVA.

Nota

Esto era de esperarse debido a que hay una variable más de agrupamiento: el periodo; es decir, que hay información del inicio del experimento un periodo intermedio y un periodo final, solo se realiza este ejemplo con fines de enseñanza.

Imaginemos que nos interesa conocer si los organismos alimentados con la dieta B tuvieron una longitud total diferente a los alimentados con la dieta C, entonces apliquemos una prueba U-Mann-Whitney.

u_data <- subset(guppys, Dieta != "A")
wilcox.test(LT~factor(Dieta), data = u_data)

    Wilcoxon rank sum test with continuity correction

data:  LT by factor(Dieta)
W = 1408, p-value = 0.223
alternative hypothesis: true location shift is not equal to 0

Este resultado corresponde bastante bien con lo visto en los gráficos de densidad en el sentido de falta de diferencias, aunque podemos representarlo utilizando un gráfico de violín que es más adecuado para propósitos de comparación:

u_mann_plot <- ggplot(data = u_data, aes(x = Dieta, y = LT, fill = Dieta)) + 
               geom_violin(trim = T,
                           alpha = 0.7,
                           show.legend = F) + 
               geom_boxplot(width = 0.1,
                            fill = "white",
                            notch = T,
                            show.legend = F) + 
               labs(title = "Comparación de LT de guppys bajo las dietas B y C",
                    subtitle = "Gráficos de violín",
                    caption = "Archivo: Datos1 2.csv",
                    x = element_blank(),
                    y = element_blank()
                    ) +
               scale_fill_manual(values = c(rgb(29,149,79,
                                                maxColorValue = 255),
                                            rgb(126,172,255,
                                                maxColorValue = 255))
                                 )
u_mann_plot

12.3.1 Comparación de dos muestras dependientes (pareadas): Prueba por rangos de Wilcoxon

Otra pregunta que podemos estar interesados en responder con esta base de datos es si las Longitudes Totales al inicio del experimento fueron iguales o diferentes al final del experimento en la dieta A y asumamos, para fines de enseñanza, a) que los datos a este nivel no siguen una distribución normal y b) que tenemos mediciones pareadas; es decir, que la medición 1 del periodo inicial corresponde al mismo individuo que la medición 1 del periodo final. Bajo estos supuestos, realizaremos la prueba de rangos de Wilcoxon.

Esta prueba es la alternativa no paramétrica a la prueba t para datos pareados y, al igual que en el caso de la prueba U de Mann-Whitney, compara si las medianas entre ambos grupos son diferentes entre sí. Por otra parte, y del mismo modo que la prueba t para datos pareados, el cálculo se realiza a partir de las diferencias entre ambos grupos, asignándoles un rango y luego sumarlas según su signo, tal que:

  • \(W^+ = \sum R^+\)
  • \(W^- = \sum R^-\)

Donde R es el rango de cada diferencia y \(+\) y \(-\) representan su signo. El valor que utilizaremos durante la prueba será el más pequeño de los dos. La prueba de hipótesis se realiza utilizando la distribución normal, tal que:

  • \(Z = \frac{min(W) - \mu}{\sigma - r}\)

Donde: \(\mu = \frac{n(n+1)}{4}\), \(\sigma = \frac{\sqrt{n(n+1)(2n+1)}}{24}\) y r es un factor de reducción cuando existen t rangos empatados: \(r = \frac{t^3-t}{48}\).

En R, por fortuna, la implementación es mucho más sencilla, solo hay que considerar que al ser una prueba pareada, el número de datos en ambos grupos debe ser igual.

# Extraemos solo los datos que nos interesan
wilcox_data <- subset(guppys, (Dieta == "A" & Periodo != "M"))
# Vemos que al final tenemos un individuo menos que al inicio,
# por lo que hay que retirarlo de la base de datos (asumamos que es el último)
aggregate(LT~Periodo, data = wilcox_data, length)
wilcox_data <- wilcox_data[-c(20),]
#Ahora tenemos ns iguales, por lo que podemos aplicar la prueba pareada
aggregate(LT~Periodo, data = wilcox_data, length)

Ahora apliquemos la prueba.

wilcox.test(LT~Periodo, data = wilcox_data, paired = T)

    Wilcoxon signed rank exact test

data:  LT by Periodo
V = 190, p-value = 3.815e-06
alternative hypothesis: true location shift is not equal to 0

Al parecer las medianas son diferentes, lo cual se sostiene al analizar los resultados de manera gráfica:

w_rank_plot <- ggplot(data = wilcox_data,
                      aes(x = Periodo, y = LT, fill = Periodo)) + 
               geom_violin(trim = T,
                           alpha = 0.7,
                           show.legend = F) + 
               geom_boxplot(width = 0.1,
                            fill = "white",
                            notch = T,
                            show.legend = F) + 
               labs(title = "Comparación de LT de guppys en dos periodos",
                    subtitle = "Gráficos de violín",
                    caption = "Archivo: Datos1 2.csv",
                    x = element_blank(),
                    y = element_blank()
                    ) +
               scale_fill_manual(values = c(rgb(29,149,79,
                                                maxColorValue = 255),
                                            rgb(126,172,255,
                                                maxColorValue = 255))
                                 )
w_rank_plot

12.4 Prueba por Rangos de Friedman y coeficiente de acuerdo de Kendall

Antes de hablar de las siguientes dos técnicas, hablemos sobre los bloques aleatorizados completos. Este caso se da cuando contamos con un factor adicional a nuestros factores de interés. En este caso realizaremos bloques para poder reducir el error experimental derivado de dichos factores. Un ejemplo es tener una serie de votantes para evaluar una serie de aspectos, o evaluar la actividad proteolítica en distintas secciones del tracto digestivo de distintos peces, o diferentes mediciones de fósforo de distintas muestras realizadas por distintas personas.

Figura 12.2: Bloques aleatorizados completos

En este sentido, la Prueba por Rangos de Friedman nos permite analizar a la vez las diferencias entre grupos, según los bloques. Para ello asignaremos rangos a las observaciones de cada bloque y compararemos el promedio de los rangos de cada grupo con respecto a la media y la varianza de los rangos de manera global. Pausa el video un momento, analiza la ecuación de \(F_R\) y trata de adivinar qué tipo de distribución sigue \(F_R\). Una pista: esa distribución ya fue vista durante el curso. La respuesta es una distribución \(\chi^2\).

Ok, una vez que realizamos esos cálculos ya evaluamos las diferencias entre los grupos tomando en consideración a los bloques, pero ¿qué tan similares son estos bloques? o en otras palabras, ¿qué tan de acuerdo estuvieron los votantes? Para ello podemos utilizar el Coeficiente de acuerdo de Kendall (W), el cuál está en términos del error cuadrático de la asignación de rangos dentro de cada bloque y el número de bloques y grupos. W indica la magnitud del acuerdo, donde 0 representa total desacuerdo y 1 total acuerdo.

Ahora vayamos a R y apliquemos estas técnicas a nuestros datos por bloques. Para ello necesitaremos cambiar el formato de la base de datos de una base de datos tabular a una base de datos codificada utilizando la función melt():

bloques_molten <- reshape2::melt(bloques,
                                 id.vars = "Pez",
                                 variable.name = "Secc",
                                 value.name = "Act")
bloques_molten

Ahora podemos aplicar la prueba de Friedman utilizando la función friedman.test(a ~ b|c), donde a representa la variable de medición, b los grupos y c los bloques.

friedman.test(Act~Secc|Pez, data = bloques_molten)

    Friedman rank sum test

data:  Act and Secc and Pez
Friedman chi-squared = 4.7586, df = 3, p-value = 0.1903

Realicemos la comparación gráficamente:

fried_plot <- ggplot(data = bloques_molten,
                     aes(x = Secc, y = Act, fill = Secc)) + 
              geom_violin(trim = T,
                           alpha = 0.7,
                           show.legend = F) + 
               geom_boxplot(width = 0.1,
                            fill = "white",
                            notch = F,
                            show.legend = F) +
              labs(title = 
                     "Actividad proteolítica entre 4 secciones del tracto
                   digestivo de peces",
                  # subtitle = "Gráficos de violín",
                  caption = "Archivo: Datos1 2.csv",
                  x = element_blank(),
                  y = element_blank()
                  )
fried_plot

De estos resultados podríamos concluir que no hubo diferencias entre la actividad proteolítica entre las distintas secciones de los tractos digestivos, pero ¿qué tan diferentes entre sí fueron los peces? Utilicemos el coeficiente de acuerdo de Kendall utilizando la función KendallW(x, correct, test, na.rm) de la librería DescTools donde x representa una matriz de k grupos (columnas) y m bloques (renglones), correct es un argumento booleano si se debe de realizar la corrección por empates en los rangos y test para realizar una prueba \(\chi^2\) con \(H_0: W = 0\); es decir, si hubo o no un desacuerdo total entre los bloques.

#Esta función sí requiere los datos en forma tabular, por lo que x será nuestra variable bloques
DescTools::KendallW(bloques, correct = T, test = T)

    Kendall's coefficient of concordance Wt

data:  bloques
Kendall chi-squared = 7.5578, df = 8, subjects = 9, raters = 5, p-value
= 0.4778
alternative hypothesis: Wt is greater 0
sample estimates:
       Wt 
0.1889447 

Y comprobémoslo gráficamente (ojo, para graficar estamos utilizando la base codificada)

kendall_plot <- ggplot(data = bloques_molten,
                       # OJO: Estamos transformando Pez a un factor
                       aes(x = as.factor(Pez),
                           y = Act,
                           fill = as.factor(Pez))) +
                geom_violin(trim = T,
                           alpha = 0.7,
                           show.legend = F) + 
               geom_boxplot(width = 0.1,
                            fill = "white",
                            notch = F,
                            show.legend = F) +
                labs(title = 
                       "Actividad proteolítica de 9 peces",
                    # subtitle = "Gráficos de violín",
                    caption = "Archivo: Datos_bloques.csv",
                    x = element_blank(),
                    y = element_blank()
                    )
kendall_plot

12.5 ANOVA no paramétrico: Kruskall-Wallis:

Esta vez nos presentan la pregunta: ¿hubo diferencias entre las longitudes finales de los guppys, según la dieta con la que fueron alimentados?

Para responder a esta pregunta podríamos utilizar un Análisis de la Varianza de una vía; sin embargo, recordemos que estamos asumiendo que se viola el supuesto de normalidad de los datos a ese nivel y, por tanto, habrá que realizar una prueba no paramétrica: ANOVA de Kruskal-Wallis. Esta prueba es una extensión de la prueba U de Mann-Whitney a más grupos, por lo tanto también trabajará con la suma de rangos de los datos y el procedimiento es muy similar:

  1. Asignar rangos a los datos unidos (sin perder el grupo)
  2. Realizar la suma de rangos de cada grupo y compararlas

Al igual que en otras pruebas que hemos visto el día de hoy, la distribución de prueba es la distribución \(\chi^2\) con \(k-1\) grados de libertad.

Una vez que realizamos la prueba solo obtendremos si hubo diferencias en alguna de las medianas, pero no sabremos cuál o cuáles fueron diferentes; por lo tanto, habrá que realizar una prueba post-hoc: la prueba de comparaciones múltiples de Dunn, dada por:

La implementación de ambas pruebas en R es sumamente sencilla. Primero, el ANOVA Kruskal-Wallis utilizando la función kruskal.test(formula, data):

H <- kruskal.test(LT~Dieta, data = subset(guppys, Periodo == "F"))
H

    Kruskal-Wallis rank sum test

data:  LT by Dieta
Kruskal-Wallis chi-squared = 3.975, df = 2, p-value = 0.137

Vemos que aparentemente no hay diferencias entre las medias; sin embargo, realicemos la prueba de Dunn para confirmar. Para ello utilizaremos la función DunnTest(formula, data) de la librería DescTools. Vemos que en efecto, pareciera no haber diferencias entre las longitudes finales de los guppys en función de las dieta con la que fueron alimentados.

DescTools::DunnTest(LT~factor(Dieta),
                    data = subset(guppys, Periodo == "F"))

 Dunn's test of multiple comparisons using rank sums : holm  

    mean.rank.diff   pval    
B-A      -7.105263 0.2578    
C-A      -9.327485 0.2578    
C-B      -2.222222 0.6851    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Corroboremos esto de manera gráfica:

kruskal_plot <- ggplot(data = subset(guppys, Periodo == "F"), 
                       aes(x = Dieta, y = LT, fill = Dieta)) +
                geom_violin(trim = T,
                           alpha = 0.7,
                           show.legend = F) + 
               geom_boxplot(width = 0.1,
                            fill = "white",
                            notch = T,
                            show.legend = F) + 
                labs(title = "Comparación de LT finales de Guppys",
                    subtitle = "Gráficos de violín",
                    caption = "Archivo: Datos1 2.csv",
                    x = element_blank(),
                    y = element_blank())
kruskal_plot

12.6 ANOVA de dos vías no paramétrico: Scheirer-Ray-Hare

La última prueba que veremos en esta sesión es la versión no paramétrica para problemas en los que tenemos dos factores y muestras independientes; i.e., la alternativa no paramétrica al ANOVA de dos vías. Para este caso cambiemos ligeramente los datos y utilicemos datos de concentraciones de Alanina (Ala) en milípedos, evaluadas por sexo y dos especies:

mili <- read.csv("datos/milipedos.csv")
head(mili)

Los datos ya se encuentran en formato largo, pero veamos primero cuántas observaciones tenemos por cada combinación de Especie:Sexo:

mili |> group_by(Especie, Sexo) |> summarise(Ala = n())
`summarise()` has grouped output by 'Especie'. You can override using the
`.groups` argument.

El diseño es balanceado, sin embargo tenemos pocas muestras por grupo. Este problema es más de representatividad y qué tan válidas son las inferencias realizadas sobre tan pocos datos porque, si comprobamos nuestros supuestos, todo parece estar en orden:

mili_aov <- aov(Ala~Especie*Sexo, data = mili)
# performance::check_model(mili_aov)
performance::check_normality(mili_aov)
OK: residuals appear as normally distributed (p = 0.240).
performance::check_homogeneity(mili_aov)
OK: There is not clear evidence for different variances across groups (Bartlett Test, p = 0.976).

Veamos entonces el resultado de la prueba:

summary(mili_aov)
             Df Sum Sq Mean Sq F value   Pr(>F)    
Especie       2  55.26   27.63  13.082  0.00031 ***
Sexo          1 138.72  138.72  65.679 2.04e-07 ***
Especie:Sexo  2   6.89    3.45   1.631  0.22331    
Residuals    18  38.02    2.11                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Concluimos que, estadísticamente y a un \(\alpha = 0.05\), parece haber diferencias significativas entre especies y entre sexos, y que la diferencia entre especies no depende del sexo (o viceversa). Podríamos confiar “ciegamente” en estos resultados y aplicar la prueba post-hoc correspondiente para las especies; sin embargo, veamos qué nos dice la prueba no paramétrica:

mili_srh <- rcompanion::scheirerRayHare(Ala~Especie+Sexo, data = mili)

DV:  Ala 
Observations:  24 
D:  0.9995652 
MS total:  50 
mili_srh

Se mantienen las diferencias entre los sexos, pero las diferencias entre especies no están en conflicto con nuestra hipótesis de nulidad a un \(\alpha = 0.05\). ¿Qué prueba utilizar? Veamos la distribución de nuestros datos:

ala_sp_plot <- ggplot(data = mili,
                      aes(x = Especie, y = Ala,
                          fill = Especie,
                          color = Especie)) +
               geom_violin(alpha = 0.2) +
               geom_boxplot(fill = NA,
                            width = 0.1)
ala_sx_plot <- ggplot(data = mili,
                      aes(x = Sexo, y = Ala,
                          fill = Sexo,
                          color = Sexo)) +
               geom_violin(alpha = 0.2) +
               geom_boxplot(fill = NA,
                            width = 0.1)
ala_int_plot <- ggplot(data = mili,
                       aes(x = Sexo, y = Ala,
                           fill = Sexo,
                           color = Sexo)) +
                geom_violin(alpha = 0.2) +
                geom_boxplot(fill =NA,
                             width = 0.1) +
                labs(y = "[Ala] (mg/L)") +
                theme(axis.text.x = element_text(angle = 90)) +
                facet_wrap(~Especie)

# Crear gráfico compuesto utilizando Patchwork:
# Gráfico de interacción a la izquierda
# Efectos principales a la derecha
(ala_int_plot|(ala_sp_plot / ala_sx_plot) &
    # Efectos principales sin etiqueta del eje y
    labs(y = element_blank())) &
  theme(legend.position = "none") &
  labs(x = element_blank())

Como era de esperarse, las diferencias entre sexos son las más marcadas: la concentración de alanina es consistentemente mayor en machos que en hembras, independientemente de la especie (no hay interacción). ¿Qué pasa entre las especies? Sin duda alguna estamos en la cuerda floja; sin embargo, al menos con los datos disponibles, la especie 2 tiende hacia los valores más bajos.

La pregunta que debemos de hacernos aquí es: ¿biológicamente tienen relevancia las diferencias observadas? La decisión dependerá totalmente de la experiencia del investigador con el problema. Fuera de este ejemplo no he tratado con datos de concentraciones de Alanina, ni tampoco he estudiado milípedos, por lo que no puedo dar una conclusión informada. ¿Entonces? Prefiero errar de precavido y decir que las diferencias entre especies no son significativas (especialmente considerando que las diferencias entre sexos son mayores) y hacer la recomendación de replicar el experimento con un tamaño de muestra más grande (OJO, no incrementar el tamaño de muestra, eso se conoce como p-hacking).

Con esto llegamos al final de esta clase. Espero que haya sido de tu agrado y nos vemos en la siguiente.