Introducción

El uso de modelos de clasificación es fundamental en el análisis predictivo y explicativo de fenómenos complejos como los tsunamis. Un modelo de clasificación se basa en un conjunto de covariables para modelar el comportamiento de una variable categórica de respuesta. El objetivo principal de estos modelos es clasificar, predecir y, sobre todo, explicar las causas subyacentes de los eventos que estudian. Mediante el análisis de datos y la aplicación de métodos matemáticos, estos modelos identifican las relaciones entre distintos factores y emplean estas relaciones para predecir el valor de una variable en función de otra. Generalmente, la predicción se limita a un número finito de resultados, como “sí” o “no”, o valores booleanos 0 y 1.

En este contexto, la regresión logística se presenta como una técnica de análisis particularmente útil. Similar a la regresión lineal, la regresión logística se utiliza cuando se desea predecir la presencia o ausencia de una característica o resultado en función de un conjunto de predictores. Sin embargo, a diferencia de la regresión lineal, en la regresión logística la variable de respuesta es categórica o dicotómica.

Este taller se centra en el análisis de datos relacionados con terremotos ocurridos en diferentes países del mundo. La base de datos df_grupo1.csv contiene 18 variables que abarcan desde características temporales y geográficas hasta impactos en infraestructura y población:

day: día
month: mes
year: año
latitude: latitud
longitude: longitud
Depth: profundidad
Magnitude: magnitud
Magnitude.Type: tipo de magnitud
intensity: intensidad
country: país
deaths: fallecidos
injuries: lesionados
damage_millions_dollars: daños en millones de dólares
houses_destroyed: casas destruidas
houses_damaged: casas dañadas
total_missing: total desaparecidos
people_affected: personas afectadas
property_damaged: propiedades dañadas

Para el propósito de este estudio, la variable dicotómica “alerta de tsunami” (flag_tsunami), que puede tomar los valores 0 o 1, se seleccionó como la variable de respuesta. Se eligieron tres de las 18 variables restantes como predictores debido a su estrecha relación con el evento del tsunami en sí mismo: magnitud (Magnitude), medida en la escala de magnitud de momento (Mw); profundidad (Depth), medida en kilómetros (km) desde la superficie terrestre hasta el hipocentro del terremoto; e intensidad (Intensity), evaluada utilizando la Escala de Intensidad de Tsunamis (ITS), que varía de 1 (muy débil) a 12 (destrucción total).

Estas tres variables fueron seleccionadas porque están directamente relacionadas con la ocurrencia del tsunami y no con las consecuencias que genera este fenómeno, permitiendo así una clasificación más precisa del evento. Finalmente, se desarrolló un modelo que permite predecir la alerta de tsunami, cumpliendo con las condiciones de bondad y ajuste. Además, se interpretó el efecto de las variables seleccionadas en la clasificación de la alerta de tsunami y se creó la ecuación matemática del modelo final.

OBJETIVOS

Desarrollar un modelo de clasificación para predecir la alerta de tsunami basado en covariables significativas.
Identificar las relaciones predictivas entre la magnitud, la profundidad y la intensidad, y su impacto en la predicción de tsunamis.
Validar el modelo de regresión logística, interpretando el efecto de las variables seleccionadas y formulando la ecuación matemática final.

Métodos y Resultados

Elección de variables significativa

En el inicio del análisis, se seleccionaron las variables más relevantes de la base de datos “datos” que podrían predecir una alerta de tsunami. Estas variables incluyen Magnitud, Profundidad e Intensidad. Se decidió descartar otras variables como Daños en millones de dólares, Muertes, Heridos, Población afectada, entre otros, la razón de esta decisión es que estas variables no son relevantes para predecir una alerta de tsunami, ya que son consecuencias del evento, no causas.

# Depurar la base de datos de variables no significativas

datos <- datos[, c("flag_tsunami", "Magnitude", "Depth", "intensity")]

str(datos)

## tibble [205 × 4] (S3: tbl_df/tbl/data.frame)
##  $ flag_tsunami: num [1:205] 0 0 1 0 0 0 1 0 0 0 ...
##  $ Magnitude   : num [1:205] 5.6 6.2 8.2 5.7 5.6 5.5 7 5.5 6.1 6.6 ...
##  $ Depth       : num [1:205] 28 12 25 15 10 10 19 10 13 212 ...
##  $ intensity   : num [1:205] 8 6 8 8 8 7 6 6 7 6 ...

