library(fields)
library(geoR)
library(akima)
library(dplyr)
library(sqldf)
library(tidyr)
library(sp)
library(sf)
library(maps)
library(ggplot2)
library(gstat)
library(knitr)
library(kableExtra)
library(mvtnorm)

Planteamiento del Problema

La calidad del aire es un factor crítico en la salud pública y el bienestar ambiental.
En particular, el ozono troposférico (O\(_3\)) y el dióxido de nitrógeno (NO\(_2\)) son contaminantes atmosféricos regulados bajo los Estándares Nacionales de Calidad del Aire Ambiental (NAAQS, por sus siglas en inglés).

En el presente informe se analizan datos oficiales de la Agencia de Protección Ambiental de Estados Unidos (EPA) correspondientes al estado de Texas en el año 2024.
El problema que se plantea es determinar si los niveles de O\(_3\) y NO\(_2\) cumplen con los límites normativos establecidos, y en qué medida podrían representar un riesgo para la salud en un momento específico del tiempo.

Descripción de las Variables

Se utilizarán dos variables provenientes de la red de monitoreo de calidad del aire en Texas:

Ozono (O\(_3\)):
Variable: Daily Max 8-hour Ozone Concentration, que corresponde al máximo promedio móvil de 8 horas de ozono reportado diariamente por cada estación de monitoreo.
- Unidad: partes por millón (ppm).
- Escala: Razón (el valor cero indica ausencia total de ozono; se dice que una estación excede los niveles de O\(_3\) cuando la concentración máxima registrada en el promedio móvil de 8 horas supera los 0.070 ppm).
Dióxido de Nitrógeno (NO\(_2\)):
Variable: Daily Max 1-hour NO2 Concentration, que corresponde a la concentración máxima de NO\(_2\) registrada en una hora dentro de cada día.
- Unidad: partes por billón (ppb).
- Escala: Razón (el valor cero indica ausencia total de NO\(_2\); se dice que una estación excede los niveles de NO\(_2\) cuando la concentración máxima registrada en una hora del día supera los 100 ppb).

Escala Temporal

Ambas variables se registran diariamente durante el año 2024.
Para efectos de este análisis preliminar, se toma como referencia un momento específico; es decir, desde el 01 de enero de 2024 hasta el 31 de diciembre de 2024.

Mapa de Texas

load("OZ_coor.RData")

load("NO_coor.RData")

texas <- st_read("Texas_State/State_Boundary.shp")

## Reading layer `State_Boundary' from data source 
##   `G:\Mi unidad\Camilo UNAL\Estadística Espacial\Proyecto\Primera Entrega\Texas_State\State_Boundary.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -11871800 ymin: 2978934 xmax: -10409240 ymax: 4369694
## Projected CRS: WGS 84 / Pseudo-Mercator

oz_sf <- st_as_sf(OZ_coor, coords = c("Site.Longitude", "Site.Latitude"), crs = 4326)
no_sf <- st_as_sf(NO_coor, coords = c("Site.Longitude", "Site.Latitude"), crs = 4326)

Se presenta el mapa de Texas con las localizaciones de las estaciones desde las cuales se tomaron los datos.

ggplot() +
  geom_sf(data = texas, fill = "gray95", color = "black") +
  geom_sf(data = oz_sf, color = "blue", shape = 16, size = 3) +
  geom_sf(data = no_sf, color = "red", shape = 17, size = 3) +
  labs(title = "Estaciones de Monitoreo de O3 y NO2 en Texas (2024)",
       caption = "Fuente: EPA") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Modelamiento con el Ozono

Análisis de Estacionariedad

Primero, se importan los datos correspondientes al Ozono en el estado de Texas, Estados Unidos.

OZ <- read.csv("OZ.csv")[,-1]

OZ$Ozone <- 1000*OZ$Ozone

OZ$OzC <- residuals(lm(Ozone ~ 1, data = OZ))

OZg <- as.geodata(OZ,
                  coords.col = 1:2,
                  data.col = 4)

head(OZ)

