Análisis de varianza de dos vías en R

Two-way ANOVA in R.

Aprende cómo hacer un ANOVA de dos vías en R. También aprenderás su objetivo, hipótesis, supuestos y cómo interpretar los resultados.

Introducción

El ANOVA de dos vías (análisis de varianza) es un método estadístico que permite evaluar el efecto simultáneo de dos variables categóricas en una variable continua cuantitativa.

El ANOVA de dos vías es una extensión del ANOVA de una vía ya que permite evaluar los efectos sobre una respuesta numérica de dos variables categóricas en lugar de una. La ventaja de un ANOVA de dos vías sobre un ANOVA de una vía es que se prueba la relación entre dos variables, teniendo en cuenta el efecto de una tercera variable. Además, también permite incluir la posible interacción de las dos variables categóricas en la respuesta.

La ventaja de un ANOVA de dos vías sobre un ANOVA de una vía es bastante similar a la ventaja de una correlación sobre una regresión lineal múltiple:

• La correlación mide la relación entre dos variables cuantitativas. La regresión lineal múltiple también mide la relación entre dos variables, pero esta vez teniendo en cuenta el posible efecto de otras covariables.
• El ANOVA de una vía prueba si una variable cuantitativa es diferente entre grupos. El ANOVA de dos vías también prueba si una variable cuantitativa es diferente entre grupos, pero esta vez teniendo en cuenta el efecto de otra variable cualitativa.

Anteriormente, hemos hablado sobre el ANOVA de una vía en R . Ahora, mostramos cuándo, por qué y cómo realizar un ANOVA de dos vías en R.

Antes de continuar, me gustaría mencionar y describir brevemente algunos métodos y pruebas estadísticas relacionados para evitar cualquier confusión:

• El software de IA en los centros de llamadas revoluciona el servici...
• Una Guía Completa sobre los Precios, Planes y Beneficios de Insight...
• Informe de McKinsey ¿Qué significa para las ventas B2B?

Un test t de Student se utiliza para evaluar el efecto de una variable categórica en una variable continua cuantitativa, cuando la variable categórica tiene exactamente 2 niveles:

• Test t de Student para muestras independientes si las observaciones son independientes (por ejemplo: si comparamos la edad entre mujeres y hombres).
• Test t de Student para muestras pareadas si las observaciones son dependientes, es decir, cuando vienen en pares (es el caso cuando se miden a los mismos sujetos dos veces, en dos momentos diferentes en el tiempo, antes y después de un tratamiento, por ejemplo).

Para evaluar el efecto de una variable categórica en una variable cuantitativa, cuando la variable categórica tiene 3 o más niveles:

• ANOVA de una vía (a menudo simplemente referido como ANOVA) si los grupos son independientes (por ejemplo, un grupo de pacientes que recibieron el tratamiento A, otro grupo de pacientes que recibieron el tratamiento B y el último grupo de pacientes que no recibieron tratamiento o un placebo).
• ANOVA de medidas repetidas si los grupos son dependientes (cuando los mismos sujetos se miden tres veces, en tres momentos diferentes en el tiempo, antes, durante y después de un tratamiento, por ejemplo).

Un ANOVA de dos vías se utiliza para evaluar los efectos de 2 variables categóricas (y su posible interacción) en una variable continua cuantitativa. Este es el tema de esta publicación.

La regresión lineal se utiliza para evaluar la relación entre una variable dependiente continua cuantitativa y una o varias variables independientes:

• Regresión lineal simple si hay sólo una variable independiente (que puede ser cuantitativa o cualitativa).
• Regresión lineal múltiple si hay al menos dos variables independientes (que pueden ser cuantitativas, cualitativas o una mezcla de ambas).

Un ANCOVA (análisis de covarianza) se utiliza para evaluar el efecto de una variable categórica en una variable cuantitativa, controlando el efecto de otra variable cuantitativa (conocida como covariable). El ANCOVA es en realidad un caso especial de regresión lineal múltiple con una mezcla de una variable independiente cualitativa y una cuantitativa.

En esta publicación, empezamos explicando cuándo y por qué es útil un ANOVA de dos vías, luego hacemos algunos análisis descriptivos preliminares y presentamos cómo realizar un ANOVA de dos vías en R. Finalmente, mostramos cómo interpretar y visualizar los resultados. También mencionamos brevemente e ilustramos cómo verificar los supuestos subyacentes.

Datos

Para ilustrar cómo realizar una ANOVA de dos vías en R, utilizamos el conjunto de datos penguins, disponible en el paquete {palmerpenguins}.

No es necesario importar el conjunto de datos, pero primero necesitamos cargar el paquete y luego llamar al conjunto de datos:

# install.packages("palmerpenguins")library(palmerpenguins)

dat <- penguins # renombrar los datosstr(dat) # estructura de los datos