Descripción de variables

En la fase inicial del estudio, se realizó un análisis descriptivo de cada variable de manera individual. Este análisis incluyó tanto las variables independientes, que se utilizaron para predecir el comportamiento, como la variable dependiente, que representa el objetivo de la predicción.Este análisis incluye la evaluación de medidas de tendencia central como la media, la mediana y la moda, así como medidas de dispersión como la varianza y la desviación estándar; además, se utilizaron gráficos como histogramas y diagramas de caja para visualizar la distribución de los datos.

Medidas de tendencia central y de dispersión de las variables

Para iniciar el análisis, se generó una tabla resumen que incluyó las principales medidas de tendencia central, como la media y la mediana, así como medidas estadísticas como el mínimo, el máximo, la desviación estándar y la desviación de error. El objetivo de este análisis fue obtener una visión general de la distribución y comportamiento de los datos iniciales, identificando los valores atípicos y centrales que podrían influir en el modelo.

La construcción de esta tabla se facilitó mediante la función get_summary_stats de la librería rstaxis, la cual generó un resumen completo de todas las medidas estadísticas de las variables numéricas. Para ahorrar tiempo en la incorporación de datos, se optó por incluir todas las variables numéricas en la función utilizando un ‘.’.

# Medidas de tendencia central y de dispersión de las variables.

resumen <- datos %>% get_summary_stats(., show = c("min","max", "median","mean","sd","se"))

resumen

## # A tibble: 4 × 8
##   variable         n   min   max median   mean     sd    se
##   <fct>        <dbl> <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>
## 1 flag_tsunami   205  0      1      0    0.317  0.466 0.033
## 2 Magnitude      205  5.5    9.1    6.6  6.64   0.75  0.052
## 3 Depth          205  0.02 664     19.7 42.5   89.5   6.25 
## 4 intensity      205  3     11      7    7.28   1.53  0.107

La tabla estadística generada proporciona una visión detallada de las diferentes medidas para cada variable. Un aspecto destacado es la alta desviación estándar en la variable Profundidad, lo que puede indicar la presencia de datos atípicos y sugerir que no sigue una distribución normal.

Por otro lado, las variables Magnitud e Intensidad muestran desviaciones estándar que podrían indicar un comportamiento de distribución normal. Esto sugiere que estos datos pueden estar bien distribuidos alrededor de la media.

Además, la desviación del error, que es una medida de la incertidumbre en las estimaciones,presenta un valor muy bajo en la mayoría de las variables, indicando una posible ausencia de datos atípicos en estas variables. Sin embargo, la variable Profundidad presenta una desviación del error alta, lo que podría sugerir la presencia de datos atípicos

Medidas de tendencia central con respecto a la variable respuesta (media y desviación)

Centrándose más en la variable respuesta, se decidió realizar un análisis específico donde se tomaron las medidas de tendencia central principales, que son la media, la desviación estándar y el coeficiente de variación. El objetivo de este análisis fue observar su comportamiento frente a la variable respuesta.

La construcción de estas tres tablas se realizó con la ayuda de la función summaryBy de la librería doBy. Esta función generó un resumen de cada variable en función de la alerta de tsunami, tanto cuando se tiene en cuenta como cuando no. Cada tabla representa una medida de tendencia.Este análisis proporcionó una visión detallada del comportamiento de la variable respuesta en relación con las variables predictoras.

# Resumen con respecto a la variable respuesta: alerta de tsunami (flag_tsunami)

coef.var=function(x){
        sd(x)/mean(x)
}
doBy::summaryBy(.~flag_tsunami, data = datos, FUN = mean)

## # A tibble: 2 × 4
##   flag_tsunami Magnitude.mean Depth.mean intensity.mean
##          <dbl>          <dbl>      <dbl>          <dbl>
## 1            0           6.36       49.6           7.25
## 2            1           7.25       27.2           7.34

doBy::summaryBy(.~flag_tsunami, data = datos, FUN = sd)

## # A tibble: 2 × 4
##   flag_tsunami Magnitude.sd Depth.sd intensity.sd
##          <dbl>        <dbl>    <dbl>        <dbl>
## 1            0        0.647    107.          1.43
## 2            1        0.584     20.2         1.73