##    Easting Northing Ozone        OzC
## 1 -1663915  3669343    34 -0.3378378
## 2 -1689317  3670665    37  2.6621622
## 3 -1822422  3509765    37  2.6621622
## 4 -1813926  3522080    39  4.6621622
## 5 -1797108  3475404    38  3.6621622
## 6 -1497689  3445674    30 -4.3378378

par(mfrow = c(2, 1), mar = c(2.75, 3, 0.5, 0.5), mgp = c(1.7, 0.7, 0))

plot(x = OZ$Ozone, y = OZ$Northing,
     xlab = "Ozono", ylab = "Norte")

plot(x = OZ$Easting, y = OZ$Ozone,
     xlab = "Este", ylab = "Ozono")

De estos gráficos es posible ver que la variable presenta estacionariedad con respecto a la media y, además, no se evidencian problemas con respecto a heteroscedasticidad. Por lo tanto, no es necesario ajustar un modelo lineal para intentar alcanzar estacionariedad.

Semivariograma

A pesar de haber ajustado el modelo lineal, el Ozono en Texas parece ser estacionario desde un comienzo.

De esta forma, se halla el semivariograma empírico para estos dos casos.

Datos originales

vargrm1 <- variog(OZg, trend = "cte")

## variog: computing omnidirectional variogram

plot(vargrm1, pch = 19, col = "springgreen3",
     xlab = "Distancia",
     ylab = "Semivarianza",
     main = "Variograma empírico de Ozono (O3) en Texas 2024")

lines(vargrm1$u, vargrm1$v, col = "springgreen3")

Como se espera, a puntos cercanos hay menor variabilidad en cuanto a la semivarianza. Se puede ver que se presenta un efecto pepita pequeño, pues la semivarianza no pasa por el orígen. Además, se puede ver que el rango ocurre en aproximadamente \(400.000\) unidades de distancia.

Resistente a datos atípicos

vg1_at <- variog(OZg, trend = "cte",
                 estimator.type = "modulus")

## variog: computing omnidirectional variogram

