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.06159927# REPRESENTACIÓN DE LA ECUACIÓN POTENCIAL
plot(1, type = "n", axes = FALSE, xlab = "", ylab = "")
text(
x = 1, y = 1.1,
labels = "Ecuación Potencial",
col = "blue",
cex = 2,
font = 2
)
text(
x = 1, y = 1,
labels = "y = a · x^b",
col = "blue",
cex = 1.8,
font = 2
)
text(
x = 1, y = 0.85,
labels = paste0(
"y = ",
round(a_pot, 4),
" · x^",
round(b_pot, 5)
),
col = "blue",
cex = 1.8,
font = 2
)
5.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 (COD)\n",
"cuando la (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
En el intervalo entre la Demanda Bioquímica de Oxígeno (BOD) y el carbono organico disuelto (COD) en los cuerpos de agua monitoreados existe una relación de tipo potencial, representada por y ̂=a⋅〖“BOD 〗^b, siendo la BOD la medida de la carga orgánica biodegradable del agua (mg/L) y y ̂la concentración estimada de COD en mg/L; dentro del intervalo seleccionado no existen restricciones, ya que el modelo potencial está definido para valores positivos de la variable independiente, y se espera que, cuando la BOD sea 4 mg/L, se obtenga una concentración de COD de aproximadamente 20 mg/L, de acuerdo con el modelo ajustado.