doBy::summaryBy(.~flag_tsunami, data = datos, FUN = coef.var)

## # A tibble: 2 × 4
##   flag_tsunami Magnitude.coef.var Depth.coef.var intensity.coef.var
##          <dbl>              <dbl>          <dbl>              <dbl>
## 1            0             0.102           2.16               0.197
## 2            1             0.0805          0.743              0.236

La tabla de la media muestra que la variable Profundidad es la más afectada cuando se tiene en cuenta la variable respuesta, en comparación con Magnitud e Intensidad, donde tener en cuenta la alerta de tsunami no afectó mucho los valores.

Por otro lado, la tabla de la desviación estándar, al igual que con la tabla de la media, afectó notablemente a la variable de Profundidad, esto sugiere que hay una menor variabilidad en los datos de Profundidad cuando se tiene en cuenta la alerta de tsunami.

Finalmente, la tabla del coeficiente de variación demuestra que la presencia de la variable respuesta disminuye el valor del coeficiente tanto en la variable Magnitud como la variable Profundidad, mientras que aumenta en la variable Intensidad.

Representación gráfica

Histogramas

Para continuar con el análisis de las variables y obtener una visión más clara, se decidió realizar histogramas, los cuales son una herramienta visual efectiva que permite mostrar y visualizar el comportamiento de las variables a trabajar. Para crear los diagramas, se utilizó la función geom_histogram de la librería ggplot2.

Para visualizar los gráficos, se decidió utilizar la función grid.arrange de la librería gridExtra. Esta función es útil para organizar múltiples gráficos en una cuadrícula, lo que facilita la comparación visual de diferentes aspectos de los datos.

#Tablero gráfico de Histogramas

grid.arrange(hist_tsunami, hist_magnitude, hist_profundiad, hist_intensity, ncol = 2)

A partir de los histogramas, se observa que los datos de la variable Intensidad tienden a tener un comportamiento con una distribución normal, significando que la mayoría de los valores de Intensidad están cerca de la media, y la distribución de los valores es simétrica alrededor de la media. Por otro lado, la variable Magnitud muestra una distribución asimétrica, donde los valores tienden a estar alrededor del valor de 6.

Boxplots

Para complementar el análisis del comportamiento de las variables, se optó por usar gráficos de boxplot, los cuales son una excelente herramienta para visualizar la distribución de los datos y detectar valores atípicos.

Para desarrollar estos gráficos, se utilizó la función geom_boxplot de la librería ggplot2. Cada boxplot muestra la mediana (la línea en el medio de la caja), el primer y tercer cuartil (los bordes de la caja), y los valores atípicos (los puntos fuera de las líneas que se extienden desde la caja).

Para visualizar los gráficos, nuevamente se decidió utilizar la función grid.arrange de la librería gridExtra.

#Tablero gráfico de Histogramas

grid.arrange(box_tsunami, box_magnitud, box_profundiad, box_intensidad, ncol = 2)

De los boxplots de cada variable, se destaca que la variable Profundidad presenta muchos datos atípicos en comparación con las otras variables. Además, la variable Intensidad no presenta ningún dato atípico.

Bloxplot en funcion de la Alerta de Tsunami

Se decidió añadir un boxplot donde se compara individualmente cada variable, evaluando si la presencia de la variable de alerta de tsunami tiene algún efecto. Para ello, se optó por usar la función boxplot de la librería lapply, esta función permite crear boxplots para cada nivel de una variable categórica, lo que facilita la comparación visual de las distribuciones en diferentes niveles de la variable de alerta de tsunami.

# Gráficos individuales de todas las variables (Cuantitativas) vs Alerta de Tsunami
color = c("#66CDAA", "#FFE4B5")

par(mfrow=c(2,4))                  
attach(datos)
tabla<-prop.table(table(flag_tsunami))
coord<-barplot(tabla, col= c("#66CDAA", "#FFE4B5"),ylim=c(0,1), main="Alerta de Tsunami")
text(coord,tabla,labels=round(tabla,2), pos=3)

lapply(names(datos)[-1],function(y){
        boxplot(datos[[y]] ~ datos$flag_tsunami, ylab= y, xlab="Alerta de Tsunami",boxwex = 0.5,col=NULL)
         stripchart(datos[[y]] ~ datos$flag_tsunami, vertical = T,
                   method = "jitter", pch = 19,
                   col = color, add = T)
})

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL

detach(datos)

De los boxplots, se puede evidenciar que la presencia de la variable Alerta de Tsunami afectó la distribución de las variables predictoras. En el caso de la Magnitud, incrementó ambos bigotes del boxplot, lo que indica una mayor dispersión de los datos, haciendo que la distribución quede más centrada.

En cuanto a la variable Profundidad, la presencia de la variable respuesta ayudó a disminuir los datos atípicos que presentaba, lo que sugiere que la alerta de tsunami puede estar asociada con eventos de mayor profundidad.

Por último, para la variable Intensidad, la presencia de la variable respuesta incrementó el tamaño de la caja del boxplot, esto indica una mayor variabilidad en los datos de Intensidad cuando se presenta una alerta de tsunami.

Para comprender mejor las relaciones entre las variables independientes y la variable respuesta, se empleó el diagrama de dispersión. Esto permitió visualizar la correlación existente entre las variables y determinar si alguna variable estaba relacionada de manera significativa con la alerta de tsunami. Se utilizó la función pairs para crear estos diagramas.

# Gráficos de dispersión 

pairs(datos[, -1],panel = panel.smooth)

Como se pudo evidenciar con los gráficos de dispersión, parece ser que no existe ninguna tendencia, lo que sugiere que ninguna variable se encuentra relacionada con otra.

Para confirmar el tipo de relación y el nivel, se calculó una matriz de correlación utilizando la función cor y se visualizó con la función corrplot. Esto proporcionó una visión más detallada de las relaciones y ayudó a identificar las variables que mostraban una correlación fuerte con la alerta de tsunami.

# Construcción de una matriz de correlaciones

m.cor <- cor(datos[,-1])

# Cambio de color de la matriz

col <- colorRampPalette(c("gray", "pink", "purple"))(50) 

# Gráfico para visualizar la matriz de correlaciones

corrplot(m.cor, method="number", type="upper", number.cex = 0.7, col = col, tl.col = "black")

Como se esperaba, ninguna variable muestra una fuerte correlación entre ellas. El valor máximo de correlación observado es de 0.27, que se encuentra entre las variables Profundidad y Magnitud. Siendo beneficioso ya que las variables predictoras altamente correlacionadas pueden causar problemas de multicolinealidad.

Modelos

Contrucción de Modelo Saturado

Para realizar el ajuste del modelo de regresión logística, se inició con la construcción de un modelo saturado. Este proceso consistió en incluir las variables predictivas previamente seleccionadas, para ello, se hizo uso de la función glm de la librería base de R, la cual es utilizada para el ajuste de modelos de regresión logística.

Mediante la función summary, se obtuvo un resumen de la información relacionada con el modelo ajustado, principalmente, y más importante, se obtuvieron los p-valores que representan si la variable es significativa o no. Estos valores permiten la descripción de la relación entre las variables predictivas y la variable de respuesta.

# Ajuste del modelo Saturado

mod_sat <- glm(formula = flag_tsunami ~., data = datos, family = "binomial")  
summary(mod_sat)

## 
## Call:
## glm(formula = flag_tsunami ~ ., family = "binomial", data = datos)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -17.787587   2.672489  -6.656 2.82e-11 ***
## Magnitude     2.908761   0.417428   6.968 3.21e-12 ***
## Depth        -0.023414   0.008118  -2.884  0.00392 ** 
## intensity    -0.265990   0.128981  -2.062  0.03918 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 256.10  on 204  degrees of freedom
## Residual deviance: 157.93  on 201  degrees of freedom
## AIC: 165.93
## 
## Number of Fisher Scoring iterations: 6

El modelo saturado muestra el p-valor, que está asociado con la prueba de hipótesis de que el coeficiente es igual a 0. En este caso, el valor p representa el nivel de significancia de las variables predictivas y es del 10%. Por lo tanto, las variables con p-valores menores a 0,10 son estadísticamente significativas, mientras que las que tengan valores mayores no lo son. Como se muestra en la tabla, todas las variables son altamente significativas.

Selección de varibles mediante la estategia Backward