plot(vg1_at, pch = 19, col = "springgreen4",
     xlab = "Distancia",
     ylab = "Semivarianza",
     main = "Variograma empírico de Ozono (O3) en Texas 2024
             (Robusto)")

lines(vg1_at$u, vg1_at$v, col = "springgreen4")

Es posible ver que es similar al semivariograma basado en el promedio. Esto se puede estar dando debido a que no hay presencia de datos atípicos que estén afectando el nivel medio de la variable.

Estimación teórica del semivariograma

Eyefit

Con base en la sección de Modelos válidos para el semivariograma, vemos que los que más se asimilan al semivariograma empírico son el Gaussiano y Matern.

Por lo tanto, con ayuda de la función eyefit, inicialmente ajustamos a “ojo” los parámetros del modelo gaussiano que se ajusta a los datos de Ozono:

# initial <- eyefit(vargrm1)
# (silla = 34.88, rango = 9e+05, nugget = 10) Guassiano

# Con el Matern
#cov.model sigmasq       phi tausq kappa kappa2   practicalRange
#   matern    46.5 1036220.5    10   0.7   <NA> 3571724.96606027

initialOZ1 <- c(c0 = 10, cs = 34.88, a = 9e+05)
initialOZ2 <- c(c0 = 10, cs = 46.5, a = 1036220.5, kappa = 0.7)
initialOZ3 <- c(c0 = 10, cs = 28.42, a = 1727035.38)

Así,

\[ \boldsymbol{\theta}^{(0)}_1=( c_0^{(0)}=10, c_s^{(0)}=34.88, a^{(0)}=900000) \quad \text{(Gaussiano)} \\ \boldsymbol{\theta}^{(0)}_2=(c_0^{(0)}=10, c_s^{(0)}=46.5, a^{(0)}=1036220.5, \kappa^{(0)} = 0.7) \quad \text{(Matern)} \]

Entonces, vemos que el semivariograma teórico ajustado se ve de la siguiente forma:

gaussian_model <- function(h, c0, cs, a){
  ifelse(h == 0,
         0,
         c0 + cs*(1 - exp(-(h/a)^2)))
}


matern_model <- function(h, c0, cs, a, kappa){
  ifelse(h == 0,
         0,
         c0 + cs * (1 - ((h / a)^kappa * besselK(h / a, nu = kappa)) /
                    (2^(kappa - 1) * gamma(kappa)))
  )
}

cubic_model <- function(h, c0, cs, a) {
  ifelse(
    h == 0, 
    0,
    ifelse(
      h <= a,  
      c0 + cs * (7*(h/a)^2 - 8.75*(h/a)^3 + 3.5*(h/a)^5 - 0.75*(h/a)^7),
      c0 + cs   
    )
  )
}

fitvar1 <- gaussian_model(h = vargrm1$u, c0 = 10, cs = 34.88, a = 9e5)

fitvar2 <- matern_model(h = vargrm1$u, c0 = 10, cs = 46.5, a = 1036220.5, kappa = 0.7)

fitvar3 <- cubic_model(h = vargrm1$u, c0 = 10, cs = 28.42, a = 1934000)


plot(vargrm1, pch = 19,
     xlab = "Distancia",
     ylab = "Semivarianza",
     main = "Estimación teórica")

lines(x = vargrm1$u,
      y = fitvar1,
      col = "orange2",
      lwd = 2)
lines(x = vargrm1$u,
      y = fitvar2,
      col = "orchid3",
      lwd = 2)
lines(x = vargrm1$u,
      y = fitvar3,
      col = "steelblue2",
      lwd = 2)

legend("topleft",
       legend = c("Gaussiano", "Matern", "Cúbico"),
       col = c("orange2", "orchid3", "steelblue2"),
       lty = 1,
       lwd = 2)

Luego, podemos calcular el Error Cuadrático Medio (MSE) con respecto al semivariograma empírico:

mse_1 <- mean((fitvar1 - vargrm1$v)^2)
mse_2 <- mean((fitvar2 - vargrm1$v)^2)
mse_3 <- mean((fitvar3 - vargrm1$v)^2)

Obteniendo los siguientes resultados:

Comparación de Modelos (MSE)
Modelo	MSE
Gaussiano	22.89337
Matern	29.72176
Cúbico	28.05368

Modelo Gaussiano

Regresión No Lineal

Ahora, con los valores iniciales escogidos anteriormente, realizamos la estimación de los parámetros mediante regresión no lineal.

Para el modelo Gaussiano,

data_OZ <- data.frame(h = vargrm1$u,
                      gamma_est = vargrm1$v,
                      n = vargrm1$n)

nlinear_gauss <- nls(gamma_est ~ c0 + cs*(1 - exp(-(h/a)^2)),
                     data = data_OZ,
                     start = list(c0 = initialOZ1[1], 
                                  cs = initialOZ1[2],
                                  a = initialOZ1[3]))

# Matriz de pesos de Cressie

W0_gauss <- (8*gaussian_model(c0 = initialOZ1[1], 
                             cs = initialOZ1[2],
                             a = initialOZ1[3],
                             h = data_OZ$h))/data_OZ$n 

nlinear_gauss0 <- nls(gamma_est ~ c0 + cs*(1 - exp(-(h/a)^2)),
                      data = data_OZ,
                      start = list(c0 = initialOZ1[1], 
                                   cs = initialOZ1[2],
                                   a = initialOZ1[3]),
                      weights = 1/W0_gauss)


# 1/n

W1_gauss <- data_OZ$n

nlinear_gauss1 <- nls(gamma_est ~ c0 + cs*(1 - exp(-(h/a)^2)),
                      data = data_OZ,
                      start = list(c0 = initialOZ1[1], 
                                   cs = initialOZ1[2],
                                   a = initialOZ1[3]),
                      weights = 1/W1_gauss)

# h/n

W2_gauss <- data_OZ$n/data_OZ$h

nlinear_gauss2 <- nls(gamma_est ~ c0 + cs*(1 - exp(-(h/a)^2)),
                      data = data_OZ,
                      start = list(c0 = initialOZ1[1], 
                                   cs = initialOZ1[2],
                                   a = initialOZ1[3]),
                      weights = 1/W2_gauss)

Y comparamos su MSE:

Regresión no Lineal (Modelo Gaussiano)
Peso	MSE
\(1\)	22.5775
Cressie	34.7773
\(n\)	25.8338
\(\frac{n}{h}\)	31.8803

Además de comparar sus gráficas:

plot(vargrm1,
     xlab = "Distancia", y = "Semivarianza",
     main = "Semivariograma Teórico con el Modelo Gaussiano",
     type = "p", pch = 19)

lines(x = data_OZ$h, y = fitted(nlinear_gauss),
      lty = "solid", lwd = 2, col = "dodgerblue2")
lines(x = data_OZ$h, y = fitted(nlinear_gauss0),
      lty = "dashed", lwd = 2, col = "tomato3")
lines(x = data_OZ$h, y = fitted(nlinear_gauss1),
      lty = "dotdash", lwd = 2, col = "springgreen3")
lines(x = data_OZ$h, y = fitted(nlinear_gauss2),
      lty = "longdash", lwd = 2, col = "darkorchid2")

legend("bottomright",
       legend = c("1", "Cressie", expression(n), expression(frac(n,h))),
       col = c("dodgerblue2", "tomato3", "springgreen3", "darkorchid2"),
       lwd = 2,
       cex = 0.9)

Tanto con la gráfica como con el MSE, nos damos cuenta que el modelo no lineal en el que no se considera una matriz de ponderación se ajusta de mejor forma al semivariograma empírico. Veamos el resúmen de este ajuste:

summary(nlinear_gauss)

## 
## Formula: gamma_est ~ c0 + cs * (1 - exp(-(h/a)^2))
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)   
## c0 9.819e+00  3.220e+00   3.049  0.01226 * 
## cs 3.451e+01  8.195e+00   4.211  0.00180 **
## a  8.577e+05  2.548e+05   3.366  0.00717 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.418 on 10 degrees of freedom
## 
## Number of iterations to convergence: 9 
## Achieved convergence tolerance: 3.892e-06

