library(dplyr)
library(gt)

1 Carga de Datos y Librerías

datos <- read.csv("~/Estudio/TERCER SEMESTRE/Estadistica/Dataset.csv",
                   sep = ";", stringsAsFactors = FALSE)

EMIS_CO2 <- as.numeric(gsub(",", ".", datos$CO2.total.emissions..tons..by.state))
EMIS_CO2 <- na.omit(EMIS_CO2)
n <- length(EMIS_CO2)
n
## [1] 2921

2 Tabla de Distribución de Frecuencias

k <- ceiling(1 + 3.322 * log10(n))
k
## [1] 13
breaks <- seq(min(EMIS_CO2), max(EMIS_CO2) + 1e-6, length.out = k + 1)

histograma_co2 <- hist(EMIS_CO2, breaks = breaks, plot = FALSE)

lis <- round(histograma_co2$breaks[1:k], 2)
lss <- round(histograma_co2$breaks[2:(k+1)], 2)
MC  <- round(histograma_co2$mids, 2)
ni  <- histograma_co2$counts
hi  <- round((ni/sum(ni))*100, 2)

TDFco2 <- data.frame(lis, lss, MC, ni, hi)

fila_total <- data.frame(
  lis = "TOTAL", lss = "", MC = "",
  ni  = sum(TDFco2$ni),
  hi  = round(sum(TDFco2$hi), 2)
)

TDFco2_total <- rbind(TDFco2, fila_total)

tabla_co2 <- TDFco2_total %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N°1**"),
    subtitle = md("Tabla de distribución de frecuencias de EMIS_CO2 (nivel-mina)")
  ) %>%
  cols_label(
    lis = "Límite inferior",
    lss = "Límite superior",
    MC  = "Marca de clase",
    ni  = "Frecuencia absoluta",
    hi  = "Frecuencia relativa (%)"
  ) %>%
  tab_source_note(
    source_note = md("Autor: Grupo 4 - Minas")
  ) %>%
  tab_options(
    table.border.top.color = "black",
    table.border.bottom.color = "black",
    table.border.top.style = "solid",
    table.border.bottom.style = "solid",
    column_labels.border.top.color = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    row.striping.include_table_body = TRUE,
    heading.border.bottom.color = "black",
    heading.border.bottom.width = px(2),
    table_body.hlines.color = "gray",
    table_body.border.bottom.color = "black"
  )

tabla_co2
Tabla N°1
Tabla de distribución de frecuencias de EMIS_CO2 (nivel-mina)
Límite inferior Límite superior Marca de clase Frecuencia absoluta Frecuencia relativa (%)
49075.86 18001114.42 9025095.14 174 5.96
18001114.42 35953152.97 26977133.7 468 16.02
35953152.97 53905191.53 44929172.25 335 11.47
53905191.53 71857230.09 62881210.81 1210 41.42
71857230.09 89809268.64 80833249.37 533 18.25
89809268.64 107761307.2 98785287.92 62 2.12
107761307.2 125713345.76 116737326.48 23 0.79
125713345.76 143665384.32 134689365.04 0 0.00
143665384.32 161617422.87 152641403.59 0 0.00
161617422.87 179569461.43 170593442.15 0 0.00
179569461.43 197521499.99 188545480.71 0 0.00
197521499.99 215473538.54 206497519.26 0 0.00
215473538.54 233425577.1 224449557.82 116 3.97
TOTAL 2921 100.00
Autor: Grupo 4 - Minas

3 Gráfica de Distribución de Frecuencias

hist(EMIS_CO2,
     breaks = breaks,
     main = "Gráfica N°1: Distribución de EMIS_CO2 (nivel-mina)",
     xlab = "Emisiones de CO2 (toneladas)",
     ylab = "Frecuencia",
     col  = "gray")

4 Conjetura del Modelo

Las emisiones de CO2 se registran a nivel estatal (no por mina), y cada una de las 2921 minas del dataset hereda el valor de su estado. Esto significa que, aunque la tabla y el histograma anteriores usan las 2921 filas (igual que en la fase Descriptiva), la muestra real e independiente para conjeturar y probar un modelo de probabilidad son los valores únicos por estado: 41 estados con dato válido de CO2. Usar n = 2921 en el test de aprobación infla artificialmente el tamaño muestral y sobre-potencia el Chi-cuadrado (mismo criterio ya documentado en la fase de Regresión para estas variables de emisiones).

Por esta razón, desde aquí en adelante se trabaja con EMIS_CO2_edo (n = 41).

EMIS_CO2_edo <- unique(EMIS_CO2)
n_edo <- length(EMIS_CO2_edo)
n_edo
## [1] 41
# Calculo de parametros Lognormal (metodo de momentos, a partir de la
# media y varianza aritmeticas de los datos originales, sin transformar)
media_edo <- mean(EMIS_CO2_edo)
var_edo   <- var(EMIS_CO2_edo)

