ANÁLISIS ESTADÍSTICO DE VOLCANES ACTIVOS A NIVEL GLOBAL
REGRESIÓN MULTIPLE LINEAL

CARRERA DE GEOLOGÍA

GRUPO N°2

ANÁLISIS ESTADÍSTICO MULTIVARIABLE

1. CARGA DE LIBRERIAS Y DATOS

2. TABLA DE PARES DE VALORES

3. DIAGRAMA DE DISPERSIÓN

4. CONJETURA DEL MODELO

5. TEST DE APROBACIÓN Y RESTRICCIONES

6. CALCULO DE PRONOSTICOS

7. CONCLUSIONES

CARGA DE LIBRERIAS Y DATOS

## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
## 
## 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
## Rows: 898 Columns: 77
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (26): volcanoLocationNum, volcano_name, location_description, country, v...
## dbl (43): noaa_event_id, year, month, day, linked_tsunami_event_id, linked_e...
## num  (1): ejecta_volume_km3_min
## lgl  (7): is_significant, publish, is_eruption, tsunami_confirmed, ring_of_f...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

JUSTIFICACIÓN

DEFINICIÓN DE VARIABLES X1, X2, Y 

# Extraemos las variables correctas de tu dataset
alt_pluma   <- as.numeric(Volcanes_Globales$est_plume_height_km)         # Y
vei         <- as.numeric(Volcanes_Globales$vei)                         # X1
conteo_erup <- as.numeric(Volcanes_Globales$volcano_eruption_count_in_dataset) # X2

# Armamos el dataset inicial eliminando vacíos y asegurando valores mayores a cero
TPV_inicial <- data.frame(y = alt_pluma, x1 = vei, x2 = conteo_erup)
TPV_inicial <- na.omit(TPV_inicial)
TPV_inicial <- TPV_inicial[TPV_inicial$y > 0 & TPV_inicial$x1 > 0 & TPV_inicial$x2 > 0, ]
total_inicial <- nrow(TPV_inicial)

USO INTERCUARTIL 

filtro_iqr <- function(v){
  Q1 <- quantile(v, 0.25, na.rm = TRUE)
  Q3 <- quantile(v, 0.75, na.rm = TRUE)
  IQRv <- Q3 - Q1
  li <- Q1 - 1.5 * IQRv   
  ls <- Q3 + 1.5 * IQRv
  return(v >= li & v <= ls)
}

# Aplicamos el filtro simultáneo
filas_buenas <- filtro_iqr(TPV_inicial$y) & filtro_iqr(TPV_inicial$x1) & filtro_iqr(TPV_inicial$x2)
TPV_sin_outliers <- TPV_inicial[filas_buenas, ]

USO DE MEDIA ARITMETICA 

# Para regresión múltiple se conservan todas las observaciones
# No se agrupan las variables porque cada volcán representa una tripleta independiente

TPV <- TPV_sin_outliers

row.names(TPV) <- NULL

TPV_final_id <- cbind(
  Nro = 1:nrow(TPV),
  TPV
)

# No existen reducciones por media aritmética
valores_consolidados_media <- 0

OUTLIERS ENCONTRADOS Y REDUCCIÓN DE FILAS REPETIDAS 

# Outliers encontrados por Intercuartil
outliers_encontrados <- total_inicial - nrow(TPV_sin_outliers)
outliers_encontrados
## [1] 40
# Reducción de filas repetidas por Media Aritmética
valores_consolidados_media <- nrow(TPV_sin_outliers) - nrow(TPV)
valores_consolidados_media
## [1] 0
# Total de tripletas de valores finales
total_tripletas <- nrow(TPV)
total_tripletas
## [1] 659

TABLA DE PARES DE VALORES 

