9 Análisis de la Varianza (ANOVA)

En esta práctica se muestra cómo realizar un contraste de análisis de la varianza (ANOVA) para comparar las medias de una variable cuantitativa entre tres o más grupos de comparación, y también cómo realizar un contraste de ANOVA de medidas repetidas para comparar las medias de una variable cuantitativa entre tres o más momentos o condiciones experimentales. En particular, se verán los siguientes contrastes de hipótesis:

ANOVA de un factor
ANOVA de medidas repetidas
ANOVA de dos factores, con y sin interacción.

9.1 Ejercicios Resueltos

Para la realización de esta práctica se requieren los siguientes paquetes:

library(tidyverse)
# Incluye los siguientes paquetes:
# - readr: para la lectura de datos.
# - dplyr: para el preprocesamiento y manipulación de datos.
# - ggplot2: para la creación de gráficos.
library(broom) # para convertir las listas con los resúmenes de los modelos de regresión a formato organizado.
library(rstatix) # para realizar ANOVAs de medidas repetidas y mixtos.
library(lme4) # para construir modelos lineales mixtos.
library(tidymodels) # para ajustar de modelos.
library(ggpubr) # para realizar gráficos de diagnóstico de los modelos.
library(knitr) # para el formateo de tablas.

Versiones de R y de los paquetes utilizados

Contexto para Copilot

Si se va a usar Copilot para resolver los ejercicios se recomienda añadir el siguiente contexto al comienzo del script o al fichero r.instructions.md.

---
name: "Reglas de programación para R"
applyTo: "**/*.[rR], **/*.qmd"
---

## Reglas de programación para R

- Usar tidyverse para todas las tareas de preprocesamiento (dplyr/tidyr/readr/stringr/lubridate/forcats/purrr).
- Preferir tuberías con |> (o %>% si se usa en el archivo).
- Evitar la manipulación de datos en base R (merge, aggregate, subset) a menos que se solicite explícitamente.
- Mantener el código legible y usar verbos explícitos: select, mutate, summarise, group_by, pivot_longer/wider.
- Cuando no haya ambigüedad no antepongas el nombre del paquete a la función.
- En los nombres de variables y valores respetar el uso de mayúsculas y minúsculas (camelCase o snake_case) y ser consistente en todo el código.
- No generar nuevos data frames, ni modificar los asistentes, salvo cuando se indique explícitamente.
- Utiliza el prefijo `df_` para los nombres de los data frames.
- Organizar las salidas con la función tidy del paquete broom y mostrar solo las columnas relevantes (estimate, std.error, statistic, p.value).
- Mostrar los data frames o las tablas con la función kable del paquete knitr.
- Para gráficos usar ggplot2 y mantener un estilo limpio y profesional.

Ejercicio 9.1 Se quiere comparar la contaminación por dióxido de nitrógeno (NO2) en tres lugares distintos de una ciudad $A$, $B$ y $C$ y para ello se han medido las concentraciones de NO2 en cada lugar en una muestra aleatoria de días, obteniendo los resultados de la siguiente tabla.

A	19.8	21.1	28.2	22.3	22.5	28.9	23.8	16.9
B	31.6	32.8	41.1	36.8	37.0	35.6	32.2	43.9	37.5	25.2
C	37.2	30.2	26.6	31.7	26.8	28.6	29.2	22.9	38.0	33.9	26.2	40.5

Crear un conjunto de datos con los datos de la muestra.

Solución

df <- data.frame(nitrogeno = c(19.8, 21.1, 28.2, 22.3, 22.5, 28.9, 23.8, 16.9, 31.6, 32.8, 41.1, 36.8, 37.0, 35.6, 32.2, 43.9, 37.5, 25.2, 37.2, 30.2, 26.6, 31.7, 26.8, 28.6, 29.2, 22.9, 38.0, 33.9, 26.2, 40.5),
lugar = rep(c("A", "B", "C"), c(8, 10, 12)))

library(tidyverse)
df <- tibble(
    lugar = factor(c(rep("A", 8), rep("B", 10), rep("C", 12))),
    nitrogeno = c(19.8, 21.1, 28.2, 22.3, 22.5, 28.9, 23.8, 16.9, 31.6, 32.8, 41.1, 36.8, 37.0, 35.6, 32.2, 43.9, 37.5, 25.2, 37.2, 30.2, 26.6, 31.7, 26.8, 28.6, 29.2, 22.9, 38.0, 33.9, 26.2, 40.5)
)

Crear un conjunto de datos df con dos columnas: lugar y nitrogeno. La columna lugar contiene los valores A, B y C, repetidos 8, 10 y 12 veces respectivamente. La columna nitrogeno contiene los valores 19.8, 21.1, 28.2, 22.3, 22.5, 28.9, 23.8, 16.9, 31.6, 32.8, 41.1, 36.8, 37.0, 35.6, 32.2, 43.9, 37.5, 25.2, 37.2, 30.2, 26.6, 31.7, 26.8, 28.6, 29.2, 22.9, 38.0, 33.9, 26.2 y 40.5.

Dibujar el diagrama de cajas con los puntos correspondientes a las mediciones de cada lugar y sus medias. A la vista del diagrama, ¿crees que existen diferencias significativas entre los niveles de NO2 de los tres lugares?

Solución

media <- mean(df$nitrogeno)
boxplot(nitrogeno ~ lugar, data = df, fill = "lightblue", alpha = 0.5)
points(df$lugar, df$nitrogeno, col = "blue", pch = 16)
points(c(1, 2, 3), tapply(df$nitrogeno, df$lugar, mean), col = "red", pch = 16, cex = 1.5)
abline(h = media, col = "red")
text(x = 0.7, y = media + 0.5, label = "Media global", col = "red")

