Perceptrón Multicapa para la Clasificación de Severidad de Verticillium sp. mediante Índices Espectrales
library(readr)
df <- read_csv(
"C:/Users/Elias Silva/Downloads/ANDREA/metodos/ntrain_median_values.csv"
)## Rows: 212 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (6): severity, ndvi_med, evi_med, ndre_med, gli_med, height_med
##
## ℹ 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.View(df)# Instalar solo la primera vez
# Cargar paquetes
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ purrr 1.2.2
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(psych)##
## Adjuntando el paquete: 'psych'
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(knitr)
# Revisar la estructura
str(df)## spc_tbl_ [212 × 6] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ severity : num [1:212] 2 1 3 1 2 1 3 1 2 1 ...
## $ ndvi_med : num [1:212] 0.878 0.865 0.877 0.879 0.893 ...
## $ evi_med : num [1:212] 0.861 0.838 0.783 0.861 0.84 ...
## $ ndre_med : num [1:212] 0.348 0.309 0.35 0.355 0.366 ...
## $ gli_med : num [1:212] 0.397 0.423 0.399 0.402 0.429 ...
## $ height_med: num [1:212] 0.497 0.681 0.847 0.759 0.713 ...
## - attr(*, "spec")=
## .. cols(
## .. severity = col_double(),
## .. ndvi_med = col_double(),
## .. evi_med = col_double(),
## .. ndre_med = col_double(),
## .. gli_med = col_double(),
## .. height_med = col_double()
## .. )
## - attr(*, "problems")=<pointer: 0x00000265aedffd70>head(df)names(df)## [1] "severity" "ndvi_med" "evi_med" "ndre_med" "gli_med"
## [6] "height_med"2. Parte I: Clasificación múltiple de severidad con perceptrón multicapa
2.1 Exploración y preprocesamiento
- Realice un análisis exploratorio de los cinco índices espectrales y la altura media de la planta. Reporte estadísticos descriptivos por categoría de severidad, incluyendo media, desviación estándar, mínimo y máximo. ¿Observa diferencias visibles entre categorías? Justifique si considera necesaria alguna transformación o estandarización de las variables antes de entrenar la red.
Tabla 1. Estadisticos descriptivos por categoria y severidad
##Estadisticos descriptivos por categoria y severidad
tabla_descriptiva <- df %>%
group_by(severity) %>%
summarise(
n = n(),
across(
.cols = c(
ndvi_med,
evi_med,
ndre_med,
gli_med,
height_med
),
.fns = list(
media = ~ mean(.x, na.rm = TRUE),
desviacion = ~ sd(.x, na.rm = TRUE),
minimo = ~ min(.x, na.rm = TRUE),
maximo = ~ max(.x, na.rm = TRUE)
),
.names = "{.col}_{.fn}"
),
.groups = "drop"
)
tabla_descriptivatabla_descriptiva_redondeada <- tabla_descriptiva %>%
mutate(
across(
where(is.numeric),
~ round(.x, 3)
)
)
tabla_descriptiva_redondeadaEl número de observaciones por categoría es desigual. Las categorías 2, 4 y 5 tienen entre 45 y 47 observaciones, mientras que la categoría 6 solo cuenta con 13. Este desbalance se debe tener en cuenta al entrenar la red, ya que la categoría 6 podría estar menos representada y, por tanto, ser clasificada con menor precisión.
## Diferencias entre categorias
df_largo <- df %>%
pivot_longer(
cols = c(
ndvi_med,
evi_med,
ndre_med,
gli_med,
height_med
),
names_to = "variable",
values_to = "valor"
)ggplot(
df_largo,
aes(
x = severity,
y = valor
)
) +
geom_boxplot(
aes(fill = severity),
alpha = 0.6,
outlier.shape = NA
) +
geom_jitter(
width = 0.15,
alpha = 0.4,
size = 1
) +
facet_wrap(
~ variable,
scales = "free_y"
) +
labs(
title = "Variación individual por categoría de severidad",
x = "Categoría de severidad",
y = "Valor"
) +
theme_minimal() +
theme(
legend.position = "none"
)## Warning: Orientation is not uniquely specified when both the x and y aesthetics are
## continuous. Picking default orientation 'x'.## Warning: Continuous x aesthetic
## ℹ did you forget `aes(group = ...)`?## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
Figura 1. Variación individual por categoria de la severidad
- El NDVI presenta una disminución progresiva y consistente conforme aumenta la severidad indicando una relación inversa: a mayor severidad, menor NDVI. Además, las diferencias se hacen más notorias a partir de las categorías 3 y 4 como se observa en la Fig. 1.La desviación estándar aumenta con la severidad (Tabla 1), desde 0,011 en la categoría 1 hasta 0,074 en la categoría 6. Esto significa que las plantas con baja severidad presentan valores de NDVI muy homogéneos, mientras que las plantas más afectadas muestran una respuesta espectral más variable.
Sin embargo, los rangos mínimo–máximo evidencian cierto solapamiento. Por ejemplo, la categoría 6 alcanza un máximo de 0,814, valor que se encuentra dentro del rango de categorías menos severas. Por tanto, aunque las medias están bien ordenadas, el NDVI por sí solo no permitiría clasificar perfectamente todos los individuos.
- El EVI también disminuye de manera general con el aumento de la severidad, la tendencia general es descendente, pero no es completamente monotónica, porque la categoría 5 presenta una media superior a la categoría 4.
La característica más importante del EVI es su elevada dispersión, especialmente en las categorías 2, 3 y 4. La desviación estándar llega a 0,316 en la categoría 3, valor muy alto respecto a su media. Los rangos también son extremadamente amplios. Por ejemplo, en la categoría 3 los valores van aproximadamente de 0,001 a 0,863. Esto indica una superposición considerable entre categorías. Aunque el EVI diferencia bien la categoría 1 frente a las categorías más severas, parece tener menor capacidad para separar con precisión las categorías intermedias.
- El NDRE muestra una disminución progresiva muy clara al igual que el NDVI. El NDRE presenta una asociación inversa consistente con la severidad. Las categorías 1 y 2 tienen valores relativamente cercanos, pero desde la categoría 3 se observa una reducción más marcada.
Las desviaciones estándar son moderadas y relativamente bajas, entre 0,020 y 0,043. En comparación con EVI, el NDRE es mucho más estable dentro de cada categoría. Esto sugiere que podría ser uno de los predictores más útiles para la red neuronal.
No obstante, también existe cierto solapamiento entre categorías contiguas. Por ejemplo, las categorías 1 y 2 presentan rangos muy similares. En cambio, la separación parece más clara entre las categorías de severidad baja y alta.
- El GLI presenta el comportamiento menos ordenado, presentando una disminución desde la categoría 1 hasta la 3, pero posteriormente los valores aumentan en las categorías 4 y 5 y vuelven a disminuir en la categoría 6. Por tanto, no se observa una relación lineal ni monotónica entre GLI y severidad.
Además, las desviaciones estándar son relativamente elevadas en las categorías 2 a 6, y los rangos muestran una amplia superposición. En consecuencia, GLI parece tener una capacidad discriminante individual limitada. Sin embargo, no necesariamente debe excluirse, porque una red neuronal puede identificar relaciones no lineales y combinaciones entre GLI y otros índices.
- La altura tiende a disminuir con el incremento de la severidad desde la categoria 1 hasta la categoría 4. Este comportamiento sugiere que el incremento inicial de la severidad estuvo asociado con una reducción del crecimiento de las plantas. Sin embargo, en las categorías 5 y 6 la altura aumentó ligeramente hasta 0,300 y 0,338, respectivamente, por lo que no se observó una relación monotónica entre ambas variables.
Adicionalmente, la altura presentó una elevada variabilidad dentro de todas las categorías, con desviaciones estándar entre 0,195 y 0,263. Los amplios intervalos entre los valores mínimos y máximos evidenciaron una fuerte superposición entre categorías. Como resultado, plantas pertenecientes a diferentes niveles de severidad pueden presentar alturas similares.
## Revisión de la distribución
ggplot(
df_largo,
aes(x = valor)
) +
geom_histogram(
bins = 25,
color = "black",
fill = "grey80"
) +
facet_wrap(
~ variable,
scales = "free"
) +
labs(
title = "Distribución de las variables predictoras",
x = "Valor",
y = "Frecuencia"
) +
theme_minimal()
Figura 2. Distribución de las variables predictorias por índice.
Los histogramas muestran que las variables predictoras no presentan una distribución normal uniforme y que cada una tiene un comportamiento diferente.
Los histogramas mostraron diferencias importantes en la forma de distribución de las variables predictoras. NDVI y GLI presentaron una asimetría negativa, caracterizada por una alta concentración de observaciones en valores elevados y una cola hacia valores bajos.
El EVI mostró una distribución irregular y con varias concentraciones de observaciones, lo que sugiere una alta heterogeneidad y una posible superposición entre categorías.
La altura media presentó una distribución amplia y aproximadamente simétrica en términos globales, aunque con acumulaciones en diferentes intervalos.
El NDRE mostró la distribución más equilibrada y sin asimetrías pronunciadas.
asimetria <- df %>%
summarise(
across(
c(
ndvi_med,
evi_med,
ndre_med,
gli_med,
height_med
),
~ psych::skew(.x, na.rm = TRUE)
)
)
asimetriaEl análisis de asimetría mostró comportamientos diferentes entre las variables predictoras. NDVI y GLI presentaron asimetrías negativas marcadas como se observó en los histogramas (Fig. 2), con valores de −1,543 y −1,462, respectivamente. Esto indica que la mayor parte de las observaciones se concentró en valores altos, mientras que un menor número de registros formó una cola hacia valores bajos.
En el caso de NDVI, este comportamiento es consistente con la disminución del índice en las plantas con mayor severidad, sin embargo para GLI sus medias no disminuyen progresivamente con la severidad y existe una amplia superposición entre categorías, por lo que la cola negativa no implica necesariamente que GLI sea pueda predecir de manera correcta la severidad.
El EVI presentó una asimetría negativa moderada de −0,564, lo que refleja una concentración menos pronunciada hacia valores altos.
NDRE y la altura media mostraron valores de asimetría cercanos a cero, 0,043 y 0,031, respectivamente, indicando distribuciones globalmente equilibradas. Estos resultados muestran que algunas variables no presentan distribuciones simétricas; sin embargo, ello no impide su utilización en una red neuronal, dado que este tipo de modelo no requiere normalidad en los predictores.
A partir del análisis exploratorio, no se consideró necesaria una transformación matemática general de las variables predictoras. Aunque NDVI y GLI presentaron asimetrías negativas marcadas, este comportamiento puede estar asociado con la concentración de observaciones en niveles bajos e intermedios de severidad y con la reducción de los índices en las plantas más afectadas. Asimismo, las redes neuronales no requieren que los predictores sigan una distribución normal. No se recomendó una transformación logarítmica, debido a la presencia de valores iguales o cercanos a cero en EVI, GLI y altura, lo cual podría generar dificultades numéricas o amplificar excesivamente las diferencias entre valores pequeños.
Si se considera necesaria la estandarización de las variables, ya que, aunque se encuentran en intervalos similares, presentan diferencias importantes en sus medias y dispersiones. Se propone aplicar una estandarización tipo Z, de modo que cada predictor tenga media cero y desviación estándar uno. Este procedimiento permitirá que todas las variables contribuyan en una escala comparable, reducirá el predominio de aquellas con mayor variabilidad y favorecerá la estabilidad y convergencia de la red neuronal. Los parámetros de estandarización deberán calcularse exclusivamente con el conjunto de entrenamiento y aplicarse posteriormente, sin recalcularlos, a los conjuntos de validación y prueba.
- Verifique si las categorías de severidad están balanceadas. En caso de desequilibrio notable, proponga y justifique al menos una estrategia para manejarlo (sobremuestreo, submuestreo, ponderación de clases u otra). Implemente la estrategia elegida y explique su efecto esperado sobre las métricas de ajuste.
## Verificación del balance de categorias
library(tidyverse)
df <- df %>%
mutate(
severity = factor(severity)
)
str(df$severity)## Factor w/ 6 levels "1","2","3","4",..: 2 1 3 1 2 1 3 1 2 1 ...##Comparación con una distribución idealmente equilibrada
numero_clases <- nlevels(df$severity)
frecuencia_ideal <- nrow(df) / numero_clases
frecuencia_ideal## [1] 35.33333# Crear la tabla de frecuencias y proporciones
balance_clases <- df %>%
count(severity, name = "frecuencia") %>%
mutate(
proporcion = frecuencia / sum(frecuencia),
porcentaje = proporcion * 100
)
# Número de categorías
numero_clases <- nlevels(df$severity)
# Frecuencia esperada si todas las categorías estuvieran equilibradas
frecuencia_ideal <- nrow(df) / numero_clases
# Agregar comparación con la frecuencia ideal
balance_clases <- balance_clases %>%
mutate(
frecuencia_ideal = frecuencia_ideal,
diferencia = frecuencia - frecuencia_ideal,
diferencia_porcentual =
((frecuencia - frecuencia_ideal) / frecuencia_ideal) * 100
)
# Mostrar resultados
balance_clasesTabla 2. Balance de categorias y comparación con distribución idealmente equilibrada
Los resultados muestran un desequilibrio notable entre las categorías de severidad (Tabla 2). En una distribución ideal cada categoría tendría aproximadamente 35 observaciones; sin embargo, la categoría 6 solo representa el 6,13 % de los datos, con 13 registros, mientras que la categoría 2 alcanza el 22,17 %, con 47 registros. La clase mayoritaria contiene 3,62 veces más observaciones que la minoritaria, por lo que una red entrenada directamente podría favorecer las categorías 2, 4 y 5 y presentar baja sensibilidad para detectar la severidad 6.
Para manejar este desequilibrio se seleccionó la ponderación de clases ya que permite asignar una mayor penalización a los errores cometidos sobre las categorías menos representadas, sin duplicar ni eliminar observaciones. Esta estrategia resulta especialmente conveniente debido al tamaño reducido de la base de datos y al bajo número de registros de la categoría 6. A diferencia del sobremuestreo aleatorio, la ponderación disminuye el riesgo de sobreajuste asociado con la repetición de los mismos casos minoritarios. Se espera que aumente la sensibilidad y el F1 de las categorías menos frecuentes, así como el F1 macro y la exactitud balanceada.
library(dplyr)
# 1. Asegurar que severity sea una variable categórica
df <- df %>%
mutate(
severity = factor(severity)
)
# 2. División estratificada 80 % entrenamiento y 20 % prueba
set.seed(123)
indices_train <- unlist(
lapply(
split(seq_len(nrow(df)), df$severity),
function(indices) {
sample(
indices,
size = floor(0.80 * length(indices))
)
}
)
)
train <- df[indices_train, ]
test <- df[-indices_train, ]
# 3. Verificar la distribución de clases
table(train$severity)##
## 1 2 3 4 5 6
## 21 37 28 36 36 10table(test$severity)##
## 1 2 3 4 5 6
## 6 10 7 9 9 3# 4. Calcular pesos inversamente proporcionales a la frecuencia
pesos_clase <- train %>%
count(severity, name = "frecuencia") %>%
mutate(
peso = nrow(train) /
(nlevels(train$severity) * frecuencia)
)
pesos_claseLos pesos de clase fueron calculados de manera inversamente proporcional a la frecuencia de cada categoría dentro del conjunto de entrenamiento. La categoría 6, con solo 10 observaciones, recibió el mayor peso, igual a 2,80, mientras que las categorías más frecuentes recibieron pesos inferiores a uno. La categoría 3, con 28 observaciones, recibió un peso igual a uno, al coincidir con la frecuencia promedio esperada. Esta ponderación permite que la contribución total de cada clase a la función de pérdida sea aproximadamente equivalente, evitando que las categorías mayoritarias dominen el proceso de aprendizaje.Se espera que la estrategia mejore la sensibilidad y el F1 de las categorías menos representadas, especialmente la severidad 6, así como el F1 macro y la exactitud balanceada del modelo.
# 5. Agregar el peso a cada fila del conjunto de entrenamiento
train_ponderado <- train %>%
left_join(
pesos_clase %>%
select(severity, peso),
by = "severity"
)
head(train_ponderado)2.2 Partición de datos
- Implemente dos esquemas de partición:
- Entrenamiento / Prueba (por ejemplo, 80%–20%).
- Entrenamiento / Validación / Prueba (por ejemplo, 70%–15%–15%).
Para cada esquema, entrene un perceptrón multicapa con la misma arquitectura inicial y compare los resultados. Explique en detalle cuál es la ventaja del esquema de tres particiones frente al de dos, en particular con respecto al riesgo de sobreajuste y a la selección de hiperparámetros.
library(dplyr)
library(nnet)## Preparación de la base
df <- df %>%
mutate(
severity = factor(
severity,
levels = c("1", "2", "3", "4", "5", "6")
)
)
predictores <- c(
"ndvi_med",
"evi_med",
"ndre_med",
"gli_med",
"height_med"
)## Partición estratificada
particion_estratificada <- function(datos, proporcion, semilla = 123) {
set.seed(semilla)
indices <- unlist(
lapply(
split(seq_len(nrow(datos)), datos$severity),
function(indice_clase) {
numero_seleccionado <- floor(
proporcion * length(indice_clase)
)
sample(
indice_clase,
size = numero_seleccionado,
replace = FALSE
)
}
)
)
list(
seleccionados = datos[indices, , drop = FALSE],
restantes = datos[-indices, , drop = FALSE]
)
}## Estandarizar sin fuga de información
estandarizar_datos <- function(train, nuevos_datos, variables) {
medias <- sapply(
train[, variables, drop = FALSE],
mean,
na.rm = TRUE
)
desviaciones <- sapply(
train[, variables, drop = FALSE],
sd,
na.rm = TRUE
)
train_std <- train
nuevos_std <- nuevos_datos
train_std[, variables] <- scale(
train[, variables, drop = FALSE],
center = medias,
scale = desviaciones
)
nuevos_std[, variables] <- scale(
nuevos_datos[, variables, drop = FALSE],
center = medias,
scale = desviaciones
)
list(
train = train_std,
nuevos = nuevos_std,
medias = medias,
desviaciones = desviaciones
)
}## Calculo de pesos por clase
calcular_pesos <- function(datos_train) {
tabla_pesos <- datos_train %>%
count(severity, name = "frecuencia") %>%
mutate(
peso = nrow(datos_train) /
(nlevels(datos_train$severity) * frecuencia)
)
datos_train %>%
left_join(
tabla_pesos %>% select(severity, peso),
by = "severity"
)
}evaluar_modelo <- function(real, predicho) {
real <- factor(real, levels = levels(df$severity))
predicho <- factor(predicho, levels = levels(df$severity))
matriz <- table(
Real = real,
Predicho = predicho
)
exactitud <- sum(diag(matriz)) / sum(matriz)
sensibilidad <- diag(matriz) / rowSums(matriz)
precision <- diag(matriz) / colSums(matriz)
f1_clase <- 2 * precision * sensibilidad /
(precision + sensibilidad)
sensibilidad[is.nan(sensibilidad)] <- NA
precision[is.nan(precision)] <- NA
f1_clase[is.nan(f1_clase)] <- NA
exactitud_balanceada <- mean(
sensibilidad,
na.rm = TRUE
)
f1_macro <- mean(
f1_clase,
na.rm = TRUE
)
list(
matriz_confusion = matriz,
exactitud = exactitud,
sensibilidad = sensibilidad,
precision = precision,
f1_clase = f1_clase,
exactitud_balanceada = exactitud_balanceada,
f1_macro = f1_macro
)
}Esquema A: entrenamiento y prueba 80 %–20 %
## Crear partición
set.seed(123)
particion_80_20 <- particion_estratificada(
datos = df,
proporcion = 0.80,
semilla = 123
)
train_80 <- particion_80_20$seleccionados
test_20 <- particion_80_20$restantes## Verificar tamaños
dim(train_80)## [1] 168 6dim(test_20)## [1] 44 6table(train_80$severity)##
## 1 2 3 4 5 6
## 21 37 28 36 36 10table(test_20$severity)##
## 1 2 3 4 5 6
## 6 10 7 9 9 3round(
prop.table(table(train_80$severity)) * 100,
2
)##
## 1 2 3 4 5 6
## 12.50 22.02 16.67 21.43 21.43 5.95round(
prop.table(table(test_20$severity)) * 100,
2
)##
## 1 2 3 4 5 6
## 13.64 22.73 15.91 20.45 20.45 6.82## Estandarizar
escalado_80_20 <- estandarizar_datos(
train = train_80,
nuevos_datos = test_20,
variables = predictores
)
train_80_std <- escalado_80_20$train
test_20_std <- escalado_80_20$nuevos## Calcular pesos
train_80_ponderado <- calcular_pesos(
train_80_std
)
train_80_ponderado %>%
distinct(severity, peso) %>%
arrange(severity)## Entrenar la red
# Asegurar que la respuesta sea categórica
train_80_ponderado$severity <- factor(
train_80_ponderado$severity,
levels = c("1", "2", "3", "4", "5", "6")
)
set.seed(123)
modelo_80_20 <- nnet::nnet(
formula = severity ~ ndvi_med + evi_med + ndre_med +
gli_med + height_med,
data = train_80_ponderado,
weights = peso,
size = 5,
decay = 0.001,
maxit = 1000,
MaxNWts = 10000,
trace = FALSE
)
modelo_80_20## a 5-5-6 network with 66 weights
## inputs: ndvi_med evi_med ndre_med gli_med height_med
## output(s): severity
## options were - softmax modelling decay=0.001## Verificación de la correcta preparación de los datos
# Verificar las categorías
levels(train_80_ponderado$severity)## [1] "1" "2" "3" "4" "5" "6"# Verificar que el peso exista
names(train_80_ponderado)## [1] "severity" "ndvi_med" "evi_med" "ndre_med" "gli_med"
## [6] "height_med" "peso"# Revisar pesos
summary(train_80_ponderado$peso)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7568 0.7778 0.7778 1.0000 1.0000 2.8000# Los pesos deben ser positivos
any(train_80_ponderado$peso <= 0)## [1] FALSE# Revisar datos faltantes
colSums(
is.na(
train_80_ponderado[, c(
"severity",
"ndvi_med",
"evi_med",
"ndre_med",
"gli_med",
"height_med",
"peso"
)]
)
)## severity ndvi_med evi_med ndre_med gli_med height_med peso
## 0 0 0 0 0 0 0library(nnet)
set.seed(123)
modelo_80_20 <- nnet::nnet(
severity ~ ndvi_med + evi_med + ndre_med +
gli_med + height_med,
data = train_80_ponderado,
weights = peso,
size = 5,
decay = 0.001,
maxit = 1000,
MaxNWts = 10000,
trace = FALSE
)
modelo_80_20## a 5-5-6 network with 66 weights
## inputs: ndvi_med evi_med ndre_med gli_med height_med
## output(s): severity
## options were - softmax modelling decay=0.001## Predecir categorias del conjunto de prueba
pred_test_80_20 <- predict(
modelo_80_20,
newdata = test_20_std,
type = "class"
)
head(pred_test_80_20)## [1] "1" "1" "1" "1" "1" "1"## Generación de matriz de confusión
matriz_80_20 <- table(
Real = test_20_std$severity,
Predicho = pred_test_80_20
)
matriz_80_20## Predicho
## Real 1 2 3 4 5 6
## 1 3 1 2 0 0 0
## 2 3 4 3 0 0 0
## 3 2 2 3 0 0 0
## 4 1 0 4 1 2 1
## 5 0 0 0 3 6 0
## 6 0 0 0 0 1 2## Calculo de exactitud global
exactitud_80_20 <- sum(diag(matriz_80_20)) / sum(matriz_80_20)
exactitud_80_20## [1] 0.4318182## Probabilidades asignadas por categoria
probabilidades_80_20 <- predict(
modelo_80_20,
newdata = test_20_std,
type = "raw"
)
head(probabilidades_80_20)## 1 2 3 4 5 6
## 1 0.7094324 0.1818660 0.1084269 0.0002599040 1.477729e-05 1.075436e-18
## 2 0.4425042 0.3924892 0.1647325 0.0002668539 7.207227e-06 2.823310e-17
## 3 0.6039694 0.2619169 0.1338395 0.0002635247 1.074671e-05 4.433639e-18
## 4 0.4376929 0.3981986 0.1638573 0.0002449973 6.265578e-06 2.964457e-17
## 5 0.6228725 0.2471851 0.1296672 0.0002639023 1.138246e-05 3.502810e-18
## 6 0.6994001 0.1892828 0.1110421 0.0002606788 1.431217e-05 1.248432e-18El modelo entrenado con el esquema 80 %–20 % alcanzó una exactitud global de 0,432, lo que indica que clasificó correctamente aproximadamente el 43,2 % de las observaciones del conjunto de prueba. Este valor es reflejo de una capacidad predictiva limitada. El análisis de las probabilidades mostró que, en las primeras observaciones evaluadas, la red concentró la mayor probabilidad en la categoría 1 y asignó probabilidades prácticamente nulas a las categorías 4, 5 y 6. Esto puede indicar una alta confianza en la clasificación de esos casos como severidad baja, pero también sugiere la necesidad de revisar la matriz de confusión y las métricas por clase para determinar si el modelo está favoreciendo determinadas categorías.
## Predicciones por categoria
table(pred_test_80_20)## pred_test_80_20
## 1 2 3 4 5 6
## 9 7 12 4 9 3resultados_80_20 <- evaluar_modelo(
real = test_20_std$severity,
predicho = pred_test_80_20
)
resultados_80_20$matriz_confusion## Predicho
## Real 1 2 3 4 5 6
## 1 3 1 2 0 0 0
## 2 3 4 3 0 0 0
## 3 2 2 3 0 0 0
## 4 1 0 4 1 2 1
## 5 0 0 0 3 6 0
## 6 0 0 0 0 1 2resultados_80_20$sensibilidad## 1 2 3 4 5 6
## 0.5000000 0.4000000 0.4285714 0.1111111 0.6666667 0.6666667resultados_80_20$precision## 1 2 3 4 5 6
## 0.3333333 0.5714286 0.2500000 0.2500000 0.6666667 0.6666667resultados_80_20$f1_clase## 1 2 3 4 5 6
## 0.4000000 0.4705882 0.3157895 0.1538462 0.6666667 0.6666667resultados_80_20$exactitud_balanceada## [1] 0.4621693resultados_80_20$f1_macro## [1] 0.4455929La matriz de confusión mostró que el desempeño del modelo fue desigual entre las categorías de severidad. Las clases 5 y 6 alcanzaron las mayores sensibilidades y valores F1, ambos cercanos a 0,667 teniendo la mejor tasa de detección, mientras que la categoría 4 presentó el menor desempeño, con una sensibilidad de 0,111 y un F1 de 0,154. La exactitud balanceada fue 0,462 y el F1 macro 0,446, lo que indica una capacidad de clasificación moderada y diferencias importantes entre clases. Los errores se concentraron principalmente entre categorías adyacentes, especialmente en los niveles intermedios, lo que sugiere una superposición considerable en sus características espectrales y morfológicas. La ponderación permitió una detección relativamente favorable de la categoría minoritaria 6, aunque el modelo todavía requiere ajuste de arquitectura e hiperparámetros para mejorar su capacidad de generalización.
Ajuste de arquitectura e hiperparámetros
library(dplyr)
library(nnet)
# Asegurar que severity sea factor
df <- df %>%
mutate(
severity = factor(
severity,
levels = c("1", "2", "3", "4", "5", "6")
)
)## Crear función de partición estratificada
particion_estratificada <- function(datos, proporcion, semilla = 123) {
set.seed(semilla)
indices <- unlist(
lapply(
split(seq_len(nrow(datos)), datos$severity),
function(indices_clase) {
cantidad <- floor(
proporcion * length(indices_clase)
)
sample(
indices_clase,
size = cantidad,
replace = FALSE
)
}
)
)
list(
seleccionados = datos[indices, , drop = FALSE],
restantes = datos[-indices, , drop = FALSE]
)
}## Partición 70 %–15 %–15 %:
# Primera partición: 70 % entrenamiento
primera_particion <- particion_estratificada(
datos = df,
proporcion = 0.70,
semilla = 123
)
train_70 <- primera_particion$seleccionados
restante_30 <- primera_particion$restantes
# Segunda partición: dividir el 30 % restante en dos partes
segunda_particion <- particion_estratificada(
datos = restante_30,
proporcion = 0.50,
semilla = 456
)
validacion_15 <- segunda_particion$seleccionados
test_15 <- segunda_particion$restantes## Definir los predictores:
predictores <- c(
"ndvi_med",
"evi_med",
"ndre_med",
"gli_med",
"height_med"
)medias_train_70 <- sapply(
train_70[, predictores],
mean,
na.rm = TRUE
)
desviaciones_train_70 <- sapply(
train_70[, predictores],
sd,
na.rm = TRUE
)train_70_std <- train_70
validacion_15_std <- validacion_15
test_15_std <- test_15
train_70_std[, predictores] <- scale(
train_70[, predictores],
center = medias_train_70,
scale = desviaciones_train_70
)
validacion_15_std[, predictores] <- scale(
validacion_15[, predictores],
center = medias_train_70,
scale = desviaciones_train_70
)
test_15_std[, predictores] <- scale(
test_15[, predictores],
center = medias_train_70,
scale = desviaciones_train_70
)dim(train_70_std)## [1] 145 6dim(validacion_15_std)## [1] 32 6dim(test_15_std)## [1] 35 6table(train_70_std$severity)##
## 1 2 3 4 5 6
## 18 32 24 31 31 9table(validacion_15_std$severity)##
## 1 2 3 4 5 6
## 4 7 5 7 7 2table(test_15_std$severity)##
## 1 2 3 4 5 6
## 5 8 6 7 7 2La partición estratificada 70 %–15 %–15 % generó 145 observaciones para entrenamiento, 32 para validación y 35 para prueba, correspondientes aproximadamente al 68,4 %, 15,1 % y 16,5 % de la base total, respectivamente.La estratificación permitió conservar la distribución relativa de los niveles de severidad en los tres subconjuntos. Las categorías 2, 4 y 5 mantuvieron la mayor representación, mientras que la severidad 6 quedó con 9 observaciones en entrenamiento y únicamente 2 en validación y prueba. En consecuencia, aunque la partición preserva el desbalance original y garantiza la presencia de todas las categorías, las métricas de la clase 6 pueden presentar alta variabilidad debido a su reducido número de casos.
library(dplyr)
library(nnet)
pesos_train_70 <- train_70_std %>%
count(severity, name = "frecuencia") %>%
mutate(
peso = nrow(train_70_std) /
(nlevels(train_70_std$severity) * frecuencia)
)
pesos_train_70train_70_ponderado <- train_70_std %>%
left_join(
pesos_train_70 %>%
select(severity, peso),
by = "severity"
)
train_70_ponderado %>%
distinct(severity, peso) %>%
arrange(severity)colSums(
is.na(
train_70_ponderado[, c(
"severity",
"ndvi_med",
"evi_med",
"ndre_med",
"gli_med",
"height_med",
"peso"
)]
)
)## severity ndvi_med evi_med ndre_med gli_med height_med peso
## 0 0 0 0 0 0 0## Crear función de evaluación
evaluar_modelo <- function(real, predicho) {
niveles <- c("1", "2", "3", "4", "5", "6")
real <- factor(real, levels = niveles)
predicho <- factor(predicho, levels = niveles)
matriz <- table(
Real = real,
Predicho = predicho
)
verdaderos_positivos <- diag(matriz)
sensibilidad <- verdaderos_positivos / rowSums(matriz)
precision <- verdaderos_positivos / colSums(matriz)
f1_clase <- 2 * precision * sensibilidad /
(precision + sensibilidad)
# Reemplazar divisiones indefinidas por cero
sensibilidad[!is.finite(sensibilidad)] <- 0
precision[!is.finite(precision)] <- 0
f1_clase[!is.finite(f1_clase)] <- 0
exactitud <- sum(verdaderos_positivos) / sum(matriz)
exactitud_balanceada <- mean(sensibilidad)
f1_macro <- mean(f1_clase)
list(
matriz_confusion = matriz,
exactitud = exactitud,
sensibilidad = sensibilidad,
precision = precision,
f1_clase = f1_clase,
exactitud_balanceada = exactitud_balanceada,
f1_macro = f1_macro
)
}Definir arquitecturas e hiperparámetros
## Definir las arquitecturas e hiperparámetros
#Probaremos cuatro tamaños de capa oculta, cuatro valores de regularización y dos números máximos de iteraciones:
rejilla_hiperparametros <- expand.grid(
size = c(3, 5, 7, 10),
decay = c(0, 0.001, 0.01, 0.1),
maxit = c(500, 1000)
)
rejilla_hiperparametrosnrow(rejilla_hiperparametros)## [1] 32## Entrenar las 32 configuraciones
resultados_ajuste <- data.frame()
for (i in seq_len(nrow(rejilla_hiperparametros))) {
size_actual <- rejilla_hiperparametros$size[i]
decay_actual <- rejilla_hiperparametros$decay[i]
maxit_actual <- rejilla_hiperparametros$maxit[i]
# Misma semilla para comparar las configuraciones
set.seed(123)
modelo_temporal <- tryCatch(
nnet::nnet(
severity ~ ndvi_med + evi_med + ndre_med +
gli_med + height_med,
data = train_70_ponderado,
weights = peso,
size = size_actual,
decay = decay_actual,
maxit = maxit_actual,
MaxNWts = 10000,
trace = FALSE
),
error = function(e) {
message(
"Error con size = ", size_actual,
", decay = ", decay_actual,
", maxit = ", maxit_actual
)
return(NULL)
}
)
if (!is.null(modelo_temporal)) {
prediccion_validacion <- predict(
modelo_temporal,
newdata = validacion_15_std,
type = "class"
)
metricas <- evaluar_modelo(
real = validacion_15_std$severity,
predicho = prediccion_validacion
)
resultados_ajuste <- rbind(
resultados_ajuste,
data.frame(
size = size_actual,
decay = decay_actual,
maxit = maxit_actual,
exactitud = metricas$exactitud,
exactitud_balanceada =
metricas$exactitud_balanceada,
f1_macro = metricas$f1_macro,
convergencia = modelo_temporal$convergence,
valor_final = modelo_temporal$value
)
)
}
}resultados_ajuste## Ordenar y seleccionar mejor resultado
resultados_ordenados <- resultados_ajuste %>%
arrange(
desc(f1_macro),
desc(exactitud_balanceada),
desc(exactitud)
)
head(resultados_ordenados, 10)resultados_ordenados %>%
mutate(
across(
c(
exactitud,
exactitud_balanceada,
f1_macro,
valor_final
),
~ round(.x, 4)
)
) %>%
head(10)mejor_configuracion <- resultados_ordenados %>%
slice(1)
mejor_configuracion## Revisar convergencia
table(resultados_ajuste$convergencia)##
## 0 1
## 27 5El resultado indica que, de las 32 configuraciones evaluadas 27 modelos tuvieron convergencia = 0, es decir, el algoritmo terminó satisfactoriamente y 5 modelos tuvieron convergencia = 1, indicando que alcanzaron el límite máximo de iteraciones antes de cumplir completamente el criterio de convergencia.
La mayoría de las arquitecturas se ajustó adecuadamente, pero cinco configuraciones necesitarían más iteraciones o una combinación distinta de hiperparámetros. Esto no significa necesariamente que los cinco modelos sean inútiles, pero sí que sus pesos podrían no haber alcanzado una solución estable. Por esa razón, es más seguro seleccionar la mejor arquitectura únicamente entre los modelos que presentaron convergencia = 0.
## Selección unicamente en modelos convergentes
resultados_convergentes <- resultados_ajuste %>%
filter(convergencia == 0) %>%
arrange(
desc(f1_macro),
desc(exactitud_balanceada),
desc(exactitud)
)
head(resultados_convergentes, 10)mejor_configuracion <- resultados_convergentes %>%
slice(1)
mejor_configuracionLa mejor configuración entre los modelos que alcanzaron convergencia estuvo compuesta por tres neuronas en la capa oculta, un parámetro de regularización decay de 0,01 y un máximo de 500 iteraciones. Esta red obtuvo una exactitud de 0,438, una exactitud balanceada de 0,449 y un F1 macro de 0,430 en el conjunto de validación. Los resultados indican que una arquitectura relativamente sencilla, acompañada de una regularización moderada, presentó el mejor desempeño entre las configuraciones evaluadas. El hecho de que una red con menor número de neuronas superara a arquitecturas más complejas sugiere que aumentar la capacidad del modelo no produjo una mejora en la generalización y pudo incrementar el riesgo de sobreajuste. Aunque esta configuración fue la mejor dentro de la rejilla considerada, las métricas evidencian un desempeño moderado, por lo que aún existen dificultades para discriminar de manera equilibrada las seis categorías de severidad.
Ventaja del esquema de tres particiones frente al de dos
El esquema de tres particiones presenta una ventaja frente al esquema de entrenamiento y prueba porque separa el aprendizaje, la selección del modelo y la evaluación final. El conjunto de entrenamiento se utiliza para ajustar la red, el conjunto de validación para comparar arquitecturas e hiperparámetros y el conjunto de prueba se reserva para evaluar el modelo seleccionado. Esta separación evita utilizar repetidamente la prueba para tomar decisiones, lo que podría producir una estimación demasiado optimista del desempeño.
La validación también permite detectar el sobreajuste. Si una red presenta buenos resultados en entrenamiento, pero un desempeño inferior en validación, significa que posiblemente aprendió características demasiado particulares del conjunto de entrenamiento. En este análisis, la validación permitió seleccionar una arquitectura sencilla, con tres neuronas ocultas, decay = 0,01 y 500 iteraciones. Esto muestra que una red más compleja no necesariamente ofrece una mejor capacidad de generalización. Aunque el esquema de tres particiones deja menos datos disponibles en cada subconjunto, es más adecuado cuando se deben seleccionar hiperparámetros y obtener una evaluación final más imparcial.
4. Dado que el tamaño muestral es mayor a 200 pero no muy grande, ¿considera que sería apropiado usar validación cruzada k-fold como alternativa o complemento? Argumente su respuesta y, si decide implementarla, reporte los resultados comparativos.
Debido a que la base contiene 212 observaciones, se considera apropiado complementar la partición convencional con validación cruzada k-fold. Aunque el tamaño muestral supera las 200 observaciones, sigue siendo limitado para dividir los datos en entrenamiento, validación y prueba sin reducir considerablemente la representación de las categorías minoritarias. En particular, la severidad 6 quedó representada por solo dos observaciones en validación y dos en prueba, por lo que sus métricas pueden variar ampliamente ante una sola clasificación correcta o incorrecta.
Se propone utilizar validación cruzada estratificada de cinco pliegues. En este procedimiento, los datos se dividen en cinco grupos y el modelo se ajusta cinco veces, utilizando en cada ocasión cuatro pliegues para entrenamiento y uno para validación. La estratificación permite conservar aproximadamente la proporción de las categorías de severidad en cada división. Como resultado, todas las observaciones participan en la validación y las métricas finales corresponden al promedio de los cinco ajustes, lo que proporciona una estimación más estable del desempeño.
La validación cruzada sería especialmente útil para seleccionar el número de neuronas ocultas, el valor de regularización y el número de iteraciones, utilizando como criterios el F1 macro y la exactitud balanceada. No obstante, se recomienda conservar un conjunto de prueba independiente para la evaluación final. Por tanto, la validación cruzada debe aplicarse dentro del conjunto de desarrollo y no sobre los datos de prueba. Esta estrategia permite aprovechar mejor la información disponible, reducir la dependencia de una sola partición y disminuir el riesgo de seleccionar hiperparámetros por una división particularmente favorable.
## Preparar datos
library(dplyr)
library(nnet)
niveles_severidad <- c("1", "2", "3", "4", "5", "6")
# Unir entrenamiento y validación
datos_desarrollo <- bind_rows(
train_70,
validacion_15
)
# Asegurar que severity sea factor
datos_desarrollo$severity <- factor(
datos_desarrollo$severity,
levels = niveles_severidad
)
test_15$severity <- factor(
test_15$severity,
levels = niveles_severidad
)
dim(datos_desarrollo)## [1] 177 6dim(test_15)## [1] 35 6table(datos_desarrollo$severity)##
## 1 2 3 4 5 6
## 22 39 29 38 38 11table(test_15$severity)##
## 1 2 3 4 5 6
## 5 8 6 7 7 2## Definir variables predictoras
predictores <- c(
"ndvi_med",
"evi_med",
"ndre_med",
"gli_med",
"height_med"
)## Crear 5 pliegues
crear_folds_estratificados <- function(datos, k = 5, semilla = 123) {
set.seed(semilla)
folds <- integer(nrow(datos))
for (clase in levels(datos$severity)) {
indices_clase <- which(datos$severity == clase)
# Mezclar aleatoriamente los registros de la clase
indices_clase <- sample(indices_clase)
# Distribuirlos entre los k pliegues
asignacion <- rep(
1:k,
length.out = length(indices_clase)
)
folds[indices_clase] <- asignacion
}
folds
}fold_id <- crear_folds_estratificados(
datos = datos_desarrollo,
k = 5,
semilla = 123
)
table(fold_id)## fold_id
## 1 2 3 4 5
## 38 37 36 34 32## Verificación de distribución de la severidad
tabla_folds <- table(
Fold = fold_id,
Severidad = datos_desarrollo$severity
)
tabla_folds## Severidad
## Fold 1 2 3 4 5 6
## 1 5 8 6 8 8 3
## 2 5 8 6 8 8 2
## 3 4 8 6 8 8 2
## 4 4 8 6 7 7 2
## 5 4 7 5 7 7 2## Función de calculo de métricas
evaluar_modelo_cv <- function(real, predicho) {
real <- factor(
real,
levels = niveles_severidad
)
predicho <- factor(
predicho,
levels = niveles_severidad
)
matriz <- table(
Real = real,
Predicho = predicho
)
verdaderos_positivos <- diag(matriz)
sensibilidad <- verdaderos_positivos / rowSums(matriz)
precision <- verdaderos_positivos / colSums(matriz)
f1_clase <- 2 * precision * sensibilidad /
(precision + sensibilidad)
# Reemplazar resultados indefinidos por cero
sensibilidad[!is.finite(sensibilidad)] <- 0
precision[!is.finite(precision)] <- 0
f1_clase[!is.finite(f1_clase)] <- 0
exactitud <- sum(verdaderos_positivos) / sum(matriz)
exactitud_balanceada <- mean(sensibilidad)
f1_macro <- mean(f1_clase)
list(
matriz = matriz,
exactitud = exactitud,
exactitud_balanceada = exactitud_balanceada,
f1_macro = f1_macro,
sensibilidad = sensibilidad,
precision = precision,
f1_clase = f1_clase
)
}## Definir rejilla de hiperparámetros
rejilla_cv <- expand.grid(
size = c(3, 5, 7, 10),
decay = c(0, 0.001, 0.01, 0.1),
maxit = c(500, 1000)
)
rejilla_cvnrow(rejilla_cv)## [1] 32evaluar_modelo_cv <- function(real, predicho) {
niveles <- c("1", "2", "3", "4", "5", "6")
real <- factor(
real,
levels = niveles
)
predicho <- factor(
predicho,
levels = niveles
)
matriz <- table(
Real = real,
Predicho = predicho
)
verdaderos_positivos <- diag(matriz)
sensibilidad <- verdaderos_positivos / rowSums(matriz)
precision <- verdaderos_positivos / colSums(matriz)
f1_clase <- 2 * precision * sensibilidad /
(precision + sensibilidad)
# Corregir divisiones indefinidas
sensibilidad[!is.finite(sensibilidad)] <- 0
precision[!is.finite(precision)] <- 0
f1_clase[!is.finite(f1_clase)] <- 0
exactitud <- sum(verdaderos_positivos) / sum(matriz)
exactitud_balanceada <- mean(sensibilidad)
f1_macro <- mean(f1_clase)
list(
matriz = matriz,
exactitud = exactitud,
exactitud_balanceada = exactitud_balanceada,
f1_macro = f1_macro,
sensibilidad = sensibilidad,
precision = precision,
f1_clase = f1_clase
)
}exists("evaluar_modelo_cv")## [1] TRUE##Ejecutar la validación cruzada. Este es el bloque principal que entrenará: 32 configuraciones×5 pliegues=160 modelos
resultados_cv <- data.frame()
for (i in seq_len(nrow(rejilla_cv))) {
size_actual <- rejilla_cv$size[i]
decay_actual <- rejilla_cv$decay[i]
maxit_actual <- rejilla_cv$maxit[i]
for (fold_actual in 1:5) {
# Datos internos de entrenamiento y validación
train_fold <- datos_desarrollo[
fold_id != fold_actual,
,
drop = FALSE
]
valid_fold <- datos_desarrollo[
fold_id == fold_actual,
,
drop = FALSE
]
# Garantizar los mismos niveles
train_fold$severity <- factor(
train_fold$severity,
levels = niveles_severidad
)
valid_fold$severity <- factor(
valid_fold$severity,
levels = niveles_severidad
)
# -----------------------------------------
# Estandarización dentro de cada pliegue
# -----------------------------------------
medias_fold <- sapply(
train_fold[, predictores, drop = FALSE],
mean,
na.rm = TRUE
)
desviaciones_fold <- sapply(
train_fold[, predictores, drop = FALSE],
sd,
na.rm = TRUE
)
# Evitar división por cero
desviaciones_fold[desviaciones_fold == 0] <- 1
train_fold_std <- train_fold
valid_fold_std <- valid_fold
train_fold_std[, predictores] <- scale(
train_fold[, predictores, drop = FALSE],
center = medias_fold,
scale = desviaciones_fold
)
valid_fold_std[, predictores] <- scale(
valid_fold[, predictores, drop = FALSE],
center = medias_fold,
scale = desviaciones_fold
)
# -----------------------------------------
# Ponderación dentro de cada pliegue
# -----------------------------------------
pesos_fold <- train_fold_std %>%
count(severity, name = "frecuencia") %>%
mutate(
peso = nrow(train_fold_std) /
(nlevels(train_fold_std$severity) * frecuencia)
)
train_fold_ponderado <- train_fold_std %>%
left_join(
pesos_fold %>% select(severity, peso),
by = "severity"
)
# -----------------------------------------
# Entrenamiento del modelo
# -----------------------------------------
set.seed(1000 + i * 10 + fold_actual)
modelo_fold <- tryCatch(
nnet::nnet(
severity ~ ndvi_med + evi_med + ndre_med +
gli_med + height_med,
data = train_fold_ponderado,
weights = peso,
size = size_actual,
decay = decay_actual,
maxit = maxit_actual,
MaxNWts = 10000,
trace = FALSE
),
error = function(e) NULL
)
# -----------------------------------------
# Evaluación
# -----------------------------------------
if (!is.null(modelo_fold)) {
pred_fold <- predict(
modelo_fold,
newdata = valid_fold_std,
type = "class"
)
metricas_fold <- evaluar_modelo_cv(
real = valid_fold_std$severity,
predicho = pred_fold
)
resultados_cv <- rbind(
resultados_cv,
data.frame(
size = size_actual,
decay = decay_actual,
maxit = maxit_actual,
fold = fold_actual,
exactitud = metricas_fold$exactitud,
exactitud_balanceada =
metricas_fold$exactitud_balanceada,
f1_macro = metricas_fold$f1_macro,
convergencia = modelo_fold$convergence
)
)
}
}
}## Revisar convergencia
table(resultados_cv$convergencia)##
## 0 1
## 132 28En la validación cruzada se ajustaron 160 modelos, correspondientes a 32 combinaciones de hiperparámetros evaluadas en cinco pliegues. De estos, 132 modelos, equivalentes al 82,5 %, alcanzaron una convergencia satisfactoria, mientras que 28 modelos, correspondientes al 17,5 %, llegaron al número máximo de iteraciones sin cumplir completamente el criterio de convergencia. Por esta razón, la selección de la arquitectura se restringió a las configuraciones que convergieron en los cinco pliegues, con el fin de comparar modelos numéricamente estables y evitar seleccionar una combinación cuyo desempeño pudiera depender de un ajuste incompleto.
resultados_cv %>%
filter(convergencia != 0) %>%
count(size, decay, maxit, name = "pliegues_sin_convergencia") %>%
arrange(desc(pliegues_sin_convergencia))El resultado muestra que los problemas de convergencia se concentraron principalmente en los modelos sin regularización.El aumento de 500 a 1000 iteraciones no resolvió completamente el problema, lo que indica que la falta de convergencia estuvo asociada principalmente con la complejidad de la red y la ausencia de penalización sobre los pesos.
## Obtener el promedio de los cinco pliegues
resumen_cv <- resultados_cv %>%
group_by(size, decay, maxit) %>%
summarise(
folds_evaluados = n(),
exactitud_media = mean(exactitud),
exactitud_sd = sd(exactitud),
exactitud_balanceada_media =
mean(exactitud_balanceada),
exactitud_balanceada_sd =
sd(exactitud_balanceada),
f1_macro_medio = mean(f1_macro),
f1_macro_sd = sd(f1_macro),
modelos_convergentes =
sum(convergencia == 0),
.groups = "drop"
)resumen_cv_ordenado <- resumen_cv %>%
arrange(
desc(f1_macro_medio),
desc(exactitud_balanceada_media),
desc(exactitud_media)
)
head(resumen_cv_ordenado, 10)resumen_cv_ordenado %>%
mutate(
across(
c(
exactitud_media,
exactitud_sd,
exactitud_balanceada_media,
exactitud_balanceada_sd,
f1_macro_medio,
f1_macro_sd
),
~ round(.x, 4)
)
) %>%
head(10)## Seleccionar la mejor configuración
mejor_configuracion_cv <- resumen_cv %>%
filter(
folds_evaluados == 5,
modelos_convergentes == 5
) %>%
arrange(
desc(f1_macro_medio),
desc(exactitud_balanceada_media),
desc(exactitud_media)
) %>%
slice(1)
mejor_configuracion_cvLa configuración seleccionada estuvo compuesta por diez neuronas ocultas, un valor de regularización decay de 0,1 y 500 iteraciones. Los cinco modelos ajustados durante la validación cruzada alcanzaron convergencia. La exactitud media fue 0,527, la exactitud balanceada media 0,569 y el F1 macro medio 0,533. Estos resultados indican un desempeño moderado y relativamente equilibrado entre las categorías. Sin embargo, las desviaciones estándar fueron 0,104, 0,117 y 0,122, respectivamente, lo que evidencia una variación apreciable entre los pliegues. Por tanto, aunque esta configuración presentó el mejor desempeño promedio, sus resultados deben interpretarse con cautela, ya que la capacidad predictiva depende parcialmente de la composición de cada subconjunto de validación.
## Grafico de resultados
library(ggplot2)
ggplot(
resumen_cv,
aes(
x = factor(size),
y = f1_macro_medio,
group = factor(decay),
linetype = factor(decay),
shape = factor(decay)
)
) +
geom_line() +
geom_point(size = 2.5) +
geom_errorbar(
aes(
ymin = f1_macro_medio - f1_macro_sd,
ymax = f1_macro_medio + f1_macro_sd
),
width = 0.1
) +
facet_wrap(~ maxit) +
labs(
title = "Validación cruzada estratificada de 5 pliegues",
subtitle = "Promedio y desviación estándar del F1 macro",
x = "Número de neuronas ocultas",
y = "F1 macro promedio",
linetype = "Decay",
shape = "Decay"
) +
theme_minimal()
Figura 3. Validación cruzada estratificada de 5 pliegues.
El gráfico mostró que el desempeño del perceptrón multicapa dependió de la interacción entre el número de neuronas ocultas y el grado de regularización (Figura 3). Las configuraciones sin regularización presentaron, en general, valores menores de F1 macro y mayores problemas de convergencia. Por el contrario, los valores positivos de decay, especialmente 0,1, produjeron mejores resultados, lo que sugiere que la regularización ayudó a controlar el sobreajuste.
## Entrenar el modelo definitivo
mejor_size_cv <- mejor_configuracion_cv$size[1]
mejor_decay_cv <- mejor_configuracion_cv$decay[1]
mejor_maxit_cv <- mejor_configuracion_cv$maxit[1]
mejor_size_cv## [1] 10mejor_decay_cv## [1] 0.1mejor_maxit_cv## [1] 500## Calcula la estandarización con todo el conjunto de desarrollo:
medias_desarrollo <- sapply(
datos_desarrollo[, predictores, drop = FALSE],
mean,
na.rm = TRUE
)
desviaciones_desarrollo <- sapply(
datos_desarrollo[, predictores, drop = FALSE],
sd,
na.rm = TRUE
)
desviaciones_desarrollo[
desviaciones_desarrollo == 0
] <- 1## Estandariza desarrollo y prueba:
desarrollo_std <- datos_desarrollo
test_final_std <- test_15
desarrollo_std[, predictores] <- scale(
datos_desarrollo[, predictores, drop = FALSE],
center = medias_desarrollo,
scale = desviaciones_desarrollo
)
test_final_std[, predictores] <- scale(
test_15[, predictores, drop = FALSE],
center = medias_desarrollo,
scale = desviaciones_desarrollo
)## Calcula los pesos con todo el conjunto de desarrollo:
pesos_desarrollo <- desarrollo_std %>%
count(severity, name = "frecuencia") %>%
mutate(
peso = nrow(desarrollo_std) /
(nlevels(desarrollo_std$severity) * frecuencia)
)
desarrollo_ponderado <- desarrollo_std %>%
left_join(
pesos_desarrollo %>% select(severity, peso),
by = "severity"
)## Entrena el modelo:
set.seed(123)
modelo_final_cv <- nnet::nnet(
severity ~ ndvi_med + evi_med + ndre_med +
gli_med + height_med,
data = desarrollo_ponderado,
weights = peso,
size = mejor_size_cv,
decay = mejor_decay_cv,
maxit = mejor_maxit_cv,
MaxNWts = 10000,
trace = FALSE
)
modelo_final_cv## a 5-10-6 network with 126 weights
## inputs: ndvi_med evi_med ndre_med gli_med height_med
## output(s): severity
## options were - softmax modelling decay=0.1## Evaluar una sola vez en prueba
prediccion_test_cv <- predict(
modelo_final_cv,
newdata = test_final_std,
type = "class"
)
resultado_test_cv <- evaluar_modelo_cv(
real = test_final_std$severity,
predicho = prediccion_test_cv
)## Resultados generales:
resultado_test_cv$matriz## Predicho
## Real 1 2 3 4 5 6
## 1 3 1 1 0 0 0
## 2 2 5 1 0 0 0
## 3 1 3 1 1 0 0
## 4 0 1 0 4 1 1
## 5 0 0 1 2 1 3
## 6 0 0 0 0 0 2resultado_test_cv$exactitud## [1] 0.4571429resultado_test_cv$exactitud_balanceada## [1] 0.5176587resultado_test_cv$f1_macro## [1] 0.4324435## Resultados por categoría:
metricas_clase_test_cv <- data.frame(
severidad = niveles_severidad,
sensibilidad = resultado_test_cv$sensibilidad,
precision = resultado_test_cv$precision,
f1 = resultado_test_cv$f1_clase
)
metricas_clase_test_cvEl modelo seleccionado mediante validación cruzada alcanzó una exactitud global de 0,457 en el conjunto de prueba, equivalente a 16 clasificaciones correctas de 35 observaciones. La exactitud balanceada fue de 0,518 y el F1 macro de 0,432, lo que evidencia un desempeño moderado y desigual entre las categorías.
Las severidades 1, 2 y 4 presentaron los resultados más consistentes, mientras que las categorías 3 y 5 mostraron las sensibilidades y valores F1 más bajos. La severidad 6 alcanzó una sensibilidad de 1,00 al clasificar correctamente sus dos casos reales; sin embargo, su precisión fue de 0,333 debido a la presencia de falsos positivos, por lo que este resultado debe interpretarse con cuidado.
En general, los errores se concentraron principalmente entre categorías próximas, especialmente entre las severidades intermedias y altas, lo que sugiere una superposición de sus características espectrales y morfológicas. El desempeño en prueba fue inferior al promedio obtenido mediante validación cruzada, aunque permaneció dentro de la variabilidad observada entre los pliegues.
- Comparación con la validación cruzada
Durante la validación cruzada, la mejor configuración había obtenido aproximadamente:
exactitud media: 52,7 %; exactitud balanceada media: 56,9 %; F1 macro medio: 53,3 %.
En la prueba independiente estos valores disminuyeron a:
exactitud: 45,7 %; exactitud balanceada: 51,8 %; F1 macro: 43,2 %.
Esta reducción indica que el desempeño estimado durante la validación cruzada fue algo más favorable que el obtenido en datos completamente independientes. No necesariamente representa un sobreajuste grave, porque las métricas de la validación cruzada presentaron desviaciones estándar relativamente amplias. El resultado de prueba se encuentra dentro de la variabilidad esperable, especialmente considerando que solo se evaluaron 35 observaciones.
2.3. Arquitectura y entrenamiento
- Entrene al menos cuatro perceptrones multicapa variando sistemáticamente:
• Número de capas ocultas: 1, 2 y 3 capas. • Número de neuronas por capa: al menos dos configuraciones por número de capas. • Tasa de aprendizaje: pruebe al menos tres valores (por ejemplo, 0.001, 0.01, 0.1).
Para cada configuración reporte la métrica de pérdida en entrenamiento y validación por época. ¿En qué configuraciones observa indicios de sobreajuste o subajuste? Identifíquelos gráficamente mediante las curvas de aprendizaje.
library(dplyr)
library(tidyr)
library(ggplot2)exists("train_70_std")## [1] TRUEexists("validacion_15_std")## [1] TRUEexists("test_15_std")## [1] TRUEpredictores <- c(
"ndvi_med",
"evi_med",
"ndre_med",
"gli_med",
"height_med"
)
x_train <- as.matrix(
train_70_std[, predictores]
)
x_val <- as.matrix(
validacion_15_std[, predictores]
)
y_train <- as.integer(
factor(train_70_std$severity, levels = 1:6)
)
y_val <- as.integer(
factor(validacion_15_std$severity, levels = 1:6)
)dim(x_train)## [1] 145 5dim(x_val)## [1] 32 5table(y_train)## y_train
## 1 2 3 4 5 6
## 18 32 24 31 31 9table(y_val)## y_val
## 1 2 3 4 5 6
## 4 7 5 7 7 2## Crear funciones básicas de la red neuronal
relu <- function(x) {
resultado <- pmax(x, 0)
dim(resultado) <- dim(x)
resultado
}
relu_derivada <- function(x) {
resultado <- 1 * (x > 0)
dim(resultado) <- dim(x)
resultado
}
softmax <- function(z) {
z_estable <- z - apply(z, 1, max)
exp_z <- exp(z_estable)
exp_z / rowSums(exp_z)
}
one_hot <- function(y, numero_clases = 6) {
matriz <- matrix(
0,
nrow = length(y),
ncol = numero_clases
)
matriz[cbind(seq_along(y), y)] <- 1
matriz
}
entropia_cruzada <- function(
probabilidades,
y_one_hot
) {
epsilon <- 1e-12
-mean(
rowSums(
y_one_hot *
log(probabilidades + epsilon)
)
)
}## Define la función que crea la arquitectura:
inicializar_red <- function(
numero_entradas,
capas_ocultas,
numero_salidas = 6,
semilla = 123
) {
set.seed(semilla)
tamanos <- c(
numero_entradas,
capas_ocultas,
numero_salidas
)
pesos <- list()
sesgos <- list()
for (i in seq_len(length(tamanos) - 1)) {
limite <- sqrt(
2 / tamanos[i]
)
pesos[[i]] <- matrix(
rnorm(
tamanos[i] * tamanos[i + 1],
mean = 0,
sd = limite
),
nrow = tamanos[i],
ncol = tamanos[i + 1]
)
sesgos[[i]] <- matrix(
0,
nrow = 1,
ncol = tamanos[i + 1]
)
}
list(
pesos = pesos,
sesgos = sesgos
)
}## Ahora agrega la función de entrenamiento:
entrenar_mlp <- function(
x_train,
y_train,
x_val,
y_val,
capas_ocultas,
tasa_aprendizaje,
epocas = 150,
semilla = 123
) {
# Asegurar que los predictores sean matrices numéricas
x_train <- as.matrix(x_train)
x_val <- as.matrix(x_val)
storage.mode(x_train) <- "double"
storage.mode(x_val) <- "double"
# Validaciones iniciales
if (ncol(x_train) == 0) {
stop("x_train no contiene columnas.")
}
if (length(capas_ocultas) == 0) {
stop("Debe especificarse al menos una capa oculta.")
}
if (any(capas_ocultas <= 0)) {
stop("Todas las capas deben tener al menos una neurona.")
}
y_train_oh <- one_hot(
y_train,
numero_clases = 6
)
y_val_oh <- one_hot(
y_val,
numero_clases = 6
)
red <- inicializar_red(
numero_entradas = ncol(x_train),
capas_ocultas = capas_ocultas,
numero_salidas = 6,
semilla = semilla
)
numero_capas_ocultas <- length(capas_ocultas)
indice_salida <- numero_capas_ocultas + 1
historial <- data.frame(
epoca = seq_len(epocas),
loss_train = NA_real_,
loss_val = NA_real_
)
for (epoca in seq_len(epocas)) {
# ==========================================================
# 1. PROPAGACIÓN HACIA ADELANTE: ENTRENAMIENTO
# ==========================================================
activaciones <- vector(
mode = "list",
length = numero_capas_ocultas + 1
)
preactivaciones <- vector(
mode = "list",
length = numero_capas_ocultas
)
activaciones[[1]] <- x_train
for (i in seq_len(numero_capas_ocultas)) {
z <- activaciones[[i]] %*%
red$pesos[[i]]
z <- sweep(
z,
MARGIN = 2,
STATS = as.numeric(red$sesgos[[i]]),
FUN = "+"
)
preactivaciones[[i]] <- z
activaciones[[i + 1]] <- relu(z)
}
activacion_final <-
activaciones[[numero_capas_ocultas + 1]]
pesos_salida <-
red$pesos[[indice_salida]]
# Verificación explícita
if (is.null(activacion_final)) {
stop("No se generó la activación de la última capa oculta.")
}
if (is.null(pesos_salida)) {
stop("No se generó la matriz de pesos de la capa de salida.")
}
if (ncol(activacion_final) != nrow(pesos_salida)) {
stop(
paste0(
"Dimensiones incompatibles: activación final = ",
paste(dim(activacion_final), collapse = " x "),
"; pesos de salida = ",
paste(dim(pesos_salida), collapse = " x ")
)
)
}
z_salida <- activacion_final %*%
pesos_salida
z_salida <- sweep(
z_salida,
MARGIN = 2,
STATS = as.numeric(
red$sesgos[[indice_salida]]
),
FUN = "+"
)
probabilidades <- softmax(z_salida)
historial$loss_train[epoca] <-
entropia_cruzada(
probabilidades,
y_train_oh
)
# ==========================================================
# 2. RETROPROPAGACIÓN
# ==========================================================
delta <- (
probabilidades - y_train_oh
) / nrow(x_train)
gradientes_pesos <- vector(
mode = "list",
length = indice_salida
)
gradientes_sesgos <- vector(
mode = "list",
length = indice_salida
)
# Gradiente de la capa de salida
gradientes_pesos[[indice_salida]] <-
t(activacion_final) %*% delta
gradientes_sesgos[[indice_salida]] <-
matrix(
colSums(delta),
nrow = 1
)
# Gradientes de las capas ocultas
for (i in seq(
from = numero_capas_ocultas,
to = 1,
by = -1
)) {
delta <- (
delta %*%
t(red$pesos[[i + 1]])
) * relu_derivada(
preactivaciones[[i]]
)
gradientes_pesos[[i]] <-
t(activaciones[[i]]) %*% delta
gradientes_sesgos[[i]] <-
matrix(
colSums(delta),
nrow = 1
)
}
# Actualización de pesos y sesgos
for (i in seq_len(indice_salida)) {
red$pesos[[i]] <-
red$pesos[[i]] -
tasa_aprendizaje *
gradientes_pesos[[i]]
red$sesgos[[i]] <-
red$sesgos[[i]] -
tasa_aprendizaje *
gradientes_sesgos[[i]]
}
# ==========================================================
# 3. PROPAGACIÓN HACIA ADELANTE: VALIDACIÓN
# ==========================================================
activacion_val <- x_val
for (i in seq_len(numero_capas_ocultas)) {
z_val <- activacion_val %*%
red$pesos[[i]]
z_val <- sweep(
z_val,
MARGIN = 2,
STATS = as.numeric(red$sesgos[[i]]),
FUN = "+"
)
activacion_val <- relu(z_val)
}
z_val_salida <- activacion_val %*%
red$pesos[[indice_salida]]
z_val_salida <- sweep(
z_val_salida,
MARGIN = 2,
STATS = as.numeric(
red$sesgos[[indice_salida]]
),
FUN = "+"
)
probabilidades_val <-
softmax(z_val_salida)
historial$loss_val[epoca] <-
entropia_cruzada(
probabilidades_val,
y_val_oh
)
}
return(
list(
red = red,
historial = historial
)
)
}dim(x_train)## [1] 145 5dim(x_val)## [1] 32 5red_revision <- inicializar_red(
numero_entradas = ncol(x_train),
capas_ocultas = c(4),
numero_salidas = 6,
semilla = 123
)
lapply(red_revision$pesos, dim)## [[1]]
## [1] 5 4
##
## [[2]]
## [1] 4 6## Ahora prueba primero una sola configuración:
modelo_prueba <- entrenar_mlp(
x_train = x_train,
y_train = y_train,
x_val = x_val,
y_val = y_val,
capas_ocultas = c(4),
tasa_aprendizaje = 0.01,
epocas = 50,
semilla = 123
)red_prueba <- inicializar_red(
numero_entradas = ncol(x_train),
capas_ocultas = c(4),
numero_salidas = 6,
semilla = 123
)
dim(x_train)## [1] 145 5lapply(red_prueba$pesos, dim)## [[1]]
## [1] 5 4
##
## [[2]]
## [1] 4 6# ============================================================
# 10. DEFINIR ARQUITECTURAS Y TASAS DE APRENDIZAJE
# ============================================================
arquitecturas <- list(
A1_1capa_4 = c(4),
A2_1capa_8 = c(8),
A3_2capas_8_4 = c(8, 4),
A4_2capas_12_6 = c(12, 6),
A5_3capas_12_8_4 = c(12, 8, 4),
A6_3capas_16_8_4 = c(16, 8, 4)
)
tasas_aprendizaje <- c(
0.001,
0.01,
0.1
)
arquitecturas## $A1_1capa_4
## [1] 4
##
## $A2_1capa_8
## [1] 8
##
## $A3_2capas_8_4
## [1] 8 4
##
## $A4_2capas_12_6
## [1] 12 6
##
## $A5_3capas_12_8_4
## [1] 12 8 4
##
## $A6_3capas_16_8_4
## [1] 16 8 4tasas_aprendizaje## [1] 0.001 0.010 0.100# ============================================================
# 11. ENTRENAR TODAS LAS CONFIGURACIONES
# ============================================================
resultados_modelos <- list()
historias_modelos <- list()
resumen_modelos <- data.frame()
contador <- 1
for (nombre_arquitectura in names(arquitecturas)) {
capas_actuales <-
arquitecturas[[nombre_arquitectura]]
for (tasa_actual in tasas_aprendizaje) {
nombre_modelo <- paste0(
nombre_arquitectura,
"_lr_",
tasa_actual
)
cat(
"\nEntrenando modelo ",
contador,
" de ",
length(arquitecturas) *
length(tasas_aprendizaje),
": ",
nombre_modelo,
"\n",
sep = ""
)
ajuste_actual <- tryCatch(
entrenar_mlp(
x_train = x_train,
y_train = y_train,
x_val = x_val,
y_val = y_val,
capas_ocultas = capas_actuales,
tasa_aprendizaje = tasa_actual,
epocas = 300,
semilla = 123
),
error = function(e) {
message(
"Error en ",
nombre_modelo,
": ",
e$message
)
NULL
}
)
if (!is.null(ajuste_actual)) {
historia_actual <-
ajuste_actual$historial
historia_actual$modelo <-
nombre_modelo
historia_actual$arquitectura <-
nombre_arquitectura
historia_actual$capas_ocultas <-
length(capas_actuales)
historia_actual$neuronas <-
paste(
capas_actuales,
collapse = "-"
)
historia_actual$learning_rate <-
tasa_actual
resultados_modelos[[nombre_modelo]] <-
ajuste_actual
historias_modelos[[nombre_modelo]] <-
historia_actual
mejor_epoca <-
which.min(
historia_actual$loss_val
)
resumen_modelos <- bind_rows(
resumen_modelos,
data.frame(
modelo = nombre_modelo,
capas_ocultas =
length(capas_actuales),
neuronas =
paste(
capas_actuales,
collapse = "-"
),
learning_rate =
tasa_actual,
mejor_epoca =
mejor_epoca,
loss_train_final =
tail(
historia_actual$loss_train,
1
),
loss_val_final =
tail(
historia_actual$loss_val,
1
),
loss_train_mejor_epoca =
historia_actual$loss_train[
mejor_epoca
],
loss_val_minima =
historia_actual$loss_val[
mejor_epoca
]
)
)
}
contador <- contador + 1
}
}##
## Entrenando modelo 1 de 18: A1_1capa_4_lr_0.001
##
## Entrenando modelo 2 de 18: A1_1capa_4_lr_0.01
##
## Entrenando modelo 3 de 18: A1_1capa_4_lr_0.1
##
## Entrenando modelo 4 de 18: A2_1capa_8_lr_0.001
##
## Entrenando modelo 5 de 18: A2_1capa_8_lr_0.01
##
## Entrenando modelo 6 de 18: A2_1capa_8_lr_0.1
##
## Entrenando modelo 7 de 18: A3_2capas_8_4_lr_0.001
##
## Entrenando modelo 8 de 18: A3_2capas_8_4_lr_0.01
##
## Entrenando modelo 9 de 18: A3_2capas_8_4_lr_0.1
##
## Entrenando modelo 10 de 18: A4_2capas_12_6_lr_0.001
##
## Entrenando modelo 11 de 18: A4_2capas_12_6_lr_0.01
##
## Entrenando modelo 12 de 18: A4_2capas_12_6_lr_0.1
##
## Entrenando modelo 13 de 18: A5_3capas_12_8_4_lr_0.001
##
## Entrenando modelo 14 de 18: A5_3capas_12_8_4_lr_0.01
##
## Entrenando modelo 15 de 18: A5_3capas_12_8_4_lr_0.1
##
## Entrenando modelo 16 de 18: A6_3capas_16_8_4_lr_0.001
##
## Entrenando modelo 17 de 18: A6_3capas_16_8_4_lr_0.01
##
## Entrenando modelo 18 de 18: A6_3capas_16_8_4_lr_0.1# ============================================================
# 12. UNIR LAS PÉRDIDAS DE TODOS LOS MODELOS
# ============================================================
curvas_completas <- bind_rows(
historias_modelos
)
dim(curvas_completas)## [1] 5400 8head(curvas_completas)# ============================================================
# 13. TABLA RESUMEN DE RESULTADOS
# ============================================================
resumen_modelos <- resumen_modelos |>
arrange(loss_val_minima) |>
mutate(
across(
c(
loss_train_final,
loss_val_final,
loss_train_mejor_epoca,
loss_val_minima
),
~ round(.x, 4)
)
)
resumen_modelosknitr::kable(
resumen_modelos,
caption = paste(
"Pérdida de entrenamiento y validación",
"de las arquitecturas evaluadas"
)
)| modelo | capas_ocultas | neuronas | learning_rate | mejor_epoca | loss_train_final | loss_val_final | loss_train_mejor_epoca | loss_val_minima |
|---|---|---|---|---|---|---|---|---|
| A2_1capa_8_lr_0.1 | 1 | 8 | 0.100 | 87 | 0.9201 | 1.6411 | 1.1181 | 1.4877 |
| A1_1capa_4_lr_0.1 | 1 | 4 | 0.100 | 177 | 0.9475 | 1.6324 | 1.0695 | 1.5554 |
| A4_2capas_12_6_lr_0.1 | 2 | 12-6 | 0.100 | 56 | 0.8267 | 1.9003 | 1.2208 | 1.5997 |
| A3_2capas_8_4_lr_0.1 | 2 | 8-4 | 0.100 | 88 | 0.9616 | 1.6852 | 1.2248 | 1.6020 |
| A2_1capa_8_lr_0.01 | 1 | 8 | 0.010 | 300 | 1.5426 | 1.6190 | 1.5426 | 1.6190 |
| A6_3capas_16_8_4_lr_0.1 | 3 | 16-8-4 | 0.100 | 36 | 0.9227 | 1.7116 | 1.3803 | 1.6355 |
| A6_3capas_16_8_4_lr_0.01 | 3 | 16-8-4 | 0.010 | 300 | 1.4116 | 1.6389 | 1.4116 | 1.6389 |
| A4_2capas_12_6_lr_0.01 | 2 | 12-6 | 0.010 | 300 | 1.4089 | 1.6629 | 1.4089 | 1.6629 |
| A5_3capas_12_8_4_lr_0.1 | 3 | 12-8-4 | 0.100 | 164 | 0.9209 | 1.8253 | 1.1658 | 1.6892 |
| A3_2capas_8_4_lr_0.01 | 2 | 8-4 | 0.010 | 300 | 1.5405 | 1.7117 | 1.5405 | 1.7117 |
| A1_1capa_4_lr_0.01 | 1 | 4 | 0.010 | 300 | 1.5432 | 1.7304 | 1.5432 | 1.7304 |
| A5_3capas_12_8_4_lr_0.01 | 3 | 12-8-4 | 0.010 | 224 | 1.6265 | 1.7961 | 1.6560 | 1.7923 |
| A6_3capas_16_8_4_lr_0.001 | 3 | 16-8-4 | 0.001 | 300 | 1.7124 | 1.8242 | 1.7124 | 1.8242 |
| A5_3capas_12_8_4_lr_0.001 | 3 | 12-8-4 | 0.001 | 300 | 1.8111 | 1.8575 | 1.8111 | 1.8575 |
| A3_2capas_8_4_lr_0.001 | 2 | 8-4 | 0.001 | 300 | 1.9770 | 1.8728 | 1.9770 | 1.8728 |
| A1_1capa_4_lr_0.001 | 1 | 4 | 0.001 | 300 | 1.8310 | 1.9032 | 1.8310 | 1.9032 |
| A4_2capas_12_6_lr_0.001 | 2 | 12-6 | 0.001 | 300 | 1.8946 | 1.9832 | 1.8946 | 1.9832 |
| A2_1capa_8_lr_0.001 | 1 | 8 | 0.001 | 300 | 2.1076 | 1.9867 | 2.1076 | 1.9867 |
# ============================================================
# 14. PREPARAR LAS CURVAS DE APRENDIZAJE
# ============================================================
curvas_loss <- curvas_completas |>
select(
epoca,
loss_train,
loss_val,
modelo,
capas_ocultas,
neuronas,
learning_rate
) |>
pivot_longer(
cols = c(
loss_train,
loss_val
),
names_to = "conjunto",
values_to = "perdida"
) |>
mutate(
conjunto = recode(
conjunto,
loss_train = "Entrenamiento",
loss_val = "Validación"
),
learning_rate = factor(
learning_rate,
levels = c(
0.001,
0.01,
0.1
)
)
)# ============================================================
# 15. CURVAS DE APRENDIZAJE DE LOS 18 MODELOS
# ============================================================
ggplot(
curvas_loss,
aes(
x = epoca,
y = perdida,
linetype = conjunto
)
) +
geom_line(
linewidth = 0.7
) +
facet_wrap(
~ modelo,
scales = "free_y",
ncol = 3
) +
labs(
title = paste(
"Curvas de aprendizaje",
"de los perceptrones multicapa"
),
subtitle = paste(
"Comparación de pérdida",
"en entrenamiento y validación"
),
x = "Época",
y = "Pérdida de entropía cruzada",
linetype = "Conjunto"
) +
theme_minimal() +
theme(
strip.text = element_text(
size = 7
)
)
# ============================================================
# 16. CURVAS ORGANIZADAS POR ARQUITECTURA Y TASA
# ============================================================
ggplot(
curvas_loss,
aes(
x = epoca,
y = perdida,
linetype = conjunto
)
) +
geom_line(
linewidth = 0.7
) +
facet_grid(
capas_ocultas + neuronas ~
learning_rate,
scales = "free_y",
labeller = label_both
) +
labs(
title = paste(
"Efecto de la arquitectura",
"y la tasa de aprendizaje"
),
x = "Época",
y = "Pérdida de entropía cruzada",
linetype = "Conjunto"
) +
theme_minimal()
# ============================================================
# 17. DIAGNÓSTICO DE LAS CURVAS
# ============================================================
diagnostico_modelos <- curvas_completas |>
group_by(
modelo,
capas_ocultas,
neuronas,
learning_rate
) |>
summarise(
mejor_epoca =
epoca[
which.min(loss_val)
],
loss_val_minima =
min(loss_val),
loss_train_mejor_epoca =
loss_train[
which.min(loss_val)
],
loss_train_final =
tail(
loss_train,
1
),
loss_val_final =
tail(
loss_val,
1
),
aumento_validacion =
loss_val_final -
loss_val_minima,
brecha_final =
loss_val_final -
loss_train_final,
reduccion_train =
first(loss_train) -
last(loss_train),
reduccion_val =
first(loss_val) -
last(loss_val),
.groups = "drop"
)diagnostico_modelos <- diagnostico_modelos |>
mutate(
diagnostico = case_when(
aumento_validacion > 0.10 &
brecha_final > 0.15 ~
"Indicios de sobreajuste",
reduccion_train < 0.10 &
reduccion_val < 0.10 ~
"Indicios de subajuste",
loss_train_final > 1.5 &
loss_val_final > 1.5 ~
"Indicios de subajuste",
TRUE ~
"Sin indicio marcado"
)
) |>
arrange(loss_val_minima)
diagnostico_modelosdiagnostico_presentacion <-
diagnostico_modelos |>
mutate(
across(
c(
loss_val_minima,
loss_train_mejor_epoca,
loss_train_final,
loss_val_final,
aumento_validacion,
brecha_final,
reduccion_train,
reduccion_val
),
~ round(.x, 4)
)
)
knitr::kable(
diagnostico_presentacion,
caption = paste(
"Diagnóstico de sobreajuste",
"y subajuste de los modelos"
)
)| modelo | capas_ocultas | neuronas | learning_rate | mejor_epoca | loss_val_minima | loss_train_mejor_epoca | loss_train_final | loss_val_final | aumento_validacion | brecha_final | reduccion_train | reduccion_val | diagnostico |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| A2_1capa_8_lr_0.1 | 1 | 8 | 0.100 | 87 | 1.4877 | 1.1181 | 0.9201 | 1.6411 | 0.1534 | 0.7210 | 1.3564 | 0.4178 | Indicios de sobreajuste |
| A1_1capa_4_lr_0.1 | 1 | 4 | 0.100 | 177 | 1.5554 | 1.0695 | 0.9475 | 1.6324 | 0.0770 | 0.6848 | 1.0469 | 0.3239 | Sin indicio marcado |
| A4_2capas_12_6_lr_0.1 | 2 | 12-6 | 0.100 | 56 | 1.5997 | 1.2208 | 0.8267 | 1.9003 | 0.3006 | 1.0736 | 1.2924 | 0.1500 | Indicios de sobreajuste |
| A3_2capas_8_4_lr_0.1 | 2 | 8-4 | 0.100 | 88 | 1.6020 | 1.2248 | 0.9616 | 1.6852 | 0.0832 | 0.7236 | 1.2485 | 0.2579 | Sin indicio marcado |
| A2_1capa_8_lr_0.01 | 1 | 8 | 0.010 | 300 | 1.6190 | 1.5426 | 1.5426 | 1.6190 | 0.0000 | 0.0763 | 0.7339 | 0.4823 | Indicios de subajuste |
| A6_3capas_16_8_4_lr_0.1 | 3 | 16-8-4 | 0.100 | 36 | 1.6355 | 1.3803 | 0.9227 | 1.7116 | 0.0760 | 0.7889 | 1.1090 | 0.2102 | Sin indicio marcado |
| A6_3capas_16_8_4_lr_0.01 | 3 | 16-8-4 | 0.010 | 300 | 1.6389 | 1.4116 | 1.4116 | 1.6389 | 0.0000 | 0.2273 | 0.6202 | 0.5813 | Sin indicio marcado |
| A4_2capas_12_6_lr_0.01 | 2 | 12-6 | 0.010 | 300 | 1.6629 | 1.4089 | 1.4089 | 1.6629 | 0.0000 | 0.2540 | 0.7103 | 0.4309 | Sin indicio marcado |
| A5_3capas_12_8_4_lr_0.1 | 3 | 12-8-4 | 0.100 | 164 | 1.6892 | 1.1658 | 0.9209 | 1.8253 | 0.1361 | 0.9044 | 1.1568 | 0.0990 | Indicios de sobreajuste |
| A3_2capas_8_4_lr_0.01 | 2 | 8-4 | 0.010 | 300 | 1.7117 | 1.5405 | 1.5405 | 1.7117 | 0.0000 | 0.1712 | 0.6695 | 0.2940 | Indicios de subajuste |
| A1_1capa_4_lr_0.01 | 1 | 4 | 0.010 | 300 | 1.7304 | 1.5432 | 1.5432 | 1.7304 | 0.0000 | 0.1872 | 0.4512 | 0.2618 | Indicios de subajuste |
| A5_3capas_12_8_4_lr_0.01 | 3 | 12-8-4 | 0.010 | 224 | 1.7923 | 1.6560 | 1.6265 | 1.7961 | 0.0037 | 0.1695 | 0.4512 | 0.2783 | Indicios de subajuste |
| A6_3capas_16_8_4_lr_0.001 | 3 | 16-8-4 | 0.001 | 300 | 1.8242 | 1.7124 | 1.7124 | 1.8242 | 0.0000 | 0.1118 | 0.3194 | 0.4369 | Indicios de subajuste |
| A5_3capas_12_8_4_lr_0.001 | 3 | 12-8-4 | 0.001 | 300 | 1.8575 | 1.8111 | 1.8111 | 1.8575 | 0.0000 | 0.0464 | 0.2666 | 0.2358 | Indicios de subajuste |
| A3_2capas_8_4_lr_0.001 | 2 | 8-4 | 0.001 | 300 | 1.8728 | 1.9770 | 1.9770 | 1.8728 | 0.0000 | -0.1042 | 0.2331 | 0.1399 | Indicios de subajuste |
| A1_1capa_4_lr_0.001 | 1 | 4 | 0.001 | 300 | 1.9032 | 1.8310 | 1.8310 | 1.9032 | 0.0000 | 0.0722 | 0.1634 | 0.0928 | Indicios de subajuste |
| A4_2capas_12_6_lr_0.001 | 2 | 12-6 | 0.001 | 300 | 1.9832 | 1.8946 | 1.8946 | 1.9832 | 0.0000 | 0.0886 | 0.2246 | 0.1155 | Indicios de subajuste |
| A2_1capa_8_lr_0.001 | 1 | 8 | 0.001 | 300 | 1.9867 | 2.1076 | 2.1076 | 1.9867 | 0.0000 | -0.1209 | 0.1689 | 0.1188 | Indicios de subajuste |
# ============================================================
# 18. MARCAR LA MEJOR ÉPOCA
# ============================================================
mejores_epocas <- curvas_completas |>
group_by(modelo) |>
slice_min(
order_by = loss_val,
n = 1,
with_ties = FALSE
) |>
ungroup()
ggplot(
curvas_loss,
aes(
x = epoca,
y = perdida,
linetype = conjunto
)
) +
geom_line(
linewidth = 0.7
) +
geom_vline(
data = mejores_epocas,
aes(
xintercept = epoca
),
linetype = "dotted"
) +
facet_wrap(
~ modelo,
scales = "free_y",
ncol = 3
) +
labs(
title = paste(
"Curvas de pérdida",
"y mejor época de validación"
),
subtitle = paste(
"La línea vertical indica",
"la menor pérdida de validación"
),
x = "Época",
y = "Pérdida",
linetype = "Conjunto"
) +
theme_minimal() +
theme(
strip.text = element_text(
size = 7
)
)
Las curvas de aprendizaje evidenciaron sobreajuste principalmente en las configuraciones con tasa de aprendizaje de 0,1. En estos modelos, la pérdida de entrenamiento continuó disminuyendo después de que la pérdida de validación alcanzó su mínimo, generando una separación progresiva entre ambas curvas. Este comportamiento fue especialmente marcado en las redes con dos y tres capas ocultas, debido a su mayor capacidad para ajustarse a las particularidades del conjunto de entrenamiento. La red con una capa oculta de ocho neuronas y tasa de aprendizaje de 0,1 obtuvo la menor pérdida de validación, 1,4877, en la época 87, pero después presentó sobreajuste. En contraste, las configuraciones con tasa de aprendizaje de 0,001 mostraron pérdidas elevadas y una reducción lenta durante las 300 épocas, lo que indicó subajuste o entrenamiento insuficiente. Las redes con tasa de 0,01 presentaron curvas más estables y una menor separación entre entrenamiento y validación, aunque varias aún se encontraban en proceso de aprendizaje al finalizar el entrenamiento. Por tanto, las curvas muestran que una tasa muy alta favoreció el sobreajuste, mientras que una tasa muy baja produjo aprendizaje insuficiente.
- Implemente y compare dos criterios para cuantificar el error durante el entrenamiento:
- Función de pérdida de entropía cruzada categórica (categorical cross-entropy).
- Error cuadrático medio (mean squared error).
¿Cuál es más apropiado para un problema de clasificación múltiple? Fundamente su respuesta tanto teórica como empíricamente a partir de los resultados obtenidos.
# ============================================================
# PUNTO 6. COMPARACIÓN DE FUNCIONES DE PÉRDIDA
# Entropía cruzada categórica vs error cuadrático medio
# ============================================================
# ============================================================
# 1. FUNCIONES DE PÉRDIDA
# ============================================================
perdida_entropia <- function(probabilidades, y_one_hot) {
epsilon <- 1e-12
-mean(
rowSums(
y_one_hot *
log(probabilidades + epsilon)
)
)
}
perdida_mse <- function(probabilidades, y_one_hot) {
mean(
(probabilidades - y_one_hot)^2
)
}
# ============================================================
# 2. FUNCIÓN PARA ENTRENAR CON CADA TIPO DE PÉRDIDA
# ============================================================
entrenar_mlp_perdida <- function(
x_train,
y_train,
x_val,
y_val,
capas_ocultas = c(8),
tasa_aprendizaje = 0.1,
epocas = 300,
tipo_perdida = c("entropia", "mse"),
semilla = 123
) {
tipo_perdida <- match.arg(tipo_perdida)
x_train <- as.matrix(x_train)
x_val <- as.matrix(x_val)
storage.mode(x_train) <- "double"
storage.mode(x_val) <- "double"
numero_clases <- 6
y_train_oh <- one_hot(
y_train,
numero_clases = numero_clases
)
y_val_oh <- one_hot(
y_val,
numero_clases = numero_clases
)
red <- inicializar_red(
numero_entradas = ncol(x_train),
capas_ocultas = capas_ocultas,
numero_salidas = numero_clases,
semilla = semilla
)
numero_capas_ocultas <- length(capas_ocultas)
indice_salida <- numero_capas_ocultas + 1
historial <- data.frame(
epoca = seq_len(epocas),
loss_train = NA_real_,
loss_val = NA_real_,
accuracy_train = NA_real_,
accuracy_val = NA_real_
)
mejor_loss_val <- Inf
mejor_epoca <- NA_integer_
mejor_red <- NULL
for (epoca in seq_len(epocas)) {
# --------------------------------------------------------
# PROPAGACIÓN HACIA ADELANTE: ENTRENAMIENTO
# --------------------------------------------------------
activaciones <- vector(
mode = "list",
length = numero_capas_ocultas + 1
)
preactivaciones <- vector(
mode = "list",
length = numero_capas_ocultas
)
activaciones[[1]] <- x_train
for (i in seq_len(numero_capas_ocultas)) {
z <- activaciones[[i]] %*%
red$pesos[[i]]
z <- sweep(
z,
MARGIN = 2,
STATS = as.numeric(red$sesgos[[i]]),
FUN = "+"
)
preactivaciones[[i]] <- z
activaciones[[i + 1]] <- relu(z)
}
activacion_final <-
activaciones[[numero_capas_ocultas + 1]]
z_salida <- activacion_final %*%
red$pesos[[indice_salida]]
z_salida <- sweep(
z_salida,
MARGIN = 2,
STATS = as.numeric(
red$sesgos[[indice_salida]]
),
FUN = "+"
)
probabilidades <- softmax(z_salida)
# --------------------------------------------------------
# CALCULAR PÉRDIDA DE ENTRENAMIENTO
# --------------------------------------------------------
if (tipo_perdida == "entropia") {
loss_train <- perdida_entropia(
probabilidades,
y_train_oh
)
} else {
loss_train <- perdida_mse(
probabilidades,
y_train_oh
)
}
historial$loss_train[epoca] <- loss_train
clases_train_pred <- max.col(
probabilidades,
ties.method = "first"
)
historial$accuracy_train[epoca] <-
mean(clases_train_pred == y_train)
# --------------------------------------------------------
# GRADIENTE DE LA CAPA DE SALIDA
# --------------------------------------------------------
if (tipo_perdida == "entropia") {
# Softmax + entropía cruzada
delta <- (
probabilidades - y_train_oh
) / nrow(x_train)
} else {
# MSE + softmax
#
# Primero se calcula dL/dp:
# 2 * (p - y) / (número de observaciones * clases)
gradiente_probabilidades <-
2 *
(probabilidades - y_train_oh) /
(
nrow(x_train) *
numero_clases
)
# Se aplica el Jacobiano de softmax:
# dL/dz = p * [g - sum(g * p)]
producto_fila <- rowSums(
gradiente_probabilidades *
probabilidades
)
delta <- probabilidades * (
gradiente_probabilidades -
producto_fila
)
}
# --------------------------------------------------------
# RETROPROPAGACIÓN
# --------------------------------------------------------
gradientes_pesos <- vector(
mode = "list",
length = indice_salida
)
gradientes_sesgos <- vector(
mode = "list",
length = indice_salida
)
gradientes_pesos[[indice_salida]] <-
t(activacion_final) %*% delta
gradientes_sesgos[[indice_salida]] <-
matrix(
colSums(delta),
nrow = 1
)
for (
i in seq(
from = numero_capas_ocultas,
to = 1,
by = -1
)
) {
delta <- (
delta %*%
t(red$pesos[[i + 1]])
) *
relu_derivada(
preactivaciones[[i]]
)
gradientes_pesos[[i]] <-
t(activaciones[[i]]) %*%
delta
gradientes_sesgos[[i]] <-
matrix(
colSums(delta),
nrow = 1
)
}
# --------------------------------------------------------
# ACTUALIZAR PESOS Y SESGOS
# --------------------------------------------------------
for (i in seq_len(indice_salida)) {
red$pesos[[i]] <-
red$pesos[[i]] -
tasa_aprendizaje *
gradientes_pesos[[i]]
red$sesgos[[i]] <-
red$sesgos[[i]] -
tasa_aprendizaje *
gradientes_sesgos[[i]]
}
# --------------------------------------------------------
# PROPAGACIÓN HACIA ADELANTE: VALIDACIÓN
# --------------------------------------------------------
activacion_val <- x_val
for (i in seq_len(numero_capas_ocultas)) {
z_val <- activacion_val %*%
red$pesos[[i]]
z_val <- sweep(
z_val,
MARGIN = 2,
STATS = as.numeric(
red$sesgos[[i]]
),
FUN = "+"
)
activacion_val <- relu(z_val)
}
z_val_salida <- activacion_val %*%
red$pesos[[indice_salida]]
z_val_salida <- sweep(
z_val_salida,
MARGIN = 2,
STATS = as.numeric(
red$sesgos[[indice_salida]]
),
FUN = "+"
)
probabilidades_val <- softmax(
z_val_salida
)
if (tipo_perdida == "entropia") {
loss_val <- perdida_entropia(
probabilidades_val,
y_val_oh
)
} else {
loss_val <- perdida_mse(
probabilidades_val,
y_val_oh
)
}
historial$loss_val[epoca] <- loss_val
clases_val_pred <- max.col(
probabilidades_val,
ties.method = "first"
)
historial$accuracy_val[epoca] <-
mean(clases_val_pred == y_val)
# --------------------------------------------------------
# GUARDAR EL MODELO DE LA MEJOR ÉPOCA
# --------------------------------------------------------
if (
is.finite(loss_val) &&
loss_val < mejor_loss_val
) {
mejor_loss_val <- loss_val
mejor_epoca <- epoca
mejor_red <- red
}
}
list(
red_final = red,
mejor_red = mejor_red,
mejor_epoca = mejor_epoca,
mejor_loss_val = mejor_loss_val,
historial = historial,
tipo_perdida = tipo_perdida
)
}
# ============================================================
# 3. ENTRENAR CON ENTROPÍA CRUZADA
# ============================================================
modelo_entropia <- entrenar_mlp_perdida(
x_train = x_train,
y_train = y_train,
x_val = x_val,
y_val = y_val,
capas_ocultas = c(8),
tasa_aprendizaje = 0.1,
epocas = 300,
tipo_perdida = "entropia",
semilla = 123
)
# ============================================================
# 4. ENTRENAR CON ERROR CUADRÁTICO MEDIO
# ============================================================
modelo_mse <- entrenar_mlp_perdida(
x_train = x_train,
y_train = y_train,
x_val = x_val,
y_val = y_val,
capas_ocultas = c(8),
tasa_aprendizaje = 0.1,
epocas = 300,
tipo_perdida = "mse",
semilla = 123
)
# ============================================================
# 5. REVISAR RESULTADOS GENERALES
# ============================================================
modelo_entropia$mejor_epoca## [1] 87modelo_entropia$mejor_loss_val## [1] 1.487721modelo_mse$mejor_epoca## [1] 300modelo_mse$mejor_loss_val## [1] 0.1383438En el modelo con entropía cruzada, el mejor resultado en validación se obtuvo en la época 87. A partir de ese punto, la pérdida de validación comenzó a aumentar, mientras que la pérdida de entrenamiento siguió disminuyendo. Esto indica que el modelo empezó a memorizar los datos de entrenamiento y a perder capacidad para clasificar correctamente datos nuevos, es decir, comenzó a sobreajustarse.
En el modelo entrenado con MSE, la pérdida de validación más baja se alcanzó en la época 300. Esto muestra que, según esta función de pérdida, el modelo continuó mejorando lentamente hasta el final del entrenamiento. Sin embargo, una pérdida menor no significa necesariamente una mejor clasificación. Por esta razón, también se deben analizar la exactitud y las matrices de confusión.
# ============================================================
# 6. FUNCIÓN DE PREDICCIÓN
# ============================================================
predecir_mlp <- function(
red,
x,
numero_capas_ocultas
) {
x <- as.matrix(x)
storage.mode(x) <- "double"
activacion <- x
for (i in seq_len(numero_capas_ocultas)) {
z <- activacion %*%
red$pesos[[i]]
z <- sweep(
z,
MARGIN = 2,
STATS = as.numeric(red$sesgos[[i]]),
FUN = "+"
)
activacion <- relu(z)
}
indice_salida <- numero_capas_ocultas + 1
z_salida <- activacion %*%
red$pesos[[indice_salida]]
z_salida <- sweep(
z_salida,
MARGIN = 2,
STATS = as.numeric(
red$sesgos[[indice_salida]]
),
FUN = "+"
)
probabilidades <- softmax(z_salida)
clases <- max.col(
probabilidades,
ties.method = "first"
)
list(
probabilidades = probabilidades,
clases = clases
)
}# ============================================================
# 7. FUNCIÓN DE MÉTRICAS MULTICLASE
# ============================================================
calcular_metricas_multiclase <- function(
observado,
predicho,
clases = 1:6
) {
observado <- factor(
observado,
levels = clases
)
predicho <- factor(
predicho,
levels = clases
)
matriz <- table(
Real = observado,
Predicho = predicho
)
exactitud <- sum(diag(matriz)) /
sum(matriz)
sensibilidad <- numeric(
length(clases)
)
precision <- numeric(
length(clases)
)
f1 <- numeric(
length(clases)
)
for (i in seq_along(clases)) {
verdaderos_positivos <- matriz[i, i]
falsos_negativos <-
sum(matriz[i, ]) -
verdaderos_positivos
falsos_positivos <-
sum(matriz[, i]) -
verdaderos_positivos
sensibilidad[i] <- ifelse(
verdaderos_positivos +
falsos_negativos == 0,
NA,
verdaderos_positivos /
(
verdaderos_positivos +
falsos_negativos
)
)
precision[i] <- ifelse(
verdaderos_positivos +
falsos_positivos == 0,
NA,
verdaderos_positivos /
(
verdaderos_positivos +
falsos_positivos
)
)
f1[i] <- ifelse(
is.na(sensibilidad[i]) ||
is.na(precision[i]) ||
sensibilidad[i] +
precision[i] == 0,
0,
2 *
sensibilidad[i] *
precision[i] /
(
sensibilidad[i] +
precision[i]
)
)
}
list(
matriz_confusion = matriz,
accuracy = exactitud,
balanced_accuracy = mean(
sensibilidad,
na.rm = TRUE
),
macro_f1 = mean(
f1,
na.rm = TRUE
),
metricas_clase = data.frame(
clase = clases,
sensibilidad = sensibilidad,
precision = precision,
f1 = f1
)
)
}# ============================================================
# 8. PREDICCIONES EN VALIDACIÓN
# ============================================================
pred_entropia <- predecir_mlp(
red = modelo_entropia$mejor_red,
x = x_val,
numero_capas_ocultas = 1
)
pred_mse <- predecir_mlp(
red = modelo_mse$mejor_red,
x = x_val,
numero_capas_ocultas = 1
)
# ============================================================
# 9. MÉTRICAS EN VALIDACIÓN
# ============================================================
metricas_entropia <- calcular_metricas_multiclase(
observado = y_val,
predicho = pred_entropia$clases
)
metricas_mse <- calcular_metricas_multiclase(
observado = y_val,
predicho = pred_mse$clases
)
metricas_entropia$matriz_confusion## Predicho
## Real 1 2 3 4 5 6
## 1 0 4 0 0 0 0
## 2 0 2 3 2 0 0
## 3 0 0 1 3 1 0
## 4 0 1 1 3 2 0
## 5 0 0 0 1 5 1
## 6 0 0 0 1 1 0metricas_mse$matriz_confusion## Predicho
## Real 1 2 3 4 5 6
## 1 0 4 0 0 0 0
## 2 0 4 0 3 0 0
## 3 1 0 1 3 0 0
## 4 1 1 0 5 0 0
## 5 2 0 0 5 0 0
## 6 0 0 0 2 0 0- Modelo de entropía cruzada
En total hubo 11 clasificaciones correctas de 32 observaciones=0,3438
El modelo de entropía cruzada reconoció mejor la clase 5. La clase 1 fue confundida completamente con la clase 2, mientras que la clase 6 fue confundida con las clases 4 y 5. También se observa que la mayoría de los errores ocurrieron entre categorías cercanas de severidad. Por ejemplo, las categorías 3 y 4 se confundieron frecuentemente entre sí. Esto puede deberse a que las clases vecinas pueden presentar características espectrales semejantes.
- Modelo de error cuadrático medio
En total fueron 10 observaciones correctas de 32=0,3125
Este modelo concentró gran parte de sus predicciones en las clases 2 y 4. Todas las observaciones reales de la clase 5 fueron clasificadas como 1 o 4, y las de la clase 6 fueron clasificadas como 4. Esto indica que el modelo produjo una clasificación menos equilibrada. Aunque reconoció relativamente bien las clases 2 y 4, ignoró por completo otras categorías, especialmente 1, 5 y 6.
# ============================================================
# 10. TABLA COMPARATIVA
# ============================================================
comparacion_perdidas <- data.frame(
funcion_perdida = c(
"Entropía cruzada",
"Error cuadrático medio"
),
mejor_epoca = c(
modelo_entropia$mejor_epoca,
modelo_mse$mejor_epoca
),
perdida_validacion_minima = c(
modelo_entropia$mejor_loss_val,
modelo_mse$mejor_loss_val
),
accuracy_validacion = c(
metricas_entropia$accuracy,
metricas_mse$accuracy
),
balanced_accuracy = c(
metricas_entropia$balanced_accuracy,
metricas_mse$balanced_accuracy
),
macro_f1 = c(
metricas_entropia$macro_f1,
metricas_mse$macro_f1
)
)
comparacion_perdidas <- comparacion_perdidas |>
mutate(
across(
where(is.numeric),
~ round(.x, 4)
)
)
comparacion_perdidasknitr::kable(
comparacion_perdidas,
caption = paste(
"Comparación de entropía cruzada",
"y error cuadrático medio"
)
)| funcion_perdida | mejor_epoca | perdida_validacion_minima | accuracy_validacion | balanced_accuracy | macro_f1 |
|---|---|---|---|---|---|
| Entropía cruzada | 87 | 1.4877 | 0.3438 | 0.2714 | 0.2439 |
| Error cuadrático medio | 300 | 0.1383 | 0.3125 | 0.2476 | 0.2056 |
La entropía cruzada alcanzó una exactitud de 34,38 %, mientras el MSE obtuvo 31,25 %. La diferencia es moderada, pero favorece a la entropía cruzada.
La exactitud balanceada calcula el promedio de la sensibilidad de las seis categorías y es especialmente importante porque las clases tienen frecuencias diferentes. La mayor exactitud balanceada de la entropía cruzada indica que, en promedio, reconoció mejor las diferentes categorías y no se concentró tanto en unas pocas clases.
F1 macro da la misma importancia a todas las categorías. El mayor valor de la entropía cruzada confirma que logró un equilibrio ligeramente mejor entre sensibilidad y precisión. Aunque ambos valores son bajos, las tres métricas coinciden en favorecer la entropía cruzada.
#CONSTRUIR CURV DE APRENDIZAJE
# ============================================================
# 11. UNIR HISTORIALES
# ============================================================
historial_entropia <-
modelo_entropia$historial |>
mutate(
funcion_perdida =
"Entropía cruzada"
)
historial_mse <-
modelo_mse$historial |>
mutate(
funcion_perdida =
"Error cuadrático medio"
)
historial_comparacion <- bind_rows(
historial_entropia,
historial_mse
)# ============================================================
# 12. CURVAS DE PÉRDIDA
# ============================================================
curvas_comparacion <-
historial_comparacion |>
select(
epoca,
loss_train,
loss_val,
funcion_perdida
) |>
pivot_longer(
cols = c(
loss_train,
loss_val
),
names_to = "conjunto",
values_to = "perdida"
) |>
mutate(
conjunto = recode(
conjunto,
loss_train = "Entrenamiento",
loss_val = "Validación"
)
)
ggplot(
curvas_comparacion,
aes(
x = epoca,
y = perdida,
linetype = conjunto
)
) +
geom_line(
linewidth = 0.8
) +
facet_wrap(
~ funcion_perdida,
scales = "free_y"
) +
labs(
title = paste(
"Comparación de funciones",
"de pérdida durante el entrenamiento"
),
subtitle = paste(
"Arquitectura: una capa oculta",
"con ocho neuronas"
),
x = "Época",
y = "Pérdida",
linetype = "Conjunto"
) +
theme_minimal()
#Comparar la exactitud por época
# ============================================================
# 13. EXACTITUD POR ÉPOCA
# ============================================================
accuracy_comparacion <-
historial_comparacion |>
select(
epoca,
accuracy_train,
accuracy_val,
funcion_perdida
) |>
pivot_longer(
cols = c(
accuracy_train,
accuracy_val
),
names_to = "conjunto",
values_to = "accuracy"
) |>
mutate(
conjunto = recode(
conjunto,
accuracy_train = "Entrenamiento",
accuracy_val = "Validación"
)
)
ggplot(
accuracy_comparacion,
aes(
x = epoca,
y = accuracy,
linetype = conjunto
)
) +
geom_line(
linewidth = 0.8
) +
facet_wrap(
~ funcion_perdida
) +
labs(
title = paste(
"Exactitud durante el entrenamiento",
"según la función de pérdida"
),
x = "Época",
y = "Exactitud",
linetype = "Conjunto"
) +
theme_minimal()
Al comparar las dos funciones de pérdida, en los dos casos se utilizó la misma red neuronal, con una capa oculta de ocho neuronas, una tasa de aprendizaje de 0,1, la misma inicialización y el mismo número de épocas.
La entropía cruzada obtuvo su menor pérdida de validación en la época 87, mientras que el MSE siguió disminuyendo hasta la época 300. Sin embargo, los valores de pérdida no se compararon directamente porque cada función utiliza una escala diferente. Por esta razón, la comparación se hizo mediante métricas de clasificación.
La entropía cruzada obtuvo mejores resultados que el MSE. Alcanzó una exactitud de 0,3438, una exactitud balanceada de 0,2714 y un F1 macro de 0,2439. En cambio, el MSE obtuvo valores de 0,3125, 0,2476 y 0,2056, respectivamente.
Las curvas mostraron que la entropía cruzada aprendió más rápido, aunque comenzó a presentar sobreajuste después de la época 87. El MSE mostró curvas más estables y cercanas entre sí, pero su exactitud fue baja y sus predicciones se concentraron en pocas categorías, lo que indica que el aprendizaje fue insuficiente.
Por lo tanto, la entropía cruzada categórica fue la opción más adecuada para este problema de clasificación múltiple. Esto se debe a que funciona mejor con la salida softmax y porque obtuvo mejores resultados de validación. Para reducir el sobreajuste, se recomienda conservar el modelo de la época 87 mediante parada temprana.
2.4. Métricas de evaluación
- Para el modelo seleccionado como mejor en la Pregunta 5, reporte las siguientes métricas calculadas sobre el conjunto de prueba:
• Exactitud global (accuracy). • Precisión, sensibilidad (recall) y F1-score por clase. • Macro-promedio y promedio ponderado de precisión, sensibilidad y F1. • Matriz de confusión completa con interpretación. • Curva ROC y área bajo la curva (AUC-ROC) para cada clase bajo el esquema unocontra-todos, si la librería lo permite.
Interprete los resultados: ¿qué categorías de severidad son más difíciles de clasificar correctamente? ¿A qué lo atribuye?
# ============================================================
# PREGUNTA 7. MÉTRICAS DE EVALUACIÓN EN PRUEBA
# ============================================================
library(dplyr)
library(tidyr)
library(ggplot2)
# ============================================================
# 1. PREPARAR EL CONJUNTO DE PRUEBA
# ============================================================
predictores <- c(
"ndvi_med",
"evi_med",
"ndre_med",
"gli_med",
"height_med"
)
x_test <- as.matrix(
test_15_std[, predictores]
)
storage.mode(x_test) <- "double"
y_test <- as.integer(
factor(
test_15_std$severity,
levels = 1:6
)
)
dim(x_test)## [1] 35 5table(y_test)## y_test
## 1 2 3 4 5 6
## 5 8 6 7 7 2# ============================================================
# 2. SELECCIONAR EL MEJOR MODELO
# ============================================================
mejor_modelo_p5 <- modelo_entropia$mejor_red
modelo_entropia$mejor_epoca## [1] 87modelo_entropia$mejor_loss_val## [1] 1.487721# ============================================================
# 3. PREDICCIONES EN EL CONJUNTO DE PRUEBA
# ============================================================
prediccion_test <- predecir_mlp(
red = mejor_modelo_p5,
x = x_test,
numero_capas_ocultas = 1
)
probabilidades_test <-
prediccion_test$probabilidades
clases_predichas_test <-
prediccion_test$clasesdim(probabilidades_test)## [1] 35 6head(probabilidades_test)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.2712215 0.4395784 0.1403482 0.02924471 0.06391375 0.05569347
## [2,] 0.2906438 0.4877297 0.1130572 0.01538613 0.04065501 0.05252820
## [3,] 0.2710579 0.5374107 0.1086576 0.01181567 0.03148892 0.03956926
## [4,] 0.2080008 0.5806604 0.1473475 0.01574559 0.02278153 0.02546421
## [5,] 0.1454605 0.5726357 0.2274655 0.02477962 0.01590118 0.01375744
## [6,] 0.1992433 0.6278977 0.1298999 0.01064111 0.01458717 0.01773076head(clases_predichas_test)## [1] 2 2 2 2 2 2head(
rowSums(probabilidades_test)
)## [1] 1 1 1 1 1 1# ============================================================
# 4. FUNCIÓN COMPLETA DE MÉTRICAS MULTICLASE
# ============================================================
calcular_metricas_completas <- function(
observado,
predicho,
clases = 1:6
) {
observado <- factor(
observado,
levels = clases
)
predicho <- factor(
predicho,
levels = clases
)
matriz <- table(
Real = observado,
Predicho = predicho
)
total <- sum(matriz)
accuracy <- sum(diag(matriz)) / total
resultados_clase <- data.frame(
clase = clases,
soporte = as.numeric(
rowSums(matriz)
),
precision = NA_real_,
sensibilidad = NA_real_,
f1 = NA_real_
)
for (i in seq_along(clases)) {
vp <- matriz[i, i]
fn <- sum(matriz[i, ]) - vp
fp <- sum(matriz[, i]) - vp
precision_i <- ifelse(
vp + fp == 0,
0,
vp / (vp + fp)
)
sensibilidad_i <- ifelse(
vp + fn == 0,
0,
vp / (vp + fn)
)
f1_i <- ifelse(
precision_i + sensibilidad_i == 0,
0,
2 *
precision_i *
sensibilidad_i /
(
precision_i +
sensibilidad_i
)
)
resultados_clase$precision[i] <-
precision_i
resultados_clase$sensibilidad[i] <-
sensibilidad_i
resultados_clase$f1[i] <-
f1_i
}
macro_precision <- mean(
resultados_clase$precision
)
macro_sensibilidad <- mean(
resultados_clase$sensibilidad
)
macro_f1 <- mean(
resultados_clase$f1
)
pesos <- resultados_clase$soporte /
sum(resultados_clase$soporte)
ponderada_precision <- sum(
resultados_clase$precision *
pesos
)
ponderada_sensibilidad <- sum(
resultados_clase$sensibilidad *
pesos
)
ponderada_f1 <- sum(
resultados_clase$f1 *
pesos
)
resumen_promedios <- data.frame(
promedio = c(
"Macro",
"Ponderado"
),
precision = c(
macro_precision,
ponderada_precision
),
sensibilidad = c(
macro_sensibilidad,
ponderada_sensibilidad
),
f1 = c(
macro_f1,
ponderada_f1
)
)
list(
matriz_confusion = matriz,
accuracy = accuracy,
metricas_clase = resultados_clase,
resumen_promedios = resumen_promedios
)
}# ============================================================
# 5. CALCULAR MÉTRICAS EN PRUEBA
# ============================================================
evaluacion_test <- calcular_metricas_completas(
observado = y_test,
predicho = clases_predichas_test,
clases = 1:6
)evaluacion_test$accuracy## [1] 0.4285714## Matriz de confusión:
evaluacion_test$matriz_confusion## Predicho
## Real 1 2 3 4 5 6
## 1 0 4 1 0 0 0
## 2 0 7 0 1 0 0
## 3 0 4 0 2 0 0
## 4 0 0 1 3 3 0
## 5 0 0 0 2 5 0
## 6 0 0 0 0 2 0## Métricas por clase:
evaluacion_test$metricas_clase## Macro-promedios y promedios ponderados:
evaluacion_test$resumen_promedio# ============================================================
# 6. TABLA DE EXACTITUD GLOBAL
# ============================================================
tabla_accuracy <- data.frame(
metrica = "Exactitud global",
valor = evaluacion_test$accuracy
) |>
mutate(
valor = round(valor, 4)
)
knitr::kable(
tabla_accuracy,
caption = "Exactitud global del modelo en el conjunto de prueba"
)| metrica | valor |
|---|---|
| Exactitud global | 0.4286 |
# ============================================================
# 7. MÉTRICAS POR CLASE
# ============================================================
metricas_por_clase <-
evaluacion_test$metricas_clase |>
mutate(
categoria = paste(
"Severidad",
clase
),
across(
c(
precision,
sensibilidad,
f1
),
~ round(.x, 4)
)
) |>
select(
categoria,
soporte,
precision,
sensibilidad,
f1
)
knitr::kable(
metricas_por_clase,
caption = paste(
"Precisión, sensibilidad y F1-score",
"por categoría de severidad"
)
)| categoria | soporte | precision | sensibilidad | f1 |
|---|---|---|---|---|
| Severidad 1 | 5 | 0.0000 | 0.0000 | 0.0000 |
| Severidad 2 | 8 | 0.4667 | 0.8750 | 0.6087 |
| Severidad 3 | 6 | 0.0000 | 0.0000 | 0.0000 |
| Severidad 4 | 7 | 0.3750 | 0.4286 | 0.4000 |
| Severidad 5 | 7 | 0.5000 | 0.7143 | 0.5882 |
| Severidad 6 | 2 | 0.0000 | 0.0000 | 0.0000 |
# ============================================================
# 8. MACRO-PROMEDIO Y PROMEDIO PONDERADO
# ============================================================
tabla_promedios <-
evaluacion_test$resumen_promedios |>
mutate(
across(
c(
precision,
sensibilidad,
f1
),
~ round(.x, 4)
)
)
knitr::kable(
tabla_promedios,
caption = paste(
"Macro-promedio y promedio ponderado",
"de las métricas de clasificación"
)
)| promedio | precision | sensibilidad | f1 |
|---|---|---|---|
| Macro | 0.2236 | 0.3363 | 0.2662 |
| Ponderado | 0.2817 | 0.4286 | 0.3368 |
# ============================================================
# 9. MATRIZ DE CONFUSIÓN
# ============================================================
matriz_confusion_df <-
as.data.frame.matrix(
evaluacion_test$matriz_confusion
)
matriz_confusion_df <-
cbind(
Clase_real = rownames(
matriz_confusion_df
),
matriz_confusion_df
)
rownames(matriz_confusion_df) <- NULL
knitr::kable(
matriz_confusion_df,
caption = "Matriz de confusión del conjunto de prueba"
)| Clase_real | 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|---|
| 1 | 0 | 4 | 1 | 0 | 0 | 0 |
| 2 | 0 | 7 | 0 | 1 | 0 | 0 |
| 3 | 0 | 4 | 0 | 2 | 0 | 0 |
| 4 | 0 | 0 | 1 | 3 | 3 | 0 |
| 5 | 0 | 0 | 0 | 2 | 5 | 0 |
| 6 | 0 | 0 | 0 | 0 | 2 | 0 |
# ============================================================
# 10. GRÁFICO DE LA MATRIZ DE CONFUSIÓN
# ============================================================
matriz_larga <- as.data.frame(
evaluacion_test$matriz_confusion
)
names(matriz_larga) <- c(
"Real",
"Predicho",
"Frecuencia"
)
ggplot(
matriz_larga,
aes(
x = Predicho,
y = Real,
fill = Frecuencia
)
) +
geom_tile() +
geom_text(
aes(
label = Frecuencia
)
) +
labs(
title = "Matriz de confusión del modelo seleccionado",
subtitle = paste(
"Filas: categoría real;",
"columnas: categoría predicha"
),
x = "Categoría predicha",
y = "Categoría real",
fill = "Frecuencia"
) +
theme_minimal()
# ============================================================
# 11. FUNCIÓN ROC Y AUC BINARIA
# ============================================================
calcular_roc_binaria <- function(
observado_binario,
probabilidad
) {
observado_binario <-
as.integer(observado_binario)
orden <- order(
probabilidad,
decreasing = TRUE
)
y_ordenado <-
observado_binario[orden]
prob_ordenada <-
probabilidad[orden]
positivos <- sum(
y_ordenado == 1
)
negativos <- sum(
y_ordenado == 0
)
if (
positivos == 0 ||
negativos == 0
) {
return(
list(
curva = data.frame(
fpr = NA_real_,
tpr = NA_real_,
umbral = NA_real_
),
auc = NA_real_
)
)
}
vp_acumulados <- cumsum(
y_ordenado == 1
)
fp_acumulados <- cumsum(
y_ordenado == 0
)
tpr <- vp_acumulados /
positivos
fpr <- fp_acumulados /
negativos
curva <- data.frame(
fpr = c(0, fpr, 1),
tpr = c(0, tpr, 1),
umbral = c(
Inf,
prob_ordenada,
-Inf
)
)
curva <- curva |>
arrange(fpr, tpr)
auc <- sum(
diff(curva$fpr) *
(
head(curva$tpr, -1) +
tail(curva$tpr, -1)
) / 2
)
list(
curva = curva,
auc = auc
)
}# ============================================================
# 12. ROC UNO-CONTRA-TODOS
# ============================================================
curvas_roc <- list()
tabla_auc <- data.frame()
for (clase_actual in 1:6) {
observado_binario <- ifelse(
y_test == clase_actual,
1,
0
)
roc_actual <- calcular_roc_binaria(
observado_binario =
observado_binario,
probabilidad =
probabilidades_test[
,
clase_actual
]
)
curva_actual <- roc_actual$curva
curva_actual$clase <- paste(
"Severidad",
clase_actual
)
curvas_roc[[clase_actual]] <-
curva_actual
tabla_auc <- bind_rows(
tabla_auc,
data.frame(
clase = paste(
"Severidad",
clase_actual
),
auc = roc_actual$auc
)
)
}
curvas_roc_df <- bind_rows(
curvas_roc
)tabla_auc_presentacion <-
tabla_auc |>
mutate(
auc = round(auc, 4)
)
knitr::kable(
tabla_auc_presentacion,
caption = paste(
"Área bajo la curva ROC",
"uno-contra-todos por clase"
)
)| clase | auc |
|---|---|
| Severidad 1 | 0.7267 |
| Severidad 2 | 0.7963 |
| Severidad 3 | 0.6724 |
| Severidad 4 | 0.8163 |
| Severidad 5 | 0.8673 |
| Severidad 6 | 1.0000 |
# ============================================================
# 13. GRÁFICO ROC POR CLASE
# ============================================================
ggplot(
curvas_roc_df,
aes(
x = fpr,
y = tpr
)
) +
geom_line(
linewidth = 0.9
) +
geom_abline(
intercept = 0,
slope = 1,
linetype = "dashed"
) +
facet_wrap(
~ clase,
ncol = 3
) +
coord_equal() +
labs(
title = "Curvas ROC uno-contra-todos",
subtitle = paste(
"Modelo seleccionado:",
"una capa oculta con ocho neuronas"
),
x = "Tasa de falsos positivos",
y = "Tasa de verdaderos positivos"
) +
theme_minimal()
# ============================================================
# 14. GRÁFICO DE MÉTRICAS POR CLASE
# ============================================================
metricas_grafico <-
evaluacion_test$metricas_clase |>
select(
clase,
precision,
sensibilidad,
f1
) |>
pivot_longer(
cols = c(
precision,
sensibilidad,
f1
),
names_to = "metrica",
values_to = "valor"
) |>
mutate(
clase = paste(
"Severidad",
clase
),
metrica = recode(
metrica,
precision = "Precisión",
sensibilidad = "Sensibilidad",
f1 = "F1-score"
)
)
ggplot(
metricas_grafico,
aes(
x = clase,
y = valor,
group = metrica,
linetype = metrica
)
) +
geom_point(
size = 2
) +
geom_line() +
scale_y_continuous(
limits = c(0, 1)
) +
labs(
title = "Desempeño del modelo por categoría de severidad",
x = "Categoría real",
y = "Valor de la métrica",
linetype = "Métrica"
) +
theme_minimal()
Las severidades 1, 3 y 6 fueron las más difíciles de clasificar, ya que
ninguna de sus observaciones fue identificada correctamente. Las clases
1 y 3 se confundieron principalmente con niveles cercanos, mientras que
los dos casos de severidad 6 fueron clasificados como severidad 5.
Aunque esta última presentó un AUC de 1,00, el resultado debe
interpretarse con cautela debido al reducido número de
observaciones.
Las severidades 4 y 5 mostraron un desempeño intermedio, mientras que la severidad 2 fue la mejor reconocida. Sin embargo, varias observaciones de otras clases también fueron asignadas a esta categoría, lo que redujo su precisión.
En general, la mayoría de los errores ocurrieron entre severidades consecutivas. Esto indica que el modelo logró reconocer parcialmente el orden del daño, pero no pudo separar con claridad las categorías cercanas.
Estas dificultades pueden explicarse porque la severidad es un proceso continuo dividido en categorías, por lo que los niveles cercanos pueden presentar valores espectrales similares. Además, los índices utilizados pueden solaparse entre clases, la altura puede estar afectada por factores distintos de la enfermedad y existe desbalance entre las categorías. El pequeño conjunto de prueba, con solo 35 observaciones, también hace que las métricas cambien considerablemente con cada predicción.
Por esta razón, la exactitud global de 0,4286 debe interpretarse con precaución. El macro F1 de 0,2662 muestra que el desempeño fue bajo y desigual entre las categorías, mientras que el promedio ponderado de 0,3368 estuvo influido por las clases más frecuentes y mejor clasificadas.
3. Parte II: Clasificación binaria de severidad
3.1. Colapso de categorías
- Proponga y justifique un criterio para colapsar las categorías de severidad en dos grupos: sano (severidad 1) versus enfermo (resto de categorías). Discuta si esta agrupación es la más apropiada desde el punto de vista agronómico o si existiría una alternativa mejor fundamentada. £Se pierde información relevante al hacer esta reducción?
Respuesta:
Se propone colapsar las seis categorías originales en dos grupos de la siguiente manera:
Sano: severidad 1. Enfermo: severidades 2, 3, 4, 5 y 6.
El criterio se basa en la presencia o ausencia de síntomas. La categoría 1 representaría plantas sin afectación visible o con una condición considerada sana, mientras que las categorías 2 a 6 indicarían algún grado de enfermedad. Esta agrupación es útil cuando el objetivo es realizar una detección temprana y sencilla, por ejemplo, separar plantas que requieren revisión de aquellas aparentemente sanas.
Desde el punto de vista estadístico, esta reducción puede mejorar la clasificación porque elimina parte de la confusión observada entre niveles consecutivos de severidad. En el modelo multiclase, varias observaciones fueron asignadas a categorías vecinas, lo que indica que los límites entre clases no estan completamente separados. Al convertir el problema en binario, esas confusiones internas entre categorías enfermas dejan de considerarse errores.
Sin embargo, desde el punto de vista agronómico, agrupar todas las plantas enfermas en una sola categoría puede ser muy general. Una planta que tiene severidad 2 no necesariamente requiere el mismo manejo que una planta con severidad 5 o 6. Por esto, una propuesta que permita la toma de desiciones agronómicas sería definir los grupos según un umbral de intervención, por ejemplo:
Sin necesidad inmediata de intervención: severidades 1. Con necesidad de seguimiento: 2. Con necesidad de intervención: severidades 3 a 6.
Este criterio sería más útil para la toma de decisiones, siempre que exista evidencia agronómica de que a partir de la severidad 3 aumenta de manera importante el riesgo de pérdida, propagación o daño productivo.
Con respecto a la pérdida de información relevante, es muy probable que se pierda información relevante. Se deja de distinguir entre niveles leves, moderados y severos, se pierde información sobre la progresión de la enfermedad y se limita la posibilidad de relacionar la severidad con cambios espectrales, fisiológicos o productivos. Por tanto, la clasificación sano versus enfermo es apropiada para detección general, pero no necesariamente para decisiones agronómicas más específicas.
- Entrene un perceptrón multicapa para el problema binario resultante, repitiendo el esquema de partición entrenamiento/validación/prueba. Compare el desempeño con el modelo de clasificación múltiple de la Parte I en términos de las métricas comunes. ¿Mejora la clasificación al reducir el número de clases? Explique por qué o por qué no.
# ============================================================
# PREGUNTA 9. PERCEPTRÓN MULTICAPA BINARIO
# Sano = severidad 1
# Enfermo = severidades 2 a 6
# ============================================================
library(dplyr)
library(tidyr)
library(ggplot2)
# ============================================================
# 1. PREPARAR LOS DATOS BINARIOS
# ============================================================
predictores <- c(
"ndvi_med",
"evi_med",
"ndre_med",
"gli_med",
"height_med"
)
x_train_bin <- as.matrix(
train_70_std[, predictores]
)
x_val_bin <- as.matrix(
validacion_15_std[, predictores]
)
x_test_bin <- as.matrix(
test_15_std[, predictores]
)
storage.mode(x_train_bin) <- "double"
storage.mode(x_val_bin) <- "double"
storage.mode(x_test_bin) <- "double"
# 0 = sano
# 1 = enfermo
y_train_bin <- ifelse(
as.integer(as.character(train_70_std$severity)) == 1,
0,
1
)
y_val_bin <- ifelse(
as.integer(as.character(validacion_15_std$severity)) == 1,
0,
1
)
y_test_bin <- ifelse(
as.integer(as.character(test_15_std$severity)) == 1,
0,
1
)
table(y_train_bin)## y_train_bin
## 0 1
## 18 127table(y_val_bin)## y_val_bin
## 0 1
## 4 28table(y_test_bin)## y_test_bin
## 0 1
## 5 30# ============================================================
# 2. CALCULAR PESOS DE CLASE
# ============================================================
frecuencias_bin <- table(y_train_bin)
peso_sano <- length(y_train_bin) /
(
2 * frecuencias_bin["0"]
)
peso_enfermo <- length(y_train_bin) /
(
2 * frecuencias_bin["1"]
)
pesos_binarios <- c(
"0" = as.numeric(peso_sano),
"1" = as.numeric(peso_enfermo)
)
pesos_binarios## 0 1
## 4.0277778 0.5708661# ============================================================
# 3. FUNCIONES DE ACTIVACIÓN Y PÉRDIDA
# ============================================================
relu_bin <- function(x) {
resultado <- pmax(x, 0)
dim(resultado) <- dim(x)
resultado
}
relu_derivada_bin <- function(x) {
resultado <- 1 * (x > 0)
dim(resultado) <- dim(x)
resultado
}
sigmoide <- function(z) {
z <- pmax(
pmin(z, 30),
-30
)
1 / (1 + exp(-z))
}
perdida_binaria_ponderada <- function(
probabilidades,
observado,
pesos_clase
) {
epsilon <- 1e-12
pesos_observacion <- ifelse(
observado == 0,
pesos_clase["0"],
pesos_clase["1"]
)
perdidas <- -(
observado *
log(probabilidades + epsilon) +
(1 - observado) *
log(1 - probabilidades + epsilon)
)
sum(
pesos_observacion *
perdidas
) /
sum(pesos_observacion)
}# ============================================================
# 4. INICIALIZAR EL MLP BINARIO
# ============================================================
inicializar_mlp_binario <- function(
numero_entradas,
neuronas_ocultas = 8,
semilla = 123
) {
set.seed(semilla)
# Inicialización He para la capa ReLU
w1 <- matrix(
rnorm(
numero_entradas *
neuronas_ocultas,
mean = 0,
sd = sqrt(2 / numero_entradas)
),
nrow = numero_entradas,
ncol = neuronas_ocultas
)
b1 <- matrix(
0,
nrow = 1,
ncol = neuronas_ocultas
)
# Pesos de la capa binaria de salida
w2 <- matrix(
rnorm(
neuronas_ocultas,
mean = 0,
sd = sqrt(1 / neuronas_ocultas)
),
nrow = neuronas_ocultas,
ncol = 1
)
b2 <- matrix(
0,
nrow = 1,
ncol = 1
)
list(
w1 = w1,
b1 = b1,
w2 = w2,
b2 = b2
)
}# ============================================================
# 5. ENTRENAR EL MLP BINARIO
# ============================================================
entrenar_mlp_binario <- function(
x_train,
y_train,
x_val,
y_val,
neuronas_ocultas = 8,
tasa_aprendizaje = 0.1,
epocas = 300,
pesos_clase,
semilla = 123
) {
x_train <- as.matrix(x_train)
x_val <- as.matrix(x_val)
storage.mode(x_train) <- "double"
storage.mode(x_val) <- "double"
y_train <- matrix(
as.numeric(y_train),
ncol = 1
)
y_val <- matrix(
as.numeric(y_val),
ncol = 1
)
red <- inicializar_mlp_binario(
numero_entradas = ncol(x_train),
neuronas_ocultas = neuronas_ocultas,
semilla = semilla
)
historial <- data.frame(
epoca = seq_len(epocas),
loss_train = NA_real_,
loss_val = NA_real_,
accuracy_train = NA_real_,
accuracy_val = NA_real_
)
mejor_loss_val <- Inf
mejor_epoca <- NA_integer_
mejor_red <- NULL
for (epoca in seq_len(epocas)) {
# --------------------------------------------------------
# PROPAGACIÓN HACIA ADELANTE: ENTRENAMIENTO
# --------------------------------------------------------
z1 <- x_train %*% red$w1
z1 <- sweep(
z1,
MARGIN = 2,
STATS = as.numeric(red$b1),
FUN = "+"
)
a1 <- relu_bin(z1)
z2 <- a1 %*% red$w2
z2 <- sweep(
z2,
MARGIN = 2,
STATS = as.numeric(red$b2),
FUN = "+"
)
probabilidades_train <- sigmoide(z2)
historial$loss_train[epoca] <-
perdida_binaria_ponderada(
probabilidades = probabilidades_train,
observado = y_train,
pesos_clase = pesos_clase
)
pred_train <- ifelse(
probabilidades_train >= 0.5,
1,
0
)
historial$accuracy_train[epoca] <-
mean(pred_train == y_train)
# --------------------------------------------------------
# RETROPROPAGACIÓN PONDERADA
# --------------------------------------------------------
pesos_observacion <- ifelse(
y_train == 0,
pesos_clase["0"],
pesos_clase["1"]
)
delta_salida <- (
probabilidades_train -
y_train
) *
pesos_observacion /
sum(pesos_observacion)
gradiente_w2 <- t(a1) %*%
delta_salida
gradiente_b2 <- matrix(
sum(delta_salida),
nrow = 1
)
delta_oculta <- (
delta_salida %*%
t(red$w2)
) *
relu_derivada_bin(z1)
gradiente_w1 <- t(x_train) %*%
delta_oculta
gradiente_b1 <- matrix(
colSums(delta_oculta),
nrow = 1
)
# --------------------------------------------------------
# ACTUALIZAR PESOS
# --------------------------------------------------------
red$w1 <- red$w1 -
tasa_aprendizaje *
gradiente_w1
red$b1 <- red$b1 -
tasa_aprendizaje *
gradiente_b1
red$w2 <- red$w2 -
tasa_aprendizaje *
gradiente_w2
red$b2 <- red$b2 -
tasa_aprendizaje *
gradiente_b2
# --------------------------------------------------------
# VALIDACIÓN
# --------------------------------------------------------
z1_val <- x_val %*% red$w1
z1_val <- sweep(
z1_val,
MARGIN = 2,
STATS = as.numeric(red$b1),
FUN = "+"
)
a1_val <- relu_bin(z1_val)
z2_val <- a1_val %*% red$w2
z2_val <- sweep(
z2_val,
MARGIN = 2,
STATS = as.numeric(red$b2),
FUN = "+"
)
probabilidades_val <- sigmoide(
z2_val
)
historial$loss_val[epoca] <-
perdida_binaria_ponderada(
probabilidades = probabilidades_val,
observado = y_val,
pesos_clase = pesos_clase
)
pred_val <- ifelse(
probabilidades_val >= 0.5,
1,
0
)
historial$accuracy_val[epoca] <-
mean(pred_val == y_val)
# Guardar la red correspondiente a la mejor época
if (
is.finite(historial$loss_val[epoca]) &&
historial$loss_val[epoca] <
mejor_loss_val
) {
mejor_loss_val <-
historial$loss_val[epoca]
mejor_epoca <- epoca
mejor_red <- red
}
}
list(
red_final = red,
mejor_red = mejor_red,
mejor_epoca = mejor_epoca,
mejor_loss_val = mejor_loss_val,
historial = historial
)
}# ============================================================
# 6. AJUSTAR EL MODELO BINARIO
# ============================================================
modelo_binario <- entrenar_mlp_binario(
x_train = x_train_bin,
y_train = y_train_bin,
x_val = x_val_bin,
y_val = y_val_bin,
neuronas_ocultas = 8,
tasa_aprendizaje = 0.1,
epocas = 300,
pesos_clase = pesos_binarios,
semilla = 123
)
modelo_binario$mejor_epoca## [1] 191modelo_binario$mejor_loss_val## [1] 0.3265802# ============================================================
# 7. CURVAS DE APRENDIZAJE
# ============================================================
curvas_binarias <-
modelo_binario$historial |>
select(
epoca,
loss_train,
loss_val
) |>
pivot_longer(
cols = c(
loss_train,
loss_val
),
names_to = "conjunto",
values_to = "perdida"
) |>
mutate(
conjunto = recode(
conjunto,
loss_train = "Entrenamiento",
loss_val = "Validación"
)
)
ggplot(
curvas_binarias,
aes(
x = epoca,
y = perdida,
linetype = conjunto
)
) +
geom_line(
linewidth = 0.8
) +
geom_vline(
xintercept =
modelo_binario$mejor_epoca,
linetype = "dotted"
) +
labs(
title = "Curvas de aprendizaje del modelo binario",
subtitle = paste(
"La línea vertical indica la mejor época:",
modelo_binario$mejor_epoca
),
x = "Época",
y = "Entropía cruzada binaria",
linetype = "Conjunto"
) +
theme_minimal()
# ============================================================
# 8. FUNCIÓN DE PREDICCIÓN
# ============================================================
predecir_mlp_binario <- function(
red,
x,
umbral = 0.5
) {
x <- as.matrix(x)
storage.mode(x) <- "double"
z1 <- x %*% red$w1
z1 <- sweep(
z1,
MARGIN = 2,
STATS = as.numeric(red$b1),
FUN = "+"
)
a1 <- relu_bin(z1)
z2 <- a1 %*% red$w2
z2 <- sweep(
z2,
MARGIN = 2,
STATS = as.numeric(red$b2),
FUN = "+"
)
probabilidades <- as.numeric(
sigmoide(z2)
)
clases <- ifelse(
probabilidades >= umbral,
1,
0
)
list(
probabilidades = probabilidades,
clases = clases
)
}
prediccion_binaria_test <-
predecir_mlp_binario(
red = modelo_binario$mejor_red,
x = x_test_bin,
umbral = 0.5
)# ============================================================
# 9. CALCULAR MÉTRICAS BINARIAS
# ============================================================
calcular_metricas_binarias <- function(
observado,
predicho
) {
matriz <- table(
Real = factor(
observado,
levels = c(0, 1),
labels = c("Sano", "Enfermo")
),
Predicho = factor(
predicho,
levels = c(0, 1),
labels = c("Sano", "Enfermo")
)
)
metricas_clase <- data.frame()
for (clase_actual in c(0, 1)) {
vp <- sum(
observado == clase_actual &
predicho == clase_actual
)
fn <- sum(
observado == clase_actual &
predicho != clase_actual
)
fp <- sum(
observado != clase_actual &
predicho == clase_actual
)
soporte <- sum(
observado == clase_actual
)
precision <- ifelse(
vp + fp == 0,
0,
vp / (vp + fp)
)
sensibilidad <- ifelse(
vp + fn == 0,
0,
vp / (vp + fn)
)
f1 <- ifelse(
precision + sensibilidad == 0,
0,
2 *
precision *
sensibilidad /
(
precision +
sensibilidad
)
)
metricas_clase <- bind_rows(
metricas_clase,
data.frame(
clase = ifelse(
clase_actual == 0,
"Sano",
"Enfermo"
),
soporte = soporte,
precision = precision,
sensibilidad = sensibilidad,
f1 = f1
)
)
}
accuracy <- mean(
observado == predicho
)
pesos <- metricas_clase$soporte /
sum(metricas_clase$soporte)
resumen <- data.frame(
promedio = c(
"Macro",
"Ponderado"
),
precision = c(
mean(metricas_clase$precision),
sum(
metricas_clase$precision *
pesos
)
),
sensibilidad = c(
mean(metricas_clase$sensibilidad),
sum(
metricas_clase$sensibilidad *
pesos
)
),
f1 = c(
mean(metricas_clase$f1),
sum(
metricas_clase$f1 *
pesos
)
)
)
list(
matriz_confusion = matriz,
accuracy = accuracy,
metricas_clase = metricas_clase,
resumen = resumen
)
}
evaluacion_binaria <- calcular_metricas_binarias(
observado = y_test_bin,
predicho =
prediccion_binaria_test$clases
)
evaluacion_binaria$matriz_confusion## Predicho
## Real Sano Enfermo
## Sano 3 2
## Enfermo 10 20evaluacion_binaria$accuracy## [1] 0.6571429evaluacion_binaria$metricas_claseevaluacion_binaria$resumen# ============================================================
# 10. TABLAS DEL MODELO BINARIO
# ============================================================
tabla_binaria_clase <-
evaluacion_binaria$metricas_clase |>
mutate(
across(
c(
precision,
sensibilidad,
f1
),
~ round(.x, 4)
)
)
knitr::kable(
tabla_binaria_clase,
caption = paste(
"Métricas por clase del modelo",
"sano versus enfermo"
)
)| clase | soporte | precision | sensibilidad | f1 |
|---|---|---|---|---|
| Sano | 5 | 0.2308 | 0.6000 | 0.3333 |
| Enfermo | 30 | 0.9091 | 0.6667 | 0.7692 |
tabla_binaria_promedios <-
evaluacion_binaria$resumen |>
mutate(
across(
c(
precision,
sensibilidad,
f1
),
~ round(.x, 4)
)
)
knitr::kable(
tabla_binaria_promedios,
caption = paste(
"Promedios macro y ponderado",
"del modelo binario"
)
)| promedio | precision | sensibilidad | f1 |
|---|---|---|---|
| Macro | 0.5699 | 0.6333 | 0.5513 |
| Ponderado | 0.8122 | 0.6571 | 0.7070 |
# ============================================================
# 11. COMPARACIÓN BINARIO VS MULTICLASE
# ============================================================
comparacion_binario_multiclase <- data.frame(
modelo = c(
"Multiclase: seis severidades",
"Binario: sano vs. enfermo"
),
accuracy = c(
0.4285714,
evaluacion_binaria$accuracy
),
precision_macro = c(
0.2236111,
evaluacion_binaria$resumen$precision[
evaluacion_binaria$resumen$promedio ==
"Macro"
]
),
sensibilidad_macro = c(
0.3363095,
evaluacion_binaria$resumen$sensibilidad[
evaluacion_binaria$resumen$promedio ==
"Macro"
]
),
f1_macro = c(
0.2661552,
evaluacion_binaria$resumen$f1[
evaluacion_binaria$resumen$promedio ==
"Macro"
]
),
precision_ponderada = c(
0.2816667,
evaluacion_binaria$resumen$precision[
evaluacion_binaria$resumen$promedio ==
"Ponderado"
]
),
sensibilidad_ponderada = c(
0.4285714,
evaluacion_binaria$resumen$sensibilidad[
evaluacion_binaria$resumen$promedio ==
"Ponderado"
]
),
f1_ponderado = c(
0.3367775,
evaluacion_binaria$resumen$f1[
evaluacion_binaria$resumen$promedio ==
"Ponderado"
]
)
)
comparacion_binario_multiclase <-
comparacion_binario_multiclase |>
mutate(
across(
where(is.numeric),
~ round(.x, 4)
)
)
knitr::kable(
comparacion_binario_multiclase,
caption = paste(
"Comparación del modelo multiclase",
"y el modelo binario"
)
)| modelo | accuracy | precision_macro | sensibilidad_macro | f1_macro | precision_ponderada | sensibilidad_ponderada | f1_ponderado |
|---|---|---|---|---|---|---|---|
| Multiclase: seis severidades | 0.4286 | 0.2236 | 0.3363 | 0.2662 | 0.2817 | 0.4286 | 0.3368 |
| Binario: sano vs. enfermo | 0.6571 | 0.5699 | 0.6333 | 0.5513 | 0.8122 | 0.6571 | 0.7070 |
La reducción de seis categorías de severidad a dos grupos mejoró el desempeño del perceptrón multicapa. La exactitud aumentó de 0,4286 en el modelo multiclase a 0,6571 en el modelo binario. También mejoraron la precisión, la sensibilidad y el F1 macro, que pasó de 0,2662 a 0,5513.
Esta mejora se debe a que el modelo ya no tuvo que diferenciar entre niveles de severidad cercanos, que pueden presentar respuestas espectrales similares. Al reunir todas las categorías enfermas en un solo grupo, muchas confusiones dejaron de considerarse errores.
Sin embargo, el desempeño todavía fue limitado. El modelo clasificó correctamente 3 de las 5 plantas sanas y 20 de las 30 enfermas, pero generó 10 falsos negativos y mostró baja precisión para la clase sana. Además, las métricas ponderadas estuvieron influenciadas por la mayor cantidad de plantas enfermas.
En conclusión, reducir el número de clases facilitó la clasificación, pero el modelo aún tuvo dificultades para separar con precisión las plantas sanas de las enfermas. Además, esta simplificación hizo que se perdiera información sobre el nivel específico de severidad.
- Reporte para el modelo binario: exactitud, precisión, sensibilidad, especificidad, F1-score, AUC-ROC y el umbral de clasificación óptimo según el criterio de Youden o el punto más cercano a la esquina superior izquierda de la curva ROC. Interprete cada métrica en el contexto del problema fitosanitario: £qué es más costoso, un falso negativo o un falso positivo?
# ============================================================
# PUNTO 10. EVALUACIÓN COMPLETA DEL MODELO BINARIO
# ============================================================
library(dplyr)
library(ggplot2)
# Probabilidades en validación
prediccion_binaria_val <- predecir_mlp_binario(
red = modelo_binario$mejor_red,
x = x_val_bin,
umbral = 0.5
)
probabilidad_val <- prediccion_binaria_val$probabilidades
# Probabilidades en prueba
prediccion_binaria_test <- predecir_mlp_binario(
red = modelo_binario$mejor_red,
x = x_test_bin,
umbral = 0.5
)
probabilidad_test <- prediccion_binaria_test$probabilidades
head(probabilidad_val)## [1] 0.3968151 0.1800170 0.5290935 0.3823159 0.9500790 0.6622206head(probabilidad_test)## [1] 0.3368069 0.2759280 0.2064852 0.1761791 0.1025023 0.1093530# ============================================================
# FUNCIÓN DE MÉTRICAS PARA UN UMBRAL DETERMINADO
# Clase positiva: Enfermo = 1
# ============================================================
calcular_metricas_umbral <- function(
observado,
probabilidad,
umbral = 0.5
) {
predicho <- ifelse(
probabilidad >= umbral,
1,
0
)
# Enfermo = clase positiva
vp <- sum(observado == 1 & predicho == 1)
fn <- sum(observado == 1 & predicho == 0)
fp <- sum(observado == 0 & predicho == 1)
vn <- sum(observado == 0 & predicho == 0)
exactitud <- (vp + vn) /
(vp + vn + fp + fn)
precision <- ifelse(
vp + fp == 0,
0,
vp / (vp + fp)
)
sensibilidad <- ifelse(
vp + fn == 0,
0,
vp / (vp + fn)
)
especificidad <- ifelse(
vn + fp == 0,
0,
vn / (vn + fp)
)
f1 <- ifelse(
precision + sensibilidad == 0,
0,
2 * precision * sensibilidad /
(precision + sensibilidad)
)
data.frame(
umbral = umbral,
verdaderos_positivos = vp,
falsos_negativos = fn,
falsos_positivos = fp,
verdaderos_negativos = vn,
exactitud = exactitud,
precision = precision,
sensibilidad = sensibilidad,
especificidad = especificidad,
f1 = f1
)
}# ============================================================
# FUNCIÓN PARA CALCULAR ROC Y AUC
# ============================================================
calcular_roc_binaria_completa <- function(
observado,
probabilidad
) {
umbrales <- sort(
unique(
c(
1,
probabilidad,
0
)
),
decreasing = TRUE
)
tabla_roc <- bind_rows(
lapply(
umbrales,
function(u) {
resultado <- calcular_metricas_umbral(
observado = observado,
probabilidad = probabilidad,
umbral = u
)
data.frame(
umbral = u,
sensibilidad = resultado$sensibilidad,
especificidad = resultado$especificidad,
tasa_falsos_positivos =
1 - resultado$especificidad,
indice_youden =
resultado$sensibilidad +
resultado$especificidad - 1,
distancia_esquina =
sqrt(
(1 - resultado$sensibilidad)^2 +
(1 - resultado$especificidad)^2
)
)
}
)
) |>
arrange(
tasa_falsos_positivos,
sensibilidad
)
# AUC por regla trapezoidal
auc <- sum(
diff(tabla_roc$tasa_falsos_positivos) *
(
head(tabla_roc$sensibilidad, -1) +
tail(tabla_roc$sensibilidad, -1)
) / 2
)
list(
tabla_roc = tabla_roc,
auc = auc
)
}# ============================================================
# SELECCIÓN DEL UMBRAL EN VALIDACIÓN
# ============================================================
roc_validacion <- calcular_roc_binaria_completa(
observado = y_val_bin,
probabilidad = probabilidad_val
)
tabla_roc_val <- roc_validacion$tabla_roc
# Umbral óptimo según Youden
resultado_youden <- tabla_roc_val |>
slice_max(
order_by = indice_youden,
n = 1,
with_ties = FALSE
)
# Umbral más cercano a la esquina superior izquierda
resultado_esquina <- tabla_roc_val |>
slice_min(
order_by = distancia_esquina,
n = 1,
with_ties = FALSE
)
resultado_youdenresultado_esquinaumbral_youden <- resultado_youden$umbral
umbral_esquina <- resultado_esquina$umbral
umbral_youden## [1] 0.6622206umbral_esquina## [1] 0.6622206# ============================================================
# MÉTRICAS EN PRUEBA
# ============================================================
metricas_umbral_05 <- calcular_metricas_umbral(
observado = y_test_bin,
probabilidad = probabilidad_test,
umbral = 0.5
)
metricas_umbral_youden <- calcular_metricas_umbral(
observado = y_test_bin,
probabilidad = probabilidad_test,
umbral = umbral_youden
)
metricas_umbral_esquina <- calcular_metricas_umbral(
observado = y_test_bin,
probabilidad = probabilidad_test,
umbral = umbral_esquina
)
comparacion_umbrales <- bind_rows(
Convencional_0.5 = metricas_umbral_05,
Youden = metricas_umbral_youden,
Esquina_superior_izquierda =
metricas_umbral_esquina,
.id = "criterio"
) |>
mutate(
across(
c(
umbral,
exactitud,
precision,
sensibilidad,
especificidad,
f1
),
~ round(.x, 4)
)
)
comparacion_umbralesknitr::kable(
comparacion_umbrales,
caption = paste(
"Comparación de métricas del modelo binario",
"según el umbral de clasificación"
)
)| criterio | umbral | verdaderos_positivos | falsos_negativos | falsos_positivos | verdaderos_negativos | exactitud | precision | sensibilidad | especificidad | f1 |
|---|---|---|---|---|---|---|---|---|---|---|
| Convencional_0.5 | 0.5000 | 20 | 10 | 2 | 3 | 0.6571 | 0.9091 | 0.6667 | 0.6 | 0.7692 |
| Youden | 0.6622 | 18 | 12 | 0 | 5 | 0.6571 | 1.0000 | 0.6000 | 1.0 | 0.7500 |
| Esquina_superior_izquierda | 0.6622 | 18 | 12 | 0 | 5 | 0.6571 | 1.0000 | 0.6000 | 1.0 | 0.7500 |
# ============================================================
# AUC-ROC EN PRUEBA
# ============================================================
roc_prueba <- calcular_roc_binaria_completa(
observado = y_test_bin,
probabilidad = probabilidad_test
)
auc_test <- roc_prueba$auc
auc_test## [1] 0.7866667# ============================================================
# REPORTE FINAL DEL MODELO BINARIO
# ============================================================
reporte_final_binario <- data.frame(
metrica = c(
"Exactitud",
"Precisión",
"Sensibilidad",
"Especificidad",
"F1-score",
"AUC-ROC",
"Umbral óptimo de Youden"
),
valor = c(
metricas_umbral_youden$exactitud,
metricas_umbral_youden$precision,
metricas_umbral_youden$sensibilidad,
metricas_umbral_youden$especificidad,
metricas_umbral_youden$f1,
auc_test,
umbral_youden
)
) |>
mutate(
valor = round(valor, 4)
)
knitr::kable(
reporte_final_binario,
caption = paste(
"Métricas finales del modelo binario",
"en el conjunto de prueba"
)
)| metrica | valor |
|---|---|
| Exactitud | 0.6571 |
| Precisión | 1.0000 |
| Sensibilidad | 0.6000 |
| Especificidad | 1.0000 |
| F1-score | 0.7500 |
| AUC-ROC | 0.7867 |
| Umbral óptimo de Youden | 0.6622 |
# ============================================================
# CURVA ROC EN VALIDACIÓN Y UMBRAL ÓPTIMO
# ============================================================
ggplot(
tabla_roc_val,
aes(
x = tasa_falsos_positivos,
y = sensibilidad
)
) +
geom_step(
linewidth = 0.9,
direction = "hv"
) +
geom_abline(
intercept = 0,
slope = 1,
linetype = "dashed"
) +
geom_point(
data = resultado_youden,
aes(
x = tasa_falsos_positivos,
y = sensibilidad
),
size = 3
) +
annotate(
"text",
x = resultado_youden$tasa_falsos_positivos,
y = resultado_youden$sensibilidad,
label = paste0(
" Umbral = ",
round(resultado_youden$umbral, 3)
),
hjust = 0,
vjust = -0.5
) +
scale_x_continuous(
limits = c(0, 1),
breaks = seq(0, 1, 0.2)
) +
scale_y_continuous(
limits = c(0, 1),
breaks = seq(0, 1, 0.2)
) +
coord_equal() +
labs(
title = "Curva ROC del modelo binario en validación",
subtitle = paste0(
"AUC = ",
round(roc_validacion$auc, 3),
"; umbral óptimo de Youden = ",
round(umbral_youden, 3)
),
x = "Tasa de falsos positivos (1 - especificidad)",
y = "Sensibilidad"
) +
theme_minimal()
Exactitud = 0,6571. El modelo clasificó correctamente 23 de las 35
plantas evaluadas. En términos fitosanitarios, acertó en aproximadamente
el 65,7 % de los casos. Sin embargo, esta métrica por sí sola puede ser
engañosa porque había muchas más plantas enfermas que sanas. Una
exactitud moderada no garantiza que el modelo esté detectando bien la
enfermedad.
Precisión = 1,0000. De todas las plantas que el modelo clasificó como enfermas, el 100 % realmente estaba enfermo. En la práctica, esto significa que el sistema no generó falsas alarmas. Por tanto, una alerta de enfermedad emitida por el modelo sería muy confiable y no llevaría, en estos datos, a inspecciones o tratamientos innecesarios sobre plantas sanas.
Sensibilidad = 0,6000. El modelo detectó 18 de las 30 plantas enfermas. En otras palabras, identificó correctamente el 60 % de los casos de enfermedad, pero dejó pasar el 40 % restante. En un contexto fitosanitario, esta es una limitación importante, porque 12 plantas enfermas fueron clasificadas como sanas y podrían quedar sin seguimiento o tratamiento.
Especificidad = 1,0000. El modelo reconoció correctamente las cinco plantas sanas. Esto significa que no confundió plantas sanas con enfermas. Desde el punto de vista operativo, reduce costos por inspecciones, aplicaciones o controles innecesarios.
F1-score = 0,7500. Esta métrica combina precisión y sensibilidad para la clase enferma. El valor es relativamente bueno, pero debe interpretarse con cuidado: está favorecido por la precisión perfecta. Aunque las alertas fueron correctas, el modelo todavía omitió una proporción considerable de plantas enfermas. Por eso, el F1 no debe hacer olvidar la sensibilidad de solo 0,60.
AUC-ROC = 0,7867. El modelo mostró una capacidad aceptable para asignar probabilidades mayores a plantas enfermas que a plantas sanas. En términos prácticos, una planta enferma tiene alrededor de 78,7 % de probabilidad de recibir una puntuación de riesgo superior a una planta sana seleccionada al azar. Esto indica una discriminación útil, pero todavía no excelente.
Umbral óptimo de Youden = 0,6622. Una planta se clasificó como enferma solo cuando su probabilidad estimada fue igual o superior a 0,6622. Este umbral aumentó la confianza de las alertas y eliminó los falsos positivos, pero redujo la sensibilidad. Declaró enfermedad solo cuando la evidencia era fuerte, pero a cambio dejó sin detectar más plantas enfermas.
En este problema, el falso negativo es generalmente más costoso que el falso positivo. Un falso negativo ocurre cuando una planta enferma es clasificada como sana. Esto puede permitir que la enfermedad avance, aumente la fuente de inóculo y se propague a otras plantas. También puede retrasar el manejo y aumentar las pérdidas productivas.
Un falso positivo implica clasificar como enferma una planta sana. Esto puede ocasionar inspecciones, análisis o tratamientos innecesarios, pero normalmente puede corregirse mediante una segunda verificación. Por ello, en vigilancia fitosanitaria suele ser preferible aceptar algunas falsas alarmas con tal de aumentar la sensibilidad y reducir el número de plantas enfermas no detectadas.
4. Parte III: Regresión logística y comparación de modelos
4.1. Modelo de regresión logística binaria
- Ajuste un modelo de regresión logística binaria (sano vs. enfermo) usando la misma partición de datos que en la Parte II. Reporte los coeficientes estimados, sus errores estándar, los valores p y los intervalos de confianza al 95%. £Qué variables son estadísticamente significativas? £Tienen el signo esperado desde el punto de vista agronómico?
# ============================================================
# PUNTO 11. REGRESIÓN LOGÍSTICA BINARIA
# Sano = 0; Enfermo = 1
# ============================================================
library(dplyr)
library(ggplot2)
# ============================================================
# 1. PREPARAR LOS DATOS
# ============================================================
datos_logit_train <- train_70_std |>
mutate(
estado_binario = ifelse(
as.integer(as.character(severity)) == 1,
0,
1
)
)
datos_logit_val <- validacion_15_std |>
mutate(
estado_binario = ifelse(
as.integer(as.character(severity)) == 1,
0,
1
)
)
datos_logit_test <- test_15_std |>
mutate(
estado_binario = ifelse(
as.integer(as.character(severity)) == 1,
0,
1
)
)
table(datos_logit_train$estado_binario)##
## 0 1
## 18 127table(datos_logit_val$estado_binario)##
## 0 1
## 4 28table(datos_logit_test$estado_binario)##
## 0 1
## 5 30# ============================================================
# 2. AJUSTAR LA REGRESIÓN LOGÍSTICA
# ============================================================
modelo_logistico <- glm(
estado_binario ~
ndvi_med +
evi_med +
ndre_med +
gli_med +
height_med,
data = datos_logit_train,
family = binomial(
link = "logit"
)
)
summary(modelo_logistico)##
## Call:
## glm(formula = estado_binario ~ ndvi_med + evi_med + ndre_med +
## gli_med + height_med, family = binomial(link = "logit"),
## data = datos_logit_train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.6322 1.7614 2.630 0.00854 **
## ndvi_med 1.2703 1.6152 0.786 0.43159
## evi_med -1.3973 1.1715 -1.193 0.23294
## ndre_med -2.4823 1.2785 -1.942 0.05218 .
## gli_med -2.2185 3.4159 -0.649 0.51604
## height_med 0.5191 0.4993 1.040 0.29850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 108.776 on 144 degrees of freedom
## Residual deviance: 66.203 on 139 degrees of freedom
## AIC: 78.203
##
## Number of Fisher Scoring iterations: 9# ============================================================
# 3. TABLA DE COEFICIENTES
# ============================================================
tabla_coeficientes <- as.data.frame(
summary(modelo_logistico)$coefficients
)
tabla_coeficientes$variable <-
rownames(tabla_coeficientes)
rownames(tabla_coeficientes) <- NULL
names(tabla_coeficientes) <- c(
"coeficiente",
"error_estandar",
"valor_z",
"valor_p",
"variable"
)
tabla_coeficientes <- tabla_coeficientes |>
select(
variable,
coeficiente,
error_estandar,
valor_z,
valor_p
)
tabla_coeficientes# ============================================================
# 4. INTERVALOS DE CONFIANZA DEL 95 %
# ============================================================
tabla_coeficientes <- tabla_coeficientes |>
mutate(
limite_inferior_95 =
coeficiente -
1.96 * error_estandar,
limite_superior_95 =
coeficiente +
1.96 * error_estandar
)
tabla_coeficientes# ============================================================
# 5. ODDS RATIOS E INTERVALOS DEL 95 %
# ============================================================
tabla_coeficientes <- tabla_coeficientes |>
mutate(
odds_ratio =
exp(coeficiente),
or_limite_inferior =
exp(limite_inferior_95),
or_limite_superior =
exp(limite_superior_95),
significancia = case_when(
valor_p < 0.001 ~ "***",
valor_p < 0.01 ~ "**",
valor_p < 0.05 ~ "*",
valor_p < 0.10 ~ ".",
TRUE ~ "No significativa"
)
)
tabla_coeficientes# ============================================================
# 6. TABLA FINAL PARA EL INFORME
# ============================================================
tabla_logistica_final <-
tabla_coeficientes |>
mutate(
across(
c(
coeficiente,
error_estandar,
valor_z,
valor_p,
limite_inferior_95,
limite_superior_95,
odds_ratio,
or_limite_inferior,
or_limite_superior
),
~ round(.x, 4)
)
)
knitr::kable(
tabla_logistica_final,
caption = paste(
"Coeficientes del modelo de regresión",
"logística binaria"
)
)| variable | coeficiente | error_estandar | valor_z | valor_p | limite_inferior_95 | limite_superior_95 | odds_ratio | or_limite_inferior | or_limite_superior | significancia |
|---|---|---|---|---|---|---|---|---|---|---|
| (Intercept) | 4.6322 | 1.7614 | 2.6298 | 0.0085 | 1.1798 | 8.0846 | 102.7377 | 3.2536 | 3244.0643 | ** |
| ndvi_med | 1.2703 | 1.6152 | 0.7865 | 0.4316 | -1.8955 | 4.4362 | 3.5620 | 0.1502 | 84.4492 | No significativa |
| evi_med | -1.3973 | 1.1715 | -1.1928 | 0.2329 | -3.6934 | 0.8987 | 0.2473 | 0.0249 | 2.4565 | No significativa |
| ndre_med | -2.4823 | 1.2785 | -1.9416 | 0.0522 | -4.9881 | 0.0235 | 0.0836 | 0.0068 | 1.0238 | . |
| gli_med | -2.2185 | 3.4159 | -0.6495 | 0.5160 | -8.9135 | 4.4766 | 0.1088 | 0.0001 | 87.9340 | No significativa |
| height_med | 0.5191 | 0.4993 | 1.0396 | 0.2985 | -0.4596 | 1.4978 | 1.6806 | 0.6316 | 4.4721 | No significativa |
# ============================================================
# 7. VARIABLES ESTADÍSTICAMENTE SIGNIFICATIVAS
# ============================================================
variables_significativas <-
tabla_coeficientes |>
filter(
variable != "(Intercept)",
valor_p < 0.05
)
variables_significativasvariables_marginales <-
tabla_coeficientes |>
filter(
variable != "(Intercept)",
valor_p >= 0.05,
valor_p < 0.10
)
variables_marginales# ============================================================
# 8. GRÁFICO DE COEFICIENTES
# ============================================================
coeficientes_grafico <-
tabla_coeficientes |>
filter(
variable != "(Intercept)"
)
ggplot(
coeficientes_grafico,
aes(
x = reorder(
variable,
coeficiente
),
y = coeficiente
)
) +
geom_point(
size = 2.5
) +
geom_errorbar(
aes(
ymin = limite_inferior_95,
ymax = limite_superior_95
),
width = 0.15
) +
geom_hline(
yintercept = 0,
linetype = "dashed"
) +
coord_flip() +
labs(
title = "Coeficientes de la regresión logística",
subtitle = "Intervalos de confianza del 95 %",
x = "Variable",
y = "Coeficiente estimado"
) +
theme_minimal()
# ============================================================
# 9. PROBABILIDADES EN VALIDACIÓN Y PRUEBA
# ============================================================
prob_logit_val <- predict(
modelo_logistico,
newdata = datos_logit_val,
type = "response"
)
prob_logit_test <- predict(
modelo_logistico,
newdata = datos_logit_test,
type = "response"
)
head(prob_logit_val)## 1 2 3 4 5 6
## 0.6843228 0.5673108 0.8926264 0.9200772 0.9999808 0.9346459head(prob_logit_test)## 1 2 3 4 5 6
## 0.8316897 0.7891204 0.6239937 0.5257479 0.2456433 0.2887990pred_logit_test <- ifelse(
prob_logit_test >= 0.5,
1,
0
)
table(
Real = datos_logit_test$estado_binario,
Predicho = pred_logit_test
)## Predicho
## Real 0 1
## 0 1 4
## 1 1 29mean(
pred_logit_test ==
datos_logit_test$estado_binario
)## [1] 0.8571429Con un nivel de significancia de 5 %, ninguna de las cinco variables predictoras fue estadísticamente significativa de manera individual.
Los valores p fueron:
Variable Coeficiente Valor p Resultado NDVI 1,2703 0,4316 No significativa EVI −1,3973 0,2329 No significativa NDRE −2,4823 0,0522 No significativa al 5 %, marginal al 10 % GLI −2,2185 0,5160 No significativa Altura 0,5191 0,2985 No significativa
El intercepto sí fue significativo, con p=0,0085, pero no corresponde a una variable agronómica. Representa el nivel basal de las posibilidades de enfermedad cuando todos los predictores estandarizados toman el valor cero, es decir, cuando se encuentran en sus valores medios.
La variable más próxima a la significancia fue NDRE, con p=0,0522. Aunque no alcanza el criterio convencional de 0,05, puede describirse como una asociación marginal o tendencia estadística al nivel de 10 %. Su intervalo de confianza fue aproximadamente de −4,99 a 0,02, por lo que todavía incluye ligeramente el cero. Su odds ratio fue 0,0836, lo que sugiere una disminución fuerte de las posibilidades de enfermedad al aumentar NDRE, aunque la incertidumbre es elevada.
EVI, NDRE y GLI presentaron signos negativos, lo cual coincide con lo esperado agronómicamente. Valores más altos de estos índices suelen asociarse con mayor verdor, contenido de clorofila, cobertura vegetal y vigor. Por tanto, sería razonable que plantas con valores mayores presentaran menor probabilidad de pertenecer al grupo enfermo.
NDRE mostró el patrón agronómico más coherente y la evidencia estadística más fuerte. Su coeficiente de −2,4823 indica que, manteniendo constantes los demás predictores, un aumento de una desviación estándar en NDRE se asoció con una reducción de las posibilidades de enfermedad. Sin embargo, al ser p=0,0522, este efecto debe presentarse como tendencia y no como resultado significativo definitivo.
NDVI presentó un coeficiente positivo de 1,2703, contrario a la expectativa inicial. En principio se esperaría que un NDVI mayor se relacionara con una menor probabilidad de enfermedad. Este signo contrario no implica necesariamente que un mayor NDVI cause enfermedad. Puede originarse porque NDVI está correlacionado con EVI, NDRE y GLI. En una regresión múltiple, el coeficiente de NDVI representa su efecto después de mantener constantes los demás índices, y la multicolinealidad puede cambiar la magnitud e incluso el signo de los coeficientes.
La altura también tuvo un signo positivo, contrario a la expectativa de que las plantas enfermas fueran más bajas. Sin embargo, la altura puede depender de otros factores, como edad, vigor inicial, estado fenológico, variedad, competencia o condiciones ambientales. Además, el coeficiente no fue significativo, por lo que no debe interpretarse como evidencia de que una mayor altura aumente la enfermedad.
- Aplique al menos un método de selección de variables sobre el modelo de regresión logística: puede usar selección paso a paso (stepwise), regularización Lasso (L1) o Ridge (L2), o criterios de información (AIC/BIC). £Cuáles son las variables más relevantes para predecir la severidad según este modelo? Compare este resultado con la importancia de variables obtenida con la red neuronal en la Pregunta 12.
# ============================================================
# PUNTO 12. SELECCIÓN DE VARIABLES EN REGRESIÓN LOGÍSTICA
# Método stepwise basado en AIC
# ============================================================
# Modelo completo previamente ajustado:
# modelo_logistico
# Modelo nulo: solamente intercepto
modelo_nulo <- glm(
estado_binario ~ 1,
data = datos_logit_train,
family = binomial(link = "logit")
)
# Modelo completo: cinco predictores
modelo_completo <- glm(
estado_binario ~
ndvi_med +
evi_med +
ndre_med +
gli_med +
height_med,
data = datos_logit_train,
family = binomial(link = "logit")
)
# Selección paso a paso en ambas direcciones
modelo_step_aic <- step(
object = modelo_nulo,
scope = list(
lower = formula(modelo_nulo),
upper = formula(modelo_completo)
),
direction = "both",
trace = TRUE
)## Start: AIC=110.78
## estado_binario ~ 1
##
## Df Deviance AIC
## + ndre_med 1 71.076 75.076
## + evi_med 1 72.555 76.555
## + ndvi_med 1 76.040 80.040
## + height_med 1 99.034 103.034
## + gli_med 1 99.308 103.308
## <none> 108.776 110.776
##
## Step: AIC=75.08
## estado_binario ~ ndre_med
##
## Df Deviance AIC
## + evi_med 1 67.715 73.715
## <none> 71.076 75.076
## + gli_med 1 69.362 75.362
## + height_med 1 70.859 76.859
## + ndvi_med 1 71.067 77.067
## - ndre_med 1 108.776 110.776
##
## Step: AIC=73.71
## estado_binario ~ ndre_med + evi_med
##
## Df Deviance AIC
## <none> 67.715 73.715
## + height_med 1 66.789 74.789
## - evi_med 1 71.076 75.076
## + ndvi_med 1 67.527 75.527
## + gli_med 1 67.684 75.684
## - ndre_med 1 72.555 76.555summary(modelo_step_aic)##
## Call:
## glm(formula = estado_binario ~ ndre_med + evi_med, family = binomial(link = "logit"),
## data = datos_logit_train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.8822 0.8481 4.578 4.7e-06 ***
## ndre_med -1.4403 0.7735 -1.862 0.0626 .
## evi_med -1.5678 1.0077 -1.556 0.1198
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 108.776 on 144 degrees of freedom
## Residual deviance: 67.715 on 142 degrees of freedom
## AIC: 73.715
##
## Number of Fisher Scoring iterations: 7formula(modelo_step_aic)## estado_binario ~ ndre_med + evi_medAIC(modelo_nulo, modelo_completo, modelo_step_aic)El procedimiento comenzó comparando diferentes combinaciones de predictores y seleccionó finalmente el modelo:
estado binario∼NDRE+EVI
Esto significa que, según el criterio AIC, las variables ndre_med y evi_med ofrecen el mejor equilibrio entre capacidad de ajuste y simplicidad. Las variables ndvi_med, gli_med y height_med fueron eliminadas porque no disminuyeron suficientemente el AIC al incorporarlas.
variables_step_aic <- attr(
terms(modelo_step_aic),
"term.labels"
)
variables_step_aic## [1] "ndre_med" "evi_med"# ============================================================
# COEFICIENTES DEL MODELO SELECCIONADO
# ============================================================
tabla_step <- as.data.frame(
summary(modelo_step_aic)$coefficients
)
tabla_step$variable <- rownames(tabla_step)
rownames(tabla_step) <- NULL
names(tabla_step)[1:4] <- c(
"coeficiente",
"error_estandar",
"valor_z",
"valor_p"
)
tabla_step <- tabla_step |>
dplyr::select(
variable,
coeficiente,
error_estandar,
valor_z,
valor_p
) |>
dplyr::mutate(
limite_inferior_95 =
coeficiente - 1.96 * error_estandar,
limite_superior_95 =
coeficiente + 1.96 * error_estandar,
odds_ratio = exp(coeficiente),
or_limite_inferior =
exp(limite_inferior_95),
or_limite_superior =
exp(limite_superior_95)
)
tabla_steptabla_step_presentacion <- tabla_step |>
dplyr::mutate(
dplyr::across(
where(is.numeric),
~ round(.x, 4)
)
)
knitr::kable(
tabla_step_presentacion,
caption = "Modelo logístico seleccionado mediante stepwise-AIC"
)| variable | coeficiente | error_estandar | valor_z | valor_p | limite_inferior_95 | limite_superior_95 | odds_ratio | or_limite_inferior | or_limite_superior |
|---|---|---|---|---|---|---|---|---|---|
| (Intercept) | 3.8822 | 0.8481 | 4.5777 | 0.0000 | 2.2200 | 5.5444 | 48.5298 | 9.2072 | 255.7936 |
| ndre_med | -1.4403 | 0.7735 | -1.8620 | 0.0626 | -2.9565 | 0.0758 | 0.2369 | 0.0520 | 1.0788 |
| evi_med | -1.5678 | 1.0077 | -1.5558 | 0.1198 | -3.5430 | 0.4073 | 0.2085 | 0.0289 | 1.5028 |
comparacion_aic <- data.frame(
modelo = c(
"Modelo nulo",
"Modelo completo",
"Modelo stepwise-AIC"
),
numero_predictores = c(
0,
length(attr(terms(modelo_completo), "term.labels")),
length(attr(terms(modelo_step_aic), "term.labels"))
),
AIC = c(
AIC(modelo_nulo),
AIC(modelo_completo),
AIC(modelo_step_aic)
),
BIC = c(
BIC(modelo_nulo),
BIC(modelo_completo),
BIC(modelo_step_aic)
),
devianza_residual = c(
deviance(modelo_nulo),
deviance(modelo_completo),
deviance(modelo_step_aic)
)
) |>
dplyr::mutate(
dplyr::across(
where(is.numeric),
~ round(.x, 4)
)
)
knitr::kable(
comparacion_aic,
caption = "Comparación de los modelos logísticos"
)| modelo | numero_predictores | AIC | BIC | devianza_residual |
|---|---|---|---|---|
| Modelo nulo | 0 | 110.7759 | 113.7526 | 108.7759 |
| Modelo completo | 5 | 78.2029 | 96.0633 | 66.2029 |
| Modelo stepwise-AIC | 2 | 73.7146 | 82.6448 | 67.7146 |
El modelo stepwise obtuvo el menor AIC. La diferencia frente al modelo completo fue:
78,2029−73,7146=4,4883
Por tanto, eliminar tres variables mejora el equilibrio entre ajuste y complejidad. El modelo completo presentó una devianza residual ligeramente menor, porque contiene más variables, pero fue penalizado por su mayor complejidad.
El BIC también favoreció claramente al modelo seleccionado:
modelo completo: 96,0633; modelo stepwise: 82,6448.
Como el BIC penaliza con mayor intensidad el número de parámetros, este resultado respalda aún más la selección de un modelo sencillo con NDRE y EVI.
# Probabilidades en prueba
prob_step_test <- predict(
modelo_step_aic,
newdata = datos_logit_test,
type = "response"
)
# Clasificación con umbral 0.5
pred_step_test <- ifelse(
prob_step_test >= 0.5,
1,
0
)
# Matriz de confusión
matriz_step <- table(
Real = datos_logit_test$estado_binario,
Predicho = pred_step_test
)
matriz_step## Predicho
## Real 0 1
## 0 2 3
## 1 1 29# Exactitud
accuracy_step <- mean(
pred_step_test ==
datos_logit_test$estado_binario
)
accuracy_step## [1] 0.8857143metricas_step_test <- calcular_metricas_umbral(
observado = datos_logit_test$estado_binario,
probabilidad = prob_step_test,
umbral = 0.5
)
metricas_step_testEste resultado indica una capacidad muy alta para detectar plantas enfermas. Solo una planta enferma fue clasificada como sana. Desde el punto de vista fitosanitario, esto es favorable porque reduce los falsos negativos.
Sin embargo, la especificidad fue baja: solo dos de las cinco plantas sanas fueron reconocidas correctamente. Tres plantas sanas fueron clasificadas como enfermas. Por tanto, el modelo privilegió la detección de enfermedad a costa de producir falsas alarmas.
La exactitud elevada también está influida por el desbalance: 30 de las 35 observaciones de prueba estaban enfermas. Por eso no debe interpretarse aisladamente.
importancia_logistica <- tabla_step |>
dplyr::filter(
variable != "(Intercept)"
) |>
dplyr::mutate(
importancia_z = abs(valor_z),
magnitud_coeficiente = abs(coeficiente)
) |>
dplyr::arrange(
dplyr::desc(importancia_z)
)
importancia_logisticaNinguna de las dos variables fue significativa al nivel convencional del 5 %. NDRE quedó cercana a la significancia y puede considerarse una asociación marginal al nivel del 10 %.
Los dos coeficientes fueron negativos. Como la categoría positiva es “enfermo”, esto indica que aumentos de NDRE o EVI se asocian con una reducción de la probabilidad de enfermedad.
Para NDRE, el odds ratio de 0,2369 indica que un aumento de una desviación estándar multiplica las posibilidades de enfermedad por aproximadamente 0,237, manteniendo constante EVI. Esto equivale a una reducción aproximada de cerca de 76,3 % en las posibilidades estimadas de enfermedad. Sin embargo, el intervalo de confianza del odds ratio incluye 1, por lo que la asociación no es estadísticamente concluyente.
NDRE está relacionado con la reflectancia en el borde rojo y suele ser sensible al contenido de clorofila y al estado fisiológico del dosel. Una planta con mayor NDRE normalmente presenta mayor actividad fotosintética o mejor condición del tejido vegetal. Por ello, resulta razonable que valores mayores se asocien con menor probabilidad de enfermedad.
EVI también representa vigor y actividad de la vegetación, con cierta corrección frente a efectos del suelo y saturación. Su signo negativo indica que plantas con mayor vigor espectral tendieron a presentar menor probabilidad de pertenecer al grupo enfermo.
La interpretación debe ser asociativa y no causal. El modelo no demuestra que aumentar estos índices reduzca directamente la enfermedad; muestra que ambos índices permiten diferenciar parcialmente plantas sanas y enfermas.
ggplot(
importancia_logistica,
aes(
x = reorder(variable, importancia_z),
y = importancia_z
)
) +
geom_col() +
coord_flip() +
labs(
title = "Importancia de variables en la regresión logística",
subtitle = "Magnitud absoluta del estadístico z",
x = "Variable",
y = "|z|"
) +
theme_minimal()
NDRE fue la variable con mayor evidencia relativa dentro del modelo
seleccionado. EVI ocupó el segundo lugar. Esto no significa que NDRE sea
significativamente superior a EVI, sino que presentó una relación más
fuerte respecto a su error estándar.
# ============================================================
# IMPORTANCIA POR PERMUTACIÓN EN LA RED NEURONAL BINARIA
# ============================================================
calcular_auc_simple <- function(
observado,
probabilidad
) {
positivos <- probabilidad[observado == 1]
negativos <- probabilidad[observado == 0]
comparaciones <- outer(
positivos,
negativos,
FUN = "-"
)
mean(comparaciones > 0) +
0.5 * mean(comparaciones == 0)
}auc_original_red <- calcular_auc_simple(
observado = y_test_bin,
probabilidad =
prediccion_binaria_test$probabilidades
)
auc_original_red## [1] 0.7866667set.seed(123)
numero_repeticiones <- 100
importancia_red <- data.frame()
for (variable_actual in predictores) {
disminuciones_auc <- numeric(
numero_repeticiones
)
for (repeticion in seq_len(
numero_repeticiones
)) {
x_permutado <- x_test_bin
columna_actual <- which(
predictores == variable_actual
)
x_permutado[, columna_actual] <-
sample(
x_permutado[, columna_actual],
replace = FALSE
)
pred_permutada <- predecir_mlp_binario(
red = modelo_binario$mejor_red,
x = x_permutado,
umbral = 0.5
)
auc_permutada <- calcular_auc_simple(
observado = y_test_bin,
probabilidad =
pred_permutada$probabilidades
)
disminuciones_auc[repeticion] <-
auc_original_red - auc_permutada
}
importancia_red <- dplyr::bind_rows(
importancia_red,
data.frame(
variable = variable_actual,
disminucion_auc_media =
mean(disminuciones_auc),
desviacion_estandar =
sd(disminuciones_auc)
)
)
}
importancia_red <- importancia_red |>
dplyr::arrange(
dplyr::desc(disminucion_auc_media)
)
importancia_redggplot(
importancia_red,
aes(
x = reorder(
variable,
disminucion_auc_media
),
y = disminucion_auc_media
)
) +
geom_col() +
geom_errorbar(
aes(
ymin =
disminucion_auc_media -
desviacion_estandar,
ymax =
disminucion_auc_media +
desviacion_estandar
),
width = 0.15
) +
coord_flip() +
labs(
title = "Importancia por permutación en la red neuronal",
subtitle = "Reducción media del AUC al desordenar cada variable",
x = "Variable",
y = "Disminución media del AUC"
) +
theme_minimal()
Cuando NDRE fue permutado, el AUC disminuyó en promedio 0,1404. Partiendo de 0,7867, esto implicaría un AUC aproximado de: 0,6463
Por tanto, destruir la información de NDRE perjudicó considerablemente la capacidad de la red para distinguir plantas sanas y enfermas. Esto confirma que fue la variable más importante. Por su parte, EVI produjo una disminución mucho menor, de 0,0245. Su contribución fue positiva, pero bastante inferior a NDRE.
La altura tuvo una reducción cercana a cero, indicando una contribución predictiva limitada.
NDVI y GLI presentaron disminuciones medias negativas. Esto no significa necesariamente que sean perjudiciales biológicamente. Significa que, al permutarlas, el AUC no disminuyó y en promedio mejoró ligeramente. Esto puede ocurrir por ruido, redundancia con otros índices, correlación entre variables y el pequeño tamaño del conjunto de prueba.
Las barras de error muestran una variabilidad amplia, especialmente para NDRE. Aunque su importancia media fue la mayor, existe incertidumbre porque solo se utilizaron 35 observaciones de prueba, incluidas únicamente cinco plantas sanas.
ranking_logistica <- importancia_logistica |>
dplyr::transmute(
variable,
importancia_logistica = importancia_z,
rango_logistica =
rank(
-importancia_z,
ties.method = "min"
)
)
ranking_red <- importancia_red |>
dplyr::transmute(
variable,
importancia_red =
disminucion_auc_media,
rango_red =
rank(
-disminucion_auc_media,
ties.method = "min"
)
)
comparacion_importancia <- dplyr::full_join(
ranking_logistica,
ranking_red,
by = "variable"
) |>
dplyr::arrange(rango_logistica)
comparacion_importanciaDentro del modelo logístico seleccionado, NDRE fue la variable más relevante, porque presentó la mayor magnitud del estadístico z, con ∣z∣=1,8620, seguida por EVI con ∣z∣=1,5558. Sus coeficientes fueron negativos, −1,4403 para NDRE y −1,5678 para EVI. Como la categoría positiva fue “enfermo”, estos signos indican que valores mayores de NDRE y EVI se asociaron con una menor probabilidad de enfermedad. Agronómicamente, este resultado es coherente, ya que ambos índices se relacionan con el vigor, el contenido de clorofila y el estado fisiológico de la vegetación.
comparacion_importancia_presentacion <-
comparacion_importancia |>
dplyr::mutate(
dplyr::across(
c(
importancia_logistica,
importancia_red
),
~ round(.x, 4)
)
)
knitr::kable(
comparacion_importancia_presentacion,
caption = paste(
"Comparación de la importancia de variables",
"entre regresión logística y red neuronal"
)
)| variable | importancia_logistica | rango_logistica | importancia_red | rango_red |
|---|---|---|---|---|
| ndre_med | 1.8620 | 1 | 0.1404 | 1 |
| evi_med | 1.5558 | 2 | 0.0245 | 2 |
| height_med | NA | NA | 0.0088 | 3 |
| ndvi_med | NA | NA | -0.0098 | 4 |
| gli_med | NA | NA | -0.0181 | 5 |
Las variables más relevantes para predecir la severidad fueron NDRE y EVI. La selección stepwise basada en AIC retuvo únicamente estos dos predictores y redujo el AIC del modelo completo de 78,2029 a 73,7146. NDRE ocupó el primer lugar según la magnitud absoluta del estadístico z, seguido por EVI. Ambos coeficientes fueron negativos, lo que indicó que valores mayores de estos índices se asociaron con una menor probabilidad de enfermedad, en concordancia con su relación con la clorofila, el vigor y la actividad fotosintética. La red neuronal produjo el mismo orden de importancia: permutar NDRE redujo el AUC en 0,1404, mientras que permutar EVI lo redujo en 0,0245. Altura, NDVI y GLI mostraron aportes bajos o nulos. En consecuencia, ambos modelos coincidieron en identificar a NDRE como el predictor principal y a EVI como una variable complementaria.