1 Carga de Librerías

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(gt)

2 Carga de Datos

datos <- read.csv("~/Estudio/TERCER SEMESTRE/Estadistica/proyecto/Dataset.csv",
                   sep = ";", stringsAsFactors = FALSE)
str(datos)
## 'data.frame':    2997 obs. of  6 variables:
##  $ EMIS_CO2 : num  130.1 318.7 446.3 71.4 301.6 ...
##  $ EMIS_NOX : num  66 151.7 369.6 30.5 262.8 ...
##  $ EMIS_SO2 : num  9.73 16.2 21.18 3.94 16.89 ...
##  $ EMIS_CO  : num  36.25 161 260.74 5.17 128.06 ...
##  $ EMIS_PM10: num  39.8 76.1 78.9 20.4 61.1 ...
##  $ EMIS_PM25: num  24.5 27.8 28.3 22.2 26.5 ...
summary(datos)
##     EMIS_CO2        EMIS_NOX         EMIS_SO2         EMIS_CO      
##  Min.   : 20.0   Min.   :  0.00   Min.   : 0.000   Min.   :  0.00  
##  1st Qu.:140.0   1st Qu.: 99.94   1st Qu.: 9.772   1st Qu.: 42.48  
##  Median :259.0   Median :187.00   Median :15.128   Median :107.51  
##  Mean   :260.4   Mean   :188.33   Mean   :13.613   Mean   :125.81  
##  3rd Qu.:380.6   3rd Qu.:275.27   3rd Qu.:18.349   3rd Qu.:199.96  
##  Max.   :499.8   Max.   :439.54   Max.   :25.148   Max.   :376.00  
##    EMIS_PM10         EMIS_PM25    
##  Min.   :  7.679   Min.   :14.45  
##  1st Qu.: 38.407   1st Qu.:24.41  
##  Median : 56.428   Median :26.94  
##  Mean   : 54.282   Mean   :26.32  
##  3rd Qu.: 70.536   3rd Qu.:28.60  
##  Max.   :104.859   Max.   :35.50
colSums(is.na(datos))
##  EMIS_CO2  EMIS_NOX  EMIS_SO2   EMIS_CO EMIS_PM10 EMIS_PM25 
##         0         0         0         0         0         0

3 Definición de Variables

Se seleccionan EMIS_NOX (variable independiente) y EMIS_CO2 (variable dependiente). Al igual que en la regresión lineal, ambos gases provienen del mismo proceso de combustión, pero en este dataset (Dataset_2_) la relación presenta mayor dispersión, lo que sugiere evaluar si una curva de segundo grado captura mejor el comportamiento que una recta.

nox <- as.numeric(datos$EMIS_NOX)
co2 <- as.numeric(datos$EMIS_CO2)

nox[is.na(nox)] <- median(nox, na.rm = TRUE)
co2[is.na(co2)] <- median(co2, na.rm = TRUE)

TPV <- data.frame(nox, co2)

dim(TPV)
## [1] 2997    2
colSums(is.na(TPV))
## nox co2 
##   0   0

4 Tabla de Pares de Valores

cat("Tamaño muestral: n =", nrow(TPV))
## Tamaño muestral: n = 2997
TPV_tabla <- head(TPV, 10)

TPV_tabla %>%
  gt() %>%
  tab_header(title = "Tabla N°1: Pares de valores EMIS_NOX - EMIS_CO2") %>%
  cols_label(nox = "EMIS_NOX", co2 = "EMIS_CO2") %>%
  tab_source_note(source_note = "Autor: Luis Cruz")
Tabla N°1: Pares de valores EMIS_NOX - EMIS_CO2
EMIS_NOX EMIS_CO2
66.010 130.080
151.677 318.677
369.579 446.328
30.482 71.436
262.797 301.645
178.235 258.165
76.423 62.459
297.944 493.452
214.639 301.545
193.373 333.784
Autor: Luis Cruz

5 Gráfica de Dispersión

x <- as.numeric(TPV$nox)
y <- as.numeric(TPV$co2)

validos <- complete.cases(x, y) &
  is.finite(x) &
  is.finite(y)

x <- x[validos]
y <- y[validos]

orden <- order(x)
x <- x[orden]
y <- y[orden]

plot(x, y,
     pch = 16,
     col = "blue",
     xlim = c(0, max(x)),
     ylim = c(0, max(y)),
     main = "Gráfica N°1: Diagrama de dispersión entre EMIS_NOX\ny EMIS_CO2 en emisiones registradas",
     xlab = "EMIS_NOX (unidades)",
     ylab = "EMIS_CO2 (unidades)")

