Analisis Estadistico sobre el cambio en el indidice de calidad de aire en Puerto Rico

Introducción

El presente estudio tiene como objetivo analizar la evolución del AQI en varios municipios de Puerto Rico durante el período 2020–2024, utilizando técnicas estadísticas abordadas en el curso. Para ello, se emplean métodos de regresión lineal y análisis de series de tiempo con el fin de evaluar diferencias entre municipios, identificar tendencias temporales y explorar la estructura dinámica del índice a nivel mensual. Los datos originales fueron procesados mediante agregación mensual e interpolación de valores faltantes, garantizando así la continuidad temporal necesaria para los análisis realizados.

Metodologia

Descripción de los datos

El estudio se realizó utilizando una base de datos de Índice de Calidad del Aire (AQI) correspondiente al período 2020–2024, con observaciones diarias para distintos municipios de Puerto Rico. Las variables principales incluyen la fecha de medición, el municipio (county Name) y el valor del AQI, que mide la calidad del aire y permite evaluar posibles riesgos para la salud.

Como primer paso, la variable de fecha fue transformada al formato adecuado (Date) y posteriormente agregada a nivel mensual, con el fin de reducir la variabilidad diaria y facilitar tanto el análisis comparativo entre municipios como el análisis de series de tiempo.

Durante la exploración inicial se identificaron valores faltantes en la variable AQI, los cuales fueron tratados mediante interpolación lineal, utilizando el método na.approx. Este procedimiento permite mantener la continuidad temporal de la serie sin introducir cambios abruptos artificiales en los datos.

Finalmente, se calculó el AQI promedio mensual por municipio, generando una base de datos agregada (data_mensual) que sirvió como insumo para los modelos de regresión y las series de tiempo.

Metodo de analisis

La metodología se dividió en dos enfoques principales:

análisis transversal mediante modelos de regresión lineal
análisis temporal mediante modelos de series de tiempo (ARIMA).

Librerías Utilizadas

library(zoo)
library(dplyr)
library(lubridate)
library(plotly)
library(ggplot2)
library(dynlm)
library(tseries)
library(forecast)
library(broom)
library(lmtest)

Preparación de Datos

data <- (read_xlsx("AQI 2020-2024.xlsx",sheet = 3))
data$Date <- as.Date(data$Date)

data <- data %>%
  mutate(YearMonth = floor_date(Date, "month"))

data$AQI <- na.approx(data$AQI, maxgap = Inf) #codigo para poder interpolar los datos faltantes

