1. Carga de Datos

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, union
library(ggplot2)

datos <- read.csv("china_water_pollution_data.csv", header = TRUE, sep = ",", dec = ".")

datos_sel <- datos %>%
  select(BOD_mg_L, COD_mg_L)

2. Gráfica

# Gráfica inicial con todos los puntos
plot(
  datos_sel$BOD_mg_L,
  datos_sel$COD_mg_L,
  pch = 19,
  col = "lightblue",
  xlab = "BOD (mg/L)",
  ylab = "COD (mg/L)",
  main = "Relación entre BOD y COD"
)

Esta gráfica muestra todas las observaciones, incluyendo valores atípicos y NA.Se observa dispersión y posibles puntos extremos que podrían distorsionar la regresión.

3. Limpieza de datos y eliminación de valores atípicos

Se eliminan valores faltantes y atípicos para evitar que distorsionen la estimación de la regresión. La limpieza mejora la representación de la tendencia general, manteniendo la variabilidad natural.

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

# Aplicar a ambas variables y eliminar NA
datos_limpios <- datos_sel[
  quitar_atipicos_IQR(datos_sel$BOD_mg_L) &
  quitar_atipicos_IQR(datos_sel$COD_mg_L),
]

4. Agrupación por intervalos (Bins) y resumen

# Crear bins de BOD
by_bin <- 0.591
datos_binned <- datos_limpios %>%
  mutate(
    BOD_BIN = cut(
      BOD_mg_L,
      breaks = seq(
        floor(min(BOD_mg_L)),
        ceiling(max(BOD_mg_L)),
        by = by_bin
      )
    )
  )

# Medianas por bin
datos_resumen <- datos_binned %>%
  group_by(BOD_BIN) %>%
  summarise(
    BOD_mediana = median(BOD_mg_L, na.rm = TRUE),
    COD_mediana = median(COD_mg_L, na.rm = TRUE),
    n = n()
  ) %>%
  filter(n >= 3)

datos_resumen
## # A tibble: 10 × 4
##    BOD_BIN     BOD_mediana COD_mediana     n
##    <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

5. Definir Variable

x <- datos_resumen$BOD_mediana
y <- datos_resumen$COD_mediana

6. Conjetura

plot(x, y,
     pch = 19,
     xlab = "BOD (mg/L)",
     ylab = "COD (mg/L)",
     main = "Gráfica Nº1: Diagrama de dispersión\nBOD y COD")

7. Parámetros del modelo

# Ajuste log-log: log(COD) = log(a) + b * log(BOD)
modelo_pot <- lm(log(y) ~ log(x))

# Parámetros del modelo
a <- exp(coef(modelo_pot)[1])
b <- coef(modelo_pot)[2]

a  # Intercepto en escala original
## (Intercept) 
##    18.36036
b  # Exponente
##     log(x) 
## 0.06159927

7.1 Ecuación del modelo

plot(1, type = "n", axes = FALSE, xlab = "", ylab = "")

text(
  x = 1, y = 1,
  labels = paste0(
    "Ecuación Potencial\n",
    "y = a · x^b\n",
    "y = ", round(a, 4), " · x^", round(b, 5)
  ),
  cex = 2,
  col = "blue",
  font = 2
)

7.2 Gráfico con recta de regresión

plot(x, y,
     col = "deepskyblue",
     pch = 16,
     xlab = "BOD (mg/L)",
     ylab = "COD (mg/L)",
     main = "Gráfica Nº2: Ajuste potencial BOD y COD")

curve(a * x^b,
      from = min(x),
      to = max(x),
      add = TRUE,
      col = "blue",
      lwd = 2)

8. Correlación y coeficiente de determinación

y_hat <- a * x^b

r <- cor(y, y_hat, method = "pearson")
r2 <- r^2 * 100

cat("Coeficiente de correlación de Pearson:", round(r, 4), "\n")
## Coeficiente de correlación de Pearson: 0.7764
cat("Coeficiente de determinación (R²) =", round(r2, 2), "%\n")
## Coeficiente de determinación (R²) = 60.28 %

9. Restricciones del modelo

cat("Restricciones del modelo potencial:\n")
## Restricciones del modelo potencial:
cat("- Dominio: x > 0 (BOD positiva)\n")
## - Dominio: x > 0 (BOD positiva)
cat("- Rango: y > 0 (COD positiva)\n")
## - Rango: y > 0 (COD positiva)
cat("- Aplicación limitada al rango observado:\n")
## - Aplicación limitada al rango observado:
cat("Rango experimental de BOD: ", range(x), "\n")
## Rango experimental de BOD:  1.465 6.425
cat("Rango experimental de COD: ", range(y), "\n")
## Rango experimental de COD:  18.695 21.345

10. Estimación

# Estimar COD para BOD = 4 mg/L
BOD_objetivo <- 4
COD_est <- a * BOD_objetivo^b

plot(1, type = "n", axes = FALSE, xlab = "", ylab = "")
text(
  1, 1,
  labels = paste0(
    "Estimación de COD para BOD = ", BOD_objetivo, " mg/L\n",
    "Resultado estimado: ", round(COD_est, 3), " mg/L"
  ),
  cex = 1.3,
  col = "blue",
  font = 2
)

11. Conclusiones

En la contaminación del agua en China existe una relación potencial entre BOD y COD. El modelo ajustado es: 𝑦 = 𝑎 ⋅ 𝑥 𝑏 y=a⋅x b . Se recomienda usarlo solo dentro del rango observado: BOD ∈ [1.46, 6.42] mg/L; COD ∈ [18.7, 21.3] mg/L. Para BOD = 4 mg/L, la concentración estimada de COD es aproximadamente 20 mg/L.