##    N° Índice VEI (X1) Conteo de Erupciones (X2) Altura de Pluma [km] (Y)
## 1   1              20                         5                       18
## 2   2              20                         5                       11
## 3   3              10                         4                        6
## 4   4               3                         3                        8
## 5   5               3                         3                        8
## 6   6              20                         5                       18
## 7   7               1                         2                        1
## 8   8              10                         4                       15
## 9   9               3                         3                       15
## 10 10               3                         3                       18
TABLA N°1
Primeras 10 tripletas de valores para la regresión múltiple lineal
N° VEI (X1) Erupciones (X2) Altura de Pluma (km) (Y)
1 20 5 18.00
2 20 5 11.00
3 10 4 6.00
4 3 3 8.00
5 3 3 8.00
6 20 5 18.00
7 1 2 1.00
8 10 4 15.00
9 3 3 15.00
10 3 3 18.00

DIAGRAMA DE DISPERSIÓN

DIAGRAMA ORIGINAL

DIAGRAMA POST FILTRADO 

#DIAGRAMA POST FILTRACIOn
x1 <- TPV$x1        # Índice VEI
x2 <- TPV$x2        # Conteo de Erupciones
y  <- TPV$y         # Altura de Pluma

# Ajustamos los márgenes para que las etiquetas de los ejes se visualicen por completo
par(mar = c(5, 6, 4, 2))

# Generamos el diagrama de dispersión 3D post-filtrado
Cobrereg <- scatterplot3d(x1, x2, y, angle = 225, pch = 16, color = "blue",
                          cex.symbols = 0.6,
                          main = "Gráfica N°2: Diagrama de dispersión de Altura de Pluma,\nÍndice VEI y Conteo de Erupciones (Post-Filtrado)",
                          xlab = "Índice VEI (x1)",
                          ylab = "Conteo de Erupciones (x2)",
                          zlab = "Altura de Pluma [km] (y)",
                          y.margin.add = 0.9,
                          las = 1)

CONJETURA DEL MODELO

DEBIDO A LA SIMILITUD VISUAL EN LA NUBE DE PUNTOS CONJETURAMOS UN MODELO MULTIVARIABLE LINEAL

CALCULO DE PARAMETROS 

# Ajustamos el modelo de regresión lineal múltiple con las nuevas variables
regresion_multiple <- lm(y ~ x1 + x2)

# Intercepto (a) - Altura de pluma estimada cuando VEI y el Conteo de Erupciones son cero
a_val <- round(coef(regresion_multiple)[1], 3)
a_val
## (Intercept) 
##      -7.563
# Pendiente b (x1) - Efecto del Índice VEI manteniendo constante el conteo de erupciones
b_val <- round(coef(regresion_multiple)[2], 5) 
b_val
##      x1 
## 4.30558
# Pendiente c (x2) - Efecto del Conteo de Erupciones manteniendo constante el VEI
c_val <- round(coef(regresion_multiple)[3], 5)
c_val
##       x2 
## -0.02025

MODELO Y REALIDAD

TESTS DE APROBACIÓN

(R²) 

#TEST DE APROBACION
# Extraemos el Coeficiente de Determinación (R-cuadrado) del nuevo modelo
r2 <- summary(regresion_multiple)$r.squared
r2*100
## [1] 80.07675

COHEFICIENTE DE CORRELACIÖN MULTIPLE 

# Cálculo del Coeficiente de Correlación Múltiple
r_multiple <- sqrt(r2)
r_multiple
## [1] 0.8948562

TABLA DE RESUMEN 

##                               Estadistico Valor_Decimal Valor_Porcentaje
## 1       Coeficiente de Determinación (R²)     0.8007675           80.08%
## 2 Coeficiente de Correlación Múltiple (R)     0.8948562           89.49%
TABLA N°2
Test de aprobación del modelo de regresión múltiple lineal
Indicador estadístico Valor decimal Resultado (%)
Coeficiente de Determinación (R²) 0.8008 80.08%
Coeficiente de Correlación Múltiple (R) 0.8949 89.49%