Y también los coeficientes estimados por este método:

Mejor modelo no lineal
Parámetro	Estimación
\(c_0\)	9.8194
\(c_s\)	34.509
\(a\)	857720.5474

Función optim

Primero, definimos la función que queremos optimizar. En este caso, queremos minimizar la suma de cuadrados entre el semivariograma empírico y el teórico:

fn_gauss <- function(par, h, gamma){
  c0 <- par[1]
  cs <- par[2]
  a <- par[3]
  
  gamma_hat <- gaussian_model(c0, cs, a, h)
  sum((gamma - gamma_hat)^2)
}

Ahora,

optim_gauss <- optim(
   par = c(10, 34.88, 900000),
   fn = fn_gauss,
   h = data_OZ$h,
   gamma = data_OZ$gamma_est,
   method = "L-BFGS-B",
   lower = c(0, 0, 0))

optim_gauss_fit <- optim_gauss$par

Y así encontramos que, con la función optim,

Función optim
Parámetro	Estimación
\(c_0\)	9.5136
\(c_s\)	25.1184
\(a\)	9e+05

Por lo tanto, graficamos el semivariograma teórico con estos parámetros y lo comparamos con el mejor modelo no lineal hallado:

plot(vargrm1,
     xlab = "Distancia", y = "Semivarianza",
     main = "Semivariograma Teórico con el Modelo Gaussiano\n (comparación optim y nls))",
     type = "p", pch = 19)

lines(x = data_OZ$h, y = fitted(nlinear_gauss),
      lty = "solid", lwd = 2, col = "dodgerblue2")

lines(x = data_OZ$h, 
      y = gaussian_model(c0 = optim_gauss_fit[1],
                         cs = optim_gauss_fit[2],
                         a = optim_gauss_fit[3],
                         h = data_OZ$h),
      lty = "dashed", lwd = 2, col = "tomato3")

