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 y de otras variables que influyen en su comportamiento, desarrollar un modelo SARIMAX 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 han seleccionado las variables de interés y los datos se han transformado a frecuencia diaria.

########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
temp <-  datos$t_agua_emb
cond   <- datos$cond_emb
ph <- datos$ph_emb


#crear base de datos

datos_limpios <- tibble(
  fecha,
  chlorophy,
  temp,
  cond,
  ph  ) 


# 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),
            temp_dia = mean(temp, na.rm = TRUE),
            cond_dia = mean(cond, na.rm = TRUE),
            ph_dia = mean(ph, na.rm = TRUE))

head(datos_diarios)
## # A tibble: 6 × 5
##   fecha_dia           chlorophy_dia temp_dia cond_dia ph_dia
##   <dttm>                      <dbl>    <dbl>    <dbl>  <dbl>
## 1 2017-10-27 00:00:00          7.72     16.2     65.2   8.11
## 2 2017-10-28 00:00:00          7.36     16.1     66.1   7.92
## 3 2017-10-29 00:00:00          9.95     16.0     66.2   7.76
## 4 2017-10-30 00:00:00         12.6      15.6     66.3   7.58
## 5 2017-10-31 00:00:00         12.5      15.4     66.5   7.67
## 6 2017-11-01 00:00:00         11.5      15.1     66.5   7.72

A continuación, se analizó el retraso temporal (lag) de las variables con mayor impacto en la concentración de clorofila-a.

# Seleccionar las variables
vars <- datos_diarios[, c("chlorophy_dia", "ph_dia", "temp_dia", "cond_dia")]

# Calcular matriz de correlación
correlacion <- cor(vars, use = "complete.obs")

# Mostrar solo las correlaciones de clorofila con las otras
correlacion["chlorophy_dia", ]
## chlorophy_dia        ph_dia      temp_dia      cond_dia 
##    1.00000000    0.46990205    0.08540501    0.05680200
# Definir máximo lag
max_lag <- 100

# Inicializar lista para resultados
cor_lags <- data.frame(
  lag = 0:max_lag,
  cor_ph = NA,
  cor_temp = NA,
  cor_cond = NA
)

# Calcular correlaciones para cada lag
for (k in 0:max_lag) {
  ph_lag   <- dplyr::lag(datos_diarios$ph_dia, k)
  temp_lag <- dplyr::lag(datos_diarios$temp_dia, k)
  cond_lag <- dplyr::lag(datos_diarios$cond_dia, k)
  
  cor_lags$cor_ph[k + 1]   <- cor(datos_diarios$chlorophy_dia, ph_lag, use = "complete.obs")
  cor_lags$cor_temp[k + 1] <- cor(datos_diarios$chlorophy_dia, temp_lag, use = "complete.obs")
  cor_lags$cor_cond[k + 1] <- cor(datos_diarios$chlorophy_dia, cond_lag, use = "complete.obs")
 
}



# Dibujar correlaciones de las variables en el tiempo con la clorofila
cor_lags_long <- cor_lags %>% 
  pivot_longer(-lag,
               names_to  = "variable",
               values_to = "cor")


ggplot(cor_lags_long, aes(x = lag, y = cor, colour = variable)) +
  geom_line(size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_colour_manual(values = c(cor_ph   = "red",
                                 cor_temp = "blue",
                                 cor_cond = "cyan"),
                      labels = c( "CE","pH","Temperatura"),
                      name   = "Variable") +
  labs(title = "Correlación de clorofila con lags",
       x     = "Lag (días)",
       y     = "Coeficiente de correlación") +
  theme_minimal(base_size = 13)

En un primer intento, utilizando los retrasos temporales identificados, los resultados no fueron satisfactorios. Por ello, se ajustaron manualmente dichos retrasos con el objetivo de encontrar una combinación que mejorara las métricas de evaluación (R², MAE y RMSE). Además, se forzó a la función auto.arima a incluir un retraso (lag) de 3 dias de la clorofila-a, debido a su comportamiento autorregresivo.

# 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)
}

# utilizar valor actual y 4 anteriores de ph (lag 0 a 4)
for (k in 0:4) {
  datos_diarios[[paste0("ph_dia_lag", k)]] <- dplyr::lag(datos_diarios$ph_dia, k)
}

# utilizar valores pasados de temperatura (lag 35 a 45)
for (k in 35:45) {
  datos_diarios[[paste0("temp_dia_lag", k)]] <- dplyr::lag(datos_diarios$temp_dia, k)
}

# utilizar valores pasados de conductividad (lag 10 a 20)
for (k in 10:20) {
  datos_diarios[[paste0("cond_dia_lag", k)]] <- lag(datos_diarios$cond_dia, k)
}



# Selección de columnas
cols_xreg <- c(
  paste0("ph_dia_lag", 0:4),
  paste0("cond_dia_lag", 10:20),
  paste0("chlorophy_dia_lag", 1:3),
  paste0("temp_dia_lag", 35:45)
)

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, frequency = 7)

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

