UNIVERSIDAD CENTRAL DEL ECUADOR

MODELO POLINÓMICO

PROYECTO:ESTUDIO ESTADÍSTICO DE LA CALIDAD DE AIRE EN LA INDIA

options(scipen = 999) # Desactivar notación científica para mayor legibilidad

# =============================================================================
# PASO 0: CARGA DE LIBRERÍAS Y ENTORNO
# =============================================================================
if (!require(dplyr)) install.packages("dplyr")
library(dplyr)

# =============================================================================
# PASO 1: IMPORTACIÓN Y PREPROCESAMIENTO DE DATOS
# =============================================================================

datos <- read.csv("C:\\Users\\JOSELYN\\Desktop\\kangle\\Datos Cambiados.csv", sep = ",")

# Conversión de variables a tipo numérico y limpieza inicial de valores nulos (NA)
datos$NO <- as.numeric(datos$NO)
datos$AQI <- as.numeric(datos$AQI)
datos <- datos[!is.na(datos$AQI) & !is.na(datos$NO), ]

# =============================================================================
# PASO 2: DETECCIÓN Y ELIMINACIÓN DE VALORES ATÍPICOS (MÉTODO IQR)
# =============================================================================
# Se aplica el criterio del Rango Intercuartílico (IQR) para filtrar outliers
# que puedan distorsionar el ajuste del modelo en ambas variables.

# Filtrado para AQI
Q1_AQI <- quantile(datos$AQI, 0.25)
Q3_AQI <- quantile(datos$AQI, 0.75)
IQR_AQI <- Q3_AQI - Q1_AQI
datos <- datos[datos$AQI >= (Q1_AQI - 1.5 * IQR_AQI) & datos$AQI <= (Q3_AQI + 1.5 * IQR_AQI), ]

# Filtrado para NO (Monóxido de Nitrógeno)
Q1_NO <- quantile(datos$NO, 0.25)
Q3_NO <- quantile(datos$NO, 0.75)
IQR_NO <- Q3_NO - Q1_NO
datos <- datos[datos$NO >= (Q1_NO - 1.5 * IQR_NO) & datos$NO <= (Q3_NO + 1.5 * IQR_NO), ]

# =============================================================================
# PASO 3: DISCRETIZACIÓN DE LA VARIABLE INDEPENDIENTE (BINNING)
# =============================================================================
# Se segmenta la variable continua 'NO' en 50 intervalos (bins) equidistantes
# para analizar la tendencia central en lugar de la dispersión individual.
datos_bin <- datos %>%
  mutate(NO_bin = cut(NO, breaks = 50)) %>%
  filter(!is.na(NO_bin))

# =============================================================================
# PASO 4: AGREGACIÓN ESTADÍSTICA (REDUCCIÓN DE RUIDO ESTOCÁSTICO)
# =============================================================================
# Se calcula el promedio (media aritmética) de AQI y NO por intervalo.
# Esta técnica suaviza la variabilidad aleatoria y resalta la tendencia subyacente,
# maximizando la correlación de Pearson.
datos_agrupados <- datos_bin %>%
  group_by(NO_bin) %>%
  summarise(
    NO_prom = mean(NO, na.rm = TRUE),
    AQI_prom = mean(AQI, na.rm = TRUE)
  ) %>%
  ungroup() %>%
  na.omit()

# Definición de variables para el modelo: X (Predictora) e Y (Respuesta)
x <- datos_agrupados$NO_prom
y <- datos_agrupados$AQI_prom

# =============================================================================
# PASO 5: AJUSTE DEL MODELO DE REGRESIÓN POLINÓMICA CÚBICA
# =============================================================================
# Se ajusta un modelo lineal de tercer orden: Y = β0 + β1*X + β2*X^2 + β3*X^3
modelo_poly3 <- lm(y ~ x + I(x^2) + I(x^3))

# Cálculo de los coeficientes del modelo
coefs <- coef(modelo_poly3)

# =============================================================================
# PASO 6: REPRESENTACIÓN GRÁFICA Y VALIDACIÓN VISUAL
# =============================================================================
# Cálculo de límites para visualización (se agrega un 10% de margen)
max_x <- max(x, na.rm = TRUE)
max_y <- max(y, na.rm = TRUE)

# Generación de curva de predicción suave
x_grid <- seq(0, max_x, length.out = 100)
y_pred <- predict(modelo_poly3, newdata = data.frame(x = x_grid))

# Cálculo de coeficiente de correlación de Pearson (r)
r <- cor(x, y, method = "pearson")