6 Tratamiento de Datos

A diferencia de las regresiones anteriores, aquí la nube de puntos es considerablemente más dispersa y compleja: no se observa una tendencia tan definida a simple vista. Por este motivo se aplica el tratamiento de datos indicado para nubes caóticas (único x, único y, omisión de outliers) antes de plantear la conjetura definitiva.

# Único x, único y: ya se garantiza en la Tabla de Pares (un par por fila)

# Omisión de outliers por rango intercuartílico (IQR)
filtro_iqr <- function(v) {
  Q1 <- quantile(v, 0.25)
  Q3 <- quantile(v, 0.75)
  IQRv <- Q3 - Q1
  li <- Q1 - 1.5 * IQRv
  ls <- Q3 + 1.5 * IQRv
  return(v >= li & v <= ls)
}

mantener <- filtro_iqr(x) & filtro_iqr(y)

x_limpio <- x[mantener]
y_limpio <- y[mantener]

cat("Observaciones originales:", length(x),
    "\nObservaciones excluidas por outliers:", length(x) - length(x_limpio),
    "\nObservaciones finales:", length(x_limpio))
## Observaciones originales: 2997 
## Observaciones excluidas por outliers: 0 
## Observaciones finales: 2997

No se identificaron valores atípicos relevantes fuera de los límites del IQR: la dispersión observada no se debe a datos anómalos, sino a una variabilidad genuina propia de este dataset (menor correlación entre EMIS_NOX y EMIS_CO2 que en el dataset original).

6.1 5.1 Tabla Pares de Valores Simplificada

No aplica: no se eliminaron observaciones adicionales por el filtro IQR, por lo que la tabla de pares de valores (Tabla N°1) ya representa los datos finales usados en el modelo.

6.2 5.2 Gráfica de Dispersión Simplificada

plot(x_limpio, y_limpio,
     pch = 16,
     col = "blue",
     xlim = c(0, max(x_limpio)),
     ylim = c(0, max(y_limpio)),
     main = "Gráfica N°1.1: Diagrama de dispersión (sin outliers)\nentre EMIS_NOX y EMIS_CO2",
     xlab = "EMIS_NOX (unidades)",
     ylab = "EMIS_CO2 (unidades)")

7 Conjetura del Modelo

La forma de la nube, aunque dispersa, sugiere una leve curvatura, por lo que se mantiene la conjetura de un modelo polinómico de grado 2.

8 Parámetros

x <- x_limpio
y <- y_limpio

modelo_poli <- lm(y ~ x + I(x^2))
modelo_poli
## 
## Call:
## lm(formula = y ~ x + I(x^2))
## 
## Coefficients:
## (Intercept)            x       I(x^2)  
##  13.2631775    1.4150485   -0.0004167
c0 <- coef(modelo_poli)[1]
c1 <- coef(modelo_poli)[2]
c2 <- coef(modelo_poli)[3]

c0
## (Intercept) 
##    13.26318
c1
##        x 
## 1.415048
c2
##        I(x^2) 
## -0.0004167489
plot(1, type = "n", axes = FALSE, xlab = "", ylab = "")

texto_ecuacion <- paste0(
  "Modelo polinómico general (grado 2)\n",
  "Y = c0 + c1*x + c2*x^2\n\n",
  "Modelo ajustado:\n",
  "Y = ", round(c0, 4), " + ", round(c1, 4), "x + (",
  round(c2, 6), ")x^2"
)

text(x = 1, y = 1, labels = texto_ecuacion, cex = 1.3, col = "blue", font = 2)

9 Comparación de la Realidad con el Modelo

x_modelo <- seq(min(x), max(x), length.out = 200)
y_modelo <- c0 + c1 * x_modelo + c2 * x_modelo^2

plot(x, y,
     pch = 16,
     col = "blue",
     xlim = c(0, max(x)),
     ylim = c(0, max(y)),
     main = "Gráfica N°2: Comparación entre la realidad y el modelo polinómico\nentre EMIS_NOX y EMIS_CO2",
     xlab = "EMIS_NOX (unidades)",
     ylab = "EMIS_CO2 (unidades)")

lines(x_modelo, y_modelo, col = "red", lwd = 3)

legend("topleft",
       legend = c("Datos reales", "Modelo polinómico"),
       col = c("blue", "red"),
       pch = c(16, NA),
       lty = c(NA, 1),
       lwd = c(NA, 3))

10 Test de Bondad

r <- cor(x, y)
r * 100
## [1] 94.99882
predicciones <- predict(modelo_poli)

