Problema de ANOVA: Frecuencia de erupciones según el tipo de volcán

Un equipo de vulcanólogos quiere determinar si la frecuencia de erupciones (en años) difiere significativamente entre tres tipos de volcán: Cono de Ceniza, Escudo y Estratovolcán. Se analizaron datos de 100 volcanes registrados en la base de datos. ¿Existen diferencias significativas en la frecuencia promedio de erupciones entre estos tipos de volcán?

Datos: A continuación se muestran las frecuencias de erupción (en años) para cada tipo de volcán:

Cono de Ceniza: 459.64, 414.10, 188.55, 242.81, 4.70, 145.87, 221.59, 194.37, 375.13, 153.44, 234.07, 204.24, 47.37, 175.63, 67.05, 266.83, 397.00, 497.67, 152.78, 356.48, 457.94, 252.13

Escudo: 447.19, 202.82, 384.74, 10.13, 228.65, 139.29, 456.18, 221.19, 484.79, 30.80, 103.65, 449.01, 14.01, 389.79, 265.44, 307.21

Estratovolcán: 117.16, 410.68, 120.53, 186.69, 228.99, 430.33, 209.74, 257.05, 328.95, 188.25, 279.34, 345.65, 358.95, 295.22, 306.90, 382.50

Preguntas:

Realiza un análisis de ANOVA para determinar si hay diferencias significativas en la frecuencia de erupciones entre los tres tipos de volcán.

Si el resultado es significativo, realiza una prueba post-hoc (como Tukey) para identificar entre qué grupos existen diferencias.

Interpreta los resultados en el contexto del estudio.

library(readxl)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(agricolae)
library(ggpubr)

# 1. Cargar la base de datos 
datos_volcanes <- read_excel("Base_Datos_Vulcanologia.xlsx", sheet = "Vulcanologia")

# Ver las primeras filas para comprobar que se cargó correctamente
head(datos_volcanes)
## # A tibble: 6 × 7
##   ID_Volcan Nombre_Volcan Tipo_Volcan    Altura_m Frecuencia_Erupcion_Anios
##       <dbl> <chr>         <chr>             <dbl>                     <dbl>
## 1         1 Cotopaxi      Cono de Ceniza    2545.                      460.
## 2         2 Cotopaxi      Escudo            5952.                      447.
## 3         3 Galeras       Escudo            1574.                      349.
## 4         4 Etna          Estratovolcán     5913.                      117.
## 5         5 Etna          Estratovolcán     1555.                      411.
## 6         6 Galeras       Escudo            5179.                      203.
## # ℹ 2 more variables: Emision_SO2_Toneladas <dbl>, Temperatura_Lava_C <dbl>
# 3. Preparación de datos
datos_analisis <- datos_volcanes %>%
  select(Tipo_Volcan, Frecuencia_Erupcion_Anios) %>%
  filter(Tipo_Volcan %in% c("Cono de Ceniza", "Escudo", "Estratovolcán")) %>%
  mutate(Tipo_Volcan = as.factor(Tipo_Volcan))

# 4. Análisis exploratorio
# Resumen estadístico
resumen_stats <- datos_analisis %>%
  group_by(Tipo_Volcan) %>%
  summarise(
    n = n(),
    Media = mean(Frecuencia_Erupcion_Anios),
    Desviacion = sd(Frecuencia_Erupcion_Anios),
    Minimo = min(Frecuencia_Erupcion_Anios),
    Maximo = max(Frecuencia_Erupcion_Anios)
  )

print(resumen_stats)
## # A tibble: 3 × 6
##   Tipo_Volcan        n Media Desviacion Minimo Maximo
##   <fct>          <int> <dbl>      <dbl>  <dbl>  <dbl>
## 1 Cono de Ceniza    30  230.       151.   4.70   498.
## 2 Escudo            34  239.       161.   3.96   485.
## 3 Estratovolcán     36  196.       135.   1.81   430.
# Visualización
boxplot <- ggplot(datos_analisis, aes(x = Tipo_Volcan, y = Frecuencia_Erupcion_Anios, fill = Tipo_Volcan)) +
  geom_boxplot() +
  geom_jitter(width = 0.1, alpha = 0.5) +
  labs(title = "Distribución de frecuencia de erupciones por tipo de volcán",
       x = "Tipo de volcán",
       y = "Frecuencia de erupciones (años)") +
  theme_minimal() +
  theme(legend.position = "none")