Para verificar si es necesario quitar alguna variable, se procedió a realizar una selección de variables mediante el método de pasos, conocido también como Stepwise selection, utilizando la función step. Esta función permitió la selección de un subconjunto óptimo de variables predictivas para el modelo en función del criterio de Akaike (AIC).

El criterio AIC tiene como objetivo hallar un equilibrio entre la capacidad de ajuste del modelo y su complejidad. Para elegir el mejor modelo, se busca aquel con un AIC más bajo, lo que indica que permite un buen ajuste de los datos observados sin ser excesivamente complejo.

# Para ello se utiliza la función Step

mod_red <- step(mod_sat,direction = "backward")

## Start:  AIC=165.93
## flag_tsunami ~ Magnitude + Depth + intensity
## 
##             Df Deviance    AIC
## <none>           157.93 165.93
## - intensity  1   162.41 168.41
## - Depth      1   187.06 193.06
## - Magnitude  1   252.15 258.15

summary(mod_sat)

## 
## Call:
## glm(formula = flag_tsunami ~ ., family = "binomial", data = datos)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -17.787587   2.672489  -6.656 2.82e-11 ***
## Magnitude     2.908761   0.417428   6.968 3.21e-12 ***
## Depth        -0.023414   0.008118  -2.884  0.00392 ** 
## intensity    -0.265990   0.128981  -2.062  0.03918 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 256.10  on 204  degrees of freedom
## Residual deviance: 157.93  on 201  degrees of freedom
## AIC: 165.93
## 
## Number of Fisher Scoring iterations: 6

Como se muestra, la funcion step no quitó niguna variable al modelo, lo que sugiere que el valor del AIC no es posible reducirlo más, lo que sugiere que el modelo saturado ya cuenta con las variables significativas para la predicción de la alerta de tsunami.

Por lo tanto, se verificó la multicolinealidad de las variables mediante el índice VIF (Variance Inflation Factor). El VIF es una medida que permite cuantificar la severidad de la multicolinealidad en un análisis de regresión ordinaria de mínimos cuadrados.

# Verificar multicolinealidad

round(vif(mod_sat),2)

## Magnitude     Depth intensity 
##      1.27      1.22      1.11

Con esta prueba se concluye que efectivamente todas las variables son independientes unas de otras y, por lo tanto, no es necesario omitir ninguna.

Por ende, se nombra al modelo saturado como “modelo”, lo que indica que es el modelo de ajuste elegido. Este modelo será el que se utilice para las predicciones futuras y el análisis posterior.

# Modelo final

modelo <- mod_sat

Construcción de Modelo Nulo

Se creó un modelo nulo, el cual se ajusta el modelo en función del intercepto y no posee variables predictoras. El propósito de este modelo es servir como una línea de base para comparar con el modelo saturado y determinar si este último es efectivamente mejor.

mod_null <- glm(formula= flag_tsunami~1, data=datos,family = "binomial")

Comparación de modelos

Mediante la prueba ANOVA (Análisis de Varianza), se realiza una comparación entre el modelo final (o saturado) y el modelo nulo.

# Modelo Final vs Modelo Nulo

