Borrador Modelo Logarítmico

#CARGA DATOS
datos <- read.csv2("china_water_pollution_data.csv", 
                   stringsAsFactors = FALSE, 
                   fileEncoding = "UTF-8")

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)

# 1. PREPARACIÓN DE VARIABLES
# Convertimos de texto con comas a números válidos
Turbidez <- as.numeric(gsub(",", ".", datos$Turbidity_NTU))
Indice_Calidad <- as.numeric(gsub(",", ".", datos$Water_Quality_Index))

datos$Turbidez <- Turbidez
datos$Indice_Calidad <- Indice_Calidad

# 2. LIMPIEZA DE DATOS
# Filtramos NAs y valores > 0 (necesario para el logaritmo)
datos_limpios <- subset(datos, !is.na(Turbidez) & !is.na(Indice_Calidad) &   
                          Turbidez > 0 & Indice_Calidad > 0)

# 3. AGRUPACIÓN POR PROMEDIOS
datos_limpios_prom <- datos_limpios %>%
  group_by(Turbidez) %>%
  summarise(Indice_Calidad_promedio = mean(Indice_Calidad, na.rm = TRUE)) %>%
  ungroup()

# 4. VERIFICACIÓN DE LONGITUDES
cat("Longitud Turbidez (Causa):", length(Turbidez), "\n")

## Longitud Turbidez (Causa): 3000

cat("Longitud Indice_Calidad (Efecto):", length(Indice_Calidad), "\n")

## Longitud Indice_Calidad (Efecto): 3000

cat("Longitud Diarrea (Validación):", length(Turbidez), "\n")

## Longitud Diarrea (Validación): 3000

# 5. AJUSTE DEL MODELO LOGARÍTMICO (lm con log)
modelo_log_prom <- lm(Indice_Calidad_promedio ~ log(Turbidez), data = datos_limpios_prom)

# Coeficientes
beta0 <- coef(modelo_log_prom)[1]
beta1 <- coef(modelo_log_prom)[2]

# Mostrar ecuación en consola
cat("Ecuación del modelo logarítmico:\n")

## Ecuación del modelo logarítmico:

cat("WQI_promedio =", round(beta0, 4), "+", round(beta1, 4), "* ln(Turbidez)\n")

## WQI_promedio = 51.045 + -0.4998 * ln(Turbidez)

# 6. FILTRADO POR RESIDUALES (Umbral del 1%)
datos_limpios_prom$pred <- predict(modelo_log_prom, newdata = datos_limpios_prom)
datos_limpios_prom <- datos_limpios_prom %>%
  mutate(residual = abs(Indice_Calidad_promedio - pred))

umbral <- 0.01 * (max(datos_limpios_prom$Indice_Calidad_promedio) - 
                    min(datos_limpios_prom$Indice_Calidad_promedio))
datos_filtrados <- filter(datos_limpios_prom, residual <= umbral)

# 7. VISUALIZACIÓN DE LA ECUACIÓN (Conjetura)
plot(1, type = "n", axes = FALSE, xlab = "", ylab = "") 
text(x = 1, y = 1,
     labels = paste0("Ecuación logarítmica\n Y = a + b*ln(x)\n Y = ", 
                     round(beta0, 4), " + (", round(beta1, 4), ")*ln(x)"),
     cex = 1.5, col = "blue", font = 2)

# 8. GRÁFICA DE REGRESIÓN FINAL
plot(datos_filtrados$Turbidez, datos_filtrados$Indice_Calidad_promedio,
     pch = 19, col = "forestgreen",
     main = "Regresión logarítmica: WQI en función de Turbidez",
     xlab = "Turbidez (NTU)", ylab = "Índice de Calidad (WQI)",
     ylim = c(0, 100), xlim = c(0, max(datos_filtrados$Turbidez) + 2))

# Agregar curva
x_seq <- seq(min(datos_filtrados$Turbidez), max(datos_filtrados$Turbidez), length.out = 300)
y_pred_plot <- beta0 + beta1 * log(x_seq)
lines(x_seq, y_pred_plot, col = "blue", lwd = 2)

# 9. CORRELACIÓN Y DOMINIO
r_con_log <- cor(log(datos_filtrados$Turbidez), datos_filtrados$Indice_Calidad_promedio)
cat("Correlación de Pearson (datos filtrados):", round(r_con_log, 4), "\n")

## Correlación de Pearson (datos filtrados): -0.6475

# Determinar rango válido (Y > 0)
indices_validos <- which(y_pred_plot > 0)
x_validos <- x_seq[indices_validos]
cat("Rango de validez: ", round(min(x_validos), 2), "< x <", round(max(x_validos), 2), "\n")

## Rango de validez:  0.04 < x < 16.6

# 10. ESTIMACIÓN PUNTUAL (Ejemplo 15 NTU)
x_val_eval <- 15
y_est <- beta0 + beta1 * log(x_val_eval)
cat("Para", x_val_eval, " de Turbidez -> WQI estimado =", round(y_est, 2), "\n")

## Para 15  de Turbidez -> WQI estimado = 49.69