1- Librerías usadas
setwd("~/Universidad en general/4to semestre/Estadistica aplicada/ea9am")
library(pacman)
p_load("DT","xfun","ggplot2", "readr", "readxl", "corrplot","scales","anytime", "tidyverse", "lubridate", "forecast", "TTR", "MLmetrics", "tseries", "fpp", "TSstudio","class","e1071","tidyr","caret","mlbench","knitr","MASS")
p_load(gridExtra,leaflet,TSstudio, dplyr, psych, GGally, corrplot)- DT: Permite mostrar matrices en tablas html.
- xfun: Permite embeber dentro del documento archivos como pueden ser los documentos de datos o del mismo proyecto.
- ggplot2: Nos permite generar gráficas.
- readr: Permite leer archivos con extensión “.csv”. *corrplot: Nos proporciona una herramienta de exploración visual en la matriz de correlación que admite el reordenamiento automático de variables para ayudar a detectar patrones ocultos entre variables.
- scales: Nos permite personalizar la apariencia de las etiquetas de ejes y leyendas. Y nos permite controlar como se generan las rupturas a partir de los límites.
- anytime: Nos permite convertir una entrada en cualquier tipo de carácter.
- tidyverse: Proporciona funciones clave de transformación de datos.
- lubridate: Nos permite manipular de manera más facil datos de fechas y horas.
- forecast: Nos ayuda a pronosticar a partir de series temporales.
- TTR: Proporciona funciones de análisis técnico más populares para los datos del mercado financiero.
- MLmetrics: Es una colección de métricas de evaluación.
- tseries: Nos ayuda con el análisis de series temporales y finanzas computacionales.
- fpp: Nos proporciona los conjuntos de datos necesarios para realizar pronósticos.
- TSstudio: Proporciona un conjunto de herramientas para el análisis descriptivo y predictivo de datos de series temporales.
- dplyr - Nos otorga la manipulación de datos que proporciona un conjunto consistente de verbos que lo ayudan a resolver los desafíos de manipulación de datos más comunes.
2.- Datos usados
datos <- read_excel("Concentracion.xlsx")
mov <- read_excel("Concentracion_Mov.xlsx")Como podemos ver en la siguiente tabla se encuentran contenidos los datos que estaremos utilizando a lo largo del análisis. Como vemos la primeras tres columnas son los datos relacionados con los tres principales gases más contaminantes los cuales son:
* Ozono (O3)
* Dióxido de azufre (SO2)
* PM10.Las demás columnas representan los datos de la movilidad que hubo de la población entre esos sectores.
datatable(datos)3.- Introducción
La contaminación del aire es un problema frecuente que podemos observar en cualquier parte del mundo. Esta tiene efectos que perjudican a las personas y al planeta. La contaminación del aire esta compuesta por partículas sólidas y gases que se encuentran en el aire. Los principales contaminantes en el aire se han generado gracias a la mano humana.
Este estudio tuvo en cuenta los siguientes contaminantes:
* Material particulado (PM10)
* Ozono (O3)
* Dióxido de azufre (SO2)
Figura 1. Calidad del aire
El ozono (O3), un gas, es un componente fundamental de la contaminación del aire en las ciudades. Cuando el ozono forma la contaminación del aire también se denomina smog.
Se denomina PM10 a pequeñas partículas sólidas o líquidas de polvo, cenizas, hollín, partículas metálicas, cemento o polen, dispersas en la atmósfera
El dióxido de azufre (SO²), es un gas presente en la atmósfera como resultado de la combinación de azufre y oxígeno y su vida media es de dos a cuatro días.
Antecedentes
El municipio de Hermosillo inicia de manera formal el Programa de Monitoreo de la Calidad del Aire en el año 2004, los aparatos de muestreo son donados al municipio por SEMARNAT y son reinstalados en sitios céntricos pese al crecimiento urbano que ha rebasado en kilómetros las zonas marginales de diez años atrás (Barajas Olvera, 2007).
En Hermosillo sonora, la contaminación del aire es elevada y la población desconoce acerca de este hecho. Existen antecedentes para esta ciudad de medición de PST (Partículas suspendidas totales similares a las PM10) de 1990 a 1995, reportándose que todos esos años se rebasó el máximo permisible anual de 75 μg/m^3 con promedios anuales que fluctuaban de 126 hasta 565 μg/m^3.
4- Objetivos
Lo que se busca con este caso de estudio es ver y analizar el comportamiento de los datos con respecto a la calidad del aire y la movilidad urbana en estos años de pandemia.
Se quiere determinar si la hipótesis de que el decremento de la movilidad urbana hizo que la calidad del aire mejorará significativamente. Para esto se requiere poder ver si existe una relación determinante entre estas variables, también poder describirlas y realizar predicciones de las mismas.
Y para explicar todo esto se utilizarán las técnicas de análisis de relación, regresión lineal múltiple, análisis de series de tiempo, análisis de correlación, y EDA.
5.- Marco teórico
La contaminación del aire representa un importante riesgo medioambiental para la salud, tanto en los países desarrollados como en los países en desarrollo.
La OMS estima que aproximadamente el 58% de las muertes prematuras relacionadas con la contaminación atmosférica que se produjeron en 2016 se debieron a cardiopatías isquémicas y accidentes cerebrovasculares, mientras que el 18% de las muertes se debieron a enfermedades pulmonares obstructivas crónicas e infecciones respiratorias agudas, y el 6% de las muertes se debieron al cáncer de pulmón. Con los diferentes cambios climáticos experimentados a lo largo de los años el tema del medio ambiente ha sido un tema que empezó hacerse más popular entre las bocas de todos los ciudadanos del mundo, además de que gracias a la propagación de esta información a través de los medios y de distintas campañas la población se ha percatado de los daños que ocasionan con distintas acciones de su día a día, creando conciencia de ello y haciendo que cambien sus hábitos.
Los principales gases contaminantes y de los cuales serán de los que informaremos dentro de este análisis son: Material particulado (Pm10), Ozono (O3) y Dioxido de azufre(SO2).
Nuestra fuente de datos para la realización de este análisis provienen de la estación de calidad del aire de la ERNO ubicado en Hermosillo. Enlace a la fuente de datos: http://www.erno.geologia.unam.mx/
6- Métodos
EDA (Exploratory Data Analysis)
El análisis exploratorio de datos es un proceso de investigación que utiliza estadísticas resumidas y herramientas gráficas para aprender sobre los datos y comprender qué se puede aprender de ellos.
Con EDA, puede encontrar incoherencias en sus datos, como valores atípicos o valores atípicos, detectar patrones, comprender las posibles relaciones entre las variables y crear patrones. Las preguntas o hipótesis interesantes se pueden probar más adelante utilizando métodos estadísticos más formales.
Analisis de correlación
El análisis de correlación es un procedimiento estadístico para determinar si dos variables están correlacionadas. El resultado del análisis es un coeficiente de correlación que puede tomar valores de -1 a 1. Indique el tipo de correlación entre dos variables. Un signo positivo indica una relación positiva entre las dos variables; Es decir, a medida que aumenta el tamaño de uno, aumenta el otro. Un signo negativo indica una relación negativa entre las dos variables. A medida que aumenta el valor de una variable, disminuye el valor de la segunda variable. Si las dos variables son independientes, entonces el coeficiente de correlación tiene una magnitud de cero. La fuerza de la relación lineal aumenta a medida que el coeficiente de correlación se aproxima a -1 o 1.
Regresión lineal multiple
La regresión lineal múltiple permite generar un modelo lineal en el que el valor de la variable dependiente o respuesta (Y) se determina a partir de un conjunto de variables independientes llamadas predictores (X1, X2, X3…).
Es una extensión de la regresión lineal simple, por lo que es fundamental comprender esta última. Los modelos de regresión múltiple pueden emplearse para predecir el valor de la variable dependiente o para evaluar la influencia que tienen los predictores sobre ella (esto último se debe que analizar con cautela para no malinterpretar causa-efecto).
Los modelos lineales múltiples siguen la siguiente ecuación:
\[ Y_{i}=(\beta_{0}+\beta_{1}X_{1i}+\beta_{2}X_{2i}+\cdots+\beta_{n}X_{ni})+e_{i} \]
$ _{0} $ : es la ordenada en el origen, el valor de la variable dependiente Y cuando todos los predictores son cero.
$ _{i} $ : es el efecto promedio que tiene el incremento en una unidad de la variable predictora Xi sobre la variable dependiente Y manteniéndose constantes el resto de variables. Se conocen como coeficientes parciales de regresión.
$ e_{i} $ : es el residuo o error, la diferencia entre el valor observado y el estimado por el modelo.
TSA (Time Series Analysis)
Una serie temporal es una serie de datos -en el sentido estadístico: variables- dispuestos en orden cronológico con algún tipo de puntero que indica cuándo corresponde esta observación. Por lo general, las series se organizan en predicados espaciados uniformemente, en este caso son intervalos iguales de series de tiempo. Ambas cosas. Las series de tiempo también pueden contener intervalos de tiempo desiguales, cuando el evento registrado no ocurre durante un período determinado. Tomemos, por ejemplo, un registro de exportación trimestral, con datos para cada trimestre, incluso si es cero, y las exportaciones se registran en años, meses, días, horas y segundos, idealmente cada mercancía vendida a agentes extranjeros.
Vecino Cercano
El análisis de series de tiempo tiene muchas aplicaciones en diversos campos de la ciencia. Por ejemplo, en la economía continuamente se esta expuesto a observaciones de los mercados financieros, indicadores de empleo, índices o indicadores del nivel de producción, índices de precios, etc. En otros campos de las ciencias sociales se emplea el análisis de series de tiempo para analizar la evolución de la población, los nacimientos, o el número de personas con matriculas escolares. Finalmente, en las ciencias exactas se pueden encontrar casos como los de un epidemiólogo que puede estar interesado en el número de casos de influenza observados en algún período de tiempo dado y si a estos se les puede asociar con algún tipo de estacionalidad. La primera aproximación que se suele tener a las series de tiempo es mediante el examen de datos puestos en una gráfica, en la cual uno de los ejes es el tiempo. No obstante, en este tipo de exa ́menes existen dos enfoques. Por un lado, existe el enfoque de la importancia del tiempo, el cual consiste en reconocer como lo que sucede hoy es afectado por lo que paso ayer – o, en general, en periodos pasados–, o como lo que pasa hoy afectara los eventos futuros. Por otro lado, existe el enfoque del análisis frecuentista o de frecuencia, mediante el cual se busca reconocer la importancia que tiene para los investigadores los ciclos: estacionales, de crisis económicas, etc.
Introducción al k Nearest Neighbors (kNN) en R
El algoritmo kNN es muy sencilla: me guardo la tabla de datos de entrenamiento y cuando me llegue un nuevo dato, encuentro los k observaciones (vecinos) más cercanos y hago la clasificación en base a esas observaciones. Al fin y al cabo, es de esperar que observaciones cercanas sean similares a la nueva observación.
Como ves, aquí vemos una gran diferencia respecto a la mayoría de algoritmos supervisados, y es que se trata de un algoritmo no paramétrico. Es decir, que el algoritmo no debe aprender el valor de ningún parámetro, por lo que no hay un proceso de entrenamiento como tal.
Así pues, la clave del algoritmo kNN que programaremos en R se basa en tres aspectos clave que debemos conocer:
Conocer las distintas medidas de distancia que existen, cómo funcionan y cuándo usar cada una de las medidas.
Entender cómo elegir la cantidad de k vecinos a los que se debe observar.
Conocer cómo hace el algortimo kNN las predicciones.
Así pues, si te parece, vamos a ir viendo cada uno de estos aspectos.
Medidas de distancia que puede usar el algoritmo kNN Dentro del algoritmo kNN las medidas de distancia más utilizadas son: distancia Euclídea, distancia de Minkowski, distancia Manhattan, distancia de Coseno y distancia Jaccard. Estas no son las únicas, el algoritmo kNN puede usar cualquier otra medida de distancia, aunque con estas cubriremos la gran mayoría de casos.
Distancia Euclídea
La distancia Euclídea es algo que ya hemos visto en este blog al programar el algoritmo K-means tanto en R como en Python. La distancia Euclídea se basa en el teorema de Pitágoras, según el cual, la hipotenusa al cuadrado es igual a la suma de catetos al cuadrado.
7- Resultados y discusión
Análisis de Relación
A continuacion se muestra la relación entre el PM10 y la movilidad de supermercados y farmacias
El análisis de relación es aquel que se encarga de indentificar una posible relación entre variables. Es una forma básica de lo que es el análisis de datos y según el tipo de variable que sea puede ser llamado de asosiación o correlación.
En este caso se mostrará la relación entre el PM10 y la movilidad de los supermercados y farmacias.
gcov <- ggplot(data = mov) +
geom_line(aes(Fecha, PM10, colour="PM10")) +
geom_line(aes(Fecha, Supermercado_Farmacia, colour="Movilidad")) +
xlab("Fecha") +
ylab("Cantidad") +
labs(colour="casos")+
ggtitle("PM10 y la movilidad de supermercados y farmacias") +
scale_y_continuous(labels = comma)
gcov
Gráfica 1. Análisis de Relación entre el PM10 y la movilidad de supermercados y farmacias
Ya viendo la gráfica nos podemos percatar de que pareciera que existe cierta relación entre el PM10 y la movilidad de farmacias, puesto que cuando la movilidad sube el PM10 también tiende a subir.
Análisis de Correlación
El análisis de correlación se utiliza para lograr determinar la relación que existe entre las variables cuantitativas. Este análisis consiste en un procedimiento para poder determinar si dos varioables están relacionadas o no. El resultado que arroja dicho análisis es un coeficiente de correlación que puede tomar valores entre -1 y +1. El signo es el indicador de que tipo de correlación es el que existe entre las variables.
Estos tipos de análisis de correlación se conocen como: * Análisis de correlación positiva. * Análisis de correlación negativa.
A continuación se realizará la correlación de los datos sobre la calidad del aire, donde se analizarán los datos de O3, SO2, PM10, Reactivacion_Comercial, Supermercado_Farmacia, Parques_Centros, Estaciones_Tramsotp, Lugares_Trabajo y Residencia, para así poder observar como se relaciona cada dato con otro y observar cuales tienen una relación mas fuerte.
aire <- as.data.frame(datos)
round(cor(x = aire, method="pearson"), 3)## O3 SO2 PM10 Reactivacion_Comercial
## O3 1.000 -0.220 -0.229 -0.193
## SO2 -0.220 1.000 0.612 0.404
## PM10 -0.229 0.612 1.000 0.281
## Reactivacion_Comercial -0.193 0.404 0.281 1.000
## Supermercado_Farmacia -0.047 0.443 0.280 0.875
## Parques_Centros -0.182 0.019 0.065 0.814
## Estaciones_Transito -0.174 0.399 0.270 0.944
## Lugares_Trabajo -0.119 0.138 0.100 0.581
## Residencia 0.130 -0.371 -0.238 -0.842
## Supermercado_Farmacia Parques_Centros
## O3 -0.047 -0.182
## SO2 0.443 0.019
## PM10 0.280 0.065
## Reactivacion_Comercial 0.875 0.814
## Supermercado_Farmacia 1.000 0.655
## Parques_Centros 0.655 1.000
## Estaciones_Transito 0.872 0.770
## Lugares_Trabajo 0.492 0.371
## Residencia -0.729 -0.617
## Estaciones_Transito Lugares_Trabajo Residencia
## O3 -0.174 -0.119 0.130
## SO2 0.399 0.138 -0.371
## PM10 0.270 0.100 -0.238
## Reactivacion_Comercial 0.944 0.581 -0.842
## Supermercado_Farmacia 0.872 0.492 -0.729
## Parques_Centros 0.770 0.371 -0.617
## Estaciones_Transito 1.000 0.474 -0.758
## Lugares_Trabajo 0.474 1.000 -0.853
## Residencia -0.758 -0.853 1.000Análisis con histogramas
Los histogramas son gráficos que indican la frecuencia de un hecho mediante una distribución de los datos. Por lo tanto para tener un mejor enfoque de las variables se procedera a realizar este analis debido a que estos utilizan variables para su elaboración.
multi.hist(x = aire, dcol = c("red", "green"), dlty = c("dotted", "solid"))
Gráfica 5. Histograma de cada dato.
Aquí se puede observar graficamente como son los datos de las variables por separado.
Análisis de Dispersión de los datos
Mediante este grafico se puede observar como estan dispersos los datos, y la correlación que tiene cada variable.
ggpairs(datos, lower = list(continuous ="smooth"),
diag = list (continuos = "barDiag"), axisLabels = "none")
Gráfica 6. Gráfico de dispersión de los datos.
Se puede observar como hay varias relaciones, como las del PM10 y SO2 con el Supermercado y Farmacias, esto es debido que proximo al lugar donde se tomaron las muestras de aire se encontraba una farmacia, lo que facilito la obtención de datos y que se le adjudicara a las farmacias.
Regresion lineal multiple
Un modelo de regresión linea multiple es un modelo estadístico que nos permite generar un modelo linea en el que el valor de la variable dependiente se determina según un conjunto de variables que son independientes las cuales son llamdas predictores. Es una extensión de lo que es la regresión lineal simple.
Modelo general de correlación
Se revisará la relación que tiene el PM10 con los supermercados y farmacias.
m <- lm(SO2 ~ Supermercado_Farmacia, data=datos)
summary(m)##
## Call:
## lm(formula = SO2 ~ Supermercado_Farmacia, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.77447 -0.16040 -0.00167 0.15632 0.98634
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.285796 0.015300 18.68 <2e-16 ***
## Supermercado_Farmacia 0.011437 0.001098 10.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2544 on 445 degrees of freedom
## Multiple R-squared: 0.1961, Adjusted R-squared: 0.1943
## F-statistic: 108.5 on 1 and 445 DF, p-value: < 2.2e-16Se puede predecir el SO2 mediante los supermercados y farmacias con un 19% de exactitud, lo cual es bajo, pero, el p-value es menor a 0.05 lo que nos da a saber que las variables estan bastante relacionadas.
Modelo de regresión lineal múltiple.
Se evaluaran las variables predictoras para dar un primer modelo para saber cuales son mas significativas para el SO2.
modelo <- lm(SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia + Parques_Centros + Estaciones_Transito + Lugares_Trabajo + Residencia, data = datos )
summary(modelo)##
## Call:
## lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Estaciones_Transito + Lugares_Trabajo +
## Residencia, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.51190 -0.09921 -0.00269 0.09870 0.52906
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0187954 0.0380324 -0.494 0.621416
## O3 -0.0057616 0.0009832 -5.860 9.10e-09 ***
## PM10 0.0056324 0.0005408 10.415 < 2e-16 ***
## Reactivacion_Comercial 0.0096304 0.0020097 4.792 2.27e-06 ***
## Supermercado_Farmacia 0.0060624 0.0015668 3.869 0.000126 ***
## Parques_Centros -0.0169489 0.0010860 -15.607 < 2e-16 ***
## Estaciones_Transito -0.0018481 0.0013556 -1.363 0.173465
## Lugares_Trabajo -0.0106817 0.0009837 -10.858 < 2e-16 ***
## Residencia -0.0352798 0.0041747 -8.451 4.32e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1582 on 438 degrees of freedom
## Multiple R-squared: 0.6938, Adjusted R-squared: 0.6882
## F-statistic: 124.1 on 8 and 438 DF, p-value: < 2.2e-16Habiendo introducido todas las variables como predictores el modelo cuenta con un R2 de 0.6882, por lo que es capaz de explicar el 68.82% de la variabilidad observada en el SO2.
Evaluación del modelo
Se realizará una evaluación del modelo anteriormente visto, esto mediante el Akaike(AIC), el cual es una medida de la calidad relativa de un modelo estadístico, para un conjunto dado de dato. Este se basa en la entropía de información y ofrece una estimación de la información perdida cuando se utiliza un modelo determinado.
step(object = modelo, direction = "both", trace = 1)## Start: AIC=-1639.29
## SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Estaciones_Transito + Lugares_Trabajo +
## Residencia
##
## Df Sum of Sq RSS AIC
## - Estaciones_Transito 1 0.0465 11.014 -1639.4
## <none> 10.968 -1639.3
## - Supermercado_Farmacia 1 0.3749 11.343 -1626.3
## - Reactivacion_Comercial 1 0.5750 11.543 -1618.5
## - O3 1 0.8600 11.828 -1607.5
## - Residencia 1 1.7883 12.756 -1573.8
## - PM10 1 2.7161 13.684 -1542.4
## - Lugares_Trabajo 1 2.9523 13.920 -1534.7
## - Parques_Centros 1 6.0997 17.067 -1443.6
##
## Step: AIC=-1639.4
## SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Lugares_Trabajo + Residencia
##
## Df Sum of Sq RSS AIC
## <none> 11.014 -1639.4
## + Estaciones_Transito 1 0.0465 10.968 -1639.3
## - Supermercado_Farmacia 1 0.3298 11.344 -1628.2
## - Reactivacion_Comercial 1 0.6091 11.623 -1617.3
## - O3 1 0.8366 11.851 -1608.7
## - Residencia 1 1.7647 12.779 -1575.0
## - PM10 1 2.7569 13.771 -1541.5
## - Lugares_Trabajo 1 2.9217 13.936 -1536.2
## - Parques_Centros 1 6.0965 17.111 -1444.5##
## Call:
## lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Lugares_Trabajo + Residencia, data = datos)
##
## Coefficients:
## (Intercept) O3 PM10
## -0.020842 -0.005669 0.005668
## Reactivacion_Comercial Supermercado_Farmacia Parques_Centros
## 0.008020 0.005439 -0.016944
## Lugares_Trabajo Residencia
## -0.010417 -0.035005Despues de realizar el proceso de AIC, nos da que el mejor modelo es
- (formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia + Parques_Centros + Lugares_Trabajo + Residencia, data = datos)
Mejor Modelo
El modelo ajustado nos elimina la variable Estaciones_Transito.
modeloAIC <- (lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
Parques_Centros + Lugares_Trabajo + Residencia, data = datos))
summary(modeloAIC)##
## Call:
## lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
## Parques_Centros + Lugares_Trabajo + Residencia, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5458 -0.1006 -0.0022 0.0964 0.5374
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0208416 0.0380399 -0.548 0.584047
## O3 -0.0056694 0.0009818 -5.775 1.46e-08 ***
## PM10 0.0056679 0.0005407 10.482 < 2e-16 ***
## Reactivacion_Comercial 0.0080203 0.0016277 4.927 1.18e-06 ***
## Supermercado_Farmacia 0.0054386 0.0015000 3.626 0.000322 ***
## Parques_Centros -0.0169444 0.0010870 -15.588 < 2e-16 ***
## Lugares_Trabajo -0.0104169 0.0009653 -10.791 < 2e-16 ***
## Residencia -0.0350045 0.0041739 -8.387 6.89e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1584 on 439 degrees of freedom
## Multiple R-squared: 0.6925, Adjusted R-squared: 0.6876
## F-statistic: 141.2 on 7 and 439 DF, p-value: < 2.2e-16El modelo ajustado nos da un %68 de exactitud.
Intervalo de Confianza de los coeficientes
En la siguiente tabla podemos observar los intervalos de confianza de los coeficientes:
confint(lm(formula = SO2 ~ O3 + PM10 + Reactivacion_Comercial + Supermercado_Farmacia +
Parques_Centros + Lugares_Trabajo + Residencia, data = datos))## 2.5 % 97.5 %
## (Intercept) -0.095604583 0.053921343
## O3 -0.007599002 -0.003739801
## PM10 0.004605203 0.006730588
## Reactivacion_Comercial 0.004821182 0.011219376
## Supermercado_Farmacia 0.002490524 0.008386774
## Parques_Centros -0.019080816 -0.014808020
## Lugares_Trabajo -0.012314119 -0.008519681
## Residencia -0.043207808 -0.026801239Se puede observar un intervalo de confianza del 97.5%.
Validación de condiciones para la regresión múltiple lineal.
plot1 <- ggplot(data = datos, aes(O3, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = datos, aes(PM10, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = datos, aes(Reactivacion_Comercial, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot4 <- ggplot(data = datos, aes(Supermercado_Farmacia, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot5 <- ggplot(data = datos, aes(Parques_Centros, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot6 <- ggplot(data = datos, aes(Lugares_Trabajo, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot7 <- ggplot(data = datos, aes(Residencia, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2, plot3, plot4, plot5, plot6, plot7)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Gráfica 7. Gráfico de residuales.
El SO2 tiene una distribucion con relacion a las demas variables anormal con un 97.5%,
Analizando la parte de Supermercado_Farmacia podemos observar que sus datos se encuentran dispersos, ademas de que no cuenta con una linea regular donde podemos observar picos por ejemplo en -20.
Variabilidad constante de los residuos
Para demostrar la variabilidad de los residuos se grafica una matriz donde se puede observar que tan relacionados estan las variables.
corrplot(cor(dplyr::select(datos, O3, PM10, Reactivacion_Comercial, Supermercado_Farmacia, Parques_Centros, Lugares_Trabajo, Residencia)),
method = "number", tl.col = "black")
Gráfica 7. Matríz de correlación.
Esta matriz nos muestra las relaciones que existen entre cada una de las variable. El color rojo nos indica que tiene un mayor impacto cuando se correlacionan.
TSA (Time Series Analysis)
Time Series Analysis es una manera especifica de analizar una secuencia de datos recolectados en un intervalo de tiempo.En TSA los analistas registran puntos de datos a intervalos constantes durante un período de tiempo determinado en lugar de simplemente registrar los puntos de datos de forma intermitente o aleatoria.
anyNA(mov)## [1] FALSEEsto nos dice que no hay ningun dato faltante
Analisís de datos y descomposición de la serie de tiempo
A continuación se realizará un analisís de la seríe de tiempo tomando como variables el PM10.
qual <- ts(mov$PM10, start = 2017, end = 2028 ,frequency = 2)
qual## Time Series:
## Start = c(2017, 1)
## End = c(2028, 1)
## Frequency = 2
## [1] 39.20583 33.75667 72.35500 59.46083 74.44708 42.33542 62.24417 30.70417
## [9] 22.72583 29.09708 37.50750 21.02875 40.74500 50.33583 38.33417 32.31750
## [17] 28.50333 12.30250 26.29667 22.18625 31.20792 33.08542 9.11250autoplot(qual)
Gráfica 8. Descomposición de serie de tiempo.
La descomposición de series de tiempo se usa para identificar componentes, reensamblando las partes para construir un pronóstico.
qualdc <- qual %>%
decompose(type="multiplicative") %>%
autoplot()
qualdc
Gráfica 9. Descomposición de serie de tiempo multiplicativa.
cada uno de los componentes representa una de las categorías subyacentes de los patrones.
Validación cruzada
La validación cruzada es una técnica para ajustar hiperparámetros y producir medidas sólidas del rendimiento del modelo.
Aquí definiremos un 20% para pruebas y un 80% para entrenamiento.
test_qual <- tail(qual, 3) #20% para pruebas
train_qual <- head(qual, length(qual)-3) #80% para entrenamientoAjuste del modelo con Holt-Winters.
Se usará el metodo Holt-Winters para el ajuste del modelo, quees una ampliación perfeccionada del enfoque de la suavización exponencial,
qual_ses <- HoltWinters(x=train_qual, seasonal = "multiplicative")
qual_ses## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = train_qual, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha: 0.868512
## beta : 0.2310532
## gamma: 1
##
## Coefficients:
## [,1]
## a 31.8137400
## b -0.4439777
## s1 0.9381583
## s2 0.6973795Forecasting
Forecasting es el proceso de hacer predicciones sobre el futuro mediante el análisis estadístico de tendencias observadas en datos históricos.
Se usará el modelo dado por Holt-Winters para asi poder hacer la predicción
qual_forecast <- forecast(object = qual_ses, h=4)Visualizando resultados
Tenemos 3 variables actualmente, los datos que ya teniamos, los de test y los entrenados.
Teniendo estos en cuenta se realiza una grafica para mostrar como se comporta cada uno.
qual %>%
autoplot(series = "actual") +
autolayer(qual_forecast$fitted, series = "train") +
autolayer(qual_forecast$mean, series = "test") +
theme_minimal()## Warning: Removed 2 row(s) containing missing values (geom_path).
Gráfica 10. Serie de tiempo de comparación entre lo actual, entrenado y el test.
Se puede observar como por momentos la variable entrenada se llega a asemejar a la actual.
Evaluar la precisión del modelo
Esto se hace para saber que tan preciso fue la predicción dada.
eval_ses <- accuracy(qual_forecast, test_qual)
eval_ses## ME RMSE MAE MPE MAPE MASE
## Training set -5.197672 19.26240 15.21923 -17.37667 43.61137 1.0504454
## Test set -2.062579 13.10816 10.92693 -57.76908 84.77694 0.7541868
## ACF1 Theil's U
## Training set -0.05535754 NA
## Test set -0.36689851 0.9559035De acuerdo con la prueba de precision el modelo tiene un error medio de: 43.4%
Utilizando el modelo SARIMA o ARIMA
Los modelos SARIMA (Seasonal Autoregressive Integrated Moving Average Model) captan el comportamiento puramente estacional de una serie.
Empezaremos haciendo una prueba al conjunto de datos
adf.test(qual)##
## Augmented Dickey-Fuller Test
##
## data: qual
## Dickey-Fuller = -2.8785, Lag order = 2, p-value = 0.2377
## alternative hypothesis: stationaryDe acuerdo con el valor de P, los datos son estacionarios porque son mayores a 0.05
Analisis de ajuste
Se realiza el analisis de ajuste y se procede a representarse de manera grafíca
qual_diff <- qual %>%
diff(lag=12) %>%
diff(lag=1)
qual_diff %>%
autoplot()
Gráfica 11. Análisis de ajuste
Analizando la serie de tiempo por medio de un modelo ARIMA
qual_auto <- auto.arima(y = train_qual, seasonal = T)
summary(qual_auto)## Series: train_qual
## ARIMA(0,1,0)
##
## sigma^2 = 312.1: log likelihood = -81.52
## AIC=165.04 AICc=165.28 BIC=165.99
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.8490189 17.21902 14.01056 -13.64256 40.69657 0.9670223
## ACF1
## Training set -0.4942984Mediante el modelo ARIMA se obtuvo un menor porcentaje de error, el cual fue del 40.2%.
Datos reales vs SARIMA
Se mostrará los resultados de ARIMA representados en una grafica
train_qual %>%
autoplot(series= "actual") +
autolayer (qual_auto$fitted, series = "SARIMA auto") +
theme_minimal()
Gráfica 12. Comparación entre datos reales y el metodo SARIMA
Se puede ver que el modelo actual y ARIMA cuentan con bastente similitud, lo que quiere decir que se podría hacer una buena estimación de los datos a futuro a pesar de su margen de error del 40%.
Analisis de los cambios diferenciales en la serie de tiempo
tsdisplay(qual_diff)
Gráfica 13. Analisis de los cambios diferenciales en la serie de tiempo
Aquí se pueden observar los cambios diferenciales de la serie de tiempo.
Vecino cercano
Inicio del proceso del tema Vecino Cercano aplicado a calidad del aire y movilidad urbana
Datos
Datos necesarios para el desarrollo de vecino cercano
datatable(mov)Filtrado
Se busca despreciar aquellas variables que no afectan al procedimiento, se definen tres variables necesarias para aplicar el procedimiento, el SO2, el PM10 y la residencia. ademas de usar la Etiqueta para clasificar si los indices de los gases presentes en la atmosfera fueron antes o durante la pandemia.
movilidadcsv <- read.csv("~/Universidad en general/4to semestre/Estadistica aplicada/ea9am/Concentracion_Mov2.csv", encoding = "UTF-8")
datos.prep <- dplyr::select(movilidadcsv, SO2, PM10, ,Residencia, Etiqueta)
datatable(datos.prep)Construir el modelo KNN
En este modelo, el modelo de vecinos mas cercanos (KNN) se construye con los mismos datos preparados y las columnas numéricas [,1:3] para posteriormente evaluar el modelo mediante los criterios de una matriz de confusión.
modeloKnn <- knn(train = datos.prep[,1:3], test = datos.prep[,1:3], k = 4, cl = datos.prep[,4] )
head(modeloKnn)## [1] Prepandemia Prepandemia Prepandemia Prepandemia Prepandemia Prepandemia
## Levels: Pandemia PrepandemiaResumen del modelo
El resumen de modelo muestra que en los datos encontramos 93 fechas pertenecientes a la pandemia y 32 a la prepandemia
summary(modeloKnn)## Pandemia Prepandemia
## 93 32Evaluando el modelo
Se construye un conjunto de datos llamado datos.r.p con valores reales y valores predichos a partir de los datos preparados incorporados en el algoritmo KNN
datos.r.p <- data.frame(reales = datos.prep$Etiqueta, prediccion = modeloKnn)
datatable(datos.r.p)datos.r.p$reales <- as.factor(datos.r.p$reales)
datos.r.p$prediccion <- as.factor(datos.r.p$prediccion)Crear matriz de confusion
La siguiente matriz de confusion indica que el modelo es capaz de predecir si los gases aumentaron o disminuyeron segun las fechas pertenecian a la pandemia o la prepandemia
matriz <- confusionMatrix(datos.r.p$reales, datos.r.p$prediccion)
matriz## Confusion Matrix and Statistics
##
## Reference
## Prediction Pandemia Prepandemia
## Pandemia 92 0
## Prepandemia 1 32
##
## Accuracy : 0.992
## 95% CI : (0.9562, 0.9998)
## No Information Rate : 0.744
## P-Value [Acc > NIR] : 3.892e-15
##
## Kappa : 0.9792
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9892
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 0.9697
## Prevalence : 0.7440
## Detection Rate : 0.7360
## Detection Prevalence : 0.7360
## Balanced Accuracy : 0.9946
##
## 'Positive' Class : Pandemia
## Calculos para las metricas
Se realizan los respectivos calculos de las metricas derivadas de la matriz y obtenemos una precision para el modelo de 99.2%, siendo esta una precision bastante alta
exactitud <- round(as.numeric(matriz$overall[1]) * 100,2)
sensibilidad <- round(as.numeric(matriz$byClass[1]) * 100,2)
especificidad <- round(as.numeric(matriz$byClass[2]) * 100,2)
precision.Pandemia <- round(as.numeric(matriz$byClass[3]) * 100,2)
precision.NOPandemia <- round(as.numeric(matriz$byClass[4]) * 100,2)
kappa <- round(as.numeric(matriz$overall[2]) * 100,2)
exactitud## [1] 99.2Hacer predicciones
Una vez mencionada la precision, creamos nuevos datos para saber si el modelo es capaz de predecir si dependiendo del aumento o disminucion de los gases, nos encontramos en prepandemia o pandemia.
SO2 <- c(.1058, -.1508,.0510)
PM10 <- c(20,30,15)
Residencia <- c(-1,15,10)
Estado <- c('?','?','?')
datos.nuevos <- data.frame(SO2, PM10, Residencia, Estado )
kable(datos.nuevos, caption = "Datos nuevos")| SO2 | PM10 | Residencia | Estado |
|---|---|---|---|
| 0.1058 | 20 | -1 | ? |
| -0.1508 | 30 | 15 | ? |
| 0.0510 | 15 | 10 | ? |
Resultados del modelo
El modelo muestra una precision correcta para todos los casos que se le adjuntaron de ejemplo
modelo <- knn(train = datos.prep[,1:3], test = datos.nuevos[,1:3], k = 4, cl = datos.prep[,4] )
modelo## [1] Prepandemia Pandemia Pandemia
## Levels: Pandemia Prepandemiadatos.nuevos <- mutate(datos.nuevos, prediccion = modelo)
kable(datos.nuevos, caption = "Predicción de datos nuevos")| SO2 | PM10 | Residencia | Estado | prediccion |
|---|---|---|---|---|
| 0.1058 | 20 | -1 | ? | Prepandemia |
| -0.1508 | 30 | 15 | ? | Pandemia |
| 0.0510 | 15 | 10 | ? | Pandemia |
8.- Conclusiones por tema
A simple vista los datos no parecen tener una relacion, pero analizandolos y teniendo los contextos de estos mismos podemos llegar a notar como ciertas cosas encajan formando asi una relacion entre datos, como por ejemplo la siguiente:
gcov Gráfica 1. Análisis de Relación entre el PM10 y la movilidad de supermercados y farmacias
Como se vio anteriormente, el PM10 y la movilidad de supermercados y farmacias se encuentran relacionadas.
matriz## Confusion Matrix and Statistics
##
## Reference
## Prediction Pandemia Prepandemia
## Pandemia 92 0
## Prepandemia 1 32
##
## Accuracy : 0.992
## 95% CI : (0.9562, 0.9998)
## No Information Rate : 0.744
## P-Value [Acc > NIR] : 3.892e-15
##
## Kappa : 0.9792
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9892
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 0.9697
## Prevalence : 0.7440
## Detection Rate : 0.7360
## Detection Prevalence : 0.7360
## Balanced Accuracy : 0.9946
##
## 'Positive' Class : Pandemia
## En el caso de Vecino cercano, se logro demostrar con precision que los gases SO2, PM10 y la residencia ademas de tener una relacion entre estas, se define claramente que los gases nocivos y la residencia disminuyeron significativamente cuando inicio la pandemia, dejando claro que la movilidad urbana se redujo considerablemente desde que la pandemia inicio, diversos motivos referente a esto existen, pero los datos son claros, se puede afirmar que la pandemia provoco una disminucion clara de la contaminacion.
Variabilidad
Esta grafica nos muestra la relación que existe entre las variables. Lo que se ve rojo nos dice que tienen una relacion mas fuerte, por ejemplo, la relacion entre reactivacion comercial y residencia esta en rojo, esto quiere decir que puede existir una relación entre las variables.
corrplot(cor(dplyr::select(datos, O3, PM10, Reactivacion_Comercial, Supermercado_Farmacia, Parques_Centros, Lugares_Trabajo, Residencia)),
method = "number", tl.col = "black") Gráfica 7. Matríz de correlación.
Entrenamiento en serie de tiempo
Al momento de realizar un entrenamiento mediante un forecast en la serie de tiempo obtuvimos lo siguiente:
qual %>%
autoplot(series = "actual") +
autolayer(qual_forecast$fitted, series = "train") +
autolayer(qual_forecast$mean, series = "test") +
theme_minimal()## Warning: Removed 2 row(s) containing missing values (geom_path). Gráfica 10. Serie de tiempo de comparación entre lo actual, entrenado y el test.
Donde se puede observar como la linea “actual” de la grafica se asemeja a la otra linea “train”, la cual fue entrenada mediante metodos antes vistos.
9- Conclusión
Mediante la realización de esta documento se logro conocer las causas y concecuencias de la contaminación, todo esto mediante un estudio de la calidad del aire, donde se encontraban 3 pincipales variables, que eran O3, PM10 Y S02.
Al realizar este estudio se pudo apreciar como se llegaban a relacionar algunas cosas que a simple vista no tenian nada que ver, o como cosas que pensaba que tenian relación al final no la tenian, por ejemplo, el como las farmacias influian tanto en los datos y relaciones, y esto se debía a que cerca del lugar de donde se tomaban las muestrar del aire, se encontraba una farmacia, lo cual le hacía ser más influyente que otras variables.
Tambien se concluye que gracias a vecino cernano logramos probar esa hipotesis de que realmente cuando inicio la pandemia, los gases contaminantes se redujeron considerablemente, pudiendo afirmar que realmente el planeta “se limpio climaticamente” durante la pandemia.
10- Bibliografía
Gobierno de México. Efectos a la salud por la contaminación del aire ambiente. gob.mx. Recuperado 7 de abril de 2022, de https://www.gob.mx/cofepris/acciones-y-programas/3-efectos-a-la-salud-por-la-contaminacion-del-aire-ambiente#
OMS. (2021, 22 septiembre). Calidad del aire ambiente (exterior) y salud. who.int. Recuperado 7 de abril de 2022, de https://www.who.int/es/news-room/fact-sheets/detail/ambient-(outdoor)-air-quality-and-health
UNAM. Instituto de Geología. Recuperado 7 de abril de 2022, de http://www.erno.geologia.unam.mx/
Análisis de correlación. (2017, 25 mayo). Conogasi. Recuperado 16 de mayo de 2022, de https://conogasi.org/articulos/analisis-de-correlacion-2/#:%7E:text=Descripci%C3%B3n-,El%20an%C3%A1lisis%20de%20correlaci%C3%B3n%20consiste%20en%20un%20procedimiento%20estad%C3%ADstico%20para,correlaci%C3%B3n%20entre%20las%20dos%20variables%20.
Análisis exploratorio de datos. (s. f.). Introducción a la estadística | JMP. Recuperado 16 de mayo de 2022, de https://www.jmp.com/es_pe/statistics-knowledge-portal/exploratory-data-analysis.html#:%7E:text=El%20an%C3%A1lisis%20exploratorio%20de%20datos%20es%20un%20proceso%20de%20investigaci%C3%B3n,se%20puede%20averiguar%20de%20ellos.
Cómo funciona Vecino más cercano promedio—Ayuda | ArcGIS Desktop. (s. f.). ArcMap. Recuperado 16 de mayo de 2022, de https://desktop.arcgis.com/es/arcmap/10.5/tools/spatial-statistics-toolbox/h-how-average-nearest-neighbor-distance-spatial-st.htm
Paladino, M. (s. f.). Series de tiempo en R. Series de tiempo. Recuperado 16 de mayo de 2022, de https://www.institutomora.edu.mx/testU/SitePages/martinpaladino/series_de_tiempo_con_R.html#:%7E:text=En%20R%20expresamos%20una%20serie,asegurarnos%20que%20se%20mantengan%20alienados.
10.- Descarga de código y datos
Descargar Proyecto:
xfun::embed_file("CalidadAire.Rmd")Descargar Datos:
xfun::embed_file("Concentracion.xlsx")xfun::embed_file("Concentracion_Mov.xlsx")Download Concentracion_Mov.xlsx
xfun::embed_file("Concentracion_Mov2.csv")