SSE <- sum((y - predicciones)^2)
SST <- sum((y - mean(y))^2)

r2 <- 1 - (SSE / SST)
r2 * 100
## [1] 90.35313

11 Restricciones

# Raíces del modelo: valores de x donde y = 0
# y = c2*x^2 + c1*x + c0  ->  fórmula general de la cuadrática
discriminante <- c1^2 - 4 * c2 * c0

raices <- sort(unname(c((-c1 - sqrt(discriminante)) / (2 * c2),
                         (-c1 + sqrt(discriminante)) / (2 * c2))))

raices
## [1]   -9.347218 3404.793611
# c2 es el coeficiente que acompaña a x^2: define hacia dónde abre la parábola
# c2 < 0 (abre hacia abajo): el modelo es negativo FUERA de las raíces
# c2 > 0 (abre hacia arriba): el modelo es negativo ENTRE las raíces
intervalo <- if (c2 < 0) {
  paste0("x < ", round(raices[1], 2), " o x > ", round(raices[2], 2))
} else {
  paste0(round(raices[1], 2), " < x < ", round(raices[2], 2))
}

ecuacion <- paste0("y = ", round(c0, 4), " + ", round(c1, 4), "x + (",
                    round(c2, 6), ")x^2")

tabla_restricciones <- data.frame(
  Modelo = c("Polinómica", "Polinómica"),
  Ecuacion = c(ecuacion, ecuacion),
  Condicion = c("y = 0 (raíces)", "y < 0 (no aplica)"),
  Valor = c(paste0("x = ", round(raices[1], 2), "  y  x = ", round(raices[2], 2)),
            intervalo)
)

tabla_restricciones %>%
  gt() %>%
  tab_header(title = "Tabla N°2: Valores donde el modelo es 0 o negativo") %>%
  cols_label(Modelo = "Modelo",
             Ecuacion = "Ecuación",
             Condicion = "Condición",
             Valor = "Valor / Intervalo de x") %>%
  cols_width(
    Modelo ~ px(110),
    Ecuacion ~ px(240),
    Condicion ~ px(140),
    Valor ~ px(220)
  ) %>%
  tab_style(
    style = cell_text(align = "center"),
    locations = cells_body()
  ) %>%
  tab_style(
    style = cell_text(align = "center"),
    locations = cells_column_labels()
  ) %>%
  tab_source_note(source_note = "Autor: Luis Cruz")
Tabla N°2: Valores donde el modelo es 0 o negativo
Modelo Ecuación Condición Valor / Intervalo de x
Polinómica y = 13.2632 + 1.415x + (-0.000417)x^2 y = 0 (raíces) x = -9.35 y x = 3404.79
Polinómica y = 13.2632 + 1.415x + (-0.000417)x^2 y < 0 (no aplica) x < -9.35 o x > 3404.79
Autor: Luis Cruz

El modelo corta el eje x en dos puntos: x = -9.35 y x = 3404.79. Como c2 es negativo, la parábola abre hacia abajo, por lo que el modelo da valores negativos de EMIS_CO2 cuando x < -9.35 o x > 3404.79.

Nota: dominio matemático de x: todos los reales (el modelo polinómico no exige ninguna restricción sobre x). Dominio de y: en la realidad EMIS_CO2 solo puede tomar valores y ≥ 0.

12 Estimación

NOX_ejemplo <- 250

CO2_esp <- c0 + c1 * NOX_ejemplo + c2 * NOX_ejemplo^2
CO2_esp
## (Intercept) 
##    340.9785
plot(1, type = "n", axes = FALSE, xlab = "", ylab = "")

texto_pronostico <- paste0(
  "¿Qué valor de EMIS_CO2 se espera\n",
  "cuando EMIS_NOX es de 250 unidades?\n\n",
  "R = ", round(CO2_esp, 2), " unidades"
)

text(x = 1, y = 1, labels = texto_pronostico, cex = 1.5, col = "blue", font = 2)

13 Conclusión

Entre EMIS_CO2 (variable dependiente) y EMIS_NOX (variable independiente) existe una relación de tipo no lineal (polinómica de grado 2) cuyo modelo matemático es:

\[Y = 13.2632 + 1.415\,X + (-4.17\times 10^{-4})\,X^2\]

siendo Y = EMIS_CO2 y X = EMIS_NOX, donde existen restricciones: el modelo produce valores negativos de EMIS_CO2 cuando x < -9.35 o x > 3404.79 (raíces en x = -9.35 y x = 3404.79).