Este documento presenta los resultados de la predicción de la clorofila-a, el principal pigmento utilizado por las cianobacterias para realizar la fotosíntesis. El objetivo es, a partir de sus valores pasados desarrollar un modelo (S)ARIMA que pueda ser útil en un sistema de alerta temprana de proliferación de las cianobacterias.

library(readr)
library(dplyr)
library(lubridate)
library(forecast)
library(tidyr)
library(ggplot2)

Se ha cargado la base de datos con registros cada 15 minutos, se ha seleccionado la clorofila-a como variable objectivo

########Cargar los datos##############

datos <- read_delim("1_bbdd_SEM.csv", delim = ";", escape_double = FALSE, trim_ws = TRUE)

#Conferir Datos
head(datos)
## # A tibble: 6 × 22
##      ID `quarter-hourly`  Hour   day month  year chl_presa chl_playa cond_emb
##   <dbl>            <dbl> <dbl> <dbl> <dbl> <dbl>     <dbl>     <dbl>    <dbl>
## 1     1                9     2    27    10  2017      4.64      3.38     65.4
## 2     2               10     2    27    10  2017      5.42      3.51     64.8
## 3     3               11     2    27    10  2017      5.6       3.43     64.9
## 4     4               12     3    27    10  2017      5.55      4.11     65.0
## 5     5               13     3    27    10  2017      5.01      3.32     65.1
## 6     6               14     3    27    10  2017      3.66      3.44     64.9
## # ℹ 13 more variables: ph_emb <dbl>, t_agua_emb <dbl>, nh4_saica <dbl>,
## #   po4_saica <dbl>, cond_saica <dbl>, do_saica <dbl>, ph_saica <dbl>,
## #   t_agua_saica <dbl>, turb_saica <dbl>, nivel_saica <dbl>, flujo_saica <dbl>,
## #   precipicacion <dbl>, t_aire <dbl>
fecha <- as.POSIXct(
  paste(datos$year, datos$month, datos$day, datos$Hour, datos$`quarter-hourly`%%4*15, sep = "-"),
  format = "%Y-%m-%d-%H-%M"
)


#Variables utilizadas

chlorophy <-  datos$chl_playa

#crear base de datos

datos_limpios <- tibble(
  fecha,
  chlorophy) 


# Convertir la base de datos en diaria
datos_limpios <- datos_limpios %>%
  mutate(fecha_dia = floor_date(fecha, unit = "day"))  # convertir fecha a solo día


# Calcular promedio diario 
datos_diarios <- datos_limpios %>%
  group_by(fecha_dia) %>%
  summarise(chlorophy_dia = mean(chlorophy, na.rm = TRUE))

head(datos_diarios)
## # A tibble: 6 × 2
##   fecha_dia           chlorophy_dia
##   <dttm>                      <dbl>
## 1 2017-10-27 00:00:00          7.72
## 2 2017-10-28 00:00:00          7.36
## 3 2017-10-29 00:00:00          9.95
## 4 2017-10-30 00:00:00         12.6 
## 5 2017-10-31 00:00:00         12.5 
## 6 2017-11-01 00:00:00         11.5

En un primer intento, utilizando la funcion auto.arima no se ha logrado resultados satisfactorios. R2 muy bajo.

df_modelo <- datos_diarios %>%
  select(chlorophy_dia) %>%
  na.omit()

##Separar datos de entrenamiento y test
set.seed(123)
n <- nrow(df_modelo)
train_idx <- 1:floor(0.8 * n)
train <- df_modelo[train_idx, ]
test <- df_modelo[-train_idx, ]

# Convertir a serie temporal diaria

clorofila_ts <- ts(train$chlorophy_dia)


# Modelo SARIMAX
modelo_sarima <- auto.arima(clorofila_ts, seasonal = TRUE)

# Predicción y evaluación
pred_sarima <- forecast(modelo_sarima, h = nrow(test))$mean
real <- test$chlorophy_dia

MAE_sarima <- mean(abs(pred_sarima - real))
RMSE_sarima <- sqrt(mean((pred_sarima - real)^2))
R2_sarima <- cor(pred_sarima, real)^2
rango <- max(real, na.rm = TRUE) - min(real, na.rm = TRUE)
MAE_norm <- MAE_sarima / rango
RMSE_norm <- RMSE_sarima / rango