sdlog_1   <- sqrt(log(1 + var_edo / media_edo^2))
meanlog_1 <- log(media_edo) - (sdlog_1^2) / 2

meanlog_1
## [1] 17.4191
sdlog_1
## [1] 0.7192199
densidad_hist <- hist(EMIS_CO2_edo, plot = FALSE)$density
densidad_curva <- dlnorm(seq(min(EMIS_CO2_edo), max(EMIS_CO2_edo), length.out = 500),
                          meanlog = meanlog_1, sdlog = sdlog_1)

ylim_max <- max(densidad_hist, densidad_curva) * 1.05

hist(EMIS_CO2_edo,
     freq = FALSE,
     ylim = c(0, ylim_max),
     main = "Gráfica N°2: Comparación de la realidad con el modelo Lognormal\nde EMIS_CO2 (nivel-estado)",
     ylab = "Densidad de probabilidad",
     xlab = "Emisiones de CO2 (toneladas)",
     col = "lightgray",
     border = "black",
     yaxt = "n")

# El eje Y (densidad) maneja numeros muy pequenos; se formatea en notacion
# cientifica para que no salga como 0.000000012 (efecto de scipen = 999)
ejeY <- axTicks(2)
axis(2, at = ejeY, labels = formatC(ejeY, format = "e", digits = 1))

x <- seq(min(EMIS_CO2_edo), max(EMIS_CO2_edo) * 1.1, length.out = 500)
curve(dlnorm(x, meanlog = meanlog_1, sdlog = sdlog_1),
      from = min(EMIS_CO2_edo),
      to = max(EMIS_CO2_edo) * 1.1,
      n = 500,
      col = "blue",
      lwd = 2,
      add = TRUE)

Nota: se conjeturó primero un modelo Normal, siguiendo el mismo tratamiento usado en el resto de variables continuas del proyecto. Sin embargo, ese modelo fue rechazado por el test de Chi-cuadrado (ver sección siguiente), debido a la fuerte asimetría de los datos —hay un estado con una emisión extrema (233,425,577 t), casi el doble del siguiente valor más alto. Por eso se reemplazó por un modelo Lognormal, adecuado para variables continuas positivas y asimétricas como ésta (de hecho, ajustó incluso mejor que un modelo Gamma probado como alternativa).

5 Test de Aprobación

k_edo <- ceiling(1 + 3.322 * log10(n_edo))
k_edo
## [1] 7
breaks_edo <- seq(min(EMIS_CO2_edo), max(EMIS_CO2_edo) + 1, length.out = k_edo + 1)

cortes_edo <- cut(EMIS_CO2_edo, breaks = breaks_edo, include.lowest = TRUE, right = FALSE)

Fo_1 <- as.vector(table(cortes_edo))
Fo_1
## [1] 17 16  6  1  0  0  1
Fe_1 <- n_edo * diff(plnorm(breaks_edo, meanlog = meanlog_1, sdlog = sdlog_1))
Fe_1
## [1] 18.3372632 14.3300235  4.9817969  1.8561254  0.7686519  0.3484000  0.1699903
# Expresar Fo y Fe en porcentaje para comparar
Fo_pct <- (Fo_1/n_edo)*100
Fe_pct <- (Fe_1/n_edo)*100

plot(Fo_pct, Fe_pct,
     xlim = c(0, max(Fo_pct, Fe_pct)),
     ylim = c(0, max(Fo_pct, Fe_pct)),
     main = "Gráfica N°3: Correlación de frecuencias observadas y esperadas\ndel modelo Lognormal de EMIS_CO2",
     xlab = "Frecuencia Observada (%)",
     ylab = "Frecuencia Esperada (%)",
     col  = "blue3",
     pch  = 19)

abline(a = 0, b = 1, col = "red", lwd = 2)

Correlacion_1 <- cor(Fo_pct, Fe_pct)*100
Correlacion_1
## [1] 98.8651
# Chi-cuadrado
grados_libertad_1 <- length(Fo_1) - 1 - 2   # 2 parametros estimados (meanlog, sdlog)
grados_libertad_1
## [1] 4
nivel_significancia <- 0.95

x2_1 <- sum((Fo_1 - Fe_1)^2 / Fe_1)
x2_1
## [1] 6.064852
umbral_aceptacion_1 <- qchisq(nivel_significancia, grados_libertad_1)
umbral_aceptacion_1
## [1] 9.487729
# Tabla resumen
Variable <- c("Emisiones de CO2 (ton)")

tabla_resumen_1 <- data.frame(
  Variable,
  round(Correlacion_1, 2),
  round(x2_1, 2),
  round(umbral_aceptacion_1, 2)
)

colnames(tabla_resumen_1) <- c("Variable", "Test Pearson (%)", "Chi Cuadrado", "Umbral de aceptación")

tabla_resumen_1 %>%
  gt() %>%
  tab_header(title = md("**Tabla N°2**"),
             subtitle = md("Resumen del test de bondad de ajuste (modelo Lognormal)")) %>%
  tab_source_note(source_note = md("Autor: Grupo 4 - Minas"))
