library(readr)
## Warning: package 'readr' was built under R version 4.6.1
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(knitr)
library(gt)
library(scatterplot3d)
datos <- read.csv("~/Estudio/TERCER SEMESTRE/Estadistica/proyecto/Dataset.csv",
sep = ";", stringsAsFactors = FALSE)
# Estructura del dataset original
str(datos)
## 'data.frame': 2997 obs. of 6 variables:
## $ EMIS_CO2 : num 74.3 106.9 325.4 129 231.5 ...
## $ EMIS_NOX : num 43.3 56.5 210.2 50.7 198.4 ...
## $ EMIS_SO2 : num 11.5 0 43.3 11.4 53.2 ...
## $ EMIS_CO : num 22.8 40.1 175.1 34.8 80.5 ...
## $ EMIS_PM10: num 46.1 43.3 56.4 45.2 75.5 ...
## $ EMIS_PM25: num 24.7 24.2 31.4 27.2 26.7 ...
# Resumen estadístico general
summary(datos)
## EMIS_CO2 EMIS_NOX EMIS_SO2 EMIS_CO
## Min. : 25.04 Min. : 0.0 Min. : 0.000 Min. : 0.00
## 1st Qu.:134.78 1st Qu.: 89.5 1st Qu.: 6.009 1st Qu.: 68.47
## Median :247.72 Median :167.4 Median : 22.010 Median :133.61
## Mean :249.87 Mean :171.2 Mean : 24.588 Mean :133.34
## 3rd Qu.:367.31 3rd Qu.:253.3 3rd Qu.: 38.786 3rd Qu.:194.60
## Max. :479.91 Max. :439.6 Max. :103.518 Max. :359.58
## EMIS_PM10 EMIS_PM25
## Min. : 5.045 Min. :14.07
## 1st Qu.: 48.875 1st Qu.:25.44
## Median : 66.432 Median :28.33
## Mean : 65.429 Mean :27.87
## 3rd Qu.: 83.030 3rd Qu.:30.55
## Max. :125.408 Max. :38.91
# Verificación de valores NA en el dataset original
colSums(is.na(datos))
## EMIS_CO2 EMIS_NOX EMIS_SO2 EMIS_CO EMIS_PM10 EMIS_PM25
## 0 0 0 0 0 0
# Selección de variables geoquímicas
co2 <- as.numeric(datos$EMIS_CO2)
nox <- as.numeric(datos$EMIS_NOX)
pm10 <- as.numeric(datos$EMIS_PM10)
# Imputación de NA con la mediana (en vez de omitir filas)
co2[is.na(co2)] <- median(co2, na.rm = TRUE)
nox[is.na(nox)] <- median(nox, na.rm = TRUE)
pm10[is.na(pm10)] <- median(pm10, na.rm = TRUE)
Geoquimica <- data.frame(
y = co2,
x1 = nox,
x2 = pm10
)
# Dimensión del dataset ya tratado
dim(Geoquimica)
## [1] 2997 3
# Verificación final: ya no existen NA
colSums(is.na(Geoquimica))
## y x1 x2
## 0 0 0
# Solo valores positivos
Geoquimica <- Geoquimica[
Geoquimica$y > 0 &
Geoquimica$x1 > 0 &
Geoquimica$x2 > 0, ]
# Función para filtrar outliers (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)
}
# Tabla sin outliers
Geoquimica <- Geoquimica[
filtro_iqr(Geoquimica$y) &
filtro_iqr(Geoquimica$x1) &
filtro_iqr(Geoquimica$x2), ]
dim(Geoquimica)
## [1] 2913 3
Geoquimica_tabla <- head(Geoquimica, 10)
Geoquimica_tabla %>%
gt() %>%
tab_header(title = "Tabla N°1: Pares de valores EMIS_NOX, EMIS_PM10 - EMIS_CO2") %>%
cols_label(y = "EMIS_CO2", x1 = "EMIS_NOX", x2 = "EMIS_PM10") %>%
tab_source_note(source_note = "Autor: Grupo 4 - Minas")
| Tabla N°1: Pares de valores EMIS_NOX, EMIS_PM10 - EMIS_CO2 | ||
| EMIS_CO2 | EMIS_NOX | EMIS_PM10 |
|---|---|---|
| 74.292 | 43.270 | 46.098 |
| 106.929 | 56.510 | 43.258 |
| 325.400 | 210.165 | 56.446 |
| 128.963 | 50.718 | 45.172 |
| 231.505 | 198.402 | 75.545 |
| 375.165 | 334.486 | 93.046 |
| 263.942 | 156.604 | 75.554 |
| 57.889 | 3.504 | 33.807 |
| 56.052 | 39.640 | 45.515 |
| 37.843 | 24.874 | 27.796 |
| Autor: Grupo 4 - Minas | ||
# Extraer variables
x1 <- Geoquimica$x1 # EMIS_NOX
x2 <- Geoquimica$x2 # EMIS_PM10
y <- Geoquimica$y # EMIS_CO2
Geoq_reg <- scatterplot3d(
x1, x2, y,
angle = 60,
pch = 16,
color = "darkred",
main = "Gráfica N°1: Diagrama de dispersión de EMIS_CO2,\nEMIS_NOX y EMIS_PM10",
xlab = "EMIS_NOX (unidades)",
ylab = "EMIS_PM10 (unidades)",
zlab = "EMIS_CO2 (unidades)"
)
regresion_multiple <- lm(y ~ x1 + x2)
summary(regresion_multiple)
##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -155.856 -28.851 -1.181 29.823 130.544
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -69.95292 2.79652 -25.01 <2e-16 ***
## x1 0.63545 0.01365 46.54 <2e-16 ***
## x2 3.21464 0.06236 51.55 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 43.76 on 2910 degrees of freedom
## Multiple R-squared: 0.8879, Adjusted R-squared: 0.8878
## F-statistic: 1.153e+04 on 2 and 2910 DF, p-value: < 2.2e-16
# Diagrama de dispersión 3D
Geoq_reg <- scatterplot3d(
x1, x2, y,
angle = 60,
pch = 16,
color = "blue",
main = "Gráfica N°2: Comparación de la realidad con el modelo\nmultivariable lineal de EMIS_CO2, EMIS_NOX y EMIS_PM10",
xlab = "EMIS_NOX (unidades)",
ylab = "EMIS_PM10 (unidades)",
zlab = "EMIS_CO2 (unidades)"
)
# Plano de regresión
Geoq_reg$plane3d(regresion_multiple)
# Ecuación del modelo
b0 <- coef(regresion_multiple)[1]
b1 <- coef(regresion_multiple)[2]
b2 <- coef(regresion_multiple)[3]
plot(1,
type = "n",
axes = FALSE,
xlab = "",
ylab = "")
texto_ecuacion <- paste0(
"Modelo lineal multivariable general\n",
"Y = b0 + b1*x1 + b2*x2\n\n",
"Modelo ajustado:\n",
"Y = ",
round(b0, 4),
" + ",
round(b1, 4),
"*x1 + ",
round(b2, 4),
"*x2"
)
text(x = 1,
y = 1,
labels = texto_ecuacion,
cex = 1.2,
col = "blue",
font = 2)
# Pearson
r <- cor(y, x1 + x2)
r * 100
## [1] 91.68799
# Determinación
r2 <- r^2 * 100
r2
## [1] 84.06688
Con dos variables independientes, el modelo y = b0 + b1x1 + b2x2 no tiene una raíz única: y = 0 define una línea en el plano (x1, x2), no un solo punto. Esa ecuación tiene infinitas combinaciones de (x1, x2) que la cumplen; no existe una única solución. Para tener valores concretos, se calcula el valor crítico de cada variable dejando la otra fija en su media.
media_x1 <- mean(x1)
media_x2 <- mean(x2)
# Valor crítico de x1 (EMIS_NOX), dejando x2 (EMIS_PM10) en su media
x1_critico <- unname((-b0 - b2 * media_x2) / b1)
# Valor crítico de x2 (EMIS_PM10), dejando x1 (EMIS_NOX) en su media
x2_critico <- unname((-b0 - b1 * media_x1) / b2)
x1_critico
## [1] -225.6581
x2_critico
## [1] -13.05044
ecuacion <- paste0("y = ", round(b0, 4), " + ", round(b1, 4),
"*x1 + ", round(b2, 4), "*x2")
# b1 y b2 positivos: el modelo es negativo por DEBAJO del valor crítico
sentido_x1 <- ifelse(b1 > 0,
paste0("x1 < ", round(x1_critico, 2)),
paste0("x1 > ", round(x1_critico, 2)))
sentido_x2 <- ifelse(b2 > 0,
paste0("x2 < ", round(x2_critico, 2)),
paste0("x2 > ", round(x2_critico, 2)))
tabla_restricciones <- data.frame(
Modelo = c("Múltiple", "Múltiple", "Múltiple"),
Ecuacion = c(ecuacion, ecuacion, ecuacion),
Condicion = c("y = 0 (plano)",
paste0("x1 crítico (x2 = ", round(media_x2, 2), ")"),
paste0("x2 crítico (x1 = ", round(media_x1, 2), ")")),
Valor = c("0 = b0 + b1*x1 + b2*x2",
paste0("x1 = ", round(x1_critico, 2), " → negativo si ", sentido_x1),
paste0("x2 = ", round(x2_critico, 2), " → negativo si ", sentido_x2))
)
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") %>%
cols_width(
Modelo ~ px(100),
Ecuacion ~ px(220),
Condicion ~ px(190),
Valor ~ px(260)
) %>%
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: Grupo 4 - Minas")
| Tabla N°2: Valores donde el modelo es 0 o negativo | |||
| Modelo | Ecuación | Condición | Valor / Intervalo |
|---|---|---|---|
| Múltiple | y = -69.9529 + 0.6354*x1 + 3.2146*x2 | y = 0 (plano) | 0 = b0 + b1*x1 + b2*x2 |
| Múltiple | y = -69.9529 + 0.6354*x1 + 3.2146*x2 | x1 crítico (x2 = 66.37) | x1 = -225.66 → negativo si x1 < -225.66 |
| Múltiple | y = -69.9529 + 0.6354*x1 + 3.2146*x2 | x2 crítico (x1 = 176.11) | x2 = -13.05 → negativo si x2 < -13.05 |
| Autor: Grupo 4 - Minas | |||
Estos valores críticos son condicionales: cada uno asume que la otra variable se mantiene en su media muestral. No representan una “raíz” absoluta del modelo, sino un corte de referencia sobre la línea de frontera y = 0.
Nota: dominio matemático de x1 y x2: todos los reales (el modelo lineal múltiple no exige ninguna restricción sobre las variables independientes). Dominio de y: en la realidad EMIS_CO2 solo puede tomar valores y ≥ 0.
"¿Qué valor de EMIS_CO2 se espera
con un EMIS_NOX de 200 y un EMIS_PM10 de 100?"
## [1] "¿Qué valor de EMIS_CO2 se espera\ncon un EMIS_NOX de 200 y un EMIS_PM10 de 100?"
CO2_Est <- regresion_multiple$coefficients[1] +
regresion_multiple$coefficients[2] * 200 +
regresion_multiple$coefficients[3] * 100
CO2_Est
## (Intercept)
## 378.6003
plot(1,
type = "n",
axes = FALSE,
xlab = "",
ylab = "")
text(x = 1,
y = 1,
labels = paste0(
"¿Qué valor de EMIS_CO2 se espera\n",
"con un EMIS_NOX de 200\n",
"y un EMIS_PM10 de 100?\n",
"R = ", round(CO2_Est, 2)
),
cex = 1.8,
col = "darkred",
font = 6)
Entre las emisiones de EMIS_NOX, EMIS_PM10 y EMIS_CO2 existe una relación positiva de tipo lineal múltiple, representada por el modelo:
\[\text{EMIS\_CO2} = -69.9529 + 0.6354\,\text{EMIS\_NOX} + 3.2146\,\text{EMIS\_PM10}\]
En esta ecuación, EMIS_NOX y EMIS_PM10 corresponden a las variables independientes y EMIS_CO2 a la variable dependiente. El coeficiente de correlación de Pearson de 91.69% evidencia una asociación positiva muy fuerte, mientras que el coeficiente de determinación de 84.07% indica la proporción de variabilidad de EMIS_CO2 explicada por el modelo.
El filtrado mediante el rango intercuartílico redujo la influencia de valores extremos y permitió observar con mayor claridad la tendencia de los datos. Dado que el modelo tiene dos variables independientes, no existe una raíz única sino una línea de frontera y = 0; por ejemplo, con EMIS_PM10 fijo en su media (66.37), el modelo se vuelve negativo cuando EMIS_NOX < -225.66.