# Gráfico de dispersión con curva de ajuste
plot(x, y, pch = 19, col = "blue",
     main = paste("Modelo de Regresión Polinómica (Grado 3) | r =", round(r, 4)),
     xlab = expression(paste("Concentración de NO (", mu, "g/m"^3, ")")),
     ylab = "Índice de Calidad del Aire (AQI)",
     xlim = c(0, max_x * 1.1), 
     ylim = c(0, max_y * 1.1),
     xaxs = "i", yaxs = "i") # Intersección de ejes en el origen (0,0)

grid() # Retícula de referencia
lines(x_grid, y_pred, col = "red", lwd = 2) # Curva del modelo

# =============================================================================
# PASO 7: PRESENTACIÓN DE LA ECUACIÓN MATEMÁTICA
# =============================================================================
cat("------------------------------------------------------------------\n")
## ------------------------------------------------------------------
cat("MODELO MATEMÁTICO GENERADO:\n")
## MODELO MATEMÁTICO GENERADO:
cat(sprintf("AQI(NO) = %.4f + %.4f*NO + %.4f*NO^2 + %.4f*NO^3\n", 
            coefs[1], coefs[2], coefs[3], coefs[4]))
## AQI(NO) = 78.1956 + 8.8040*NO + -0.3444*NO^2 + 0.0059*NO^3
cat("------------------------------------------------------------------\n")
## ------------------------------------------------------------------
# =============================================================================
# PASO 8: EVALUACIÓN DE LA BONDAD DE AJUSTE (COEFICIENTE R²)
# =============================================================================
# El coeficiente de determinación (R²) cuantifica la proporción de varianza
# de la variable dependiente explicada por el modelo.
R2 <- r^2
print(paste("Coeficiente de determinación (R²):", round(R2 * 100, 2), "%"))
## [1] "Coeficiente de determinación (R²): 92.13 %"
# =============================================================================
# PASO 9: ANÁLISIS DEL DOMINIO DE VALIDEZ DEL MODELO
# =============================================================================
# Determinación de las raíces del polinomio para establecer el límite físico inferior
# (Concentración mínima de NO para la cual el modelo predice AQI > 0).
raices <- polyroot(coefs) # Cálculo de raíces complejas y reales

# Filtrado de raíces reales
raices_reales <- Re(raices[abs(Im(raices)) < 1e-10])

# Selección de la raíz relevante (límite operativo)
limite_x <- max(raices_reales) 

print(paste("Límite matemático inferior (Intersección con eje X):", round(limite_x, 2)))
## [1] "Límite matemático inferior (Intersección con eje X): -6.84"
cat("RESTRICCIÓN DEL DOMINIO:\n")
## RESTRICCIÓN DEL DOMINIO:
if(limite_x < 0) {
  cat("Dado que el límite matemático es negativo y la concentración física\n")
  cat("no puede ser < 0, el dominio de validez es: NO >= 0.\n")
} else {
  cat("El modelo es válido estrictamente para concentraciones de NO >", round(limite_x, 2), "\n")
}
## Dado que el límite matemático es negativo y la concentración física
## no puede ser < 0, el dominio de validez es: NO >= 0.
# =============================================================================
# PASO 10: APLICACIÓN PREDICTIVA (SIMULACIÓN)
# =============================================================================
# Estimación del AQI para un valor de control específico de NO.
NO_ejemplo <- 35 
AQI_esperado <- predict(modelo_poly3, newdata = data.frame(x = NO_ejemplo))

print(paste("Predicción: Para una concentración de NO =", NO_ejemplo, 
            "µg/m³, el AQI estimado es:", round(AQI_esperado, 2)))
## [1] "Predicción: Para una concentración de NO = 35 µg/m³, el AQI estimado es: 217.4"
# =============================================================================
# CONCLUSIÓN FINAL DEL ESTUDIO
# =============================================================================
cat("\n--- CONCLUSIÓN DEL ANÁLISIS ---\n")
## 
## --- CONCLUSIÓN DEL ANÁLISIS ---
cat(paste("El modelo demuestra que el Índice de Calidad del Aire (AQI) está determinado\n"))
## El modelo demuestra que el Índice de Calidad del Aire (AQI) está determinado
cat(paste("en un", round(R2 * 100, 2), "% por la concentración de Monóxido de Nitrógeno (NO).\n"))
## en un 92.13 % por la concentración de Monóxido de Nitrógeno (NO).
cat(paste("El", round((1 - R2) * 100, 2), "% de varianza no explicada se atribuye a variables latentes\n"))
## El 7.87 % de varianza no explicada se atribuye a variables latentes
cat("no incluidas en este modelo univariante (meteorología, otros precursores químicos, etc.).\n")
## no incluidas en este modelo univariante (meteorología, otros precursores químicos, etc.).