Este tutorial describe el proceso de Machine Learning para la clasificación de uso de suelo y estimación de índice de quemado a partir de datos de sensores remotos en diferentes ubicaciones geográficas en Ecuador.
Etapa 1: Dominio del problema
Los datos proporcionados representan mediciones obtenidas mediante sensores remotos en diferentes ubicaciones geográficas en Ecuador. Estos datos son útiles para una variedad de aplicaciones, desde la monitorización ambiental hasta la gestión de recursos naturales y la planificación urbana.
El objetivo principal de utilizar estos datos es estimar el índice de áreas quemadas en función de variables disponibles. El objetivo secundario es clasificar la clase de uso del suelo (bosque, cultivo, urbana, pastizal).
El dataset proporcionado incluye la siguiente descripción de variables:
longitud: Coordenada geográfica de longitud donde se encuentra la parcela.
latitud: Coordenada geográfica de latitud donde se encuentra la parcela.
altitud: Altitud medida en metros sobre el nivel del mar en la ubicación de la parcela. Puede contener algunos datos faltantes.
temp.media: Temperatura media registrada en grados Celsius en la ubicación de la parcela. Algunos valores pueden contener errores.
temp.max: Temperatura máxima registrada en grados Celsius en la ubicación de la parcela.
temp.min: Temperatura mínima registrada en grados Celsius en la ubicación de la parcela.
precipitacion: Cantidad de precipitación medida en mm en la ubicación de la parcela. Puede contener algunos datos faltantes.
ndvi: Índice de Vegetación de Diferencia Normalizada (NDVI), que indica la densidad de vegetación en la ubicación de la parcela.
-1, indica una superficie completamente no vegetada o cubierta por agua.
1, indica una superficie con vegetación muy densa.
ndwi: Índice de Agua de Diferencia Normalizada (NDWI), que indica la presencia de agua en la ubicación de la parcela.
-1, indica la presencia de agua pura.
1, indica la ausencia total de agua.
ndbi: Índice de Quema Normalizada (NDBI), que indica la presencia de áreas quemadas en la ubicación de la parcela.
reflectancia.rojo: Reflectancia en la banda roja en la ubicación de la parcela.
reflectancia.infrarrojo: Reflectancia en la banda infrarroja en la ubicación de la parcela.
clase.uso.suelo: Clase de uso del suelo asignada a la ubicación de la parcela.
indice.quemado: Índice de áreas quemadas en la ubicación de la parcela, esto es, propensión de un área a sufrir incendios o la presencia de áreas quemadas.
A partir de la información proporcionada, se responden las siguientes preguntas:
¿Qué exactamente deseamos hacer?
Clasificar el uso de suelo y estimar índice de quemado a partir de datos de sensores remotos en diferentes ubicaciones geográficas en Ecuador.
¿Es factible alcanzar lo que buscamos con los datos disponibles?
Si, las variables proporcionadas aparentemente parecen ser suficientes para alcanzar la clasificación y la regresión.
¿Cómo podemos lograrlo?
Aplicar técnicas predictivas de Machine Learning.
¿Qué tipo de problema se va a resolver?
Problema de clasificación (uso de suelo). Problema de regresión (estimación de índice de quemado).
¿El objetivo es?
Predecir. A nivel numérico (índice de quemado) y a nivel categórico (uso de suelo).
Precisando:
Variable a Predecir:
Uso del suelo (clase.uso.suelo): clasificación de la cobertura del suelo en función de las características proporcionadas por los sensores remotos. Clases disponibles:
1: Bosque o Vegetación Natural: Esta clase representa áreas cubiertas principalmente por vegetación natural, como bosques, selvas tropicales, manglares u otras formaciones vegetales no modificadas significativamente por la actividad humana.
2: Agricultura o Cultivos: Esta clase incluye áreas utilizadas para la agricultura, como campos de cultivo, plantaciones, pastizales agrícolas u otros terrenos destinados a la producción agrícola.
3: Áreas Urbanas o de Desarrollo: Esta clase abarca áreas urbanizadas, incluyendo ciudades, pueblos, zonas industriales, infraestructuras urbanas y áreas desarrolladas destinadas a usos residenciales, comerciales o industriales.
4: Pastizales o Tierras de Pastoreo: Esta clase incluye áreas de pastizales naturales o áreas de pastoreo gestionadas, donde se crían animales para la producción de carne, leche u otros productos.
Índice de áreas quemadas (indice.quemado): regresión mediante estimación del nivel de áreas quemadas en una determinada ubicación.
Valores bajos de índice.quemado indican una menor probabilidad de incendios o la ausencia de áreas quemadas en la región. Esto podría significar que el área ha experimentado pocos incendios en el pasado reciente o que ha habido una efectiva gestión de incendios para prevenir la propagación del fuego.
Valores altos de índice.quemado sugieren una mayor probabilidad de incendios o la presencia de áreas quemadas en la región. Esto podría indicar que el área ha sido afectada por incendios recientes o que hay condiciones propicias para que ocurran incendios, como sequías, altas temperaturas, o una acumulación de material combustible.
En esta etapa, se procederá a la adquisición de los datos necesarios para el análisis. Los datos serán leídos desde un archivo CSV y se inspeccionarán para asegurar que se hayan cargado correctamente.
Para comenzar, utilizamos la función read_delim de la biblioteca readr para leer los datos desde un archivo CSV. El archivo se encuentra en la carpeta “Datos” y se llama datos_sensores_remotos.csv.
Code
# Leer los datos desde un archivo CSVdatos_originales <- readr::read_delim("Datos/datos_sensores_remotos.csv", delim =",")# Mostrar una vista preliminar de los datosdatos_originales %>%glimpse()
La función glimpse nos ofrece una descripción rápida del conjunto de datos, y permite contestar la pregunta:
¿Cuántas y cuáles son las instancias y variables del conjunto de datos?
Son 2900 instancias y 14 variables. Se observa además que clase.uso.suelo es una variable cualitativa y todas las demás variables son de tipo cuantitativo.
Etapa 3. Preprocesamiento de datos
En esta etapa, realizaremos el preprocesamiento de los datos para asegurarnos de que estén listos para el análisis posterior. Esto incluye: a. limpieza, b. transformación de los datos y c. análisis exploratorio de datos (EDA)
Para comenzar, contestamos las siguientes preguntas:
¿Existen valores ausentes o datos faltantes en las variables y cómo manejarlos?
Para aquello usamos la función glimpse para inspeccionamos los datos
Code
# Inspeccionar las variables del conjunto de datosdatos_originales %>%glimpse()
Se eliminan cerca de 300 instancias. Lo que representa un 10% de pérdida aproximadamente. Un porcentaje dentro de lo permitido.
¿Cuál es la distribución de cada variable en el conjunto de datos?
Distribución de variables continuas
Para explorar la distribución de las variables cuantitativas las convertimos en un formato largo y generamos histogramas.
Code
# Seleccionar solo las variables numéricasdatos_numericos <- datos_limpios %>%select(-clase.uso.suelo) # Excluir la variable cualitativa nominal# Convertir el conjunto de datos en formato largodatos_long <-pivot_longer(datos_numericos, cols =everything())# Trazar histogramas facetas para todas las variables numéricasggplot(datos_long, aes(x = value)) +geom_histogram(binwidth =5, fill ="skyblue", color ="black") +facet_wrap(~name, scales ="free") +labs(title ="Distribución de variables cuantitativas", x ="Valor", y ="Frecuencia") +theme_minimal()
Distribución de variables cualitativas
Para la variable cualitativa clase.uso.suelo, utilizamos un gráfico de barras para visualizar su distribución.
Code
# Gráfico de barras para la variable 'clase.uso.suelo'ggplot(datos_limpios, aes(x = clase.uso.suelo)) +geom_bar(fill ="lightgreen", color ="black") +labs(title ="Distribución de clases en variable uso del suelo", x ="Clases", y ="Frecuencia") +theme_minimal()
A partir de la distribución de las variables, respondemos la pregunta: ¿Existen valores fuera de rango o inconsistentes en las variables y cómo manejarlos?
Si, hemos observado estos problemas. Identificamos valores fuera de rango o extremos en las variables de precipitación y temperatura media.
Code
# Corregir precipitación eliminando los casos negativosdatos_limpios <- datos_limpios %>%filter(precipitacion >=0)# Corregir temp.media eliminando temperaturas superiores a 50datos_limpios <- datos_limpios %>%filter(temp.media <=37) # Verificar la estructura de los datos limpiosdatos_limpios %>%glimpse()
Para las demás variables, aparentemente los datos tiene un comportamiento más adecuado.
En R, es común convertir variables cualitativas (categóricas) en factores, especialmente cuando se trabaja con modelos de machine learning. Los factores son útiles porque permiten que R trate las variables categóricas adecuadamente, ya sea para hacer análisis descriptivos, gráficos o modelado.
En nuestro caso, la variable clase.uso.suelo es una variable categórica que representa diferentes clases de uso del suelo (Bosque, Cultivos, Urbano, Pastizales). Para asegurar que esta variable se maneje correctamente en análisis posteriores, la convertimos en un factor.
Code
# Convertir la variable cualitativa 'clase.uso.suelo' en factordatos_limpios <- datos_limpios %>%mutate(clase.uso.suelo =factor(clase.uso.suelo, levels =c("Bosque", "Cultivos", "Urbano", "Pastizales"),labels =c(1, 2, 3, 4)))# Contar los registros que existen de cada clase de uso del suelodatos_limpios %>%select_if(is.factor) %>%count(clase.uso.suelo)
El análisis de correlación es fundamental para comprender cómo se relacionan las diferentes variables entre sí dentro de un conjunto de datos. Por ello se plantea la interrogante ¿Existen correlaciones entre pares de variables independientes en el conjunto de datos?
Primero, calculamos la matriz de correlación para las variables numéricas. La matriz de correlación muestra cómo cada variable se relaciona con las demás en términos de una escala que va de -1 a 1, donde 1 indica una correlación positiva perfecta, -1 indica una correlación negativa perfecta y 0 indica que no hay correlación.
Code
# Calcular la matriz de correlación excluyendo la variable categóricacorrelaciones <- datos_limpios %>%select(-clase.uso.suelo) %>%# Excluir la variable clasecor()# Convertir la matriz de correlación a un data frame para su visualizacióncorrelaciones_df <-as.data.frame(correlaciones) %>%rownames_to_column(var ="variable1") %>%pivot_longer(cols =-variable1, names_to ="variable2", values_to ="correlacion")# Trazar la matriz de correlación mediante un heatmap (mapa de calor)ggplot(correlaciones_df, aes(x = variable1, y = variable2, fill = correlacion)) +geom_tile() +scale_fill_gradient(low ="blue", high ="green") +labs(title ="Matriz de Correlación",x ="Variable 1", y ="Variable 2", fill ="Correlación") +theme_minimal() +theme(axis.text.x =element_text(angle =45, hjust =1),plot.title =element_text(hjust =0.5))
En el mapa de calor se aprecia que parece existir una correlación entre las variables reflectancia.rojo y reflectancia.infrarojo. Entonces se procede a verificar el par de variables correlacionadas considerando aquellas con correlación absoluta mayor que el umbral de 0.50.
Code
# Convertir la matriz de correlación a un data framecor_df <-as.data.frame(as.table(correlaciones))# Filtrar los pares de variables con una correlación absoluta mayor que un umbral (por ejemplo, 0.5)correlaciones_significativas <- cor_df %>%filter(abs(Freq) >0.5& Var1 != Var2)# Mostrar las correlaciones significativascorrelaciones_significativas
Como identificamos pares de variables con alta correlación, procedemos a eliminar una de ellas para reducir la multicolinealidad en nuestros modelos.
Code
# Eliminar la variables correlacionaddatos_limpios <- datos_limpios %>%select(-reflectancia.rojo)# Verificar la estructura de los datos limpiosdatos_limpios %>%glimpse()
Además de entender las correlaciones entre las variables independientes, es crucial analizar cómo cada una de ellas se relaciona con la variable objetivo (indice.quemado para regresión y clase.uso.suelo para clasificación).
¿Hay correlación entre cada variable independiente y la variable objetivo?
La idea es identificar qué variables independientes tienen la mayor influencia en la variable objetivo (indice.quemado) y (clase.uso.suelo). Esto es crucial para modelos predictivos y para entender los factores más relevantes que afectan el fenómeno estudiado. Variables con alta correlación con la variable objetivo son generalmente más importantes en los modelos predictivos.
Code
# Calcular la correlación de cada variable numérica con la variable objetivo 'indice.quemado'correlaciones_objetivo <- datos_limpios %>%select(-c(clase.uso.suelo, indice.quemado)) %>%# Excluir la variable categórica y la variable objetivocor(datos_limpios$indice.quemado) %>%# Calcular la correlación con la variable objetivoas.data.frame() %>%rownames_to_column(var ="Variable") %>%rename(correlacion_con_objetivo =2)# Mostrar las correlacionescorrelaciones_objetivo
Las correlaciones calculadas muestran coeficientes muy cercanos a cero, lo cual indica que no existe una relación lineal significativa entre las variables independientes analizadas y la variable indice.quemado. Este resultado sugiere que será difícil encontrar un modelo de regresión que ajuste adecuadamente los datos y produzca predicciones precisas.
Para las variables independientes y clase.uso.suelo se realiza el comportamiento para cada una de las clases:
Code
# correlacion entre variables independientes y la dependiente datos_limpios %>%select( c(1:11, 12)) %>%ggpairs(., legend =1, mapping = ggplot2::aes(colour=clase.uso.suelo), lower =list(continuous =wrap("smooth", alpha =0.3, size=0.1))) +theme(legend.position ="bottom")
De la gráfica se observa que ciertas variables parecen tener alguna relación (se muestrán con asteriscos) con ciertas clases de la variable dependiente.
Ahora se calcularán algunas medidas de tendencia central. Para ello se carga un archivo funciones.R que contiene funciones personalizadas para el cálculo de medidas de tendencia central. Posteriormente, se copian los datos limpios (datos_limpios) a datos_finales para realizar el análisis detallado.
Code
# Análisis Exploratorio de Datos (EDA)# Cargar funciones personalizadas desde el archivo funciones.Rsource("Scripts/funciones.R")# Copiar los datos limpios a un nuevo objeto para el análisis finaldatos_finales <- datos_limpios
Se calculan medidas de tendencia central:
Code
# Medidas de tendencia central y dispersión# Calcular media, mediana, varianza y desviación estándar para las variables numéricasestadisticas <- datos_finales %>%select_if(is.numeric) %>%summarize_all(list(media = mean,mediana = median,varianza = var,p.desviacion = sd)) %>%pivot_longer(cols =everything(),names_to =c("variable", "medida"),names_sep ="_",values_to ="valor")# Se almacena el resultado en medidas.central.disper medidas.central.disper <-spread(estadisticas, key="variable", value="valor" )
Se calculan las medidas de posición, con sus cuartiles (Q1, Q2, Q3, Q4) para cada variable numérica en datos_finales. Los cuartiles son medidas de posición que dividen el conjunto de datos en cuatro partes iguales, proporcionando información sobre la distribución de los datos y la presencia de valores atípicos (outliers). Se genera un diagrama de caja (boxplot) para visualizar las distribuciones de las variables numéricas (df_num) y detectar la presencia de outliers.
Code
# Medidas de posición (cuartiles)# Calcular los cuartiles (Q1, Q2, Q3, Q4) para las variables numéricascuartiles <- datos_finales %>%summarise(across(where(is.numeric),~quibble(.x, c(0.25, 0.50, 0.75, 1), dropNA =TRUE)),.groups ='drop') %>% tidyr::unnest(names_repair ='unique', names_sep ="_") %>%select(ends_with("_x")) %>%rotate_df() %>%rename(Q1 = V1, Q2 = V2, Q3 = V3, Q4 = V4)# Seleccionar solo las columnas que contienen variables numéricasdf_num <- datos_finales %>%select(where(is.numeric))# Convertir los datos a formato largo para usar ggplotdf_num <- reshape2::melt(df_num)# Diagrama de caja para visualizar distribuciones y outliersdf_num %>%ggplot(aes(x = variable, y = value)) +geom_boxplot(fill ="dodgerblue", color ="black", outlier.colour ="red") +labs(x ="Variable", y ="Valor") +facet_grid(. ~ variable, scales ="free_x")
A partir de lo anterior se observan algunos outliers (puntos rojos). Eliminaremos eso para la variable temp.max.
Code
# Eliminar outliers en la variable temp.maxdatos_finales <- datos_finales %>%mutate(ruido =es_outlier(temp.max)) %>%filter(ruido ==FALSE) %>%select(-ruido)
Finalmente realizamos la normalización de las variables mediante el escalado por z-score. Esto se realiza con el objetivo de estandarizar las variables numéricas, asegurando que todas las variables estén en la misma escala, lo que facilita la comparación y el aprendizaje de los modelos de machine learning.
Code
# Transformación o escalado de variables por z-score# Paso 1: Selección de variables numéricas del conjunto de datosdatos_escalados <- datos_finales %>%select(where(is.numeric)) %>%# Paso 2: Aplicación del escalado z-score a todas las columnas numéricasmutate(across(c(1:11), ~escalado.z(.x), .names ='nueva_{.col}'))# Paso 3: Selección de variables escaladas por z-score y combinación con las variables originalesdatos_finales_escalados <- datos_escalados %>%select(starts_with("nueva", ignore.case =TRUE, vars =NULL)) %>%bind_cols(datos_finales %>%select(clase.uso.suelo, indice.quemado))# Mostrar el resultado final del conjunto de datos escaladosdatos_finales_escalados %>%glimpse()
Para concluir las tres primeras etapas del proceso de Machine Learning, guardamos los datos procesados en un archivo .csv.
Code
# Guardar el resultado en un archivo CSVwrite_csv(datos_finales_escalados, "Datos/datos_finales_escalados.csv")
Etapa 4 y 5. Modelado y Validación
Antes del modelo vamos a cargar los datos procesados y limpios:
Code
# Carga de datos preprocesados# Leer los datos desde un archivo CSVdatos_procesados <- readr::read_delim("Datos/datos_finales_escalados.csv", delim =",")# Mostrar una vista preliminar de los datosdatos_procesados %>%glimpse()
En estas etapas inicialmente se realizará la separación del conjunto de datos en entrenamiento y prueba, para luego aplicar procesos de selección de características, modelado y validación.
Separación de datos
Con esta separación se asegura que podamos evaluar el rendimiento del modelo en datos no vistos, proporcionando una estimación realista de su capacidad de generalización del algoritmo de clasificación.
Code
# Se separa en 70% entrenamiento y 30% prueba entrenamiento_ids <-sample(1:nrow(datos_procesados), nrow(datos_procesados)*0.7, replace =FALSE)# Crea dataset de entrenamiento y prueba ds_entrenamiento <- datos_procesados[entrenamiento_ids, ] ds_prueba <- datos_procesados[-entrenamiento_ids, ] # Convertir a factor la variable clase en ambos conjuntos de datos ds_entrenamiento$clase.uso.suelo <-factor(ds_entrenamiento$clase.uso.suelo) ds_prueba$clase.uso.suelo <-factor(ds_prueba$clase.uso.suelo)
Mostramos una vista preliminar de los datos separados:
Este proceso consiste en seleccionar las características (variables) más importantes de un dataset para facilitar y mejorar el aprendizaje del modelo. En nuestro caso, el conjunto de datos contiene 11 variables independientes y una variable objetivo clase.uso.suelo con 4 posibles valores. El objetivo del modelo es predecir la clase uso de suelo (bosque, cultivo, urbano, pastizales) de la parcela en función de las 11 características. Para ello identificamos, a través de 2 algoritmos, las variables que parecen ser más importantes para clasificación de la variable clase.uso.suelo.
Algoritmo 1. Boruta
La idea central de Boruta es comparar la importancia de las variables observadas con la importancia de variables aleatorias, también conocidas como “variables sombra” o “shadow features”, que se generan a partir de permutaciones aleatorias de los valores de las variables originales.
Code
# Dejar solo variables independientes (sin indice quemado) + clase.uso.suelodatos_boruta <- ds_entrenamiento %>%select(-indice.quemado)# Variables independientesx <- datos_boruta %>%select(-clase.uso.suelo)# Variable dependientey <- datos_boruta$clase.uso.suelo# Garantizamos que las variables categóricas estén como factoresy <-as.factor(y)# Aplicar Borutamodelo_boruta <-Boruta(x, y)
Los resultados obtenidos por Boruta son los siguientes:
Code
# Resultados de Borutamodelo_boruta
Boruta performed 46 iterations in 7.88006 secs.
1 attributes confirmed important: nueva_precipitacion;
10 attributes confirmed unimportant: nueva_altitud, nueva_latitud,
nueva_longitud, nueva_ndbi, nueva_ndvi and 5 more;
La información nos indica que Boruta llevó a cabo 44 iteraciones para evaluar la importancia de las características. En cada iteración, el algoritmo compara la importancia de las variables originales con las variables sombra (shadow features). Boruta identificó nueva_precipitacion como una característica importante. Esto significa que esta variable es útil para predecir la variable objetivo (clase.uso.suelo) y debe ser retenida en el modelo.
Podemos generar un gráfico de resultados:
Code
# Gráfico de importancia de características Borutaplot(modelo_boruta, sort =TRUE, las =2, cex.axis =0.5) #Muestra ordenadas las variables por importancia #Ejes deben estar en una orientación vertical.
Code
# Verde: Características confirmadas como importantes.# Rojo: Características confirmadas como no importantes.# Amarillo: Características de importancia incierta (tentativas).
Se observa claramente que la variable seleccionada por Boruta como importante (color verde) es nueva_precipitacion.
Algoritmo 2. Recursive Feature Elimination (RFE)
Es un método que entrena un modelo con todas las variables independientes y luego elimina las variables menos importantes una a una, hasta que se alcance el número deseado de variables.
Code
# Preparar los datos para RFE. Usa validación cruzada con k=5ctrl <-rfeControl(functions = caretFuncs, method ="cv", number =5)# Dejar solo variables independientes (sin indice quemado) + clase.uso.suelodatos_rfe <- ds_entrenamiento %>%select(-indice.quemado)# Variables independientesx <- datos_rfe %>%select(-clase.uso.suelo)# Variable dependientey <- datos_rfe$clase.uso.suelo# Garantizamos que las variables categóricas estén como factoresy <-as.factor(y)# Aplicar RFEmodelo_rfe <-rfe(x = x, y = y, sizes =c(1:11), # Un vector numérico de enteros que indica la cantidad de características que se deben conservar.rfeControl = ctrl) #Parámetros de RFE libreria caret
note: only 2 unique complexity parameters in default grid. Truncating the grid to 2 .
note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .
note: only 2 unique complexity parameters in default grid. Truncating the grid to 2 .
note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .
note: only 2 unique complexity parameters in default grid. Truncating the grid to 2 .
note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .
note: only 2 unique complexity parameters in default grid. Truncating the grid to 2 .
note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .
note: only 2 unique complexity parameters in default grid. Truncating the grid to 2 .
note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .
Los resultados obtenidos por RFE son los siguientes:
Code
# Imprimir los resultados de la selección de característicasprint(modelo_rfe)
La información nos indica que RFE utilizó RFE con validación cruzada de 5 particiones (5-fold cross-validation) para evaluar el rendimiento del modelo al seleccionar diferentes tamaños de subconjuntos de características. Además, se muestra las 5 principales variables seleccionadas entre un total de 10. Estas son las características que RFE ha identificado como las más importantes para predecir la variable objetivo (clase.uso.suelo)
Podemos generar un gráfico de resultados:
Code
# Gráfico de características seleccionadasplot(modelo_rfe, type =c("g", "o"))
Code
# Eje X (Número de Características): muestra el número de características utilizadas en el modelo. # En este caso, puede variar desde 1 hasta 11, según el parámetro sizes = c(1:11).# Eje Y (Rendimiento del Modelo): muestra una métrica de rendimiento del modelo, # como la precisión (accuracy), la cual se calcula mediante validación cruzada.
El gráfico muestra cómo cambia el rendimiento del modelo a medida que se utilizan diferentes cantidades de características. A partir del gráfico, se puede determinar cuántas características son necesarias para lograr un buen rendimiento del modelo sin incluir características irrelevantes.
Finalmenre RFE indica que las variables relevantes son 10, esto es todas las independientes excepto altitud. Tal como se puede observar con el siguiente código:
Code
# Obtener las variables relevantes seleccionadas por RFErfe_vars <- modelo_rfe$optVariables rfe_vars
Se utilizarán dos técnicas habituales para procesos de clasificación: árboles de decisión y máquinas de soporte vectorial (SVM).
Modelado con árboles de decisión
Utilizaremos árboles de decisión para el problema de clasificación multiclase donde clase.uso.suelo tiene 4 niveles. En este caso, se usará la función rpart del paquete rpart de R, especificando la medida de selección de atributos “ganancia de información”. Se considerarán tres escenarios: 1) usando todas las variables independientes, 2) variables independientes importantes según Boruta y 3) variables importantes según RFE:
Escenario 1. Aplicar árboles de decisión usando todas las variables + clase.uso.suelo
Code
# Modelo con todas las variablestodas_vars <-colnames(ds_entrenamiento)[colnames(ds_entrenamiento) !="clase.uso.suelo"&colnames(ds_entrenamiento) !="indice.quemado"]formula_todas <-as.formula(paste("clase.uso.suelo ~", paste(todas_vars, collapse =" + ")))modelo_todas <-rpart(formula_todas, data = ds_entrenamiento, method ="class", #indica que se realizará clasificaciónparms =list(split ="information")) #information (otro es gini) utiliza la ganancia de información como criterio de división.# Evaluar el modelopredicciones_todas <-predict(modelo_todas, ds_prueba, type ="class")confusionMatrix(predicciones_todas, ds_prueba$clase.uso.suelo)
Confusion Matrix and Statistics
Reference
Prediction 1 2 3 4
1 138 97 0 96
2 0 0 0 0
3 0 2 150 1
4 38 70 4 97
Overall Statistics
Accuracy : 0.5556
95% CI : (0.5177, 0.593)
No Information Rate : 0.2799
P-Value [Acc > NIR] : < 0.00000000000000022
Kappa : 0.4036
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: 1 Class: 2 Class: 3 Class: 4
Sensitivity 0.7841 0.0000 0.9740 0.5000
Specificity 0.6267 1.0000 0.9944 0.7756
Pos Pred Value 0.4169 NaN 0.9804 0.4641
Neg Pred Value 0.8950 0.7561 0.9926 0.7996
Prevalence 0.2540 0.2439 0.2222 0.2799
Detection Rate 0.1991 0.0000 0.2165 0.1400
Detection Prevalence 0.4776 0.0000 0.2208 0.3016
Balanced Accuracy 0.7054 0.5000 0.9842 0.6378
La exactitud (Accuracy) del modelo indica que ha clasificado correctamente el 55.56% de las instancias totales y la exactitud está entre el 51.77% y el 59.3%. El No Information Rate (NIR) es la proporción de la clase más grande en el conjunto de datos de prueba. En este caso, la clase mayoritaria representa el 27.99% de las instancias. Es un punto de referencia para evaluar si el modelo está haciendo mejor que una clasificación al azar que siempre predice la clase mayoritaria. El valor p extremadamente bajo (< 0.00000000000000022) indica que podemos rechazar la hipótesis nula y concluir que la precisión del modelo es significativamente mejor que el NIR.
La estadística Kappa mide la concordancia entre las predicciones del modelo y las clases verdaderas, ajustada por la concordancia que ocurriría por azar. Un valor de 0.403 indica una concordancia justa. Se suelen considerar las siguientes equivalencias:
Rango de Kappa
Interpretación
0
Concordancia equivalente al azar
0.01 ≤ Kappa ≤ 0.20
Sin concordancia representativa
0.21 ≤ Kappa ≤ 0.39
Concordancia mínima
0.40 ≤ Kappa ≤ 0.59
Concordancia justa
0.60 ≤ Kappa ≤ 0.79
Concordancia moderada
0.80 ≤ Kappa ≤ 0.90
Concordancia fuerte
Sobre 0.90
Concordancia casi perfecta
El valor p de la prueba de McNemar se usa para evaluar si hay una diferencia significativa en las proporciones de errores cometidos por dos clasificadores. En este contexto, no está disponible (NA), porque no es aplicable para el análisis actual.
Ahora mostraremos el árbol de decisión:
Code
# Visualización del árbol de decisión para el modelo con todas las variablesrpart.plot(modelo_todas, type=4)
Cada rectángulo representa un nodo de nuestro árbol de decisión y contiene su regla de clasificación. El color del nodo refleja la clase mayoritaria entre los datos que agrupa, es decir, la clase predicha por el modelo para ese grupo. Dentro de cada nodo, se muestra la proporción de instancias pertenecientes a cada clase y la proporción del total de datos agrupados en ese nodo.
Por ejemplo, el rectángulo anaranjado de la izquierda (un nodo terminal) tiene 44% de casos en la clase 1, 31% en la clase 2, 0% en la clase 3, y 24% en la clase 4. Este nodo representa el 48% de todos los datos. Estas proporciones nos brindan una idea de la precisión del modelo al hacer predicciones.
Escenario 2. Aplicar árboles de decisión usando las variables importantes Boruta + clase.uso.suelo
Code
# Modelo con variables Boruta formula_boruta <-as.formula(paste("clase.uso.suelo ~", paste(boruta_vars, collapse =" + ")))modelo_boruta <-rpart(formula_boruta, data = ds_entrenamiento,method ="class", parms =list(split ="information"))# Evaluar el modelopredicciones_boruta <-predict(modelo_boruta, ds_prueba, type ="class")confusionMatrix(predicciones_boruta, ds_prueba$clase.uso.suelo)
Confusion Matrix and Statistics
Reference
Prediction 1 2 3 4
1 138 97 0 96
2 0 0 0 0
3 0 2 150 1
4 38 70 4 97
Overall Statistics
Accuracy : 0.5556
95% CI : (0.5177, 0.593)
No Information Rate : 0.2799
P-Value [Acc > NIR] : < 0.00000000000000022
Kappa : 0.4036
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: 1 Class: 2 Class: 3 Class: 4
Sensitivity 0.7841 0.0000 0.9740 0.5000
Specificity 0.6267 1.0000 0.9944 0.7756
Pos Pred Value 0.4169 NaN 0.9804 0.4641
Neg Pred Value 0.8950 0.7561 0.9926 0.7996
Prevalence 0.2540 0.2439 0.2222 0.2799
Detection Rate 0.1991 0.0000 0.2165 0.1400
Detection Prevalence 0.4776 0.0000 0.2208 0.3016
Balanced Accuracy 0.7054 0.5000 0.9842 0.6378
Como se observa, los resultados son similares al escenario 1, esto debido a que justamente la variable más relevante encontrada por Boruta, coincide con el atributo que mejor separa las clases según el árbol de decisión. De hecho, esto ocasiona que se genere un árbol similar:
Code
# Visualización del árbol de decisión para el modelo con todas las variablesrpart.plot(modelo_boruta, type=1)
Escenario 3. Aplicar árboles de decisión usando las variables importantes RFE + clase.uso.suelo
Code
# Modelo con variables RFE formula_rfe <-as.formula(paste("clase.uso.suelo ~", paste(rfe_vars, collapse =" + ")))modelo_rfe <-rpart(formula_rfe, data = ds_entrenamiento, method ="class", parms =list(split ="information"))# Evaluar el modelopredicciones_rfe <-predict(modelo_rfe, ds_prueba, type ="class")confusionMatrix(predicciones_rfe, ds_prueba$clase.uso.suelo)
Confusion Matrix and Statistics
Reference
Prediction 1 2 3 4
1 138 97 0 96
2 0 0 0 0
3 0 2 150 1
4 38 70 4 97
Overall Statistics
Accuracy : 0.5556
95% CI : (0.5177, 0.593)
No Information Rate : 0.2799
P-Value [Acc > NIR] : < 0.00000000000000022
Kappa : 0.4036
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: 1 Class: 2 Class: 3 Class: 4
Sensitivity 0.7841 0.0000 0.9740 0.5000
Specificity 0.6267 1.0000 0.9944 0.7756
Pos Pred Value 0.4169 NaN 0.9804 0.4641
Neg Pred Value 0.8950 0.7561 0.9926 0.7996
Prevalence 0.2540 0.2439 0.2222 0.2799
Detection Rate 0.1991 0.0000 0.2165 0.1400
Detection Prevalence 0.4776 0.0000 0.2208 0.3016
Balanced Accuracy 0.7054 0.5000 0.9842 0.6378
El rendimiento del árbol de decisión en este escenario es similar nuevamente. Lo que parece ser un indicativo que la única variable necesaria para el proceso de clasificación es nueva_precipitacion con el siguiente árbol:
Code
# Visualización del árbol de decisión para el modelo con todas las variablesrpart.plot(modelo_rfe, type=1)
Modelado con SVM
Ahora se aplicará SVM para el problema de clasificación multiclase donde clase.uso.suelo tiene 4 niveles. En este caso se usará la función svm del paquete e1071. Se considerarán tres escenarios: 1) usando todas las variables independientes, 2) variables independientes importantes según Boruta y 3) variables importantes según RFE:
Escenario 1. Aplicar SVM usando todas las variables + clase.uso.suelo
Code
# Modelo con todas las variablescaso_todas_vars <-colnames(ds_entrenamiento)[colnames(ds_entrenamiento) !="clase.uso.suelo"&colnames(ds_entrenamiento) !="indice.quemado"]resultados_todas <-procesa_svm(ds_entrenamiento, ds_prueba, caso_todas_vars)print("Modelo con todas las variables:")
La exactitud (Accuracy) del modelo indica que ha clasificado correctamente el 57.0% de las instancias totales y la exactitud está entre el 53.22% y el 60.72%. El No Information Rate (NIR) es la proporción de la clase más grande en el conjunto de datos de prueba. En este caso, la clase mayoritaria representa el 26.7% de las instancias. El valor p extremadamente bajo (< 0.00000000000000022) indica que podemos rechazar la hipótesis nula y concluir que la precisión del modelo es significativamente mejor que el NIR.
La estadística Kappa con un valor de 0.4266 indica una concordancia justa. El valor p de la prueba de McNemar no está disponible (NA), porque no es aplicable para el análisis actual.
Escenario 2. Aplicar SVM usando las variables importantes Boruta + clase.uso.suelo
Code
# Modelo con variables Boruta resultados_boruta <-procesa_svm(ds_entrenamiento, ds_prueba, boruta_vars) print(resultados_boruta$train_confusion)
Para el escenario 2, se ha clasificado correctamente el 58.01% de las instancias totales con un intervalo de entre el 54.23% y el 61.71%. El valor p extremadamente bajo (< 0.00000000000000022) indica que podemos rechazar la hipótesis nula y concluir que la precisión del modelo es significativamente mejor que el NIR. La estadística Kappa de 0.4399 indica una concordancia justa. Al igual que en el primer escenario, la prueba de McNemar no es aplicable.
Escenario 3. Aplicar SVM usando las variables importantes RFE + clase.uso.suelo
Code
# Modelo con variables RFE resultados_RFE <-procesa_svm(ds_entrenamiento, ds_prueba, rfe_vars) print(resultados_RFE$train_confusion)
Finalmente, para el escenario 3, se ha clasificado correctamente el 55.7% de las instancias totales con un intervalo de entre el 51.91% y el 59.44%. El valor p extremadamente bajo (< 0.00000000000000022) indica que podemos rechazar la hipótesis nula y concluir que la precisión del modelo es significativamente mejor que el NIR. La estadística Kappa de 0.409 indica una concordancia justa. Al igual que en el primer escenario, la prueba de McNemar no es aplicable.
Al comparar los tres escenarios, observamos que en general el rendimiento de la clasificación no es el mejor. Partiendo de este hecho, el modelo que utiliza las variables seleccionadas por Boruta (Escenario 2) muestra el mejor rendimiento en términos de exactitud y concordancia, seguido por el modelo con todas las variables (Escenario 1), y finalmente el modelo con las variables seleccionadas por RFE (Escenario 3). Si se tendría que escoger una opción, sería el modelo SVM aplicado en el escenario 2 (Boruta).
Evaluación del mejor modelo
Del mejor modelo de SVM se obtienen los hiperparámetros ajustados y se genera las clasificaciones para el conjunto de prueba (o validación en caso de existir). Los hiperparámetros de una SVM son:
Parámetro
Descripción
C
Parámetro de regularización. Controla el equilibrio entre el margen grande y clasificaciones correctas. Se utiliza para kernel lineal, polinómico y radial. El valor predeterminado es 1. Rango: [0.1, 100].
gamma
Coeficiente del kernel. Controla el alcance de influencia de una única muestra de entrenamiento. Se utiliza para el kernel radial. El valor predeterminado es 1/(dimensión de los datos). Rango típico: [0.001, 10].
degree
Grado del kernel polinómico (solo para kernel polinómico). El valor predeterminado es 3. Rango: [2, 5].
coef0
Término independiente en la función kernel (solo para kernel polinómico). El valor predeterminado es 0. Rango: [-1, 1].
Code
# Obtener los hiperparámetros del mejor modelo SVMmejor_modelo <- resultados_boruta$modelprint(mejor_modelo)
Call:
svm(formula = clase.uso.suelo ~ ., data = train_datos[, c(vars, "clase.uso.suelo")],
type = "C-classification", kernel = "linear")
Parameters:
SVM-Type: C-classification
SVM-Kernel: linear
cost: 1
Number of Support Vectors: 1194
El mejor modelo corresponde a una SVM con función de Kernel de base radial (Gaussiano). Para obtener los hiperparámetros usamos:
Code
C = mejor_modelo$costgamma = mejor_modelo$gammaprint(C)
[1] 1
Code
print(gamma)
[1] 1
Por lo tanto, si se usa una SVM para clasificar estos datos y obtener el rendimiento analizado previamente se deben configurar los siguientes hiperparámetros con C=1, gamma=1, con función de Kernel radial.
Ninguna de las variables independientes tiene un p-valor significativo a un nivel de significancia común (por ejemplo, 0.05), lo que sugiere que ninguna de las variables independientes está significativamente asociada con la variable indice.quemado en este modelo. El R² negativo (-0.01145242) indica que el modelo de regresión lineal no es adecuado para explicar la variabilidad en el indice.quemado utilizando las variables predictoras proporcionadas. Aunque el RMSE (0.2886069) no es extremadamente alto, la falta de significancia en la mayoría de las variables y el bajo R² sugieren que el modelo no está capturando bien las relaciones entre las variables predictoras y la variable dependiente.
Por ello, se exploran otros modelos de predicción: Lasso, XGBoost.
Entrenamiento del modelo de regresión Lasso
Code
# Modelo Lasso # Preparar los datosx_train <- ds_entrenamiento %>%select(-c(indice.quemado, clase.uso.suelo)) %>%as.matrix()y_train <- ds_entrenamiento$indice.quemadox_test <- ds_prueba %>%select(-c(indice.quemado, clase.uso.suelo))%>%as.matrix()y_test <- ds_prueba$indice.quemado# Ajustar el modelo Lasso utilizando validación cruzada para seleccionar el mejor lambdamodelo_lasso <-cv.glmnet(x_train, y_train, alpha =1, lambda =10^seq(10, -2, length =100))modelo_lasso
Call: cv.glmnet(x = x_train, y = y_train, lambda = 10^seq(10, -2, length = 100), alpha = 1)
Measure: Mean-Squared Error
Lambda Index Measure SE Nonzero
min 0 100 0.08217 0.001389 3
1se 10000000000 1 0.08235 0.001398 0
Code
# Obtener las predicciones en el conjunto de pruebapredicciones_lasso <-predict(modelo_lasso, s ="lambda.min", newx = x_test)# Evaluación del modelomse_lasso <-mean((predicciones_lasso - y_test)^2)rmse_lasso <-sqrt(mse_lasso)r2_lasso <-1-sum((predicciones_lasso - y_test)^2) /sum((y_test -mean(y_test))^2)cat("MSE (Lasso):", mse_lasso, "\n")
MSE (Lasso): 0.08737396
Code
cat("RMSE (Lasso):", rmse_lasso, "\n")
RMSE (Lasso): 0.2955909
Code
cat("R² (Lasso):", r2_lasso, "\n")
R² (Lasso): -0.004957665
La regresión Lasso tampoco es capaz de encontrar una relación fuerte con la variable de respuesta.
Entrenamiento del modelo
XGBoost (Extreme Gradient Boosting) es un algoritmo de aprendizaje automático basado en árboles de decisión que funciona mediante la construcción secuencial de una serie de árboles, donde cada árbol intenta corregir los errores de los árboles anteriores.
Utiliza una técnica llamada boosting, que ajusta iterativamente nuevos modelos sobre los residuos de los modelos anteriores, con el objetivo de minimizar una función específica (como el error cuadrático medio para problemas de regresión.)
Code
# Preparar los datosx_train <- ds_entrenamiento %>%select(-c(indice.quemado, clase.uso.suelo)) %>%as.matrix()y_train <- ds_entrenamiento$indice.quemadox_test <- ds_prueba %>%select(-c(indice.quemado, clase.uso.suelo)) %>%as.matrix()y_test <- ds_prueba$indice.quemado# Crear matrices DMatrix para XGBoost (Formato requerido por el algoritmo)dtrain <-xgb.DMatrix(data = x_train, label = y_train)dtest <-xgb.DMatrix(data = x_test, label = y_test)# Parámetros del modeloparams <-list(booster ="gbtree", # Usar árboles de decisiónobjective ="reg:squarederror", # Objetivo de regresióneta =0.1, # Tasa de aprendizajemax_depth =6, # Profundidad máxima de los árbolessubsample =0.8, # Submuestreo de los datos de entrenamientocolsample_bytree =0.8# Submuestreo de las columnas)# Ajustar el modelomodelo_xgb <-xgb.train(params = params,data = dtrain,nrounds =100, # Número de iteraciones de boostingwatchlist =list(train = dtrain, test = dtest), # Para monitorear el rendimientoearly_stopping_rounds =10, # Detener el entrenamiento temprano si no mejoraprint_every_n =10# Imprimir el error cada 10 iteraciones)
[1] train-rmse:0.283991 test-rmse:0.295568
Multiple eval metrics are present. Will use test_rmse for early stopping.
Will train until test_rmse hasn't improved in 10 rounds.
[11] train-rmse:0.256821 test-rmse:0.297780
Stopping. Best iteration:
[2] train-rmse:0.281756 test-rmse:0.295406
Code
# Predicciones en el conjunto de pruebapredicciones_xgb <-predict(modelo_xgb, newdata = dtest)# Evaluación del modelomse_xgb <-mean((predicciones_xgb - y_test)^2)rmse_xgb <-sqrt(mse_xgb)r2_xgb <-1-sum((predicciones_xgb - y_test)^2) /sum((y_test -mean(y_test))^2)cat("MSE (XGBoost):", mse_xgb, "\n")
MSE (XGBoost): 0.08726466
Code
cat("RMSE (XGBoost):", rmse_xgb, "\n")
RMSE (XGBoost): 0.2954059
Code
cat("R² (XGBoost):", r2_xgb, "\n")
R² (XGBoost): -0.003700474
XGBoost tampoco logra obtener buen rendimiento, de hecho el R² negativo (-0.01228988) es un signo de alerta. Significa que el modelo no está explicando bien la variabilidad en los datos.