## tibble [344 × 8] (S3: tbl_df/tbl/data.frame)##  $ species          : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...##  $ island           : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...##  $ bill_length_mm   : num [1:344] 39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...##  $ bill_depth_mm    : num [1:344] 18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...##  $ flipper_length_mm: int [1:344] 181 186 195 NA 193 190 181 195 193 190 ...##  $ body_mass_g      : int [1:344] 3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...##  $ sex              : Factor w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ...##  $ year             : int [1:344] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...

El conjunto de datos contiene 8 variables para 344 pingüinos, resumidos a continuación:

summary(dat)

##       species          island    bill_length_mm  bill_depth_mm  ##  Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  ##  Chinstrap: 68   Dream    :124   1st Qu.:39.23   1st Qu.:15.60  ##  Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  ##                                  Mean   :43.92   Mean   :17.15  ##                                  3rd Qu.:48.50   3rd Qu.:18.70  ##                                  Max.   :59.60   Max.   :21.50  ##                                  NA's   :2       NA's   :2      ##  flipper_length_mm  body_mass_g       sex           year     ##  Min.   :172.0     Min.   :2700   female:165   Min.   :2007  ##  1st Qu.:190.0     1st Qu.:3550   male  :168   1st Qu.:2007  ##  Median :197.0     Median :4050   NA's  : 11   Median :2008  ##  Mean   :200.9     Mean   :4202                Mean   :2008  ##  3rd Qu.:213.0     3rd Qu.:4750                3rd Qu.:2009  ##  Max.   :231.0     Max.   :6300                Max.   :2009  ##  NA's   :2         NA's   :2

En esta publicación, nos centraremos en las siguientes tres variables:

• species: la especie del pingüino (Adelie, Chinstrap o Gentoo)
• sex: sexo del pingüino (hembra y macho)
• body_mass_g: masa corporal del pingüino (en gramos)

Si es necesario, se puede encontrar más información sobre este conjunto de datos ejecutando ?penguins en R.

body_mass_g es la variable cuantitativa continua y será la variable dependiente, mientras que species y sex son ambas variables cualitativas.

Estas dos últimas variables serán nuestras variables independientes, también conocidas como factores. Asegúrese de que R las lea como factores. Si no es así, deberán ser transformadas en factores.

Objetivos e hipótesis de un ANOVA de dos vías

Como se mencionó anteriormente, un ANOVA de dos vías se utiliza para evaluar simultáneamente el efecto de dos variables categóricas en una variable continua cuantitativa.

Se le llama ANOVA de dos vías porque estamos comparando grupos que están formados por dos variables categóricas independientes.

Aquí, nos gustaría saber si la masa corporal depende de la especie y/o del sexo. En particular, estamos interesados en:

1. medir y probar la relación entre la especie y la masa corporal,
2. medir y probar la relación entre el sexo y la masa corporal y
3. potencialmente comprobar si la relación entre la especie y la masa corporal es diferente para las hembras y los machos (lo que equivale a comprobar si la relación entre el sexo y la masa corporal depende de la especie)

Las dos primeras relaciones se denominan efectos principales, mientras que el tercer punto se conoce como el efecto de interacción.

Los efectos principales prueban si al menos un grupo es diferente de otro (mientras se controla la otra variable independiente). Por otro lado, el efecto de interacción tiene como objetivo probar si la relación entre dos variables difiere dependiendo del nivel de una tercera variable.

Cuando se realiza un ANOVA de dos vías, no es obligatorio probar el efecto de interacción. Sin embargo, omitir un efecto de interacción puede conducir a conclusiones erróneas si el efecto de interacción está presente.

Si volvemos a nuestro ejemplo, tenemos las siguientes pruebas de hipótesis :

Efecto principal del sexo en la masa corporal:

• H0: la media de la masa corporal es igual entre hembras y machos
• H1: la media de la masa corporal es diferente entre hembras y machos

Efecto principal de la especie en la masa corporal:

• H0: la media de la masa corporal es igual entre las 3 especies
• H1: la media de la masa corporal es diferente para al menos una especie

Interacción entre sexo y especie:

• H0: no hay interacción entre sexo y especie, lo que significa que la relación entre la especie y la masa corporal es la misma para hembras y machos (de manera similar, la relación entre el sexo y la masa corporal es la misma para las 3 especies)
• H1: hay una interacción entre sexo y especie, lo que significa que la relación entre la especie y la masa corporal es diferente para hembras que para machos (de manera similar, la relación entre el sexo y la masa corporal depende de la especie)

Supuestos de un ANOVA de dos vías

La mayoría de las pruebas estadísticas requieren algunos supuestos para que los resultados sean válidos, y un ANOVA de dos vías no es una excepción.

Los supuestos de un ANOVA de dos vías son similares a los de un ANOVA de una vía. Para resumir:

• Tipo de variable: la variable dependiente debe ser continua cuantitativa, mientras que las dos variables independientes deben ser categóricas (con al menos dos niveles).
• Independencia: las observaciones deben ser independientes entre grupos y dentro de cada grupo.
• Normalidad:
• Para muestras pequeñas, los datos deben seguir aproximadamente una distribución normal
• Para muestras grandes (generalmente n ≥ 30 en cada grupo/muestra), no se requiere normalidad (gracias al teorema del límite central)
• Igualdad de varianzas: las varianzas deben ser iguales en todos los grupos.
• Valores atípicos: no debe haber valores atípicos significativos en ningún grupo.

Se pueden encontrar más detalles sobre estos supuestos en los supuestos de un ANOVA de una vía.

Ahora que hemos visto los supuestos subyacentes del ANOVA de dos vías, los revisamos específicamente para nuestro conjunto de datos antes de aplicar la prueba e interpretar los resultados.

Tipo de variable

La variable dependiente, la masa corporal, es continua cuantitativa, mientras que ambas variables independientes, el sexo y la especie, son variables cualitativas (con al menos 2 niveles).

Por lo tanto, se cumple este supuesto.

Independencia

La independencia se verifica generalmente en función del diseño del experimento y de cómo se han recopilado los datos.

Para simplificar, las observaciones suelen ser:

• independientes si cada unidad experimental (en este caso, un pingüino) ha sido medida solo una vez y las observaciones se recogen de una porción representativa y seleccionada al azar de la población, o
• dependientes si cada unidad experimental ha sido medida al menos dos veces (como suele ser el caso en el campo médico, por ejemplo, con dos mediciones en los mismos sujetos; antes y después del tratamiento).

En nuestro caso, la masa corporal se ha medido solo una vez en cada pingüino, y en una muestra representativa y aleatoria de la población, por lo que se cumple la suposición de independencia.

Normalidad

Tenemos una muestra grande en todos los subgrupos (cada combinación de los niveles de los dos factores, llamada celda):

table(dat$species, dat$sex)

##            ##             female male##   Adelie        73   73##   Chinstrap     34   34##   Gentoo        58   61

por lo que no es necesario comprobar la normalidad.

Para completitud, aún mostramos cómo verificar la normalidad, como si tuviéramos una muestra pequeña.

Existen varios métodos para probar la suposición de normalidad. Los métodos más comunes son:

• un gráfico QQ por grupo o en los residuos, y/o
• un histograma por grupo o en los residuos, y/o
• una prueba de normalidad (por ejemplo, la prueba de Shapiro-Wilk) por grupo o en los residuos.

La forma más fácil/corta es verificar la normalidad con un gráfico QQ en los residuos. Para dibujar este gráfico, primero debemos guardar el modelo:

# guardar modelo
mod <- aov(body_mass_g ~ sex * species,  data = dat)

Este fragmento de código se explicará más adelante.

Ahora podemos dibujar el gráfico QQ en los residuos. Mostramos dos formas de hacerlo, primero con la función plot() y segundo con la función qqPlot() del paquete {car}:

# método 1plot(mod, which = 2)

# método 2library(car)