Tabla N°2
Resumen del test de bondad de ajuste (modelo Lognormal)
Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Emisiones de CO2 (ton) 98.87 6.06 9.49
Autor: Grupo 4 - Minas

Con Chi² = 6.06 por debajo del umbral de aceptación 9.49 (95% de confianza, gl = 4), el modelo Lognormal se acepta.

6 Cálculo de Probabilidades

¿Cuál es la probabilidad de que un estado productor tenga emisiones de CO2 entre el primer y el tercer invervalo de la distribución observada?

q1 <- round(quantile(EMIS_CO2_edo, 0.25), 0)
q3 <- round(quantile(EMIS_CO2_edo, 0.75), 0)

Probabilidad_1 <- (plnorm(q3, meanlog_1, sdlog_1) - plnorm(q1, meanlog_1, sdlog_1)) * 100
Probabilidad_1
##      75% 
## 41.96187
x <- seq(min(EMIS_CO2_edo), max(EMIS_CO2_edo), length.out = 300)
plot(x,
     dlnorm(x, meanlog_1, sdlog_1),
     col = "skyblue3",
     lwd = 1,
     type = "l",
     main = "Gráfica N°4: Cálculo de probabilidades",
     ylab = "Densidad de probabilidad",
     xlab = "Emisiones de CO2 (toneladas)",
     yaxt = "n")

ejeY <- axTicks(2)
axis(2, at = ejeY, labels = formatC(ejeY, format = "e", digits = 1))

x_section <- seq(q1, q3, length.out = 200)
y_section <- dlnorm(x_section, meanlog_1, sdlog_1)

lines(x_section, y_section, col = "red", lwd = 2)
polygon(c(x_section, rev(x_section)),
        c(y_section, rep(0, length(y_section))),
        col = rgb(1, 0, 0, 0.6))

legend("topright",
       legend = c("Modelo Lognormal", "Área de Probabilidad"),
       col = c("skyblue3", "red"),
       lwd = 2,
       cex = 0.7)

texto_prob <- paste0("Probabilidad = ", round(Probabilidad_1, 2), " %")
text(x = median(EMIS_CO2_edo), y = max(dlnorm(x, meanlog_1, sdlog_1))*0.9,
     labels = texto_prob,
     col = "black",
     cex = 0.8,
     font = 2)

# De 41 estados productores, cuantos tendrian emisiones en ese rango
cantidad_1 <- (plnorm(q3, meanlog_1, sdlog_1) - plnorm(q1, meanlog_1, sdlog_1)) * n_edo
cantidad_1
##      75% 
## 17.20437

7 Intervalos de Confianza

x_media <- mean(EMIS_CO2_edo)
x_media
## [1] 47571283
sigma <- sd(EMIS_CO2_edo)
sigma
## [1] 39154783
n_ic <- length(EMIS_CO2_edo)
n_ic
## [1] 41
e <- (sigma/sqrt(n_ic)) * qt(0.975, n_ic - 1)
e
## [1] 12358775
li <- x_media - e
li
## [1] 35212508
ls <- x_media + e
ls
## [1] 59930057
tabla_media <- data.frame(
  round(li, 2),
  Variable,
  round(ls, 2),
  round(e, 2)
)

colnames(tabla_media) <- c("Límite inferior", "Media poblacional", "Límite superior", "Error estándar")

tabla_media %>%
  gt() %>%
  tab_header(title = md("**Tabla N°3**"),
             subtitle = md("Intervalo de confianza de la media poblacional (95%)")) %>%
  tab_source_note(source_note = md("Autor: Grupo 4 - Minas"))
Tabla N°3
Intervalo de confianza de la media poblacional (95%)
Límite inferior Media poblacional Límite superior Error estándar
35212508 Emisiones de CO2 (ton) 59930057 12358775
Autor: Grupo 4 - Minas

8 Conclusiones

Las emisiones de CO2 se registran a nivel estatal, por lo que la muestra real e independiente corresponde a n = 41 estados (no a las 2921 minas que heredan ese valor). Un primer modelo Normal fue descartado tras rechazar el test de Chi-cuadrado, debido a la fuerte asimetria generada por un estado con emision extrema. Bajo el criterio de muestra por estado, el modelo Lognormal(meanlog = 17.42, sdlog = 0.72) ajusta adecuadamente a la distribucion de EMIS_CO2, con una correlacion de Pearson de 98.87% y un Chi-cuadrado de 6.06, por debajo del umbral de aceptacion de 9.49 (95% de confianza). La probabilidad de que un estado productor emita entre 27,632,119 y 61,869,277 toneladas de CO2 es de 41.96%, equivalente a aproximadamente 17 de los 41 estados productores. Finalmente, la media poblacional de emisiones se estima, con 95% de confianza, entre 35,212,508 y 59,930,057 toneladas.