data_mensual <- data %>%
  group_by(`county Name`, YearMonth) %>%
  summarise(AQI = mean(AQI, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    Year  = year(YearMonth),
    Month = month(YearMonth)
  )

## `summarise()` has grouped output by 'county Name'. You can override using the
## `.groups` argument.

Modelo RLM

modelorlm<- lm(
  AQI ~ Year + Month + `county Name`,
  data = data_mensual
)

summary(modelorlm)

## 
## Call:
## lm(formula = AQI ~ Year + Month + `county Name`, data = data_mensual)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -35.558  -8.015  -1.463   7.228  57.774 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            663.02791 1185.39587   0.559 0.576376    
## Year                    -0.31666    0.58625  -0.540 0.589523    
## Month                   -0.01077    0.23953  -0.045 0.964156    
## `county Name`Caguas     -6.36678    2.63095  -2.420 0.016148 *  
## `county Name`Catano     23.41974    2.60791   8.980  < 2e-16 ***
## `county Name`Mayagnez    8.70337    2.60776   3.337 0.000958 ***
## `county Name`Ponce       3.10488    2.63095   1.180 0.238930    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.22 on 285 degrees of freedom
## Multiple R-squared:  0.3404, Adjusted R-squared:  0.3265 
## F-statistic: 24.51 on 6 and 285 DF,  p-value: < 2.2e-16

El modelo de regresión lineal múltiple con todos los predictores, nos indica que no existe evidencia estadísticamente significativa de una tendencia temporal anual ni de estacionalidad mensual en el AQI durante el periodo analizado. Sin embargo, se observan diferencias significativas entre municipios. En particular, Cataño y Mayagüez presentan niveles de AQI significativamente mayores que Bayamón, mientras que Caguas presenta valores significativamente menores. El modelo explica aproximadamente el 33% de la variabilidad del AQI y resulta globalmente significativo.

#Modelos de nulo y saturado
formulalm <- AQI~Year + Month + `county Name`

m0_log <- lm(AQI~ 1, data = data_mensual)   # sin predictores
mF_log <- lm(formulalm, data = data_mensual)        # todos los predictores

# Hacia adelante 
m1_step_log <- step(m0_log,
                    scope     = list(lower = ~1, upper = formula(mF_log)),
                    direction = "forward",
                    k         = log(nrow(data_mensual)),
                    trace     = TRUE)

## Start:  AIC=1670.54
## AQI ~ 1
## 
##                 Df Sum of Sq   RSS    AIC
## + `county Name`  4   29695.0 57713 1572.0
## <none>                       87408 1670.5
## + Year           1      67.3 87341 1676.0
## + Month          1      11.5 87397 1676.2
## 
## Step:  AIC=1572.04
## AQI ~ `county Name`
## 
##         Df Sum of Sq   RSS    AIC
## <none>               57713 1572.0
## + Year   1    59.035 57654 1577.4
## + Month  1     0.425 57713 1577.7

#Modelo Final
data_mensual$`county Name` <- as.factor(data_mensual$`county Name`)
data_mensual$`county Name` <- relevel(
  data_mensual$`county Name`,
  ref = "Catano"
)
mod.final <-lm(AQI~`county Name`,data= data_mensual)
summary(mod.final)

## 
## Call:
## lm(formula = AQI ~ `county Name`, data = data_mensual)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -35.172  -8.341  -1.297   7.039  57.178 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             46.089      1.846  24.965  < 2e-16 ***
## `county Name`Bayamon   -23.410      2.600  -9.004  < 2e-16 ***
## `county Name`Caguas    -29.794      2.634 -11.313  < 2e-16 ***
## `county Name`Mayagnez  -14.701      2.611  -5.631 4.27e-08 ***
## `county Name`Ponce     -20.323      2.634  -7.716 1.99e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.18 on 287 degrees of freedom
## Multiple R-squared:  0.3397, Adjusted R-squared:  0.3305 
## F-statistic: 36.92 on 4 and 287 DF,  p-value: < 2.2e-16

Modelo final luego de seleccion de variables hacia adelanete:El modelo final de regresión, obtenido luego de seleccionar las variables significativas, indica que existen diferencias en el AQI promedio mensual entre los municipios analizados. Cataño presenta el AQI promedio más alto, mientras que Bayamón, Caguas, Mayagüez y Ponce presentan valores significativamente menores. El municipio explica aproximadamente el 34% de la variabilidad del AQI, lo que sugiere que la localización geográfica es un factor clave en la calidad del aire.

Analisis de los Supuestos del modelo

df  <- data.frame(
  yhat = fitted.values(mod.final),
  res  = rstandard(mod.final))

ggplot(df, aes(sample = res)) +
  stat_qq(color = "blue") +
  stat_qq_line(linewidth = 1) +  
  labs(x = "Cuantiles teóricos", y = "Cuantiles muestrales") +
  theme_minimal(base_size = 14)

shapiro.test(df$res)

## 
##  Shapiro-Wilk normality test
## 
## data:  df$res
## W = 0.96068, p-value = 4.181e-07

Prueba de normalidad: Aunque la prueba de Shapiro–Wilk sugiere que los residuos no son perfectamente normales, el análisis gráfico indica que la desviación no es sustancial ya que a simple vista, se puede ver como los datos siguen la línea bastante bien. Esto se puede deber a el tamaño de la muestra y la robustez del modelo de regresión lineal.

#Varianza Constante
ggplot(df, aes(x = yhat, y = res)) +
  geom_point(alpha = 0.6, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey40") +
  labs(x = "Valores ajustados", y = "Residuales estandarizados") +
  theme_minimal(base_size = 14)

bptest(mod.final)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod.final
## BP = 27.43, df = 4, p-value = 1.627e-05

Prueba de varianza: El gráfico de residuos estandarizados frente a los valores ajustados muestra que los errores se distribuyen de forma aproximadamente aleatoria alrededor de cero, sin patrones sistemáticos evidentes, lo que sugiere que el supuesto de linealidad se cumple. No obstante, se observa una ligera mayor dispersión de los residuos para valores ajustados más altos, lo que podría indicar una leve heterocedasticidad, aunque no severa.

#Independencia 
df1 <- data.frame(
  res   =  rstandard(mod.final)) %>%
  mutate(orden = 1:length(res))   


ggplot(df1, aes(x = orden, y = res)) +
  geom_point(alpha = 0.6, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey40") +
  labs(x = "Orden/tiempo", y = "Residuales estandarizados") +
  theme_minimal(base_size = 14)

dwtest(mod.final)

## 
##  Durbin-Watson test
## 
## data:  mod.final
## DW = 0.80827, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

Prueba de independencia: El test de Durbin–Watson indicó la presencia de autocorrelación positiva en los residuos (DW = 0.81, p < 0.001), lo que sugiere que los errores no son independientes. Esta violación del supuesto es consistente con la naturaleza temporal de los datos mensuales de calidad del aire.

Analisis de Series de Tiempo (Dividido por Municipio)

Caguas

caguas_m <- data %>%
  filter(`county Name` == "Caguas") %>%
  mutate(YearMonth = floor_date(Date, "month")) %>%
  group_by(YearMonth) %>%
  summarise(AQI = mean(AQI, na.rm = TRUE))

Serie de Tiempo

ts_caguas <- ts(caguas_m$AQI, start = c(2020,1), end = c(2024,12), frequency = 12)
plot(ts_caguas)

La serie presenta fluctuaciones marcadas a lo largo del periodo 2020–2025, con picos elevados alrededor de 2021 y 2023. A partir de mediados de 2023 se observa una disminución en los valores y un comportamiento más estable hacia 2024–2025. No se identifica una estacionalidad clara, y el patrón sugiere posibles cambios en el nivel de la serie.

Descomposición

d_caguas <- decompose(ts_caguas)
plot(d_caguas)

La tendencia muestra un empeoramiento inicial de la calidad del aire hasta alrededor de 2021 (valores más altos de AQI), seguido de una mejoría gradual desde 2022 hasta 2024, evidenciada por la disminución de los valores. La componente estacional presenta un patrón recurrente y estable, indicando variaciones estacionales consistentes en la calidad del aire. El componente aleatorio no presenta una estructura clara, lo que sugiere que la tendencia y la estacionalidad explican adecuadamente el comportamiento sistemático del AQI.

Prueba de Estacionariedad

adf.test(ts_caguas)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_caguas
## Dickey-Fuller = -3.5917, Lag order = 3, p-value = 0.0415
## alternative hypothesis: stationary

El resultado de la prueba ADF arroja un p-value = 0.0415, menor que 0.05, por lo que se rechaza la hipótesis nula de no estacionariedad. Esto indica que la serie de tiempo del AQI en Caguas es estacionaria, lo que la hace adecuada para la aplicación de modelos ARMA sin necesidad de diferenciación adicional.

ACF y PACF

acf(ts_caguas, lag.max = 36, main = "ACF - Serie original")

pacf(ts_caguas, lag.max = 36, main = "PACF - Serie Caguas")

La ACF muestra autocorrelaciones significativas en los primeros rezagos, lo que indica dependencia temporal en los valores del AQI. Posteriormente, las autocorrelaciones disminuyen y oscilan alrededor de cero, sugiriendo que la dependencia se disipa con el tiempo. Este patrón es consistente con una serie estacionaria y apoya el uso de modelos autorregresivos y/o de media móvil para modelar el comportamiento del AQI.
La PACF muestra un rezago significativo en el primer periodo, mientras que los rezagos posteriores no presentan valores estadísticamente significativos. Este patrón sugiere que la dependencia temporal del AQI puede explicarse principalmente por un componente autorregresivo de orden 1 (AR(1)), siendo una especificación razonable para el modelaje de la serie.

Pruebas de Autoregresión

AR1.caguas <- dynlm(ts_caguas ~ L(ts_caguas, 1))
AR2.caguas <- dynlm(ts_caguas ~ L(ts_caguas, 1:2))

summary(AR1.caguas); summary(AR2.caguas)

## 
## Time series regression with "ts" data:
## Start = 2020(2), End = 2024(12)
## 
## Call:
## dynlm(formula = ts_caguas ~ L(ts_caguas, 1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.5320  -3.6080  -0.6528   1.3234  20.6816 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       5.6929     1.8793   3.029  0.00368 ** 
## L(ts_caguas, 1)   0.6466     0.1014   6.374 3.48e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.929 on 57 degrees of freedom
## Multiple R-squared:  0.4162, Adjusted R-squared:  0.4059 
## F-statistic: 40.63 on 1 and 57 DF,  p-value: 3.484e-08

## 
## Time series regression with "ts" data:
## Start = 2020(3), End = 2024(12)
## 
## Call:
## dynlm(formula = ts_caguas ~ L(ts_caguas, 1:2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.8333  -3.5493  -0.5448   2.1534  20.2125 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          7.2364     1.9907   3.635 0.000613 ***
## L(ts_caguas, 1:2)1   0.8259     0.1302   6.342 4.48e-08 ***
## L(ts_caguas, 1:2)2  -0.2760     0.1301  -2.121 0.038449 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.779 on 55 degrees of freedom
## Multiple R-squared:  0.4607, Adjusted R-squared:  0.441 
## F-statistic: 23.49 on 2 and 55 DF,  p-value: 4.231e-08

El coeficiente del rezago 1 es positivo y estadísticamente significativo (β = 0.6466, p < 0.001), lo que indica que el AQI actual depende de manera importante del valor del mes anterior. Esto sugiere persistencia temporal en la calidad del aire: periodos con AQI alto (peor calidad) tienden a ser seguidos por valores también elevados. El modelo explica aproximadamente 41.6% de la variabilidad del AQI.
El modelo AR(2) muestra que el primer rezago tiene un efecto positivo y significativo (β₁ = 0.8259), mientras que el segundo rezago tiene un efecto negativo y significativo (β₂ = −0.2760), indicando un ajuste correctivo después de dos periodos. Este modelo presenta un mejor ajuste que el AR(1), al explicar cerca del 46.1% de la variabilidad del AQI y reducir ligeramente el error estándar de los residuos.
Ambos modelos confirman la existencia de dependencia temporal en el AQI de Caguas, pero el AR(2) captura mejor la dinámica de corto plazo, sugiriendo que los cambios en la calidad del aire no solo dependen del mes inmediato anterior, sino también de ajustes provenientes de periodos previos.

Seleccion de Modelo

ma.caguas   <- arima(ts_caguas, order = c(0,0,1))
ar.caguas  <- arima(ts_caguas, order = c(1,0,0))
arma.caguas <- arima(ts_caguas, order = c(1,0,1))
ma2.caguas   <- arima(ts_caguas, order = c(0,0,2))
ar2.caguas  <- arima(ts_caguas, order = c(2,0,0))
arma2.caguas <- arima(ts_caguas, order = c(2,0,2))
AIC(ma.caguas, ar.caguas, arma.caguas,ma2.caguas,ar2.caguas,arma2.caguas)

##              df      AIC
## ma.caguas     3 409.8418
## ar.caguas     3 406.0589
## arma.caguas   4 404.3293
## ma2.caguas    4 401.2185
## ar2.caguas    4 403.5017
## arma2.caguas  6 404.8518

De acuerdo con el criterio de información de Akaike (AIC), el modelo con menor valor es MA(2) (AIC = 401.22), por lo que se considera el mejor modelo entre los evaluados, al lograr el mejor balance entre ajuste y complejidad. Aunque modelos como AR(2) y ARMA presentan un buen desempeño, sus valores de AIC son mayores, lo que indica un ajuste relativamente inferior para explicar la dinámica del AQI en Caguas.

Prediccion

prediccion.caguas <- forecast(ma2.caguas, h=12) 
autoplot(prediccion.caguas)

El modelo ARIMA(0,0,2) proyecta que el AQI se mantendrá relativamente estable en el corto plazo, sin cambios drásticos en la tendencia. El valor central del pronóstico se mantiene en niveles similares a los observados recientemente, lo que sugiere una continuidad en la calidad del aire en Caguas. Los intervalos de confianza se amplían a medida que avanza el horizonte de pronóstico, reflejando una mayor incertidumbre en las predicciones futuras.

Bayamón

bayamon_m <- data %>%
  filter(`county Name` == "Bayamon") %>%
  mutate(YearMonth = floor_date(Date, "month")) %>%
  group_by(YearMonth) %>%
  summarise(AQI = mean(AQI, na.rm = TRUE))

Serie de Tiempo

ts_bayamon <- ts(bayamon_m$AQI,
                 start = c(2020,1),end = c(2024,12),
                 frequency = 12)
plot(ts_bayamon)

La serie muestra fluctuaciones considerables a lo largo del periodo 2020–2025, con un aumento marcado del AQI alrededor de 2023, lo que indica un episodio de deterioro temporal de la calidad del aire. Posteriormente, los valores disminuyen y se estabilizan parcialmente durante 2024–2025, aunque permanecen en niveles relativamente más altos que en años anteriores, sugiriendo una recuperación incompleta hacia mejores condiciones de calidad del aire.

Descomposicion de la serie

d_bayamon <- decompose(ts_bayamon)
plot(d_bayamon)

La tendencia muestra un empeoramiento progresivo de la calidad del aire desde aproximadamente 2022 hasta 2023, reflejado en el aumento de los valores del AQI, seguido de una leve estabilización hacia 2024. La componente estacional presenta un patrón recurrente y consistente en el tiempo, indicando variaciones estacionales regulares en la calidad del aire. El componente aleatorio no exhibe una estructura marcada, lo que sugiere que la tendencia y la estacionalidad capturan adecuadamente el comportamiento sistemático del AQI.

Prueba de Estacionariedad

adf.test(ts_bayamon)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_bayamon
## Dickey-Fuller = -3.6814, Lag order = 3, p-value = 0.03394
## alternative hypothesis: stationary

El p-value obtenido (0.0339) es menor que 0.05, por lo que se rechaza la hipótesis nula de no estacionariedad. Esto indica que la serie del AQI en Bayamón es estacionaria, lo que permite aplicar modelos ARMA sin necesidad de diferenciación adicional.

ACF y PACF

acf(ts_bayamon, lag.max = 36, main = "ACF - Serie original")

pacf(ts_bayamon, lag.max = 36, main = "PACF - Serie Caguas")

La ACF muestra autocorrelaciones significativas en los primeros rezagos, lo que indica dependencia temporal en los valores del AQI. Esta dependencia se reduce progresivamente y las autocorrelaciones oscilan alrededor de cero en rezagos posteriores, lo cual es consistente con una serie estacionaria y sugiere la presencia de componentes AR y/o MA de bajo orden.
La PACF presenta un rezago inicial significativo, mientras que los rezagos siguientes no muestran efectos estadísticamente relevantes. Este patrón sugiere que la dinámica del AQI en Bayamón puede modelarse adecuadamente con un componente autorregresivo de orden 1 (AR(1)), capturando la dependencia de corto plazo en la calidad del aire.

Pruebas de Autoregresion

AR1.bayamon <- dynlm(ts_bayamon ~ L(ts_bayamon, 1))
AR2.bayamon <- dynlm(ts_bayamon ~ L(ts_bayamon, 1:2))
summary(AR1.bayamon); summary(AR2.bayamon)

## 
## Time series regression with "ts" data:
## Start = 2020(2), End = 2024(12)
## 
## Call:
## dynlm(formula = ts_bayamon ~ L(ts_bayamon, 1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.0841  -4.2185  -0.0818   4.5026  26.0417 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        9.8566     2.6415   3.731 0.000441 ***
## L(ts_bayamon, 1)   0.5594     0.1083   5.164  3.2e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.738 on 57 degrees of freedom
## Multiple R-squared:  0.3188, Adjusted R-squared:  0.3068 
## F-statistic: 26.67 on 1 and 57 DF,  p-value: 3.199e-06

## 
## Time series regression with "ts" data:
## Start = 2020(3), End = 2024(12)
## 
## Call:
## dynlm(formula = ts_bayamon ~ L(ts_bayamon, 1:2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.0058  -4.5903   0.1594   4.1758  25.5248 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          9.08012    2.99150   3.035 0.003667 ** 
## L(ts_bayamon, 1:2)1  0.52195    0.13551   3.852 0.000308 ***
## L(ts_bayamon, 1:2)2  0.07295    0.13337   0.547 0.586588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.852 on 55 degrees of freedom
## Multiple R-squared:  0.3191, Adjusted R-squared:  0.2944 
## F-statistic: 12.89 on 2 and 55 DF,  p-value: 2.566e-05

El coeficiente del rezago 1 es positivo y estadísticamente significativo (β = 0.5594, p < 0.001), lo que indica persistencia temporal en el AQI: meses con peor calidad del aire (AQI más alto) tienden a ser seguidos por valores también elevados. El modelo explica aproximadamente 31.9% de la variabilidad del AQI en Bayamón.
En el modelo AR(2), el primer rezago continúa siendo significativo, mientras que el segundo rezago no es estadísticamente significativo (p = 0.587). Además, el ajuste del modelo no mejora respecto al AR(1), ya que el R² ajustado disminuye ligeramente, lo que sugiere que el segundo rezago no aporta información adicional relevante.
El modelo AR(1) resulta más adecuado para describir la dinámica del AQI en Bayamón, confirmando una dependencia de corto plazo sin evidencia clara de efectos más prolongados en el tiempo.

Seleccion de Modelo

ma.bayamon   <- arima(ts_bayamon, order = c(0,0,1))
ar.bayamon   <- arima(ts_bayamon, order = c(1,0,0))
arma.bayamon <- arima(ts_bayamon, order = c(1,0,1))
ma2.bayamon   <- arima(ts_bayamon, order = c(0,0,2))
ar2.bayamon   <- arima(ts_bayamon, order = c(2,0,0))
arma2.bayamon <- arima(ts_bayamon, order = c(2,0,2))
AIC(ma.bayamon, ar.bayamon, arma.bayamon,ma2.bayamon,ar2.bayamon,arma2.bayamon)

##               df      AIC
## ma.bayamon     3 425.8962
## ar.bayamon     3 420.9369
## arma.bayamon   4 422.4684
## ma2.bayamon    4 426.4465
## ar2.bayamon    4 422.6160
## arma2.bayamon  6 425.0508

Según el criterio de información de Akaike (AIC), el modelo AR(1) presenta el menor valor de AIC (420.94) entre los modelos evaluados, por lo que se considera el modelo más adecuado para describir la dinámica del AQI en Bayamón. Modelos más complejos, como AR(2) o ARMA, no ofrecen una mejora suficiente en el ajuste como para justificar su mayor complejidad.

Prediccion

prediccion.bayamon <- forecast(ar.bayamon, h=12) 
autoplot(prediccion.bayamon)

El modelo ARIMA(1,0,0) proyecta una leve disminución del AQI en el corto plazo, lo que sugiere una mejoría moderada en la calidad del aire en Bayamón. No obstante, los intervalos de confianza son amplios, indicando una incertidumbre considerable en las proyecciones futuras, especialmente a medida que aumenta el horizonte de pronóstico.

Cataño

catano_m <- data %>%
  filter(`county Name` == "Catano") %>%
  mutate(YearMonth = floor_date(Date, "month")) %>%
  group_by(YearMonth) %>%
  summarise(AQI = mean(AQI, na.rm = TRUE))

Serie de Tiempo

ts_catano <- ts(catano_m$AQI, start = c(2020,1), end = c(2024,12), frequency = 12)
plot_catano <- plot(ts_catano)

La serie muestra alta variabilidad a lo largo del periodo 2020–2025, con picos elevados de AQI especialmente entre 2023 y 2024, lo que indica episodios de deterioro significativo de la calidad del aire. Hacia finales del periodo se observa una disminución de los valores, sugiriendo una posible mejoría reciente en la calidad del aire.

Descomposicion

d_catano <- decompose(ts_catano)
plot(d_catano)

La tendencia refleja un empeoramiento progresivo del AQI hasta alrededor de 2024, seguido de una caída hacia finales del periodo. La componente estacional es marcada y recurrente, lo que indica variaciones estacionales importantes en la calidad del aire. El componente aleatorio presenta fluctuaciones, pero sin una estructura dominante adicional.

Prueba de Estacionariedad

adf.test(ts_catano)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_catano
## Dickey-Fuller = -3.3516, Lag order = 3, p-value = 0.07202
## alternative hypothesis: stationary

log.datos.ts.catano  <- log(ts_catano)
dlog.datos.ts.catano <- diff(log.datos.ts.catano)

El p-value obtenido (0.0720) es mayor que 0.05, por lo que no se rechaza la hipótesis nula de no estacionariedad al nivel de significancia del 5%. Esto sugiere que la serie del AQI en Cataño no es claramente estacionaria, por lo que podría requerir transformación o diferenciación antes de aplicar modelos ARIMA.

Prueba de Estacionariedad con serie transformada

adf.test(dlog.datos.ts.catano)

## Warning in adf.test(dlog.datos.ts.catano): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlog.datos.ts.catano
## Dickey-Fuller = -5.0429, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

El p-value obtenido (0.01) es menor que 0.05, por lo que se rechaza la hipótesis nula de no estacionariedad. Esto indica que la serie del AQI en Cataño, luego de aplicar la transformación (log-diferenciación/diferenciación), se vuelve estacionaria, lo que justifica el uso de modelos ARIMA con diferenciación para su análisis y pronóstico.

ADF y PACF

acf(dlog.datos.ts.catano, lag.max = 36, main = "ACF - Serie Cataño")

pacf(dlog.datos.ts.catano, lag.max = 36, main = "PACF - Serie Cataño ")

La ACF muestra que, tras la transformación, la mayoría de las autocorrelaciones se mantienen dentro de los límites de significancia, lo que indica que la dependencia temporal ha sido adecuadamente reducida. Esto sugiere que la serie transformada del AQI en Cataño es aproximadamente estacionaria y que puede modelarse.
La PACF no muestra rezagos claramente significativos, ya que la mayoría de los valores se encuentran dentro de los límites de confianza. Esto sugiere que, una vez transformada la serie, no existe una estructura autorregresiva fuerte, apoyando el uso de modelos ARIMA simples.

Pruebas de Autoregresión

AR1.catano <- dynlm(dlog.datos.ts.catano ~ L(dlog.datos.ts.catano, 1))
AR2.catano <- dynlm(dlog.datos.ts.catano ~ L(dlog.datos.ts.catano, 1:2))
summary(AR1.catano); summary(AR2.catano)

## 
## Time series regression with "ts" data:
## Start = 2020(3), End = 2024(12)
## 
## Call:
## dynlm(formula = dlog.datos.ts.catano ~ L(dlog.datos.ts.catano, 
##     1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.04153 -0.19141  0.01888  0.24977  0.83514 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                -0.00561    0.05059  -0.111    0.912
## L(dlog.datos.ts.catano, 1) -0.02788    0.13962  -0.200    0.842
## 
## Residual standard error: 0.3849 on 56 degrees of freedom
## Multiple R-squared:  0.0007114,  Adjusted R-squared:  -0.01713 
## F-statistic: 0.03987 on 1 and 56 DF,  p-value: 0.8425

## 
## Time series regression with "ts" data:
## Start = 2020(4), End = 2024(12)
## 
## Call:
## dynlm(formula = dlog.datos.ts.catano ~ L(dlog.datos.ts.catano, 
##     1:2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.02675 -0.19443  0.01985  0.26504  0.86203 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)
## (Intercept)                   -0.00731    0.05183  -0.141    0.888
## L(dlog.datos.ts.catano, 1:2)1 -0.03358    0.14291  -0.235    0.815
## L(dlog.datos.ts.catano, 1:2)2 -0.05591    0.14262  -0.392    0.697
## 
## Residual standard error: 0.391 on 54 degrees of freedom
## Multiple R-squared:  0.003876,   Adjusted R-squared:  -0.03302 
## F-statistic: 0.1051 on 2 and 54 DF,  p-value: 0.9004

Ni el intercepto ni el coeficiente del rezago 1 resultan estadísticamente significativos (p > 0.80), lo que indica ausencia de dependencia temporal en la serie transformada del AQI. El valor de R^2 es prácticamente nulo, sugiriendo que el modelo no explica variabilidad relevante.
En el modelo AR(2), ninguno de los rezagos es significativo y el ajuste no mejora respecto al AR(1), con un R^2 ajustado negativo. Esto confirma que no existe una estructura autorregresiva relevante tras la transformación aplicada.

Seleccion de Modelo

ima.catano  <- arima(dlog.datos.ts.catano, order = c(0,1,1))
ari.catano   <- arima(dlog.datos.ts.catano, order = c(1,1,0))
arima.catano <- arima(dlog.datos.ts.catano, order = c(1,1,1))
ima2.catano  <- arima(dlog.datos.ts.catano, order = c(0,1,2))
ari2.catano   <- arima(dlog.datos.ts.catano, order = c(2,1,0))
arima2.catano <- arima(dlog.datos.ts.catano, order = c(2,1,2))
AIC(ima.catano, ari.catano, arima.catano,ima2.catano,ari2.catano,arima2.catano)

##               df      AIC
## ima.catano     2 60.80490
## ari.catano     2 81.21995
## arima.catano   3 62.79980
## ima2.catano    3 62.79937
## ari2.catano    3 77.38081
## arima2.catano  5 64.50649

El modelo IMA(1) presenta el menor valor de AIC (60.80) entre los modelos evaluados, por lo que se considera el modelo más adecuado para describir la dinámica del AQI en Cataño luego de la diferenciación.

Prediccion

prediccion.catano<- forecast(ima.catano, h=12) 
autoplot(prediccion.catano)

El modelo ARIMA(0,1,1) proyecta que las variaciones del AQI se mantendrán cercanas a cero en el corto plazo, lo que indica ausencia de cambios sistemáticos fuertes en la calidad del aire. Los intervalos de confianza amplios reflejan alta incertidumbre, por lo que no se anticipa una tendencia clara de mejora o deterioro sostenido del AQI en el periodo pronosticado.

Mayaguez

Mayaguez_m <- data %>%
  filter(`county Name` == "Mayagnez") %>%
  mutate(YearMonth = floor_date(Date, "month")) %>%
  group_by(YearMonth) %>%
  summarise(AQI = mean(AQI, na.rm = TRUE))

Serie de Tiempo

ts_mayaguez <- ts(Mayaguez_m$AQI, start = c(2020,1), end = c(2024,12), frequency = 12)
plot_mayaguez <- plot(ts_mayaguez)

La serie muestra valores elevados de AQI entre 2020 y 2022, indicando periodos de peor calidad del aire, seguidos por una disminución marcada a partir de 2023, lo que sugiere una mejoría sostenida en la calidad del aire hacia el final del periodo.

Descomposicion

d_mayaguez <- decompose(ts_mayaguez)
plot(d_mayaguez)

La tendencia refleja un deterioro inicial de la calidad del aire hasta alrededor de 2022, seguido de una mejoría clara y continua desde 2023 en adelante. La componente estacional presenta un patrón recurrente, indicando variaciones estacionales consistentes en el AQI. El componente aleatorio no muestra una estructura dominante adicional.

Prueba de Estacionariedad

adf.test(ts_mayaguez)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_mayaguez
## Dickey-Fuller = -2.9852, Lag order = 3, p-value = 0.1766
## alternative hypothesis: stationary

log.datos.ts.mayaguez  <- log(ts_mayaguez)
dlog.datos.ts.mayaguez <- diff(log.datos.ts.mayaguez)

El p-value obtenido (0.1766) es mayor que 0.05, por lo que no se rechaza la hipótesis nula de no estacionariedad. Esto indica que la serie del AQI en Mayagüez no es estacionaria, por lo que será necesario aplicar diferenciación o transformación antes de utilizar modelos ARIMA.

Prueba de estacionariedad con serie transformada

adf.test(dlog.datos.ts.mayaguez)

## Warning in adf.test(dlog.datos.ts.mayaguez): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlog.datos.ts.mayaguez
## Dickey-Fuller = -4.8527, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

El p-value obtenido (0.01) es menor que 0.05, por lo que se rechaza la hipótesis nula de no estacionariedad. Esto indica que, tras aplicar la transformación (log-diferenciación/diferenciación), la serie del AQI en Mayagüez se vuelve estacionaria, permitiendo el uso adecuado de modelos ARIMA para su análisis y pronóstico.

ACF y PACF

acf(ts_mayaguez, lag.max = 36, main = "ACF - Serie Mayaguez")

pacf(ts_mayaguez, lag.max = 36, main = "PACF - Serie Mayaguez")

La ACF presenta autocorrelaciones significativas en los primeros rezagos que disminuyen gradualmente, lo que indica dependencia temporal de corto plazo.
La PACF muestra un rezago inicial claramente significativo, mientras que los rezagos posteriores no presentan efectos dominantes. Esto sugiere que, aunque existe dependencia inicial, no hay una estructura autorregresiva fuerte.

Pruebas de Autoregresion

AR1.mayaguez <- dynlm(dlog.datos.ts.mayaguez ~ L(dlog.datos.ts.mayaguez, 1))
AR2.mayaguez <- dynlm(dlog.datos.ts.mayaguez ~ L(dlog.datos.ts.mayaguez, 1:2))
summary(AR1.mayaguez); summary(AR2.mayaguez)

## 
## Time series regression with "ts" data:
## Start = 2020(3), End = 2024(12)
## 
## Call:
## dynlm(formula = dlog.datos.ts.mayaguez ~ L(dlog.datos.ts.mayaguez, 
##     1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.04929 -0.24513  0.03042  0.21840  1.38498 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                  -0.00413    0.06136  -0.067    0.947
## L(dlog.datos.ts.mayaguez, 1) -0.05642    0.13360  -0.422    0.674
## 
## Residual standard error: 0.4673 on 56 degrees of freedom
## Multiple R-squared:  0.003174,   Adjusted R-squared:  -0.01463 
## F-statistic: 0.1783 on 1 and 56 DF,  p-value: 0.6744

## 
## Time series regression with "ts" data:
## Start = 2020(4), End = 2024(12)
## 
## Call:
## dynlm(formula = dlog.datos.ts.mayaguez ~ L(dlog.datos.ts.mayaguez, 
##     1:2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.02636 -0.27426  0.03899  0.23721  1.14723 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                     -0.01370    0.06023  -0.227   0.8209  
## L(dlog.datos.ts.mayaguez, 1:2)1 -0.08024    0.13048  -0.615   0.5411  
## L(dlog.datos.ts.mayaguez, 1:2)2 -0.26563    0.13020  -2.040   0.0462 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4546 on 54 degrees of freedom
## Multiple R-squared:  0.07552,    Adjusted R-squared:  0.04128 
## F-statistic: 2.206 on 2 and 54 DF,  p-value: 0.12

Ni el intercepto ni el coeficiente del rezago 1 resultan estadísticamente significativos (p = 0.674), lo que indica que no existe dependencia temporal relevante de corto plazo en las variaciones del AQI tras la transformación. El R2 es prácticamente nulo, sugiriendo un ajuste muy limitado del modelo.
En el modelo AR(2), el segundo rezago resulta negativo y marginalmente significativo (p = 0.046), lo que sugiere un efecto correctivo a dos periodos en las variaciones del AQI. No obstante, el poder explicativo del modelo sigue siendo bajo (R2 ajustado ≈ 0.04), por lo que la estructura autorregresiva es débil.
Luego de la transformación, la dinámica del AQI en Mayagüez muestra poca dependencia autorregresiva, con evidencia limitada de ajustes de corto plazo.

Seleccion de Modelo

ima.mayaguez  <- arima(dlog.datos.ts.mayaguez, order = c(0,1,1))
ari.mayaguez  <- arima(dlog.datos.ts.mayaguez, order = c(1,1,0))
arima.mayaguez <- arima(dlog.datos.ts.mayaguez, order = c(1,1,1))
ima2.mayaguez  <- arima(dlog.datos.ts.mayaguez, order = c(0,1,2))
ari2.mayaguez  <- arima(dlog.datos.ts.mayaguez, order = c(2,1,0))
arima2.mayaguez <- arima(dlog.datos.ts.mayaguez, order = c(2,1,2))
AIC(ima.mayaguez, ari.mayaguez, arima.mayaguez,ima2.mayaguez,ari2.mayaguez,arima2.mayaguez)

##                 df       AIC
## ima.mayaguez     2  82.80429
## ari.mayaguez     2 111.53462
## arima.mayaguez   3  84.71496
## ima2.mayaguez    3  84.62231
## ari2.mayaguez    3 101.85716
## arima2.mayaguez  5  79.48979

el modelo ARIMA(2) es el que presenta el menor valor de AIC (79.49) entre los modelos evaluados, por lo que se considera el modelo más adecuado para describir la dinámica del AQI en Mayagüez luego de la transformación.

Prediccion

prediccion.mayaguez<- forecast(arima2.mayaguez, h=12) 
autoplot(prediccion.mayaguez)

El modelo ARIMA(2,1,2) proyecta que las variaciones del AQI se mantendrán cercanas a cero en el corto plazo, sin evidencia de una tendencia clara de deterioro o mejoría sostenida en la calidad del aire.

Ponce

Ponce_m <- data %>%
  filter(`county Name` == "Ponce") %>%
  mutate(YearMonth = floor_date(Date, "month")) %>%
  group_by(YearMonth) %>%
  summarise(AQI = mean(AQI, na.rm = TRUE))

serie de tiempo

ts_ponce <- ts(Ponce_m$AQI, start = c(2020,1), end = c(2024,12), frequency = 12)
plot_ponce <- plot(ts_ponce)

La serie presenta variabilidad moderada con picos de AQI más altos entre 2020 y 2022, indicando episodios de peor calidad del aire. A partir de 2023 se observa una disminución gradual de los valores, sugiriendo una mejoría reciente en la calidad del aire en Ponce.

descomposicion

d_ponce <- decompose(ts_ponce)
plot(d_ponce)

La tendencia muestra un ligero empeoramiento inicial seguido de una mejoría progresiva hacia el final del periodo. La componente estacional es clara y recurrente, lo que indica variaciones estacionales consistentes en el AQI. El componente aleatorio no presenta una estructura dominante adicional.

Prueba de Estacionariedad

adf.test(ts_ponce)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_ponce
## Dickey-Fuller = -3.9734, Lag order = 3, p-value = 0.01686
## alternative hypothesis: stationary

El p-value obtenido (0.0169) es menor que 0.05, por lo que se rechaza la hipótesis nula de no estacionariedad. Esto indica que la serie del AQI en Ponce es estacionaria.

ACF y PACF

acf(ts_ponce, lag.max = 36, main = "ACF - Serie Ponce")

pacf(ts_ponce, lag.max = 36, main = "PACF - Serie Ponce")

La ACF muestra autocorrelaciones significativas en los primeros rezagos que disminuyen gradualmente, lo que indica dependencia temporal de corto plazo en el AQI.
La PACF presenta un rezago inicial significativo, mientras que los rezagos posteriores no muestran efectos dominantes.

Prueba de Autoregresion

AR1.ponce <- dynlm(ts_ponce ~ L(ts_ponce, 1))
AR2.ponce <- dynlm(ts_ponce ~ L(ts_ponce, 1:2))
summary(AR1.ponce); summary(AR2.ponce)

## 
## Time series regression with "ts" data:
## Start = 2020(2), End = 2024(12)
## 
## Call:
## dynlm(formula = ts_ponce ~ L(ts_ponce, 1))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.905  -8.322  -1.915   5.128  45.166 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     14.9138     3.4341   4.343 5.84e-05 ***
## L(ts_ponce, 1)   0.4162     0.1212   3.433  0.00112 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.05 on 57 degrees of freedom
## Multiple R-squared:  0.1713, Adjusted R-squared:  0.1568 
## F-statistic: 11.78 on 1 and 57 DF,  p-value: 0.001119

## 
## Time series regression with "ts" data:
## Start = 2020(3), End = 2024(12)
## 
## Call:
## dynlm(formula = ts_ponce ~ L(ts_ponce, 1:2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.800  -7.671  -1.016   4.951  45.394 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)        10.9539     3.9134   2.799  0.00705 **
## L(ts_ponce, 1:2)1   0.3103     0.1315   2.360  0.02183 * 
## L(ts_ponce, 1:2)2   0.2578     0.1313   1.964  0.05464 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.86 on 55 degrees of freedom
## Multiple R-squared:  0.2271, Adjusted R-squared:  0.199 
## F-statistic: 8.079 on 2 and 55 DF,  p-value: 0.0008389

El coeficiente del rezago 1 es positivo y estadísticamente significativo (β = 0.4162, p = 0.0011), lo que indica persistencia temporal en el AQI: meses con peor calidad del aire tienden a ser seguidos por valores también elevados. El modelo explica alrededor del 17% de la variabilidad del AQI en Ponce.
En el modelo AR(2), el primer rezago es significativo y el segundo rezago es marginalmente significativo, lo que sugiere una dependencia que se extiende hasta dos periodos. Este modelo mejora el ajuste respecto al AR(1), explicando aproximadamente 23% de la variabilidad del AQI y reduciendo el error estándar de los residuos.

Seleccion de Modelo

ma.ponce  <- arima(ts_ponce, order = c(0,0,1))
ar.ponce  <- arima(ts_ponce, order = c(1,0,0))
arma.ponce <- arima(ts_ponce, order = c(1,0,1))
ma2.ponce  <- arima(ts_ponce, order = c(0,0,2))
ar2.ponce  <- arima(ts_ponce, order = c(2,0,0))
arma2.ponce <- arima(ts_ponce, order = c(2,0,2))
AIC(ma.ponce, ar.ponce, arma.ponce,ma2.ponce,ar2.ponce,arma2.ponce)

##             df      AIC
## ma.ponce     3 465.9470
## ar.ponce     3 461.8253
## arma.ponce   4 460.6637
## ma2.ponce    4 461.9600
## ar2.ponce    4 459.8888
## arma2.ponce  6 463.6903

el modelo AR(2) presenta el menor valor de AIC (459.89) entre los modelos evaluados, por lo que se considera el modelo más adecuado para describir la dinámica del AQI en Ponce.

Prediccion

prediccion.ponce<- forecast(ar2.ponce, h=12) 
autoplot(prediccion.ponce)

El modelo ARIMA(2,0,0) proyecta que el AQI se mantendrá relativamente estable en el corto plazo, con una leve tendencia al aumento respecto a los valores más recientes, lo que podría indicar un ligero deterioro de la calidad del aire.

Conclusion

El análisis del AQI en Puerto Rico durante 2020–2024 muestra que la calidad del aire varía significativamente entre municipios, siendo Cataño y Mayagüez los de mayor exposición, mientras que Caguas y Bayamón presentan valores menores. Los modelos de regresión lineal indicaron que la localización geográfica explica alrededor del 34% de la variabilidad del AQI. Por su parte, el análisis de series de tiempo reveló dependencia temporal en la mayoría de los municipios, con modelos ARIMA capaces de proyectar tendencias estables a corto plazo. Estos hallazgos sugieren que tanto factores espaciales como temporales influyen en la calidad del aire, proporcionando información útil para políticas ambientales y estrategias de mitigación específicas por municipio.

Evaluación estadística y modelaje predictivo del Índice de Calidad del Aire (AQI) en Puerto Rico (2020–2024)

Julio Clemente Vélez, Eleazar Ramos Vélez

2025-12-16

Analisis Estadistico sobre el cambio en el indidice de calidad de aire en Puerto Rico

Introducción

Metodologia

Metodo de analisis

Librerías Utilizadas

Preparación de Datos

Modelo RLM

Analisis de los Supuestos del modelo

Analisis de Series de Tiempo (Dividido por Municipio)

Caguas

Serie de Tiempo

Descomposición

Prueba de Estacionariedad

ACF y PACF

Pruebas de Autoregresión

Seleccion de Modelo

Prediccion

Bayamón

Serie de Tiempo

Descomposicion de la serie

Prueba de Estacionariedad

ACF y PACF

Pruebas de Autoregresion

Seleccion de Modelo

Prediccion

Cataño

Serie de Tiempo

Descomposicion

Prueba de Estacionariedad

Prueba de Estacionariedad con serie transformada

ADF y PACF

Pruebas de Autoregresión

Seleccion de Modelo

Prediccion

Mayaguez

Serie de Tiempo

Descomposicion

Prueba de Estacionariedad

Prueba de estacionariedad con serie transformada

ACF y PACF

Pruebas de Autoregresion

Seleccion de Modelo

Prediccion

Ponce

serie de tiempo

descomposicion

Prueba de Estacionariedad

ACF y PACF

Prueba de Autoregresion

Seleccion de Modelo

Prediccion

Conclusion