legend("topleft",
       legend = c("nls", "optim"),
       col = c("dodgerblue2", "tomato3"),
       lwd = 2,
       cex = 0.9)

Comparación nls y optim
Método	MSE
nls	22.5775210133
optim	57.8922389274

Como nos podemos dar cuenta, ambas estimaciones son muy distintas. Las obtenidas con ayuda de la función optim no se ajustaron muy bien al semivariograma empírico.

Por lo tanto, los semivariogramas teóricos que mejor se ajustaron al empírico fueron los de la regresión no lineal (sin considerar ponderaciones y con los pesos \(1/n\)).

MV

Para este método, es necesario obtener la matriz de distancias y también la matriz de covarianza \(\Sigma(\theta)\).

dist_matrix_OZ <- as.matrix(dist(OZ[, c("Easting", "Northing")]))

sigma2_OZ <- sum(initialOZ1[1:2])

Suponiendo que las observaciones provienen de una distribución normal, se quiere maximizar

\[ l(\theta \, |\,Z) \] Realizamos esto con la ayuda de la función optim:

loglik_gaussOZ <- function(par, Z, coords) {
  c0 <- par[1]  
  cs <- par[2]   
  a  <- par[3]   
  
  
  h <- as.matrix(dist(coords))
  
  
  Sigma <- cs * exp(-(h / a)^2)
  diag(Sigma) <- diag(Sigma) + c0
  
  n <- length(Z)
  L <- try(chol(Sigma), silent = TRUE)
  if (inherits(L, "try-error")) return(1e6)
  
  
  SinvZ <- backsolve(L, forwardsolve(t(L), Z))
  
  logdet <- 2 * sum(log(diag(L)))
  nll <- 0.5 * (n * log(2*pi) + logdet + sum(Z * SinvZ))
  
  return(nll)  
}

optim_gaussOZ <- optim(
  par = initialOZ1,
  fn = loglik_gaussOZ,
  Z = OZ$OzC,
  coords = OZ[, c("Easting", "Northing")],
  method = "L-BFGS-B",
  lower = c(0, 0, 1)
)

optim_gauss <- optim_gaussOZ$par

fit_ml_gaussOZ <- gaussian_model(h = data_OZ$h,
                                 c0 = optim_gauss[1],
                                 cs = optim_gauss[2],
                                 a = optim_gauss[3])
plot(vargrm1)
lines(x = data_OZ$h, y = fit_ml_gaussOZ)

mse_4 <- mean((data_OZ$gamma_est - fit_ml_gaussOZ)^2)

Modelo Cúbico

Regresión No Lineal

Se ajustan modelos no lineales con ayuda de la función nls().

nlinear_cubic <- nls(gamma_est ~  c0 + cs * (7*(h/a)^2 - 8.75*(h/a)^3 
                                            + 3.5*(h/a)^5 - 0.75*(h/a)^7),
                     data = data_OZ,
                     start = list(c0 = initialOZ3[1], 
                                  cs = initialOZ3[2],
                                  a = initialOZ3[3]))

# Matriz de pesos de Cressie

W0_cubic <- (8*cubic_model(c0 = initialOZ3[1], 
                           cs = initialOZ3[2],
                           a = initialOZ3[3],
                           h = data_OZ$h))/data_OZ$n 

nlinear_cubic0 <- nls(gamma_est ~ c0 + cs * (7*(h/a)^2 - 8.75*(h/a)^3 
                                            + 3.5*(h/a)^5 - 0.75*(h/a)^7),
                      data = data_OZ,
                      start = list(c0 = initialOZ3[1], 
                                  cs = initialOZ3[2],
                                  a = initialOZ3[3]),
                      weights = 1/W0_cubic)


# 1/n

W1_cubic <- data_OZ$n

nlinear_cubic1 <- nls(gamma_est ~ c0 + cs * (7*(h/a)^2 - 8.75*(h/a)^3 
                                            + 3.5*(h/a)^5 - 0.75*(h/a)^7),
                      data = data_OZ,
                      start = list(c0 = initialOZ3[1], 
                                   cs = initialOZ3[2],
                                   a = initialOZ3[3]),
                      weights = 1/W1_cubic)

