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)
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
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
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 | |
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)")
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).
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.
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)")
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.
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)
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))
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
# 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.
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)
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 sí 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).