RESTRICCIONES DEL MODELO 

# Parámetros del modelo

intercepto <- -7.563
b1 <- 4.30558
b2 <- -0.02025


# ------------------------------------------------------------------------------
# Y = 0
# ------------------------------------------------------------------------------

ecuacion_raiz <- paste0(
  "x2 = ",
  round(b1/abs(b2), 3),
  "x1 - ",
  round(abs(intercepto)/abs(b2), 3)
)


# ------------------------------------------------------------------------------
# Y < 0
# ------------------------------------------------------------------------------

intervalo_negativo <- paste0(
  "x2 > ",
  round(b1/abs(b2), 3),
  "x1 - ",
  round(abs(intercepto)/abs(b2), 3)
)


# ------------------------------------------------------------------------------
# Dominios
# ------------------------------------------------------------------------------

dominio_x1 <- "(0,+∞)"
dominio_x2 <- "(0,+∞)"
dominio_y  <- "(0,+∞)"


# ==============================================================================
# TABLA
# ==============================================================================

tabla_restricciones_multiple <- data.frame(
  
  Punto = c(
    "Y = 0",
    "Y < 0",
    "Dom. x1",
    "Dom. x2",
    "Dom. Y"
  ),
  
  Calculo_matematico = c(
    "-7.563 + 4.30558x1 - 0.02025x2 = 0",
    "-7.563 + 4.30558x1 - 0.02025x2 < 0",
    "x1 ∈ (0,+∞)",
    "x2 ∈ (0,+∞)",
    "Y ∈ (0,+∞)"
  ),
  
  Resultado = c(
    ecuacion_raiz,
    intervalo_negativo,
    dominio_x1,
    dominio_x2,
    dominio_y
  ),
  
  Conclusion = c(
    
    "El modelo presenta infinitas soluciones que forman una recta donde la altura de la pluma es igual a cero.",
    
    "Los pares (x1,x2) que satisfacen esta desigualdad generan alturas negativas, por lo que el modelo deja de tener significado físico.",
    
    "El Índice de Explosividad Volcánica solo puede tomar valores reales positivos.",
    
    "El número de erupciones solo puede tomar valores reales positivos.",
    
    "La altura de la pluma solo puede tomar valores reales positivos."
    
  )
  
)

RESTRICCIONES DEL MODELO

Restricciones matemáticas del modelo de regresión lineal múltiple
Modelo: Y = -7.563 + 4.30558·x1 - 0.02025·x2
Punto evaluado Procedimiento matemático Resultado Conclusión
Y = 0 -7.563 + 4.30558x1 - 0.02025x2 = 0 x2 = 212.621x1 - 373.481 El modelo presenta infinitas soluciones que forman una recta donde la altura de la pluma es igual a cero.
Y < 0 -7.563 + 4.30558x1 - 0.02025x2 < 0 x2 > 212.621x1 - 373.481 Los pares (x1,x2) que satisfacen esta desigualdad generan alturas negativas, por lo que el modelo deja de tener significado físico.
Dom. x1 x1 ∈ (0,+∞) (0,+∞) El Índice de Explosividad Volcánica solo puede tomar valores reales positivos.
Dom. x2 x2 ∈ (0,+∞) (0,+∞) El número de erupciones solo puede tomar valores reales positivos.
Dom. Y Y ∈ (0,+∞) (0,+∞) La altura de la pluma solo puede tomar valores reales positivos.

CALCULO DE PRONOSTICOS

CALCULOS 

 # 1. Definimos los valores de prueba para el pronóstico volcánico
  x1_test <- 4     # Valor de Índice VEI (X1)
x2_test <- 15    # Valor de Conteo de Erupciones (X2)

# 2. Evaluamos en la ecuación del modelo para obtener la respuesta (Y)
y_pred <- a_val + (b_val * x1_test) + (c_val * x2_test)
y_pred_round <- round(y_pred, 1)

CONCLUSIÓN