# Modelo SARIMAX
modelo_sarimax <- auto.arima(clorofila_ts, xreg = xreg_train, seasonal = TRUE)
summary(modelo_sarimax)
## Series: clorofila_ts 
## Regression with ARIMA(0,0,1)(1,0,2)[7] errors 
## 
## Coefficients:
##           ma1     sar1    sma1    sma2  ph_dia_lag0  ph_dia_lag1  ph_dia_lag2
##       -0.7038  -0.6035  0.6156  0.1180       2.8611      -5.8119       3.3885
## s.e.   0.0902   0.1446  0.1432  0.0397       0.2469       0.5985       0.7389
##       ph_dia_lag3  ph_dia_lag4  cond_dia_lag10  cond_dia_lag11  cond_dia_lag12
##           -0.7534       0.3081         -0.0169          0.0411         -0.0437
## s.e.       0.5778       0.2501          0.0155          0.0333          0.0376
##       cond_dia_lag13  cond_dia_lag14  cond_dia_lag15  cond_dia_lag16
##               0.0204          0.0227         -0.0048         -0.0445
## s.e.          0.0376          0.0376          0.0375          0.0373
##       cond_dia_lag17  cond_dia_lag18  cond_dia_lag19  cond_dia_lag20
##               0.0510         -0.0176         -0.0308          0.0241
## s.e.          0.0375          0.0374          0.0331          0.0153
##       chlorophy_dia_lag1  chlorophy_dia_lag2  chlorophy_dia_lag3
##                   1.6466             -0.5056             -0.1525
## s.e.              0.0932              0.1174              0.0368
##       temp_dia_lag35  temp_dia_lag36  temp_dia_lag37  temp_dia_lag38
##              -0.1152          0.0374          0.0852          0.1315
## s.e.          0.1345          0.3077          0.3545          0.3566
##       temp_dia_lag39  temp_dia_lag40  temp_dia_lag41  temp_dia_lag42
##               0.0960         -0.0743         -0.2519         -0.2296
## s.e.          0.3582          0.3583          0.3572          0.3597
##       temp_dia_lag43  temp_dia_lag44  temp_dia_lag45
##               0.4914         -0.0778         -0.0850
## s.e.          0.3594          0.3200          0.1437
## 
## sigma^2 = 1.893:  log likelihood = -1538.71
## AIC=3147.43   AICc=3150.36   BIC=3315.32
## 
## Training set error measures:
##                        ME     RMSE       MAE       MPE    MAPE      MASE
## Training set 0.0005267796 1.349503 0.8974238 -1.788738 11.1783 0.2752012
##                      ACF1
## Training set -0.004908915

Inicialmente, la función auto.arima identificó como modelo óptimo el SARIMA (0,0,1)(1,0,2)[7]. Sin embargo, al incorporar como variable externa la clorofila-a con un retraso de tres días (lags 1, 2 y 3), lo cual refleja su naturaleza autorregresiva, se consideró más adecuado el modelo SARIMAX (3,0,1)(1,0,2)[7]. En este caso, el componente autorregresivo de orden 3 en la parte no estacional representa explícitamente la influencia de los valores pasados de clorofila-a en la predicción actual. La validez del modelo se ha comprobado mediante el análisis de los residuos utilizando el conjunto de entrenamiento, así como a través de las métricas R², MAE y RMSE evaluadas sobre el conjunto de test.

residuos <- residuals(modelo_sarimax)
# Gráfico de residuos a lo largo del tiempo que comproba que no hay tendencias
plot(residuos, main = "Residuos del modelo SARIMAX", ylab = "Residuos", col = "gray")
abline(h = 0, col = "red", lty = 2)

acf(residuos, main = "ACF de residuos") #comproba que no hay correlacion de los residuos

pacf(residuos, main = "PACF de residuos") #comproba que no hay correlacion de los residuos

qqnorm(residuos); qqline(residuos, col = "red") #comproba que los residuos siguen una distribucion normal

# Predicción y evaluación
pred_sarimax <- forecast(modelo_sarimax, xreg = xreg_test, h = nrow(test))$mean
real <- test$chlorophy_dia

MAE_sarimax <- mean(abs(pred_sarimax - real))
RMSE_sarimax <- sqrt(mean((pred_sarimax - real)^2))
R2_sarimax <- cor(pred_sarimax, real)^2
rango <- max(real, na.rm = TRUE) - min(real, na.rm = TRUE)
MAE_norm <- MAE_sarimax / rango
RMSE_norm <- RMSE_sarimax / rango

cat("SARIMAX:\n")
## SARIMAX:
cat("MAE:", MAE_sarimax, "\n")
## MAE: 0.9405274
cat("RMSE:", RMSE_sarimax, "\n")
## RMSE: 1.95553
cat("R²:", R2_sarimax, "\n")  
## R²: 0.8324701
cat("MAE normalizado :", round(MAE_norm, 4), "\n")
## MAE normalizado : 0.0289
cat("RMSE normalizado :", round(RMSE_norm, 4), "\n")
## RMSE normalizado : 0.0601

Una vez validado el modelo, 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 = "SARIMAX vs Real",ylim = c(0, 40))
lines(tiempo, pred_sarimax, 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_sarimax)
pred_extendida[(length(clorofila_real) - n_pred + 1):length(clorofila_real)] <- pred_sarimax

# 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_sarimax > 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.8636
cat("Sensibilidad: ", round(sensibilidad, 4), "\n")
## Sensibilidad:  0.95
cat("Exactitud:    ", round(exactitud, 4), "\n")
## Exactitud:     0.9821
cat("Puntuación F1:", round(f1_score, 4), "\n")
## Puntuación F1: 0.9048

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,9.

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.