qqPlot(mod$residuals,  id = FALSE # quitar identificación de puntos)

El código del método 1 es ligeramente más corto, pero no muestra el intervalo de confianza alrededor de la línea de referencia.

Si los puntos siguen la línea recta (llamada línea de Henry) y caen dentro de la banda de confianza, podemos asumir la normalidad. Este es el caso aquí.

Si prefiere verificar la normalidad en función del histograma de los residuos, aquí está el código:

# histogramahist(mod$residuals)

El histograma de los residuos muestra una distribución gaussiana, lo que concuerda con la conclusión del gráfico QQ.

Aunque el gráfico QQ y el histograma son suficientes para verificar la normalidad, si desea probarlo más formalmente con una prueba estadística, la prueba de Shapiro-Wilk también se puede aplicar a los residuos:

# prueba de normalidadshapiro.test(mod$residuals)

## ##  Shapiro-Wilk normality test## ## data:  mod$residuals## W = 0.99776, p-value = 0.9367

⇒ No rechazamos la hipótesis nula de que los residuos siguen una distribución normal (valor p = 0,937).

A partir del QQ-plot, el histograma y la prueba de Shapiro-Wilk, concluimos que no rechazamos la hipótesis nula de normalidad de los residuos.

La suposición de normalidad está así verificada, ahora podemos verificar la igualdad de las varianzas. 2

Homogeneidad de varianzas

La igualdad de varianzas, también conocida como homogeneidad de varianzas o homocedasticidad, puede ser verificada visualmente con la función plot():

plot(mod, which = 3)

Dado que la dispersión de los residuos es constante, la línea suave roja es horizontal y plana, por lo que parece que se cumple la suposición de varianza constante aquí.

El gráfico de diagnóstico anterior es suficiente, pero si lo prefiere, también se puede probar de manera más formal con la prueba de Levene (también del paquete {car}): 3

leveneTest(mod)

## Levene's Test for Homogeneity of Variance (center = median)##        Df F value Pr(>F)## group   5  1.3908 0.2272##       327

⇒ No rechazamos la hipótesis nula de que las varianzas son iguales (valor-p = 0.227).

Ambos enfoques, visual y formal, llegan a la misma conclusión; no rechazamos la hipótesis de homogeneidad de varianzas.

Valores atípicos

La forma más fácil y común de detectar valores atípicos es visualmente gracias a los diagramas de caja por grupos.

Para hembras y machos:

library(ggplot2)

# diagramas de caja por sexo
ggplot(dat) + aes(x = sex, y = body_mass_g) + geom_boxplot()

Para las tres especies:

# diagramas de caja por especie
ggplot(dat) + aes(x = species, y = body_mass_g) + geom_boxplot()

Hay, según lo definido por el criterio del rango intercuartil, dos valores atípicos para la especie Chinstrap. Sin embargo, estos puntos no son lo suficientemente extremos como para sesgar los resultados.

Por lo tanto, consideramos que se cumple la suposición de no tener valores atípicos significativos.

ANOVA de dos vías

Hemos demostrado que se cumplen todas las suposiciones, por lo que ahora podemos proceder a la implementación de la ANOVA de dos vías en R.

Esto nos permitirá responder las siguientes preguntas de investigación:

• Controlando la especie, ¿la masa corporal es significativamente diferente entre los dos sexos?
• Controlando el sexo, ¿la masa corporal es significativamente diferente para al menos una especie?
• ¿Es diferente la relación entre la especie y la masa corporal entre pingüinos hembra y macho?

Análisis preliminares

Antes de realizar cualquier prueba estadística, es una buena práctica hacer algunas estadísticas descriptivas para tener una primera visión general de los datos, y quizás tener una idea de los resultados que se esperan.

Esto se puede hacer a través de estadísticas descriptivas o gráficos.

Estadísticas descriptivas

Si queremos mantenerlo simple, podemos calcular solo la media para cada subgrupo:

# media por grupo
aggregate(body_mass_g ~ species + sex,  data = dat,  FUN = mean)

##     species    sex body_mass_g## 1    Adelie female    3368.836## 2 Chinstrap female    3527.206## 3    Gentoo female    4679.741## 4    Adelie   male    4043.493## 5 Chinstrap   male    3938.971## 6    Gentoo   male    5484.836

O, eventualmente, la media y la desviación estándar para cada subgrupo utilizando el paquete {dplyr}:

# media y desviación estándar por grupobiblioteca(dplyr)

group_by(dat, sexo, especie) %>%  summarise(    mean = round(mean(body_mass_g, na.rm = TRUE)),    sd = round(sd(body_mass_g, na.rm = TRUE))  )

## # A tibble: 8 × 4## # Grupos:   sexo [3]##   sexo    especie    mean    sd##   <fct>  <fct>     <dbl> <dbl>## 1 hembra Adelie     3369   269## 2 hembra Chinstrap  3527   285## 3 hembra Gentoo     4680   282## 4 macho  Adelie     4043   347## 5 macho  Chinstrap  3939   362## 6 macho  Gentoo     5485   313## 7 <NA>   Adelie     3540   477## 8 <NA>   Gentoo     4588   338

Gráficos

Si eres un lector frecuente del blog, sabes que me gusta dibujar gráficos para visualizar los datos antes de interpretar los resultados de una prueba.

El gráfico más apropiado cuando tenemos una variable cuantitativa y dos variables cualitativas es un diagrama de caja por grupo. Esto se puede hacer fácilmente con el paquete {ggplot2}:

# diagrama de caja por grupobiblioteca(ggplot2)

ggplot(dat) +  aes(x = especie, y = body_mass_g, fill = sexo) +  geom_boxplot()

Algunas observaciones están ausentes para el sexo, podemos eliminarlas para tener un gráfico más conciso:

dat %>%  filter(!is.na(sexo)) %>%  ggplot() +  aes(x = especie, y = body_mass_g, fill = sexo) +  geom_boxplot()

Note que también podríamos haber hecho el siguiente gráfico:

dat %>%  filter(!is.na(sexo)) %>%  ggplot() +  aes(x = sexo, y = body_mass_g, fill = especie) +  geom_boxplot()

Pero para un gráfico más legible, tiendo a preferir poner la variable con el menor número de niveles como color (que es de hecho el argumento fill en la capa aes()) y la variable con el mayor número de categorías en el eje x (es decir, el argumento x en la capa aes()).

A partir de las medias y los diagramas de caja por subgrupo, ya podemos ver que, en nuestra muestra:

• las pingüinas hembras tienden a tener una masa corporal más baja que los machos, y eso es cierto para todas las especies consideradas, y
• la masa corporal es mayor para los pingüinos Gentoo que para las otras dos especies.

Tenga en cuenta que estas conclusiones solo son válidas dentro de nuestra muestra. Para generalizar estas conclusiones a la población, necesitamos realizar el ANOVA de dos vías y verificar la significación de las variables explicativas. Este es el objetivo de la siguiente sección.

ANOVA de dos vías en R

Como se mencionó anteriormente, incluir un efecto de interacción en un ANOVA de dos vías no es obligatorio. Sin embargo, para evitar conclusiones erróneas, se recomienda primero verificar si la interacción es significativa o no, y en función de los resultados, incluirla o no.

Si la interacción no es significativa, es seguro eliminarla del modelo final. Por el contrario, si la interacción es significativa, debe incluirse en el modelo final que se utilizará para interpretar los resultados.

Por lo tanto, comenzamos con un modelo que incluye los dos efectos principales (es decir, sexo y especie) y la interacción:

# ANOVA de dos vías con interacción# guardar modelo mod <- aov(body_mass_g ~ sex * species,  data = dat)

# imprimir resultadosresumen(mod)

##              Df    Sum Sq  Mean Sq F value   Pr(>F)    ## sex           1  38878897 38878897 406.145  < 2e-16 ***## species       2 143401584 71700792 749.016  < 2e-16 ***## sex:species   2   1676557   838278   8.757 0.000197 ***## Residuals   327  31302628    95727                     ## ---## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 11 observaciones eliminadas debido a faltantes

La suma de cuadrados (columna Sum Sq) muestra que las especies explican una gran parte de la variabilidad de la masa corporal. Es el factor más importante en la explicación de esta variabilidad.

Los valores de p se muestran en la última columna de la salida anterior (Pr(>F)). A partir de estos valores de p, concluimos que, en el nivel de significancia del 5%:

• controlando la especie, la masa corporal es significativamente diferente entre los dos sexos,
• controlando el sexo, la masa corporal es significativamente diferente para al menos una especie, y
• la interacción entre el sexo y la especie (mostrada en la línea sex:species en la salida anterior) es significativa.

Entonces, a partir del efecto de interacción significativo que acabamos de ver, hemos visto que la relación entre la masa corporal y la especie es diferente entre machos y hembras. Dado que es significativo, tenemos que mantenerlo en el modelo y debemos interpretar los resultados de ese modelo.

Si, por el contrario, la interacción no fuera significativa (es decir, si el valor de p ≥ 0,05), habríamos eliminado este efecto de interacción del modelo. A modo ilustrativo, a continuación se muestra el código para una ANOVA de dos vías sin interacción, referida como un modelo aditivo:

# ANOVA de dos vías sin interacciónaov(body_mass_g ~ sex + species,  data = dat)

Para los lectores acostumbrados a realizar regresiones lineales en R, notarán que la estructura del código para una ANOVA de dos vías es de hecho similar:

• la fórmula es variable dependiente ~ variables independientes
• se utiliza el signo + para incluir variables independientes sin interacción
• se utiliza el signo * para incluir variables independientes con interacción

La semejanza con una regresión lineal no es una sorpresa porque una ANOVA de dos vías, como todas las ANOVA, es en realidad un modelo lineal.

Tenga en cuenta que el siguiente código también funciona y da los mismos resultados:

# método 2mod2 <- lm(body_mass_g ~ sex * species,  data = dat)

Anova(mod2)

## Tabla de Anova (Pruebas de tipo II)## ## Respuesta: body_mass_g##                Sum Sq  Df F value    Pr(>F)    ## sex          37090262   1 387.460 < 2.2e-16 ***## species     143401584   2 749.016 < 2.2e-16 ***## sex:species   1676557   2   8.757 0.0001973 ***## Residuals    31302628 327                      ## ---## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tenga en cuenta que la función aov() asume un diseño balanceado, lo que significa que tenemos tamaños de muestra iguales dentro de los niveles de nuestras variables de agrupamiento independientes.

Para diseños no balanceados, es decir, números desiguales de sujetos en cada subgrupo, se recomiendan los siguientes métodos:

• la ANOVA de tipo II cuando no hay ninguna interacción significativa, que se puede hacer en R con Anova(mod, type = "II"), 5 y
• la ANOVA de tipo III cuando hay una interacción significativa, que se puede hacer en R con Anova(mod, type = "III").

Esto está más allá del alcance del post y asumimos un diseño equilibrado aquí. Para el lector interesado, vea esta discusión detallada sobre ANOVA de tipo I, tipo II y tipo III.

Comparaciones por pares

A través de los dos efectos principales significativos, concluimos que:

• controlando la especie, la masa corporal es diferente entre hembras y machos, y
• controlando el sexo, la masa corporal es diferente para al menos una especie.

Si la masa corporal es diferente entre los dos sexos, dado que hay exactamente dos sexos, debe ser porque la masa corporal es significativamente diferente entre hembras y machos.

Si uno quiere saber qué sexo tiene la masa corporal más alta, se puede deducir de las medias y / o diagramas de caja por subgrupo. Aquí, está claro que los machos tienen una masa corporal significativamente más alta que las hembras.

Sin embargo, no es tan sencillo para la especie. Permítanme explicar por qué no es tan fácil como para los sexos.

Hay tres especies (Adelia, Barbijo y Papúa), por lo que hay 3 pares de especies:

1. Adelia y Barbijo
2. Adelia y Papúa
3. Barbijo y Papúa

Si la masa corporal es significativamente diferente para al menos una especie, podría ser que:

• la masa corporal es significativamente diferente entre Adelia y Barbijo pero no significativamente diferente entre Adelia y Papúa, y no significativamente diferente entre Barbijo y Papúa, o
• la masa corporal es significativamente diferente entre Adelia y Papúa pero no significativamente diferente entre Adelia y Barbijo, y no significativamente diferente entre Barbijo y Papúa, o
• la masa corporal es significativamente diferente entre Barbijo y Papúa pero no significativamente diferente entre Adelia y Barbijo, y no significativamente diferente entre Adelia y Papúa.

O también podría ser que:

• la masa corporal es significativamente diferente entre Adelia y Barbijo, y entre Adelia y Papúa, pero no significativamente diferente entre Barbijo y Papúa, o
• la masa corporal es significativamente diferente entre Adelia y Barbijo, y entre Barbijo y Papúa, pero no significativamente diferente entre Adelia y Papúa, o
• la masa corporal es significativamente diferente entre Barbijo y Papúa, y entre Adelia y Papúa, pero no significativamente diferente entre Adelia y Barbijo.

Por último, también podría ser que la masa corporal sea significativamente diferente entre todas las especies.

Como en una ANOVA de un factor, no podemos, en esta etapa, saber exactamente qué especie es diferente de cuál en términos de masa corporal. Para saber esto, necesitamos comparar cada especie de dos en dos gracias a pruebas post-hoc (también conocidas como comparaciones por pares).

Existen varias pruebas post-hoc, siendo la más común la Tukey HSD, que prueba todos los pares posibles de grupos. Como se mencionó anteriormente, esta prueba solo necesita hacerse en la variable de especies porque hay solo dos niveles para el sexo.

Al igual que en la ANOVA de un factor, la Tukey HSD se puede realizar en R de la siguiente manera:

# método 1TukeyHSD(mod, which = "species")

##   Tukey multiple comparisons of means##     95% family-wise confidence level## ## Fit: aov(formula = body_mass_g ~ sex * species, data = dat)## ## $species##                        diff       lwr       upr     p adj## Chinstrap-Adelie   26.92385  -80.0258  133.8735 0.8241288## Gentoo-Adelie    1377.65816 1287.6926 1467.6237 0.0000000## Gentoo-Chinstrap 1350.73431 1239.9964 1461.4722 0.0000000

O usando el paquete {multcomp}:

# método 2library(multcomp)

summary(glht(  aov(body_mass_g ~ sex + species,    data = dat  ),  linfct = mcp(species = "Tukey")))

## ##   Pruebas simultáneas para hipótesis lineales generales## ## Comparaciones múltiples de medias: Contrastes de Tukey## ## ## Ajustado a: aov(formula = body_mass_g ~ sex + species, data = dat)## Hipótesis lineales:##                         Estimación Error estándar  t valor   Pr (> | t |)  ## Chinstrap - Adelie == 0    26.92      46.48   0.579     0.83    ## Gentoo - Adelie == 0     1377.86      39.10  35.236   <1e-05 ***## Gentoo - Chinstrap == 0  1350.93      48.13  28.067   <1e-05 ***## ---## Códigos significativos:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## (Se informan valores p ajustados -- método de un solo paso)

O usando la función pairwise.t.test() utilizando el método de ajuste de valor p de su elección: 6

# método 3pairwise.t.test(dat$body_mass_g, dat$species,  p.adjust.method = "BH")

## ##  Comparaciones por pares utilizando pruebas t con SD agrupada ## ## datos:  dat$body_mass_g y dat$species ## ##           Adelie Chinstrap## Chinstrap 0.63   -        ## Gentoo    <2e-16 <2e-16   ## ## Método de ajuste de valor p: BH

Tenga en cuenta que al utilizar el segundo método, es el modelo sin interacción el que debe especificarse en la función glht(), incluso si la interacción es significativa. Además, no olvide reemplazar mod y species en mi código con el nombre de su modelo y el nombre de su variable independiente.

Ambos métodos dan los mismos resultados, es decir:

• la masa corporal no es significativamente diferente entre Chinstrap y Adelie (valor p ajustado = 0.83),
• la masa corporal es significativamente diferente entre Gentoo y Adelie (valor p ajustado < 0.001), y
• la masa corporal es significativamente diferente entre Gentoo y Chinstrap (valor p ajustado < 0.001).

Recuerde que son los valores p ajustados los que se informan, para evitar el problema de la prueba múltiple que ocurre al comparar varios pares de grupos.

Si desea comparar todas las combinaciones de grupos, puede hacerlo con la función TukeyHSD() y especificando la interacción en el argumento which:

# todas las combinaciones de sexo y especieTukeyHSD(mod,  which = "sex:species")

##   Comparaciones múltiples de Tukey de medias##     Nivel de confianza familiar del 95%## ## Ajustado a: aov(formula = body_mass_g ~ sex * species, data = dat)## ## $`sex:species`##                                      diff       lwr       upr     p adj## male:Adelie-female:Adelie        674.6575  527.8486  821.4664 0.0000000## female:Chinstrap-female:Adelie   158.3703  -25.7874  342.5279 0.1376213## male:Chinstrap-female:Adelie     570.1350  385.9773  754.2926 0.0000000## female:Gentoo-female:Adelie     1310.9058 1154.8934 1466.9181 0.0000000## male:Gentoo-female:Adelie       2116.0004 1962.1408 2269.8601 0.0000000## female:Chinstrap-male:Adelie    -516.2873 -700.4449 -332.1296 0.0000000## male:Chinstrap-male:Adelie      -104.5226 -288.6802   79.6351 0.5812048## female:Gentoo-male:Adelie        636.2482  480.2359  792.2606 0.0000000## male:Gentoo-male:Adelie         1441.3429 1287.4832 1595.2026 0.0000000## male:Chinstrap-female:Chinstrap  411.7647  196.6479  626.8815 0.0000012## female:Gentoo-female:Chinstrap  1152.5355  960.9603 1344.1107 0.0000000## male:Gentoo-female:Chinstrap    1957.6302 1767.8040 2147.4564 0.0000000## female:Gentoo-male:Chinstrap     740.7708  549.1956  932.3460 0.0000000## male:Gentoo-male:Chinstrap      1545.8655 1356.0392 1735.6917 0.0000000## male:Gentoo-female:Gentoo        805.0947  642.4300  967.7594 0.0000000

O bien con la función HSD.test() del paquete {agricolae}, que indica los subgrupos que no son significativamente diferentes entre sí con la misma letra:

library(agricolae)

HSD.test(mod,  trt = c("sex", "species"),  console = TRUE # print results)

## ## Estudio: mod ~ c("sex", "species")## ## Prueba HSD para body_mass_g ## ## Error cuadrático medio:  95726.69 ## ## sex:species,  medias## ##                  body_mass_g      std  r  Min  Max## female:Adelie       3368.836 269.3801 73 2850 3900## female:Chinstrap    3527.206 285.3339 34 2700 4150## female:Gentoo       4679.741 281.5783 58 3950 5200## male:Adelie         4043.493 346.8116 73 3325 4775## male:Chinstrap      3938.971 362.1376 34 3250 4800## male:Gentoo         5484.836 313.1586 61 4750 6300## ## Alpha: 0.05 ; DF Error: 327 ## Valor crítico del rango estudiantizado: 4.054126 ## ## Grupos según la probabilidad de diferencias de medias y nivel alfa (0.05)## ## Tratamientos con la misma letra no son significativamente diferentes.## ##                  body_mass_g grupos## male:Gentoo         5484.836      a## female:Gentoo       4679.741      b## male:Adelie         4043.493      c## male:Chinstrap      3938.971      c## female:Chinstrap    3527.206      d## female:Adelie       3368.836      d

Si tiene muchos grupos para comparar, puede ser más fácil interpretarlos mediante gráficos:

# establecer márgenes de eje para que las etiquetas no se cortenpar(mar = c(4.1, 13.5, 4.1, 2.1))

# crear intervalo de confianza para cada comparaciónplot(TukeyHSD(mod, which = "sex:species"),  las = 2 # rotar las etiquetas del eje x)

A partir de los resultados y el gráfico anterior, se concluye que todas las combinaciones de sexo y especie son significativamente diferentes, excepto entre la Chinstrap hembra y la Adelie hembra (valor-p = 0.138) y entre la Chinstrap macho y la Adelie macho (valor-p = 0.581).

Estos resultados, que por cierto están en línea con los diagramas de caja mostrados anteriormente y que serán confirmados con las visualizaciones a continuación, concluyen el ANOVA de dos vías en R.

Visualizaciones

Si desea visualizar los resultados de una manera diferente a lo que ya se ha presentado en los análisis preliminares, a continuación se presentan algunas ideas de gráficos útiles.

En primer lugar, con la media y el error estándar de la media por subgrupo utilizando la función allEffects() del paquete {effects}:

# método 1library(effects)

plot(allEffects(mod))

O utilizando el paquete {ggpubr}:

# método 2library(ggpubr)

ggline(subset(dat, !is.na(sex)), # eliminar nivel NA para sexo  x = "species",  y = "body_mass_g",  color = "sex",  add = c("mean_se") # agregar media y error estándar) +  labs(y = "Media de masa corporal (g)")

Alternativamente, utilizando {Rmisc} y {ggplot2}:

library(Rmisc)

# calcular la media y el error estándar de la media por subgrupo
summary_stat <- summarySE(dat,  measurevar = "body_mass_g",  groupvars = c("species", "sex"))

# trazar la media y el error estándar de la media
ggplot(  subset(summary_stat, !is.na(sex)), # eliminar el nivel NA para el sexo  
aes(x = species, y = body_mass_g, colour = sex)) +  geom_errorbar(aes(ymin = body_mass_g - se, ymax = body_mass_g + se), # agregar barras de error    
width = 0.1 # ancho de las barras de error  
) +  geom_point() +  labs(y = "Media de la masa corporal (g)")

En segundo lugar, si prefiere dibujar solo la media por subgrupo:

with(dat,  interaction.plot(species, sex, body_mass_g))

Por último, pero no menos importante, para aquellos de ustedes que están familiarizados con GraphPad, es probable que estén familiarizados con la representación de la media y las barras de error de la siguiente manera:

# trazar la media y el error estándar de la media como gráfico de barras
ggplot(  subset(summary_stat, !is.na(sex)), # eliminar el nivel NA para el sexo  
aes(x = species, y = body_mass_g, fill = sex)) +  geom_bar(position = position_dodge(), stat = "identity") +  geom_errorbar(aes(ymin = body_mass_g - se, ymax = body_mass_g + se), # agregar barras de error    
width = 0.25, # ancho de las barras de error    
position = position_dodge(.9)  ) +  labs(y = "Media de la masa corporal (g)")

Conclusión

En este post, comenzamos con algunos recordatorios sobre las diferentes pruebas que existen para comparar una variable cuantitativa entre grupos. Luego nos centramos en el ANOVA de dos vías, comenzando desde su objetivo e hipótesis hasta su implementación en R, junto con las interpretaciones y algunas visualizaciones. También mencionamos brevemente sus supuestos subyacentes y una prueba post-hoc para comparar todos los subgrupos.

Todo esto se ilustró con el conjunto de datos penguins disponible en el paquete {palmerpenguins}.

Gracias por leer.

Espero que este artículo te ayude a realizar un ANOVA de dos vías con tus datos.

Como siempre, si tienes una pregunta o sugerencia relacionada con el tema tratado en este artículo, por favor agrégala como un comentario para que otros lectores puedan beneficiarse de la discusión.

1. En teoría, un ANOVA de una vía también se puede utilizar para comparar 2 grupos, y no solo 3 o más. Sin embargo, en la práctica, a menudo se realiza una prueba t de Student para comparar 2 grupos, y un ANOVA de una vía para comparar 3 o más grupos. Las conclusiones obtenidas a través de una prueba t de Student para muestras independientes y un ANOVA de una vía con 2 grupos serán similares. ↩︎
2. Tenga en cuenta que si la suposición de normalidad no se cumple, se pueden aplicar muchas transformaciones para mejorarla, siendo la más común la transformación logarítmica (función log() en R). ↩︎
3. Tenga en cuenta que la prueba de Bartlett también es apropiada para probar la suposición de varianzas iguales. ↩︎
4. Un modelo aditivo supone que las 2 variables explicativas son independientes; no interactúan entre sí. ↩︎
5. Donde mod es el nombre de su modelo guardado. ↩︎
6. Aquí, utilizamos la corrección de Benjamini & Hochberg (1995), pero puede elegir entre varios métodos. Consulte ?p.adjust para más detalles. ↩︎

Artículos relacionados

• ANOVA en R
• Prueba t de Student en R y a mano: cómo comparar dos grupos en diferentes escenarios?
• Prueba de Wilcoxon de una muestra en R
• Prueba de Kruskal-Wallis, o la versión no paramétrica del ANOVA
• Prueba de hipótesis a mano

Publicado originalmente en https://statsandr.com el 19 de junio de 2023.

---

• Education
• programming
• Science
• Statistics
• Technology