cat("SARIMAX:\n")
## SARIMAX:
cat("MAE:", MAE_sarima, "\n")
## MAE: 4.106056
cat("RMSE:", RMSE_sarima, "\n")
## RMSE: 4.882003
cat("R²:", R2_sarima, "\n")  
## R²: 0.0001438328
cat("MAE normalizado :", round(MAE_norm, 4), "\n")
## MAE normalizado : 0.1262
cat("RMSE normalizado :", round(RMSE_norm, 4), "\n")
## RMSE normalizado : 0.15

A continuación, se incluyen los valores previos de clorofila-a como variables exógenas para alimentar el modelo con el propósito de optimizar el rendimiento predictivo del modelo.

# utilizar valores pasados de clorofila (lag 1 a 3)
for (k in 1:3) {
  datos_diarios[[paste0("chlorophy_dia_lag", k)]] <- dplyr::lag(datos_diarios$chlorophy_dia, k)
}



# Selección de columnas
cols_xreg <- c(
  paste0("chlorophy_dia_lag", 1:3)
 
)

df_modelo <- datos_diarios %>%
  select(chlorophy_dia, all_of(cols_xreg)) %>%
  na.omit()

##Separar datos de entrenamiento y test
set.seed(123)
n <- nrow(df_modelo)
train_idx <- 1:floor(0.8 * n)
train <- df_modelo[train_idx, ]
test <- df_modelo[-train_idx, ]

# Convertir a serie temporal diaria

clorofila_ts <- ts(train$chlorophy_dia)

# Regresores externos
xreg_train <- as.matrix(train[, -1])
xreg_test  <- as.matrix(test[, -1])

# Modelo (S)ARIMA
modelo_sarima <- auto.arima(clorofila_ts,xreg = xreg_train, seasonal = TRUE)
summary(modelo_sarima)
## Series: clorofila_ts 
## Regression with ARIMA(1,0,3) errors 
## 
## Coefficients:
##           ar1     ma1     ma2     ma3  intercept  chlorophy_dia_lag1
##       -0.2948  0.3297  0.4107  0.2694     0.4023              1.0223
## s.e.   0.3145  0.1633  0.1400  0.1044     0.1613              0.3055
##       chlorophy_dia_lag2  chlorophy_dia_lag3
##                  -0.3314              0.2684
## s.e.              0.3752              0.1449
## 
## sigma^2 = 2.275:  log likelihood = -1694.36
## AIC=3406.72   AICc=3406.92   BIC=3450.22
## 
## Training set error measures:
##                        ME     RMSE       MAE       MPE     MAPE      MASE
## Training set -0.001248426 1.501754 0.9604489 -2.145457 10.58895 0.9794633
##                       ACF1
## Training set -0.0009660661
# Predicción y evaluación
pred_sarima <- forecast(modelo_sarima, xreg = xreg_test, h = nrow(test))$mean
real <- test$chlorophy_dia

MAE_sarima <- mean(abs(pred_sarima - real))
RMSE_sarima <- sqrt(mean((pred_sarima - real)^2))
R2_sarima <- cor(pred_sarima, real)^2
rango <- max(real, na.rm = TRUE) - min(real, na.rm = TRUE)
MAE_norm <- MAE_sarima / rango
RMSE_norm <- RMSE_sarima / rango

cat("SARIMAX:\n")
## SARIMAX:
cat("MAE:", MAE_sarima, "\n")
## MAE: 0.9273332
cat("RMSE:", RMSE_sarima, "\n")
## RMSE: 1.862999
cat("R²:", R2_sarima, "\n")  
## R²: 0.8204025
cat("MAE normalizado :", round(MAE_norm, 4), "\n")
## MAE normalizado : 0.0285
cat("RMSE normalizado :", round(RMSE_norm, 4), "\n")
## RMSE normalizado : 0.0573

Inicialmente, la función auto.arima identificó como modelo óptimo un SARIMA de tipo ARIMA(1,0,3). No obstante, al incorporar como variables externas los valores rezagados de clorofila-a (lags 1, 2 y 3), que reflejan su carácter autorregresivo, el modelo más adecuado resultó ser un ARIMA(3,0,3). En este último, el componente autorregresivo de orden 3 en la parte no estacional captura explícitamente la influencia de los valores pasados de clorofila-a sobre la predicción actual. El modelo obtenido presenta un coeficiente de determinación (R²) de 0,82.

Se analizó su capacidad predictiva de dos formas: cualitativamente, mediante gráficos comparativos entre los valores observados y predichos para una evaluación visual; y cuantitativamente, mediante una matriz de confusión, con el fin de evaluar la eficacia del modelo en el contexto de un sistema de alertas tempranas.