library(tidyverse)
media <- mean(df$nitrogeno)
df |> 
    ggplot(aes(x = lugar, y = nitrogeno, fill = lugar)) +
    geom_boxplot(alpha = 0.5) +
    geom_point() + 
    stat_summary(fun = mean, geom = "point", size = 3, color = "red") +
    geom_hline(yintercept = media, color = "red") +
    geom_text(aes(x = 0.7, y = media + 0.5, label = "Media global") )

Dibujar un diagrama de cajas con los puntos correspondientes a las mediciones de cada lugar y sus medias. Dibujar también la media global con una línea horizontal.

Como se puede apreciar en el diagrama existen diferencias muy claras entre los niveles de NO2 de los tres lugares, especialmente entre $B$ y $A$.

Realizar un contraste ANOVA para ver si hay diferencias estadísticamente significativas entre las concentraciones medias de NO2 de los tres lugares.

NotaAyuda

Antes de realizar el contraste de ANOVA hay que comprobar que se cumplen los supuestos del modelo ANOVA: normalidad, homocedasticidad e independencia de las observaciones.
TipSolución
Comprobemos en primer lugar si se cumplen los supuestos del modelo ANOVA.

9.1.0.1 Normalidad
Stats

Copilot

Para comprobar la normalidad de la variable dependiente se puede utilizar el test de normalidad de Shapiro-Wilk mediante la función shapiro.test.
Parámetros:

x: un vector numérico de datos a los que se les aplicará el test de normalidad.

library(knitr) library(broom) shapiro.test(df$nitrogeno) |> tidy() |> kable()

statistic p.value method

0.982475 0.8868395 Shapiro-Wilk normality test

Realizar un test de normalidad de Shapiro-Wilk para la columna nitrogeno del data frame df.
Como el p-valor es mayor que el nivel de significación $\alpha=0.05$, no podemos rechazar la hipótesis nula de normalidad de los datos.

9.1.0.2 Homocedasticidad
Stats

Copilot

Para comprobar la homocedasticidad de las varianzas de los grupos de comparación se puede utilizar el test de Barlett de homogeneidad de varianzas mediante la función bartlett.test.
Parámetros:

formula: una fórmula que especifica el modelo a ajustar, con la variable dependiente a la izquierda del símbolo ~ y la variable independiente a la derecha.

data: el data frame que contiene las variables utilizadas en la fórmula.

bartlett.test(nitrogeno ~ lugar, data = df) |> tidy() |> kable()

statistic p.value parameter method

0.7074 0.7020855 2 Bartlett test of homogeneity of variances

Realizar un test de homocedasticidad de Bartlett para la variable nitrogeno agrupada por la variable lugar del data frame df.
Como el p-valor es mayor que el nivel de significación $\alpha=0.05$, no podemos rechazar la hipótesis nula de homocedasticidad de las varianzas.

9.1.0.3 Independencia
Car

Copilot

Para comprobar la independencia de las observaciones se puede utilizar el test de independencia de Durbin-Watson mediante la función durbinWatsonTest del paquete car.
Parámetros:

model: un modelo de regresión lineal ajustado con la función lm o aov.

data: el data frame que contiene las variables utilizadas en el modelo.

library(car) durbinWatsonTest(aov(nitrogeno ~ lugar, data = df)) |> tidy() |> kable()

statistic p.value autocorrelation method alternative

2.291765 0.672 -0.2191414 Durbin-Watson Test two.sided

Realizar un test de independencia de Durbin-Watson para el modelo ANOVA con la variable dependiente nitrogeno y la variable independiente lugar del data frame df.
Como el p-valor es mayor que el nivel de significación $\alpha=0.05$, no podemos rechazar la hipótesis nula de independencia de las observaciones.

Así pues, todas las condiciones del modelo ANOVA se cumplen, por lo que podemos realizar el contraste de comparación de medias que es
\[\begin{align*} H_0 &: \mu_A = \mu_B = \mu_C \\ H_1 &: \mbox{Existen diferencias entre al menos dos medias} \end{align*}\]
Stats

Rstatix

Tidymodels

Copilot

Para realizar el contraste de ANOVA se puede utilizar la función aov del paquete stats.
Parámetros:

formula: una fórmula que especifica el modelo a ajustar, con la variable dependiente a la izquierda del símbolo ~ y la variable independiente a la derecha.

data: el data frame que contiene las variables utilizadas en la fórmula.

aov1 <- aov(nitrogeno ~ lugar, data = df) aov1 |> tidy() |> kable()

term df sumsq meansq statistic p.value

lugar 2 696.3036 348.15179 13.71644 7.75e-05

Residuals 27 685.3164 25.38209 NA NA

También podemos utilizar la función lm del paquete stats y luego aplicar la función anova al modelo ajustado.
Parámetros:

formula: una fórmula que especifica el modelo a ajustar, con la variable dependiente a la izquierda del símbolo ~ y la variable independiente a la derecha.

data: el data frame que contiene las variables utilizadas en la fórmula.

lm(nitrogeno ~ lugar, data = df) |> anova() |> tidy() |> kable()

term df sumsq meansq statistic p.value

lugar 2 696.3036 348.15179 13.71644 7.75e-05

Residuals 27 685.3164 25.38209 NA NA

Para realizar el contraste de ANOVA se puede utilizar la función anova_test del paquete rstatix.
Parámetros:

formula: una fórmula que especifica el modelo a ajustar, con la variable dependiente a la izquierda del símbolo ~ y la variable independiente a la derecha.