# Crear histogramas separados para cada tipo de volcán
histogramas <- ggarrange(
  ggplot(data = datos_analisis %>% filter(Tipo_Volcan == "Cono de Ceniza"), 
         aes(x = Frecuencia_Erupcion_Anios)) + 
    geom_histogram(bins = 15, fill = "#F8766D", color = "white") + 
    labs(title = "Cono de Ceniza", x = "Frecuencia de erupciones", y = "Frecuencia") +
    theme_minimal(),
  
  ggplot(data = datos_analisis %>% filter(Tipo_Volcan == "Escudo"), 
         aes(x = Frecuencia_Erupcion_Anios)) + 
    geom_histogram(bins = 15, fill = "#00BA38", color = "white") + 
    labs(title = "Escudo", x = "Frecuencia de erupciones", y = "Frecuencia") +
    theme_minimal(),
  
  ggplot(data = datos_analisis %>% filter(Tipo_Volcan == "Estratovolcán"), 
         aes(x = Frecuencia_Erupcion_Anios)) + 
    geom_histogram(bins = 15, fill = "#619CFF", color = "white") + 
    labs(title = "Estratovolcán", x = "Frecuencia de erupciones", y = "Frecuencia") +
    theme_minimal(),
  
  ncol = 3,  # Organizar en 3 columnas
  common.legend = FALSE
)

# Mostrar los gráficos
print(boxplot)

print(histogramas)

# 5. Verificación de supuestos del ANOVA

# Homogeneidad de varianzas (Levene's Test)
levene_test <- leveneTest(Frecuencia_Erupcion_Anios ~ Tipo_Volcan, data = datos_analisis)
cat("\nPrueba de Levene para homogeneidad de varianzas:\n")
## 
## Prueba de Levene para homogeneidad de varianzas:
print(levene_test)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.4375 0.6469
##       97
# Normalidad por grupo (Shapiro-Wilk)
normalidad <- datos_analisis %>%
  group_by(Tipo_Volcan) %>%
  summarise(
    p_valor = shapiro.test(Frecuencia_Erupcion_Anios)$p.value
  )
cat("\nPruebas de normalidad por grupo:\n")
## 
## Pruebas de normalidad por grupo:
print(normalidad)
## # A tibble: 3 × 2
##   Tipo_Volcan    p_valor
##   <fct>            <dbl>
## 1 Cono de Ceniza  0.0486
## 2 Escudo          0.0159
## 3 Estratovolcán   0.0364
# 6. Análisis de ANOVA
modelo_anova <- aov(Frecuencia_Erupcion_Anios ~ Tipo_Volcan, data = datos_analisis)
anova_summary <- summary(modelo_anova)

cat("\nResultados del ANOVA:\n")
## 
## Resultados del ANOVA:
print(anova_summary)
##             Df  Sum Sq Mean Sq F value Pr(>F)
## Tipo_Volcan  2   35482   17741   0.799  0.453
## Residuals   97 2152949   22195
# 7. Prueba post-hoc si el ANOVA es significativo
p_valor <- anova_summary[[1]]$`Pr(>F)`[1]
if (!is.na(p_valor) && p_valor < 0.05) {
  cat("\nANOVA significativo (p =", p_valor, "). Realizando prueba post-hoc de Tukey...\n")
  
  tukey_result <- HSD.test(modelo_anova, "Tipo_Volcan", group = TRUE, console = TRUE)
  
  # Visualización resultados post-hoc
  grupos_tukey <- as.data.frame(tukey_result$groups)
  grupos_tukey$Tipo_Volcan <- rownames(grupos_tukey)
  
  tukey_plot <- ggplot(grupos_tukey, aes(x = reorder(Tipo_Volcan, -Frecuencia_Erupcion_Anios), 
                                     y = Frecuencia_Erupcion_Anios,
                                     fill = Tipo_Volcan)) +
    geom_bar(stat = "identity", alpha = 0.7) +
    geom_errorbar(aes(ymin = Frecuencia_Erupcion_Anios - sd, 
                     ymax = Frecuencia_Erupcion_Anios + sd), 
                  width = 0.2) +
    geom_text(aes(label = groups), vjust = -0.5) +
    labs(title = "Comparación de medias con prueba de Tukey",
         subtitle = paste("ANOVA p-valor =", round(p_valor, 4)),  # PARÉNTESIS CERRADO CORRECTAMENTE
         x = "Tipo de volcán",
         y = "Frecuencia media de erupciones (años)") +
    theme_minimal() +
    theme(legend.position = "none")
  
  print(tukey_plot)
  
} else if (!is.na(p_valor)) {
  cat("\nEl ANOVA no resultó significativo (p =", p_valor, "). No se requieren pruebas post-hoc.\n")
} else {
  cat("\nNo se pudo calcular el p-valor del ANOVA.\n")
}
## 
## El ANOVA no resultó significativo (p = 0.4525771 ). No se requieren pruebas post-hoc.
# 8. Análisis de tamaño del efecto
if (!is.na(p_valor) && p_valor < 0.05) {
  eta_squared <- etaSquared(modelo_anova, type = 1, anova = TRUE)
  cat("\nTamaño del efecto (Eta cuadrado):\n")
  print(eta_squared)
}