tiempo <- seq_along(real)

plot(tiempo, real, type = "l", col = "blue", lwd = 2,
     ylab = "Clorofila-a (µg/L)", xlab = "Tiempo (días)", main = "ARIMA vs Real",ylim = c(0, 40))
lines(tiempo, pred_sarima, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("Real", "Predicción"), col = c("blue", "red"), lty = c(1,2), lwd = 2, cex = 0.7)

fechas <- datos_diarios$fecha_dia

# Crear un vector con todas las observaciones reales
clorofila_real <- datos_diarios$chlorophy_dia

# Extender predicciones SARIMAX para ubicarlas en la parte final
pred_extendida <- rep(NA, length(clorofila_real))
n_pred <- length(pred_sarima)
pred_extendida[(length(clorofila_real) - n_pred + 1):length(clorofila_real)] <- pred_sarima

# Graficar
plot(fechas, clorofila_real, type = "l", col = "blue", lwd = 1,
     ylab = "Clorofila-a (µg/L)", xlab = "Año", main = "SARIMAX vs Real", ylim = c(0, 50))
lines(fechas, pred_extendida, col = "red", lwd = 1, lty = 2)

legend("topleft", legend = c("Real", "Predicción"),
       col = c("blue", "red"), lty = c(1, 2), lwd = 2, cex = 0.7)
abline(h = 10, col = "black", lty = 1)

##################################MATRIZ DE CONFUSION###########################

# Clasificación binaria
pred_bin <- ifelse(pred_sarima > 10, 1, 0)
real_bin <- ifelse(real > 10, 1, 0)

# Matriz de confusión
TP <- sum(pred_bin == 1 & real_bin == 1)
FP <- sum(pred_bin == 1 & real_bin == 0)
TN <- sum(pred_bin == 0 & real_bin == 0)
FN <- sum(pred_bin == 0 & real_bin == 1)

# Abrir nueva ventana gráfica
plot(NA, xlim = c(0, 2), ylim = c(0, 2), type = "n", axes = FALSE, xlab = "VALORES REALES", ylab = "VALORES PREDICCIÓN", main = "Matriz de Confusión")

# Dibujar rectángulos con colores
rect(0, 1, 1, 2, col = "olivedrab", border = "black")  # TP
rect(1, 1, 2, 2, col = "lightcoral", border = "black")    # FP
rect(0, 0, 1, 1, col = "lightcoral", border = "black")    # FN
rect(1, 0, 2, 1, col = "olivedrab", border = "black")  # TN

# Añadir textos en cada celda
text(0.5, 1.5, paste("TP =", TP), cex = 1)
text(1.5, 1.5, paste("FP =", FP), cex = 1)
text(0.5, 0.5, paste("FN =", FN), cex = 1)
text(1.5, 0.5, paste("TN =", TN), cex = 1)

# Etiquetas personalizadas del eje X (en el medio debajo de las columnas)
text(x = c(0.5, 1.5), y = -0.1, labels = c("Positivo", "Negativo"), xpd = TRUE)

# Etiquetas personalizadas del eje Y (a la izquierda de las filas, rotadas)
text(x = -0.1, y = c(0.5, 1.5), labels = c("Negativo", "Positivo"), xpd = TRUE, srt = 90) 

# Métricas
precision    <- TP / (TP + FP)
sensibilidad <- TP / (TP + FN)
exactitud    <- (TP + TN) / (TP + TN + FP + FN)
f1_score     <- 2 * (precision * sensibilidad) / (precision + sensibilidad)


cat("Precisión:    ", round(precision, 4), "\n")
## Precisión:     0.75
cat("Sensibilidad: ", round(sensibilidad, 4), "\n")
## Sensibilidad:  0.75
cat("Exactitud:    ", round(exactitud, 4), "\n")
## Exactitud:     0.9571
cat("Puntuación F1:", round(f1_score, 4), "\n")
## Puntuación F1: 0.75

Entre las métricas derivadas de la matriz de confusión, la puntuación F1 es la considerada de mayor robustez y relevancia en este trabajo, por su capacidad de balancear la precisión y la sensibilidad, con un valor de 0,75.

El modelo predictivo permite la implementación de un sistema de alarmas automatizado, capaz de activar de forma inmediata protocolos de intervención ante el riesgo de floraciones de cianobacterias.

Cabe señalar que su capacidad predictiva podría mejorarse incorporando otras variables exógenas que influyen en la dinámica de las cianobacterias.