data: el data frame que contiene las variables utilizadas en la fórmula.

Después podemos obtener la tabla de ANOVA con la función get_anova_table.

library(rstatix) df |> anova_test(nitrogeno ~ lugar) |> get_anova_table() |> kable()

Effect DFn DFd F p p<.05 ges

lugar 2 27 13.716 7.75e-05 * 0.504

Para realizar un contraste de ANOVA podemos realizar un ajuste de modelo con la función lm del paquete tidymodels y luego aplicar la función anova al modelo ajustado.

library(tidymodels) lm(nitrogeno ~ lugar, data = df) |> anova() |> tidy() |> kable()

term df sumsq meansq statistic p.value

lugar 2 696.3036 348.15179 13.71644 7.75e-05

Residuals 27 685.3164 25.38209 NA NA

Realizar un contraste de ANOVA para comparar las medias de la variable nitrogeno entre los grupos definidos por la variable lugar del data frame df.
Como el p-valor del contraste es 0.00007 que es mucho menor que el nivel de significación $\alpha=0.05$, rechazamos la hipótesis nula y se concluye que existen diferencias estadísticamente significativas entre las concentraciones medias de NO2 de al menos dos lugares.
Analizar los residuos del modelo ANOVA.
TipSolución
Base

Copilot

El análisis de los residuos se suele realizar mediante la función plot, pasándole como argumento el modelo ANOVA.

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

Realizar un análisis de residuos del modelo ANOVA y mostrar los cuatro diagramas de diagnóstico: Residuals vs Fitted, Normal Q-Q, Scale-Location y Residuals vs Leverage.
El primer diagrama “Residuals vs Fitted” muestra si los residuos tienen tendencia lineal o no. En este caso, la línea roja que representa las medias es prácticamente horizontal por lo que se puede asumir que la tendencia es lineal.

El segundo diagrama “Normal Q-Q” muestra si los residuos siguen una distribución normal. En este caso los puntos se ajustan bastante bien a la línea recta, por lo que se puede asumir que los residuos siguen una distribución normal.

El tercer diagrama “Scale-Location” muestra si los residuos tienen una varianza constante (homocedasticidad). En este caso, la línea roja que representa la media de los residuos es prácticamente horizontal, por lo que se puede asumir que la varianza de los residuos es constante.

El cuarto diagrama “Residuals vs Leverage” muestra si hay observaciones influyentes en el modelo. En este caso, no hay observaciones que se salgan de la línea discontinua roja que representa el límite de influencia (distancia de Cook), por lo que no hay datos atípicos que sesguen el modelo.

statistic	p.value	method
0.982475	0.8868395	Shapiro-Wilk normality test

statistic	p.value	parameter	method
0.7074	0.7020855	2	Bartlett test of homogeneity of variances

statistic	p.value	autocorrelation	method	alternative
2.291765	0.672	-0.2191414	Durbin-Watson Test	two.sided

term	df	sumsq	meansq	statistic	p.value
lugar	2	696.3036	348.15179	13.71644	7.75e-05
Residuals	27	685.3164	25.38209	NA	NA

term	df	sumsq	meansq	statistic	p.value
lugar	2	696.3036	348.15179	13.71644	7.75e-05
Residuals	27	685.3164	25.38209	NA	NA

Effect	DFn	DFd	F	p	p<.05	ges
lugar	2	27	13.716	7.75e-05	*	0.504

term	df	sumsq	meansq	statistic	p.value
lugar	2	696.3036	348.15179	13.71644	7.75e-05
Residuals	27	685.3164	25.38209	NA	NA

Realizar un contraste post-hoc de comparación de las medias de N02 por pares. ¿Entre qué lugares existe una diferencia estadísticamente significativa en la concentración media de NO2?

Solución

Para realizar un contraste post-hoc de comparación de medias por pares podemos usar la función TukeyHSD del paquete stats.
Parámetros:

model: un modelo de ANOVA ajustado con la función aov o con lm.

TukeyHSD(aov(nitrogeno ~ lugar, data = df)) |> 
    tidy() |> 
    kable()

term	contrast	estimate	conf.low	conf.high	adj.p.value
lugar	B-A	12.432500	6.507278	18.3577222	0.0000513
lugar	C-A	8.045833	2.344286	13.7473810	0.0045194
lugar	C-B	-4.386667	-9.735193	0.9618592	0.1234729

Otra opción es utilizar la función pairwise.t.test del paquete stats que aplica correcciones a los p-valores.

Parámetros:

x: un vector numérico de datos a los que se les aplicará el test de comparación de medias por pares.
g: un factor que define los grupos de comparación.
p.adjust.method: un método de corrección de p-valores para los tests múltiples. En este caso, se puede utilizar el método “bonferroni”.

pairwise.t.test(df$nitrogeno, df$lugar, p.adjust.method = "bonferroni", conf.level = 0.95) |> 
    tidy() |> 
    kable()

group1	group2	p.value
B	A	0.0000531
C	A	0.0049147
C	B	0.1558055

Para realizar un contraste post-hoc de comparación de medias por pares se puede utilizar la función tukey_hsd del paquete rstatix.
Parámetros:

model: un modelo de ANOVA ajustado con la función aov o con lm.

df |> 
    tukey_hsd(nitrogeno ~ lugar) |>  
    kable()

term	group1	group2	estimate	conf.low	conf.high	p.adj	p.adj.signif
lugar	A	B	12.432500	6.507278	18.3577222	5.13e-05	****
lugar	A	C	8.045833	2.344286	13.7473810	4.52e-03	**
lugar	B	C	-4.386667	-9.735193	0.9618592	1.23e-01	ns

