Regresión potencial entre (COD) y (BOD) en cuerpos de agua de China (2023)
Edgar Saavedra/ Grupo 1
2026-02-08
# UNIVERSIDAD CENTRAL DEL ECUADOR
# Facultad de Ingeniería en Geología, Minas, Petroleos y Ambiental
# Ingeniería Ambiental
# Autor: Grupo 1
# Fecha: 07/02/2026
library(readr)
library(dplyr)##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(tidyr)1. Carga de datos
datos <- read.csv(
"china_water_pollution_data.csv",
header = TRUE,
sep = ",",
dec = ".",
fileEncoding = "UTF-8"
)
str(datos)## 'data.frame': 3000 obs. of 25 variables:
## $ Province : chr "Zhejiang" "Sichuan" "Zhejiang" "Beijing" ...
## $ City : chr "Ningbo" "Mianyang" "Ningbo" "Beijing" ...
## $ Monitoring_Station : chr "Ningbo_Station_2" "Mianyang_Station_1" "Ningbo_Station_8" "Beijing_Station_10" ...
## $ Latitude : num 25.5 32.2 30 30 43.5 ...
## $ Longitude : num 123 113 125 118 122 ...
## $ Date : chr "2023-06-01" "2023-03-05" "2023-07-13" "2023-02-17" ...
## $ Water_Temperature_C : num 22.5 27.3 21 16.6 21.8 ...
## $ pH : num 6.93 6.89 6.02 7.31 7.77 7.16 7.2 7.38 7.77 7.17 ...
## $ Dissolved_Oxygen_mg_L : num 9.3 8.14 5.34 10.06 7.93 ...
## $ Conductivity_uS_cm : num 652 358 520 593 656 ...
## $ Turbidity_NTU : num 0.85 4.49 17.46 7.38 3.7 ...
## $ Nitrate_mg_L : num 2.14 2.06 2.11 1.9 1.8 2.69 1.99 2.7 1.82 2.11 ...
## $ Nitrite_mg_L : num 0.03 0.015 0.029 0.014 0.019 0.022 0.003 0.034 0.026 0.028 ...
## $ Ammonia_N_mg_L : num 0.38 0.38 0.3 0.2 0.22 0.44 0.29 0.32 0.71 0.53 ...
## $ Total_Phosphorus_mg_L : num 0.074 0.147 0.021 0.155 0.152 0.134 0.05 0.104 0.126 0.039 ...
## $ Total_Nitrogen_mg_L : num 2.71 3.15 3.39 2.91 3.45 2.87 3.05 2.75 3.68 3.52 ...
## $ COD_mg_L : num 15.4 16.8 17.3 17.9 20.1 ...
## $ BOD_mg_L : num 1.39 2.98 2.65 5.18 3.47 5.27 3.34 5.71 3.68 4.73 ...
## $ Heavy_Metals_Pb_ug_L : num 6.9 4.68 3.24 3.2 2.01 4.42 5.01 6.75 4.59 3.09 ...
## $ Heavy_Metals_Cd_ug_L : num 0.66 0.39 0.27 0.67 0.34 -0.03 0.6 0.5 0.21 0.35 ...
## $ Heavy_Metals_Hg_ug_L : num 0.02 0.1 0.11 0.11 0.14 0.12 0.09 0.08 0.16 0.04 ...
## $ Coliform_Count_CFU_100mL: int 87 116 110 99 82 91 92 91 95 108 ...
## $ Water_Quality_Index : num 36.6 66.2 98.7 71.3 16.1 ...
## $ Pollution_Level : chr "Very Poor" "Excellent" "Poor" "Poor" ...
## $ Remarks : chr "High pollution spike detected" "High pollution spike detected" "High pollution spike detected" "Monitoring recommended" ...2. Depuración de datos: eliminación de atípicos (IQR) y agrupación por bins.
datos$COD_mg_L <- as.numeric(datos$COD_mg_L)
datos$BOD_mg_L <- as.numeric(datos$BOD_mg_L)
datos_depu <- datos %>%
select(BOD_mg_L, COD_mg_L) %>%
na.omit() %>%
filter(BOD_mg_L > 0, COD_mg_L > 0)
cat("Filtrado NA/positivos:", nrow(datos), "→", nrow(datos_depu), "filas\n")## Filtrado NA/positivos: 3000 → 3000 filasremove_iqr_outliers <- function(x) {
Q1 <- quantile(x, 0.25, na.rm = TRUE)
Q3 <- quantile(x, 0.75, na.rm = TRUE)
IQR <- Q3 - Q1
x >= Q1 - 1.5 * IQR & x <= Q3 + 1.5 * IQR
}
mask_bod <- remove_iqr_outliers(datos_depu$BOD_mg_L)
mask_cod <- remove_iqr_outliers(datos_depu$COD_mg_L)
datos_depu <- datos_depu[mask_bod & mask_cod, ]
cat("Filtrado IQR:", nrow(datos_depu), "filas limpias tras IQR\n")## Filtrado IQR: 2966 filas limpias tras IQRby_bin <- 0.591
datos_depu <- datos_depu %>%
mutate(BOD_BIN = cut(
BOD_mg_L,
breaks = seq(floor(min(BOD_mg_L)),
ceiling(max(BOD_mg_L)),
by = by_bin),
include.lowest = TRUE
))
# Medianas por bin
med_bin <- datos_depu %>%
group_by(BOD_BIN) %>%
summarise(
BOD_median = median(BOD_mg_L, na.rm = TRUE),
COD_median = median(COD_mg_L, na.rm = TRUE),
n_bin = n()
) %>%
drop_na()
med_bin## # A tibble: 10 × 4
## BOD_BIN BOD_median COD_median n_bin
## <fct> <dbl> <dbl> <int>
## 1 [1,1.59] 1.46 18.7 16
## 2 (1.59,2.18] 2.00 18.9 70
## 3 (2.18,2.77] 2.54 19.8 246
## 4 (2.77,3.36] 3.13 20.0 447
## 5 (3.36,3.96] 3.69 20.4 621
## 6 (3.96,4.55] 4.22 20.0 662
## 7 (4.55,5.14] 4.82 19.8 495
## 8 (5.14,5.73] 5.34 19.5 276
## 9 (5.73,6.32] 5.94 20.3 107
## 10 (6.32,6.91] 6.42 21.3 26# Mantener positivos
med_bin_pot <- med_bin %>%
filter(BOD_median > 0, COD_median > 0)
cat("Bins usados:", nrow(med_bin_pot), "\n")## Bins usados: 10by_bin <- 0.591
datos_depu <- datos_depu %>%
mutate(BOD_BIN = cut(
BOD_mg_L,
breaks = seq(floor(min(BOD_mg_L)),
ceiling(max(BOD_mg_L)),
by = by_bin),
include.lowest = TRUE
))
med_bin <- datos_depu %>%
group_by(BOD_BIN) %>%
summarise(
BOD_median = median(BOD_mg_L, na.rm = TRUE),
COD_median = median(COD_mg_L, na.rm = TRUE),
n_bin = n()
) %>%
drop_na()
med_bin_pot <- med_bin %>% filter(BOD_median > 0, COD_median > 0)
cat("Bins usados:", nrow(med_bin_pot), "\n")## Bins usados: 10med_bin_pot## # A tibble: 10 × 4
## BOD_BIN BOD_median COD_median n_bin
## <fct> <dbl> <dbl> <int>
## 1 [1,1.59] 1.46 18.7 16
## 2 (1.59,2.18] 2.00 18.9 70
## 3 (2.18,2.77] 2.54 19.8 246
## 4 (2.77,3.36] 3.13 20.0 447
## 5 (3.36,3.96] 3.69 20.4 621
## 6 (3.96,4.55] 4.22 20.0 662
## 7 (4.55,5.14] 4.82 19.8 495
## 8 (5.14,5.73] 5.34 19.5 276
## 9 (5.73,6.32] 5.94 20.3 107
## 10 (6.32,6.91] 6.42 21.3 263. Definir variables de análisis
x <- med_bin_pot$BOD_median
y <- med_bin_pot$COD_median4. Conjetura
# A partir del diagrama de dispersión se observa una
# relación no lineal entre la Demanda Bioquímica de Oxígeno (BOD)
# y la Demanda Química de Oxígeno (COD).
# La tendencia muestra un crecimiento suave de la COD conforme
# aumenta la BOD, lo que sugiere el uso de un modelo potencial
# del tipo y = a * x^b.
# Variables para la conjetura (medianas por bin)
x <- med_bin_pot$BOD_median
y <- med_bin_pot$COD_median
plot(
x, y,
pch = 19,
col = "steelblue",
xlab = "BOD (mediana por bin) [mg/L]",
ylab = "COD (mediana por bin) [mg/L]",
main = "Gráfica N°1: Relación entre la BOD y la COD",
font.main = 2
)
5. Parámetros
# Ajuste potencial en escala log-log:
# log(COD) = log(a) + b*log(BOD)
modelo_pot <- lm(log(COD_median) ~ log(BOD_median), data = med_bin_pot)
summary(modelo_pot)##
## Call:
## lm(formula = log(COD_median) ~ log(BOD_median), data = med_bin_pot)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.043940 -0.013812 -0.003557 0.016481 0.036037
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.91019 0.02354 123.633 2.05e-14 ***
## log(BOD_median) 0.06160 0.01731 3.558 0.00742 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02535 on 8 degrees of freedom
## Multiple R-squared: 0.6128, Adjusted R-squared: 0.5644
## F-statistic: 12.66 on 1 and 8 DF, p-value: 0.007418# Parámetros del modelo potencial: COD = a * BOD^b
coef_pot <- coef(modelo_pot)
a_pot <- exp(coef_pot[1]) # a = e^(intercepto)
b_pot <- coef_pot[2] # b = pendiente
a_pot## (Intercept)
## 18.36036b_pot## log(BOD_median)
## 0.061599275.1 Gráfica del modelo ajustado
# Curva ajustada
x_vals <- seq(min(med_bin_pot$BOD_median),
max(med_bin_pot$BOD_median),
length.out = 200)
y_pred <- a_pot * x_vals^b_pot
# Gráfica (puntos + curva)
plot(med_bin_pot$BOD_median, med_bin_pot$COD_median,
pch = 19, col = "steelblue",
xlab = "BOD (mediana por bin) [mg/L]",
ylab = "COD (mediana por bin) [mg/L]",
main = "Gráfica N°2: Modelo potencial ajustado entre BOD y COD\nChina, 2023",
font.main = 2)
lines(x_vals, y_pred, col = "darkgreen", lwd = 2)
legend("topleft",
legend = c("Medianas por bin",
paste0("COD = ", round(a_pot,4), " * BOD^", round(b_pot,4))),
col = c("steelblue", "darkgreen"),
pch = c(19, NA),
lwd = c(NA, 2),
bty = "n")
6. TEST
6.1 Correlación de Pearson
# Correlación de Pearson sobre variables transformadas (log-log)
r <- cor(
log(med_bin_pot$BOD_median),
log(med_bin_pot$COD_median),
method = "pearson"
)
r## [1] 0.7828295cat("Correlación de Pearson (%):", round(r * 100, 2), "%\n")## Correlación de Pearson (%): 78.28 %6.2 Coeficiente de determinación
r2 <- summary(modelo_pot)$r.squared * 100
cat("R² potencial =", round(r2, 2), "%\n")## R² potencial = 61.28 %7. Restricciones del modelo
## 7. Restricciones del modelo# Restricción 1: el modelo potencial requiere valores positivos (por el log)
# (en tu caso ya filtraste BOD_median > 0 y COD_median > 0, pero lo verificamos)
sum(med_bin_pot$BOD_median <= 0)## [1] 0sum(med_bin_pot$COD_median <= 0)## [1] 0# Restricción 2: el modelo es válido solo dentro del rango observado de BOD
rango_BOD <- range(med_bin_pot$BOD_median)
rango_COD <- range(med_bin_pot$COD_median)
cat("Rango válido de BOD (mg/L): [", round(rango_BOD[1], 2), ",", round(rango_BOD[2], 2), "]\n")## Rango válido de BOD (mg/L): [ 1.46 , 6.42 ]cat("Rango observado de COD (mg/L): [", round(rango_COD[1], 2), ",", round(rango_COD[2], 2), "]\n")## Rango observado de COD (mg/L): [ 18.7 , 21.34 ]## numeric(0)
# El modelo potencial COD = a * BOD^b solo está definido para valores positivos
# de la Demanda Bioquímica de Oxígeno (BOD > 0), debido a la transformación
# logarítmica utilizada en el ajuste del modelo.
# Asimismo, el modelo es válido únicamente dentro del rango observado de la BOD
# en los datos analizados y no debe extrapolarse fuera de dicho intervalo.8. Estimación
# Valor objetivo (puedes cambiarlo)
BOD_objetivo <- 4 # mg/L
# Estimación usando el modelo potencial: COD = a * BOD^b
COD_est <- a_pot * BOD_objetivo^b_pot
# Mostrar como en el ejemplo (texto centrado grande)
plot(1, type = "n", axes = FALSE, xlab = "", ylab = "")
text(
1, 1,
labels = paste0(
"¿Cuál es la concentración esperada de DQO (COD)\n",
"cuando la DBO (BOD) del agua es ", BOD_objetivo, " mg/L?\n\n",
"Resultado estimado: ", round(COD_est, 3), " mg/L"
),
cex = 1.3,
col = "blue",
font = 2
)
9. Conclusión
Utilizando las medianas agrupadas por intervalos de la Demanda Bioquímica de Oxígeno (BOD), y transformando ambas variables a escala logarítmica (log-log), se modeló la relación no lineal entre la BOD y la Demanda Química de Oxígeno (COD) mediante un modelo de tipo potencial, común en el análisis de procesos de contaminación orgánica en cuerpos de agua.
El ajuste obtuvo un coeficiente de determinación R² ≈ 61 % y una correlación de Pearson log-log de 78 %, lo que indica una relación positiva y moderada entre la BOD y la COD. En el contexto del estudio de la contaminación del agua en China en 2023, el modelo potencial evidencia una relación escalable entre ambas variables; sin embargo, la COD no depende exclusivamente de la BOD, ya que intervienen otros factores fisicoquímicos y aportes externos.