# h/n

W2_cubic <- data_OZ$n/data_OZ$h

nlinear_cubic2 <- nls(gamma_est ~ c0 + cs * (7*(h/a)^2 - 8.75*(h/a)^3 
                                            + 3.5*(h/a)^5 - 0.75*(h/a)^7),
                      data = data_OZ,
                      start = list(c0 = initialOZ3[1], 
                                   cs = initialOZ3[2],
                                   a = initialOZ3[3]),
                      weights = 1/W2_cubic)

Regresión no Lineal (Modelo Cúbico)
Peso	MSE
\(1\)	22.9361
Cressie	36.6435
\(n\)	26.1399
\(\frac{n}{h}\)	32.1966

De esta forma y, al igual que con el modelo gaussiano, se tiene que el ajuste con un menor Error Cuadrático Medio (MSE) fue el que no considera ponderaciones.

De todas formas, si lo comparamos con el obtenido con modelo gaussiano, vemos que este último sigue siendo mejor.

plot(vargrm1,
     xlab = "Distancia", y = "Semivarianza",
     main = "Comparación de Semivariogramas Teóricos",
     type = "p", pch = 19)

lines(x = data_OZ$h, y = fitted(nlinear_gauss),
      lty = "solid", lwd = 2, col = "dodgerblue2")
lines(x = data_OZ$h, y = fitted(nlinear_cubic),
      lty = "solid", lwd = 2, col = "tomato3")

legend("bottomright",
       legend = c(
         paste0("Gaussiano (MSE = ", round(temp_nlinear_gauss[1], 4), ")"),
         paste0("Cúbico (MSE = ", round(temp_nlinear_cubic[1], 4), ")")
       ),
       col = c("dodgerblue2", "tomato3"),
       lwd = 2,
       cex = 0.9)

Función optim

Se define la función a optimizar:

fn_cubic <- function(par, h, gamma_est) {
  c0 <- par[1]
  cs <- par[2]
  a  <- par[3]
  
  gamma <- cubic_model(h, c0, cs, a)
  
  mean((gamma_est - gamma)^2)
}

fit_cubic <- optim(
  par = initialOZ3,
  fn = fn_cubic,
  h = data_OZ$h,
  gamma_est = data_OZ$gamma_est,
  method = "L-BFGS-B",  # permite restricciones
  lower = c(0, 0, 0),   # no negativos
  upper = c(Inf, Inf, Inf)
)

optim_cubic <- fit_cubic$par

optim_cubic_svm <- cubic_model(c0 = optim_cubic[1],
                               cs = optim_cubic[2],
                               a = optim_cubic[3],
                               h = data_OZ$h)


temp_gauss <- c(mean((data_OZ$gamma_est - fitted(nlinear_cubic))^2),
                mean((data_OZ$gamma_est - optim_cubic_svm)^2))

temp_peso <- c("nls", "optim")


kable(cbind(temp_peso, round(temp_gauss, 4)),
      col.names = c("Método", "MSE"),
      format = "html",
      align = "c",
      booktabs = TRUE,
      escape = FALSE,
      caption = "Comparación nls y optim (Modelo Cúbico)") %>%
      kable_styling(position = "center",
                    full_width = FALSE)

Comparación nls y optim (Modelo Cúbico)
Método	MSE
nls	22.9361
optim	23.5749

Así, podemos ver que, aunque se aproximan más los MSE, el obtenido por nls es mejor.

MV

loglik_cubicOZ <- function(par, Z, coords) {
  c0 <- par[1]  
  cs <- par[2]   
  a  <- par[3]   
  
  
  h <- as.matrix(dist(coords))
  
  
  Sigma <- cs * (7*(h/a)^2 - 8.75*(h/a)^3 + 3.5*(h/a)^5 - 0.75*(h/a)^7)
  diag(Sigma) <- diag(Sigma) + c0
  
  n <- length(Z)
  L <- try(chol(Sigma), silent = TRUE)
  if (inherits(L, "try-error")) return(1e6)
  
  
  SinvZ <- backsolve(L, forwardsolve(t(L), Z))
  
  logdet <- 2 * sum(log(diag(L)))
  nll <- 0.5 * (n * log(2*pi) + logdet + sum(Z * SinvZ))
  
  return(nll)  
}