Realizar un contraste post-hoc de comparación de medias por pares para la variable nitrogeno agrupada por la variable lugar del data frame df, aplicando la corrección de Bonferroni o de Tukey a los p-valores.

Existe una diferencia muy significativa entre la concentración media de NO2 de los lugares $A$ y $B$ (p-valor 0.00005), y también entre los lugares $A$ y $C$ (p-valor 0.0049), pero no entre los lugares $B$ y $C$ (p-valor 0.1558).

Ejercicio 9.2 El fichero trigo.csv contiene información sobre la cosecha de trigo (en toneladas) obtenida para dos tipos de semillas y tres tipos de abonos ($A$, $B$ y $C$).

Crear conjunto de datos con los datos de la muestra a partir del fichero trigo.csv.

Solución

df <- read.csv("https://aprendeconalf.es/estadistica-practicas-r/datos/trigo.csv")
head(df)

semilla	parcela	abono	cosecha
1	1	A	4.73
2	2	A	4.76
1	3	A	4.17
2	4	A	5.23
1	1	A	4.39
2	2	A	5.15

library(tidyverse)
df <- read_csv("https://aprendeconalf.es/estadistica-practicas-r/datos/trigo.csv") |> 
    mutate(semilla = factor(semilla), abono = factor(abono))
head(df)

semilla	parcela	abono	cosecha
1	1	A	4.73
2	2	A	4.76
1	3	A	4.17
2	4	A	5.23
1	1	A	4.39
2	2	A	5.15

Crear un conjunto de datos df a partir del fichero csv con la siguiente URL: https://aprendeconalf.es/estadistica-practicas-r/datos/trigo.csv. La variable semilla y la variable abono deben ser factores.

Realizar un contraste ANOVA para ver si la cosecha de trigo depende del abono. ¿Cuál es la suma de cuadrados de los residuos del modelo?
TipSolución
Tenemos que realizar el contraste

\[\begin{align*} H_0 &: \mu_A = \mu_B = \mu_C \\ H_1 &: \mbox{Existen diferencias entre al menos dos medias} \end{align*}\]
Stats

Copilot

Con la función aov

library(knitr) library(broom) aov(cosecha ~ abono, data = df) |> tidy() |> kable()

term df sumsq meansq statistic p.value

abono 2 0.9097896 0.4548948 4.758202 0.0107787

Residuals 93 8.8910094 0.0956023 NA NA

Con la función lm

lm(cosecha ~ abono, data = df) |> anova() |> tidy() |> kable()

term df sumsq meansq statistic p.value

abono 2 0.9097896 0.4548948 4.758202 0.0107787

Residuals 93 8.8910094 0.0956023 NA NA

Realizar un contraste de ANOVA para comparar las medias de la variable cosecha entre los grupos definidos por la variable abono del data frame df.
Como el p-valor del contraste es 0.0108 que es menor que el nivel de significación $\alpha=0.05$, rechazamos la hipótesis nula y se concluye que existen diferencias estadísticamente significativas entre las medias de las cosechas de trigo de los distintos abonos, y por tanto, la cosecha depende del tipo de abono.

La suma de cuadrados de los residuos del modelo es 8.891.
Realizar un contraste ANOVA para ver si la cosecha de trigo depende del abono y del tipo de semilla. ¿Cuánto se reduce la suma de cuadrados de los residuos al incluir el tipo de semilla en el modelo?
TipSolución
Tenemos que realizar el contraste ANOVA de dos factores sin interacción. Para realizar un contraste ANOVA de dos factores sin interacción se puede utilizar tanto la función aov como la función lm pero incluyendo la fórmula del modelo vd ~ f1 + f2, donde vd es la variable dependiente, f1 es el primer factor y f2 el segundo.
Stats

Copilot

Para realizar un contraste ANOVA de dos factores sin interacción se puede utilizar tanto la función aov como la función lm pero incluyendo la fórmula del modelo vd ~ f1 + f2, donde vd es la variable dependiente, f1 es el primer factor y f2 el segundo.

Con la función aov

library(knitr) library(broom) aov(cosecha ~ semilla + abono, data = df) |> tidy() |> kable()

term df sumsq meansq statistic p.value

semilla 1 2.5382510 2.5382510 36.75869 0.0000000

abono 2 0.9097896 0.4548948 6.58774 0.0021192

Residuals 92 6.3527583 0.0690517 NA NA

Con la función lm

lm(cosecha ~ semilla + abono, data = df) |> anova() |> tidy() |> kable()

term df sumsq meansq statistic p.value

semilla 1 2.5382510 2.5382510 36.75869 0.0000000

abono 2 0.9097896 0.4548948 6.58774 0.0021192

Residuals 92 6.3527583 0.0690517 NA NA

Realizar un contraste de ANOVA de dos factores sin interacción para comparar las medias de la variable cosecha entre los grupos definidos por la variable semilla y la variable abono del data frame df.
Como tanto el p-valor correspondiente a la semilla, que es prácticamente 0, como el correspondiente al abono, que es 0.002, son menores que el nivel de significación $\alpha=0.05$, rechazamos la hipótesis nula y se concluye que la cosecha depende tanto del tipo de semilla como del tipo de abono.

Al incluir el tipo de semilla en el modelo, la suma de cuadrados de los residuos se reduce de 8.8689 a 6.3528, por lo que este modelo explica mejor la variabilidad de la cosecha.

term	df	sumsq	meansq	statistic	p.value
abono	2	0.9097896	0.4548948	4.758202	0.0107787
Residuals	93	8.8910094	0.0956023	NA	NA