anova(mod_null,mod_sat,test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: flag_tsunami ~ 1
## Model 2: flag_tsunami ~ Magnitude + Depth + intensity
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       204     256.10                          
## 2       201     157.93  3   98.171 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ho = ambos modelos presentan similar ajuste H1 = el modelo con más términos presenta mejor ajuste

Como resultado del análisis, se encontró que el p-valor de la comparación entre el modelo nulo y el modelo saturado tiende a 0. Esto indica que hay una diferencia significativa entre estos dos modelos. Por lo tanto, se acepta la hipótesis alternativa, lo que sugiere que el modelo saturado proporciona un mejor ajuste a los datos en comparación con el modelo nulo.

Analisis del Modelo Final

Después de obtener el modelo final, se comienza a realizar un análisis de este y verificar su efectividad, con el objetivo de obtener la fórmula que se utilizará para predecir una alerta de tsunami en base a las variables definidas.

Interpretación de los parametros

Primero, se analizaron las variables predictivas y su comportamiento en el modelo, este análisis es crucial para entender cómo cada variable contribuye a la predicción de una alerta de tsunami. Para visualizar el comportamiento de estas variables, se utilizó la función cbind de la librería base de R.

# Creación de función para la interpretación
Int.par <- function(Modelo){
  Tabla_coef <- cbind(Coef=round(coef(Modelo),2),IC=round(confint(Modelo,level=0.95),2))
  Tabla_odds <- round(exp(Tabla_coef),2)
  colnames(Tabla_odds)[1]="e-beta"
  list(Coeficientes= Tabla_coef,OR=Tabla_odds)
}

Int.par(modelo)

## Waiting for profiling to be done...

## $Coeficientes
##               Coef  2.5 % 97.5 %
## (Intercept) -17.79 -23.48 -12.94
## Magnitude     2.91   2.15   3.80
## Depth        -0.02  -0.04  -0.01
## intensity    -0.27  -0.53  -0.02
## 
## $OR
##             e-beta 2.5 % 97.5 %
## (Intercept)   0.00  0.00   0.00
## Magnitude    18.36  8.58  44.70
## Depth         0.98  0.96   0.99
## intensity     0.76  0.59   0.98

Como se observa en la columna e-beta, las variables Magnitud, Profundidad e Intensidad tienen un impacto significativo en la probabilidad de una alerta de tsunami:

En el caso de la Magnitud, por cada incremento de 18 unidades, la probabilidad de que haya una alerta de tsunami aumenta un 36%. Esto sugiere que los eventos de mayor magnitud tienen una mayor probabilidad de generar una alerta de tsunami.
Con la variable Profundidad, por cada incremento de una unidad, la probabilidad de que ocurra una alerta de tsunami disminuye un 0,02%. Esto indica que los eventos de menor profundidad tienen una mayor probabilidad de generar una alerta de tsunami.
Finalmente, en la variable Intensidad, por cada incremento de una unidad, la probabilidad de que haya una alerta de tsunami disminuye un 31%. Esto sugiere que los eventos de menor intensidad tienen una mayor probabilidad de generar una alerta de tsunami.

Análisis de la Curva ROC (PC = 0.5)

Después de obtener el modelo final de regresión logística, es fundamental evaluar su efectividad y comprender cómo cada variable contribuye a la predicción de una alerta de tsunami. Uno de los métodos más utilizados para evaluar la capacidad predictiva de un modelo es la Curva ROC (Receiver Operating Characteristic).

La Curva ROC representa la relación entre la tasa de verdaderos positivos (TPR) y la tasa de falsos positivos (FPR) a medida que se varía el umbral de clasificación. En otras palabras, muestra cómo el modelo se desempeña en términos de sensibilidad y especificidad en diferentes puntos de corte.

Para crear la Curva ROC en R, se utilizó la función roc de la librería pROC, en este caso se opto por tener un punto de corte en 0.5, con el objetivo de ajustarlo despues en el punto de corte ideal. Esta función calculó la TPR y la FPR para diferentes umbrales y trazó la curva, además, el área bajo la curva (AUC) proporcionó una medida global de la capacidad predictiva del modelo.

# Creación de la curva ROC
prob_pred <- predict(modelo, type = "response")
roc_curva <- roc(datos$flag_tsunami, prob_pred, plot=TRUE, main="Curva ROC del Modelo")

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

Evalución del Poder de Clasificación

Después de obtener el modelo final de regresión logística, se procedió a evaluar su capacidad de clasificación. El objetivo era comprender cómo cada variable contribuía a la predicción de una alerta de tsunami y verificar la efectividad general del modelo.

Para evaluar la capacidad de clasificación de nuestro modelo de regresión logística, se llevaron a cabo tres operaciones clave:

Obtención de las probabilidades pronosticadas: Utilizamos la función predict con el argumento type = “response” para obtener las probabilidades pronosticadas por nuestro modelo.
Predicciones del modelo: Establecimos un punto de corte (pc) en 0.5 para clasificar las probabilidades en 0 o 1. Si la probabilidad pronosticada es mayor que el punto de corte, asignamos 1 (alerta de tsunami); de lo contrario, asignamos 0 (no alerta de tsunami).
Indicadores de correcta clasificación: Utilizamos la función confusionMatrix() de la librería caret para evaluar la clasificación del modelo

# Obtención de las probabilidades pronosticadas

prob_pred <- predict(modelo, type = "response")

# Predicciones del modelo

pc = 0.5
predicciones <- as.factor(ifelse(prob_pred > pc, 1, 0))

# Indicadores de correcta clasificación

caret::confusionMatrix(predicciones, as.factor(datos$flag_tsunami), dnn = c("Alertas de Tsunami observadas", "Alertas de Tsunamis pronosticadas" ), positive = "0")

## Confusion Matrix and Statistics
## 
##                              Alertas de Tsunamis pronosticadas
## Alertas de Tsunami observadas   0   1
##                             0 126  22
##                             1  14  43
##                                           
##                Accuracy : 0.8244          
##                  95% CI : (0.7653, 0.8739)
##     No Information Rate : 0.6829          
##     P-Value [Acc > NIR] : 3.465e-06       
##                                           
##                   Kappa : 0.5807          
##                                           
##  Mcnemar's Test P-Value : 0.2433          
##                                           
##             Sensitivity : 0.9000          
##             Specificity : 0.6615          
##          Pos Pred Value : 0.8514          
##          Neg Pred Value : 0.7544          
##              Prevalence : 0.6829          
##          Detection Rate : 0.6146          
##    Detection Prevalence : 0.7220          
##       Balanced Accuracy : 0.7808          
##                                           
##        'Positive' Class : 0               
##

Según los datos observados se puede ver que el modelo tiene una buena capacidad para predecir correctamente las alertas de tsunami, al poseer una precisión del 85%, entre proporción de predicciones positivas correctas con respecto a todas las predicciones positivas.

También, los resultados indican que la concordancia entre las predicciones del modelo y las obsevaciones reales no es muy alta, esto lo indico el indice kappa, que tiene un valor de 0.5.

Además, al modelo tener una sensibilidad de un 0.9 indica que el modelo tiene alta capacidad de detectar alertas de tsunami reales, por el contrario el modelo presenta una especificidad de un 0.6, indicando que el propio modelo tiene una capaciad moderada para clasificar correctamente los casos de no alerta.

Optimización del Punto de Corte

Para mejorar la especificidad del modelo y, por ende, el índice kappa, se procedió a optimizar el punto de corte, para lograrlo, se realizó una visualización de las probabilidades predichas utilizando diferentes puntos de corte. El objetivo principal de esta estrategia es encontrar el umbral óptimo que maximice la precisión y minimice los errores de clasificación.

Visualización de las probabilidades predichas

Para visualizar las probabilidades predichas y obtener el punto de corte óptimo, se opto por utilizar la siguiente funcion:

with(datos,{
        dD=density(prob_pred[flag_tsunami == 1])
        dS=density(prob_pred[flag_tsunami == 0])
        plot(NULL, xlim=range(prob_pred), ylim=range(c(dD$y,dS$y)), 
             flag_tsunami = "n", ylab = "Densidad",xlab = "Probabilidad Predicha")
        lines(dD,col = "#104E8B"); lines(dS, col = "#8B3A62")
        text( 0.08,3, "Sin Alerta de Tsunami"); text(0.82,1, "Alerta de Tsunami")
        PC.option=c(0.5,0.25,0.75)
        abline(v=PC.option,col="#828282",lty=2)
        text(PC.option,3,paste("PC = ",PC.option),cex=0.7,adj=0)})

Como se observa, la gráfica muestra que el punto de corte se encuentra cerca de 0.25, específicamente alrededor de 0.27 o 0.28.

Curva ROC especificando el punto de corte óptimo

Despues de conocer el punto de corte optimo, se procedio a realizar nuevamente la curva Una vez identificado el punto de corte óptimo, se procedió a generar nuevamente la curva ROC. El objetivo principal era visualizar la eficiencia de esta curva al comparar la sensibilidad y la especificidad del modelo, para facilitar el gráfico, se utilizó la función plot.roc de la librería pROC.

pROC::plot.roc(roc_curva, legacy.axes = TRUE, print.thres = "best", print.auc = TRUE,
               auc.polygon = FALSE, max.auc.polygon = FALSE, auc.polygon.col = "gainsboro",
               col = 2, grid = TRUE)

Evalución del Poder de Clasificación con PC óptimo

Para evaluar el desempeño estadístico del modelo, se realizó una nueva evaluación, con la diferencia de que en esta ocasión fue la selección del punto de corte óptimo.Realizar la evaluación, tiene como objetivo de evaluar la efectividad general del modelo.

# Predicciones del modelo

pc = 0.270
predicciones <- as.factor(ifelse(prob_pred > pc, 1, 0))

# Indicadores de correcta clasificación

caret::confusionMatrix(predicciones, as.factor(datos$flag_tsunami), dnn = c("Alertas de Tsunami observadas", "Alertas de Tsunamis pronosticadas" ), positive = "0")

## Confusion Matrix and Statistics
## 
##                              Alertas de Tsunamis pronosticadas
## Alertas de Tsunami observadas   0   1
##                             0 112   6
##                             1  28  59
##                                          
##                Accuracy : 0.8341         
##                  95% CI : (0.776, 0.8823)
##     No Information Rate : 0.6829         
##     P-Value [Acc > NIR] : 6.624e-07      
##                                          
##                   Kappa : 0.6489         
##                                          
##  Mcnemar's Test P-Value : 0.0003164      
##                                          
##             Sensitivity : 0.8000         
##             Specificity : 0.9077         
##          Pos Pred Value : 0.9492         
##          Neg Pred Value : 0.6782         
##              Prevalence : 0.6829         
##          Detection Rate : 0.5463         
##    Detection Prevalence : 0.5756         
##       Balanced Accuracy : 0.8538         
##                                          
##        'Positive' Class : 0              
##

Según los datos observados y en comparacion, se puede ver que el modelo tiene una buena capacidad para predecir correctamente las alertas de tsunami, al poseer una precisión del 83%, aunque disminuyo un 2% que con el punto de corte en 0.5, aun sigue teniendo una buena capacidad para predecir.

También, los resultados indican que la concordancia entre las predicciones del modelo y las obsevaciones reales aumento con este nuevo punto de corte, llegando a ser de 0.6 en el indice kappa.

Además, al modelo nuevo, la sensibilidad disminuyo a 0.8, aun indica que el modelo tiene alta capacidad de detectar alertas de tsunami reales, por el contrario al modelo con el nuevo punto de corte aumento la especificidad de un 0.9, indicando que el propio modelo aumento la capacidad para clasificar correctamente los casos de no alerta.

Ecuación del Modelo Final

summary(modelo)

## 
## Call:
## glm(formula = flag_tsunami ~ ., family = "binomial", data = datos)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -17.787587   2.672489  -6.656 2.82e-11 ***
## Magnitude     2.908761   0.417428   6.968 3.21e-12 ***
## Depth        -0.023414   0.008118  -2.884  0.00392 ** 
## intensity    -0.265990   0.128981  -2.062  0.03918 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 256.10  on 204  degrees of freedom
## Residual deviance: 157.93  on 201  degrees of freedom
## AIC: 165.93
## 
## Number of Fisher Scoring iterations: 6

Después de completar todo el proceso de construcción y ajuste del modelo de regresión logística, se llegó a la ecuación final que permite predecir la probabilidad de una alerta de tsunami.

\[P(Y = 1 | \vec{X} = \vec{x}) = \frac{1}{1 + e^{-(-17.787587+2.908761\times \text{Magnitude}-0.023414 \times \text{Depth} -0.265990\times\text{intensity})}}\]

Conclusiones

Las variables seleccionadas para el modelo saturado demostraron ser estadísticamente significativas. Inicialmente, la clasificación del modelo se realizó con un punto de corte de 0.5. La elección de este punto de corte se basa en la generación de una mayor cantidad de registros de Verdaderos Positivos y Verdaderos Negativos en la matriz de confusión 2x2. La sensibilidad alcanzó un valor del 90%, lo que indica una alta capacidad para clasificar los Verdaderos Positivos (posición 1,1), mientras que la especificidad fue del 65%, reflejando la capacidad para clasificar los Verdaderos Negativos (posición 2,2).
Se identificó un hiperparámetro que permite mejorar tanto la sensibilidad como la especificidad, según el gráfico de Distribución de Densidad de Probabilidades para la variable respuesta flag_tsunami. Este gráfico mostró que un punto de corte cercano a 0.25 es el más óptimo para lograr un equilibrio adecuado entre ambas métricas. Nuevamente, con un punto de corte del 0.27 según la grafica ROC se obtuvo que la sensibilidad alcanzó un valor de 0.8, en contraste con la especificidad de 0.9.

Taller 3: Análisis de Regresión Logística

Arias Eva, Durán Mariana, Restrepo Nicolás

2024-05-10