optim_cubicOZ <- optim(
  par = initialOZ3,
  fn = loglik_cubicOZ,
  Z = OZ$OzC,
  coords = OZ[, c("Easting", "Northing")],
  method = "L-BFGS-B",
  lower = c(0, 0, 1)
)

ml_cubic <- optim_cubicOZ$par

fit_ml_cubicOZ <- cubic_model(h = data_OZ$h,
                              c0 = ml_cubic[1],
                              cs = ml_cubic[2],
                              a = ml_cubic[3])
plot(vargrm1)
lines(x = data_OZ$h, y = fit_ml_cubicOZ)

mean((data_OZ$gamma_est - fit_ml_cubicOZ)^2)

## [1] 24.53791

Modelamiento con el Dióxido de Carbono

Análisis de Estacionariedad

Se procede a importar los datos relacionados al dióxido de nitrógeno en el estado de Texas.

NO <- read.csv("NO.csv")[,-1]

NOg <- as.geodata(NO)


head(NO)

##    Easting Northing   NO
## 1 -1689317  3670665  5.9
## 2 -1822422  3509765 10.4
## 3 -1813926  3522080  1.0
## 4 -1797108  3475404  8.3
## 5 -1798947  3506464 19.9
## 6 -1497689  3445674 20.1

A forma de exploración, analizamos la relación que hay entre el dióxido de carbono y las coordenadas.

par(mfrow = c(2, 1), mar = c(2.75, 3, 0.5, 0.5), mgp = c(1.7, 0.7, 0))

plot(x = NO$NO, y = NO$Northing,
     xlab = "Dióxido de Nitrógeno", ylab = "Norte")
plot(x = NO$Easting, y = NO$NO,
     xlab = "Este", ylab = "Dióxido de Nitrógeno")

Parece ser estacionario con respecto a la media pues no tiene una tendencia que dependa del espacio.

Semivariograma

Al igual que con la variable de ozono, el dióxido de nitrógeno parece tener un comportamiento estacionario a través del espacio desde un principio.

Datos originales

vgOZ_1 <- variog(NOg, trend = "cte")

## variog: computing omnidirectional variogram

plot(vgOZ_1, pch = 19, col = "#CD4F39",
     xlab = "Distancia",
     ylab = "Semivarianza",
     main = "Variograma empírico de NO2 en Texas")

lines(vgOZ_1$u, vgOZ_1$v, col = "#CD4F39")

Resistente a datos atípicos

vgOZ_at <- variog(NOg, trend = "cte",
                  estimator.type = "modulus")

## variog: computing omnidirectional variogram

plot(vgOZ_at, pch = 19, col = "#912CEE",
     xlab = "Distancia",
     ylab = "Semivarianza",
     main = "Variograma empírico de NO2 en Texas")

lines(vgOZ_at$u, vgOZ_at$v, col = "#912CEE")

Evaluación de la Calidad del Aire en Texas 2024: Ozono (O\(_3\)) y Dióxido de Nitrógeno (NO\(_2\))

Camilo Gómez - Sebastián Olarte

2025-10-21

Planteamiento del Problema

Descripción de las Variables

Escala Temporal

Mapa de Texas

Modelamiento con el Ozono

Análisis de Estacionariedad

Semivariograma

Datos originales

Resistente a datos atípicos

Estimación teórica del semivariograma

Eyefit

Modelo Gaussiano

Regresión No Lineal

Función optim

MV

Modelo Cúbico

Regresión No Lineal

Función optim

MV

Modelamiento con el Dióxido de Carbono

Análisis de Estacionariedad

Semivariograma

Datos originales

Resistente a datos atípicos