term	df	sumsq	meansq	statistic	p.value
abono	2	0.9097896	0.4548948	4.758202	0.0107787
Residuals	93	8.8910094	0.0956023	NA	NA

term	df	sumsq	meansq	statistic	p.value
semilla	1	2.5382510	2.5382510	36.75869	0.0000000
abono	2	0.9097896	0.4548948	6.58774	0.0021192
Residuals	92	6.3527583	0.0690517	NA	NA

term	df	sumsq	meansq	statistic	p.value
semilla	1	2.5382510	2.5382510	36.75869	0.0000000
abono	2	0.9097896	0.4548948	6.58774	0.0021192
Residuals	92	6.3527583	0.0690517	NA	NA

Incluir en el modelo anterior también la interacción entre el tipo de semilla y el tipo de abono. ¿Es significativa la interacción entre los dos factores? ¿Cuánto se reduce la suma de cuadrados de los residuos al incluir la interacción en el modelo?

Solución

Tenemos que realizar el contraste ANOVA de dos factores con interacción. Para realizar un contraste ANOVA de dos factores con interacción se puede utilizar tanto la función aov como la función lm pero incluyendo la fórmula del modelo vd ~ f1 * f2, donde vd es la variable dependiente, f1 es el primer factor y f2 el segundo.

Con la función aov

library(knitr)
library(broom)
aov(cosecha ~ semilla * abono, data = df) |> 
    tidy() |> 
    kable()

term	df	sumsq	meansq	statistic	p.value
semilla	1	2.5382510	2.5382510	38.464056	0.0000000
abono	2	0.9097896	0.4548948	6.893368	0.0016388
semilla:abono	2	0.4136396	0.2068198	3.134098	0.0483265
Residuals	90	5.9391188	0.0659902	NA	NA

Con la función lm

lm(cosecha ~ semilla * abono, data = df) |> 
    anova() |> 
    tidy() |> 
    kable()

term	df	sumsq	meansq	statistic	p.value
semilla	1	2.5382510	2.5382510	38.464056	0.0000000
abono	2	0.9097896	0.4548948	6.893368	0.0016388
semilla:abono	2	0.4136396	0.2068198	3.134098	0.0483265
Residuals	90	5.9391188	0.0659902	NA	NA

Realizar un contraste de ANOVA de dos factores con interacción para comparar las medias de la variable cosecha entre los grupos definidos por la variable semilla y la variable abono del data frame df, incluyendo también la interacción entre ambos factores.

Al igual que antes, tanto el tipo de semilla como el tipo de abono son significativos en el modelo, pero además lo es su interacción ya que el p-valor correspondiente es 0.0483 que es menor que el nivel de significación $\alpha = 0.05$.

Al incluir la interacción entre el tipo de semilla y el tipo de abono en el modelo, la suma de cuadrados de los residuos se reduce de 6.3528 a 5.9391, por lo que este modelo explica todavía mejor la variabilidad de la cosecha.

Realizar un diagrama de interacción entre los factores con las medias de los distintos grupos experimentales.
TipSolución
Base

Tidyverse

Copilot

Para realizar un diagrama de interacción entre los factores con las medias de los distintos grupos experimentales se puede utilizar la función interaction.plot del paquete stats.
Parámetros:

x.factor: un factor que define los grupos en el eje x.

trace.factor: un factor que define los grupos que se representan con diferentes colores o símbolos.

response: la variable a representar.

fun: una función que se aplica a los grupos de la variable respuesta para obtener el estadístico resumen. Por defecto es mean.

interaction.plot(x.factor = df$semilla, trace.factor = df$abono, response = df$cosecha, fun = mean, type = "b", col = c("red", "blue", "green"), pch = c(16, 17, 18), xlab = "Tipo de semilla", ylab = "Media de la cosecha", legend = TRUE)

df |> group_by(semilla, abono) |> summarise(media = mean(cosecha)) |> ggplot(aes(x = semilla, y = media, color = abono)) + geom_line(aes(group = abono)) + geom_point()

Realizar un diagrama de interacción entre los factores semilla y abono con las medias de la variable cosecha para cada combinación de semilla y abono del data frame df.
Como se aprecia en el diagrama, las líneas que unen las medias no son paralelas, lo que significa que hay interacción entre el tipo de semilla y el tipo de abono en la cosecha de trigo. ```

Realizar un contraste post-hoc para comparar las medias de los distintos grupos experimentales por pares. ¿Qué combinación de semilla y abono ofrece el mayor rendimiento de la cosecha?

Solución

Primero comparamos las medias de los grupos según el tipo de abono.

TukeyHSD(aov(cosecha ~ abono, data = df)) |> 
    tidy() |> 
    kable()

term	contrast	estimate	conf.low	conf.high	adj.p.value
abono	B-A	0.0946875	-0.0894244	0.2787994	0.4415513
abono	C-A	0.2368750	0.0527631	0.4209869	0.0079604
abono	C-B	0.1421875	-0.0419244	0.3262994	0.1625203

df |>
    tukey_hsd(cosecha ~ abono) |>  
    kable()

term	group1	group2	estimate	conf.low	conf.high	p.adj	p.adj.signif
abono	A	B	0.0946875	-0.0894244	0.2787994	0.44200	ns
abono	A	C	0.2368750	0.0527631	0.4209869	0.00796	**
abono	B	C	0.1421875	-0.0419244	0.3262994	0.16300	ns

Realizar un contraste post-hoc de comparación de medias por pares para la variable cosecha agrupada por la variable abono del data frame df, aplicando la corrección de Tukey a los p-valores.

Se observa que solo hay una diferencia estadísticamente significativa entre las medias de las cosechas con los tipos de abono $A$ y $C$.

Ahora comparamos las medias de los grupos según el tipo de semilla.

TukeyHSD(aov(cosecha ~ abono, data = df)) |> 
    tidy() |> 
    kable()

term	contrast	estimate	conf.low	conf.high	adj.p.value
abono	B-A	0.0946875	-0.0894244	0.2787994	0.4415513
abono	C-A	0.2368750	0.0527631	0.4209869	0.0079604
abono	C-B	0.1421875	-0.0419244	0.3262994	0.1625203

df |>
    tukey_hsd(cosecha ~ semilla) |>  
    kable()

term	group1	group2	null.value	estimate	conf.low	conf.high	p.adj	p.adj.signif
semilla	1	2	0	0.3252083	0.2125535	0.4378632	1e-07	****

Realizar un contraste post-hoc de comparación de medias por pares para la variable cosecha agrupada por la variable semilla del data frame df, aplicando la corrección de Tukey a los p-valores.

Como solo hay dos grupos, solo hay una comparación posible, y se concluye hay diferencias significativas entre las medias de las cosechas con los tipos de semillas 1 y 2.

Finalmente, comparamos las medias de los grupos que surgen de la interacción del tipo de semilla con el tipo de abono.

TukeyHSD(aov(cosecha ~ semilla : abono, data = df)) |> 
    tidy() |> 
    arrange(adj.p.value) |> 
    kable()

term	contrast	estimate	conf.low	conf.high	adj.p.value
semilla:abono	2:C-1:A	0.566250	0.3017713	0.8307287	0.0000002
semilla:abono	2:B-1:A	0.498125	0.2336463	0.7626037	0.0000055
semilla:abono	2:A-1:A	0.490000	0.2255213	0.7544787	0.0000080
semilla:abono	1:C-1:A	0.397500	0.1330213	0.6619787	0.0004539
semilla:abono	2:C-1:B	0.385000	0.1205213	0.6494787	0.0007519
semilla:abono	2:B-1:B	0.316875	0.0523963	0.5813537	0.0095145
semilla:abono	1:B-2:A	-0.308750	-0.5732287	-0.0442713	0.0125345
semilla:abono	1:C-1:B	0.216250	-0.0482287	0.4807287	0.1740364
semilla:abono	1:B-1:A	0.181250	-0.0832287	0.4457287	0.3529759
semilla:abono	2:C-1:C	0.168750	-0.0957287	0.4332287	0.4346043
semilla:abono	1:C-2:B	-0.100625	-0.3651037	0.1638537	0.8769068
semilla:abono	1:C-2:A	-0.092500	-0.3569787	0.1719787	0.9106159
semilla:abono	2:C-2:A	0.076250	-0.1882287	0.3407287	0.9592452
semilla:abono	2:C-2:B	0.068125	-0.1963537	0.3326037	0.9748943
semilla:abono	2:B-2:A	0.008125	-0.2563537	0.2726037	0.9999992

df |>
    tukey_hsd(cosecha ~ semilla : abono) |> 
    arrange(p.adj) |>
    kable()

term	group1	group2	estimate	conf.low	conf.high	p.adj	p.adj.signif
semilla:abono	1:A	2:C	0.566250	0.3017713	0.8307287	2.00e-07	****
semilla:abono	1:A	2:B	0.498125	0.2336463	0.7626037	5.50e-06	****
semilla:abono	1:A	2:A	0.490000	0.2255213	0.7544787	8.00e-06	****
semilla:abono	1:A	1:C	0.397500	0.1330213	0.6619787	4.54e-04	***
semilla:abono	1:B	2:C	0.385000	0.1205213	0.6494787	7.52e-04	***
semilla:abono	1:B	2:B	0.316875	0.0523963	0.5813537	9.51e-03	**
semilla:abono	2:A	1:B	-0.308750	-0.5732287	-0.0442713	1.25e-02	*
semilla:abono	1:B	1:C	0.216250	-0.0482287	0.4807287	1.74e-01	ns
semilla:abono	1:A	1:B	0.181250	-0.0832287	0.4457287	3.53e-01	ns
semilla:abono	1:C	2:C	0.168750	-0.0957287	0.4332287	4.35e-01	ns
semilla:abono	2:B	1:C	-0.100625	-0.3651037	0.1638537	8.77e-01	ns
semilla:abono	2:A	1:C	-0.092500	-0.3569787	0.1719787	9.11e-01	ns
semilla:abono	2:A	2:C	0.076250	-0.1882287	0.3407287	9.59e-01	ns
semilla:abono	2:B	2:C	0.068125	-0.1963537	0.3326037	9.75e-01	ns
semilla:abono	2:A	2:B	0.008125	-0.2563537	0.2726037	1.00e+00	ns

Realizar un contraste post-hoc de comparación de medias por pares para la variable cosecha agrupada por la interacción de las variables semilla y abono del data frame df, aplicando la corrección de Tukey a los p-valores, y ordenando las comparaciones por el p-valor ajustado de menor a mayor.

De más a menos significativas, se observa que hay diferencias estadísticamente significativas entre las medias de las cosechas de los grupos 2:C (semilla 2 y abono C) y 1:A (semilla 1 y abono A), entre los grupos 2:B y 1:A, entre los grupos 2:A y 1:A, entre los grupos 2:C y 1:B, entre los grupos 1:C y 1:A, entre los grupos 2:B y 1:B y ente los grupos 1:B y 2:A. Entre el resto de grupos no hay diferencias estadísticamente significativas.

A la vista del diagrama del apartado anterior, se observa que el mayor rendimiento se obtiene con el tipo de semilla 2 y el abono $C$, seguido del tipo de semilla 2 y el abono $B$, aunque la diferencia entre las medias de estos dos grupos no es significativa.

Ejercicio 9.3 El fichero cata-vinos.csv contiene información la puntuación obtenida por unos vinos en una cata realizada por 15 expertos.

Crear conjunto de datos con los datos de la muestra a partir del fichero cata-vinos.csv.

Solución

library(tidyverse)
df <- read_csv("datos/cata-vinos.csv")
head(df)

vinoA	vinoB	vinoC
7.8	8.0	7.8
7.9	7.6	7.6
8.6	6.9	8.0
8.0	7.7	8.0
8.1	7.4	8.0
8.7	7.2	7.9

library(tidyverse)
df <- read_csv("datos/cata-vinos.csv")
head(df)

vinoA	vinoB	vinoC
7.8	8.0	7.8
7.9	7.6	7.6
8.6	6.9	8.0
8.0	7.7	8.0
8.1	7.4	8.0
8.7	7.2	7.9

Crear un conjunto de datos df a partir del fichero csv con la siguiente URL: https://aprendeconalf.es/estadistica-practicas-r/datos/cata-vinos.csv.

Dibujar un diagrama de líneas de las puntuaciones de cada catador.

Solución

Para dibujar un diagrama de líneas de las puntuaciones de cada catador para cada vino, se puede utilizar la función matplot del paquete base.
Parámetros:

y: un vector o matriz de valores para el eje y. En este caso, se puede utilizar la transpuesta del data frame df para que cada fila corresponda a un catador y cada columna a un vino.

matplot(t(df), type = "l", xaxt = "n", xlab = "Vino", ylab = "Puntuación",
col = seq_len(nrow(df)), lty = 1)
axis(1, at = seq_along(names(df)), labels = names(df))

df_largo <- df |> 
    # Generamos un identificador para cada individuo.
    mutate(id = 1:nrow(df)) |>
    # Convertimos el data frame a formato largo.
    pivot_longer(cols = -id, names_to = "Vino", values_to = "Puntuación") |> 
    # Convertimos el id y el momento en factores.
    mutate_at(c("id", "Vino"), as.factor)

# Dibujamos el diagrama de líneas.
df_largo |> ggplot(aes(x = Vino, y = Puntuación, group = id, color = id)) +
    geom_line()

Dibujar un diagrama de líneas de las puntuaciones de cada catador para cada vino. El eje x debe representar los vinos, el eje y las puntuaciones, y cada línea debe representar a un catador.

A primera vista, parece que las puntuaciones del vino B son más bajas.

Realizar un contraste ANOVA para ver si hay diferencias estadísticamente significativas entre las medias del peso en los tres momentos.

Solución

Tenemos que realizar el contraste ANOVA de medidas repetidas ya que los tres vinos se catan por las mismas personas.

Para realizar un contraste ANOVA de medidas repetidas se puede utilizar la función aov del paquete Stats introduciendo en el modelo el termino de error Error(id/Vino).

modelo <- aov(Puntuación ~ Vino + Error(id/Vino), data = df_largo) 
modelo  |>  
    tidy() |> 
    kable()

stratum	term	df	sumsq	meansq	statistic	p.value
id	Residuals	14	0.9831111	0.0702222	NA	NA
id:Vino	Vino	2	3.0484444	1.5242222	12.88767	0.0001077
id:Vino	Residuals	28	3.3115556	0.1182698	NA	NA

Para realizar un contraste de ANOVA de medidas repetidas se puede utilizar la función anova_test del paquete rstatix, introduciendo la variable dependiente en el argumento dv, la variable identificadora en el argumento wid y la variable de medidas repetidas en el argumento within.

library(knitr)
library(broom)
library(rstatix)
anova <- anova_test(data = df_largo, dv = Puntuación, wid = id, within = Vino)
get_anova_table(anova) |>
    kable()

Effect	DFn	DFd	F	p	p<.05	ges
Vino	2	28	12.888	0.000108	*	0.415

Para realizar un contraste de ANOVA de medidas repetidas se puede ajustar un modelo lineal mixto con la función lmer del paquete lme4 introduciendo la variable dependiente, la variable de medidas repetidas como variable independiente y el identificador como efecto aleatorio (1|id), y luego aplicar la función anova al modelo ajustado.

library(lme4)
lmer(Puntuación ~ Vino + (1|id), data = df_largo) |> 
anova(data.frame(ames = NULL, check.rows = FALSE, check.names = TRUE, stringsAsFactors = default.stringsAsFactors())) |>
    tidy() |> 
    kable()

term	npar	sumsq	meansq	statistic
Vino	2	3.048444	1.524222	14.90624

Realizar un contraste de ANOVA de medidas repetidas para comparar las medias de la variable Puntuación entre los grupos definidos por la variable Vino del data frame df_largo.

Como el p-valor del contraste es 0.0001 < 0.05 podemos rechazar la hipótesis nula de igualdad de medias y concluimos que existen diferencias significativas entre las puntuaciones medias de al menos dos vinos.

Realizar un contraste de comparación de medias por pares para ver entre qué vinos existen diferencias estadísticamente significativas en la puntuación media.

Solución

Para realizar un contraste post-hoc de comparación de medias por pares para poblaciones pareadas podemos usar la función pairwise.t.test del paquete stats que aplica la corrección de Bonferroni a los p-valores.
Parámetros:

x: un vector numérico de datos a los que se les aplicará el test de comparación de medias por pares.
g: un factor que define los grupos de comparación.
p.adjust.method: un método de corrección de p-valores para los tests múltiples. En este caso, se puede utilizar el método “bonferroni”.
conf.level: el nivel de confianza para los intervalos de confianza de las diferencias de medias. Por defecto es 0.95.
paired: un valor lógico que indica si se trata de poblaciones pareadas. En este caso, se debe establecer a TRUE.

pairwise.t.test(df_largo$Puntuación, df_largo$Vino, paired = TRUE, p.adjust.method = "bonferroni", conf.level = 0.95) |> 
    tidy() |> 
    kable()

group1	group2	p.value
vinoB	vinoA	0.0020015
vinoC	vinoA	0.1286003
vinoC	vinoB	0.0154065

Para realizar un contraste post-hoc de comparación de medias por pares para poblaciones pareadas se puede utilizar la función pairwise_t_test del paquete Rstatix.
Parámetros:

formula: una fórmula que especifica el modelo a ajustar, con la variable dependiente a la izquierda del símbolo ~ y la variable independiente a la derecha.
data: el data frame que contiene las variables utilizadas en la fórmula.
p.adjust.method: un método de corrección de p-valores para los tests múltiples. En este caso, se puede utilizar el método “bonferroni”.
paired: un valor lógico que indica si se trata de poblaciones pareadas. En este caso, se debe establecer a TRUE.

pwc <- df_largo |> 
    pairwise_t_test(Puntuación ~ Vino, paired = TRUE, p.adjust.method = "bonferroni")
pwc |> 
    kable()

.y.	group1	group2	n1	n2	statistic	df	p	p.adj	p.adj.signif
Puntuación	vinoA	vinoB	15	15	4.349100	14	0.000667	0.002	**
Puntuación	vinoA	vinoC	15	15	2.227053	14	0.043000	0.129	ns
Puntuación	vinoB	vinoC	15	15	-3.312258	14	0.005000	0.015	*

Realizar un contraste post-hoc de comparación de medias por pares para la variable Puntuación agrupada por la variable Vino del data frame df_largo, aplicando la corrección de Bonferroni a los p-valores, y teniendo en cuenta que se trata de poblaciones pareadas.

A la vista de los p-valores de la comparación por pares existe una diferencia estadísticamente muy significativa entre el la puntuación media de los vinos A y B (p-valor 0.002 < 0.01), y una diferencia estadísticamente significativa entre las puntuaciones medias de los vinos B y C (p-valor 0.015 < 0.05), pero no existe una diferencia estadísticamente significativa entre la puntuación media de los vinos A y C (p-valor 0.129 > 0.05).

Dibujar los diagramas de cajas de las puntuaciones de los tres vinos, mostrando las diferencias significativas entre los vinos.
TipSolución
Ggpubr

Copilot

Para dibujar los diagramas de cajas de las puntuaciones de los tres vinos, mostrando las diferencias significativas entre los vinos, se puede utilizar la función ggboxplot del paquete ggpubr junto con la función stat_pvalue_manual para añadir las anotaciones de las comparaciones por pares.
Parámetros de ggboxplot:

data: el data frame que contiene las variables a representar.

x: la variable que se representa en el eje x. En este caso, es la variable Vino.

y: la variable que se representa en el eje y. En este caso, es la variable Puntuación.

fill: la variable que se utiliza para rellenar las cajas. En este caso, es la variable Vino.

library(ggpubr) pwc <- pwc |> add_xy_position(x = "Vino") ggboxplot(df_largo, x = "Vino", y = "Puntuación", fill = "Vino", alpha = 0.5) + stat_pvalue_manual(pwc) + labs( subtitle = get_test_label(anova, detailed = TRUE), caption = get_pwc_label(pwc) )

Dibujar los diagramas de cajas de las puntuaciones de los tres vinos, mostrando las diferencias significativas entre los vinos, utilizando la función ggboxplot del paquete ggpubr y añadiendo las anotaciones de las comparaciones por pares con la función stat_pvalue_manual.

9.2 Ejercicios propuestos

Ejercicio 9.4 La tabla siguiente contiene las notas medias en la prueba de acceso a la universidad (EVAU) de una muestra de alumnos de los cinco institutos de una ciudad.

\[ \begin{array}{ccccc} \hline \mbox{Instituto 1} & \mbox{Instituto 2} & \mbox{Instituto 3} & \mbox{Instituto 4} & \mbox{Instituto 5} \\ 5.5 & 6.1 & 4.9 & 3.2 & 6.7 \\ 5.2 & 7.2 & 5.5 & 3.3 & 5.8 \\ 5.9 & 5.5 & 6.1 & 5.5 & 5.4 \\ 7.1 & 6.7 & 6.1 & 5.7 & 5.5 \\ 6.2 & 7.6 & 6.2 & 6.0 & 4.9 \\ 5.9 & 5.9 & 6.4 & 6.1 & 6.2 \\ 5.3 & 8.1 & 6.9 & 4.7 & 6.1 \\ 6.2 & 8.3 & 4.5 & 5.1 & 7.0 \\ \hline \end{array} \]

Dibujar el diagrama de cajas y puntos de las notas de cada centro con sus respectivas medias. ¿Se observan diferencias entre los centros en el diagrama?
Realizar un contraste ANOVA para ver si hay diferencias estadísticamente significativas entre las notas medias de los cinco centros.
Realizar un contraste post-hoc de comparación de las medias de las notas de los cinco centros por pares. ¿Entre qué centros existen diferencias estadísticamente significativas en la nota media?

9.1 Ejercicios Resueltos

9.1.0.1 Normalidad

9.1.0.2 Homocedasticidad

9.1.0.3 Independencia

9.2 Ejercicios propuestos