Cargamos la libreria sp, disponible para el tratamiento de datos espaciales y la libreria readr para la carga de la base de datos con extención .CSV
library(sp)
library(readr)
data_atlantic <- read.csv("C:/Users/asus/Documents/data_atlantic_1998_2012.csv")
View(data_atlantic)# Observar la base de datos
dim(data_atlantic) #Observar las dimensiones de la base de datos
## [1] 666 9
names(data_atlantic) # Da el nombre de las variables de la base de datos
## [1] "FIPS" "x" "y" "Rate" "POV" "SMOK" "PM25" "NO2" "SO2"
class(data_atlantic) #Tipo de datos de la base de datos ("data.frame")
## [1] "data.frame"
Observamos los primeros seis datos y los ultimos dos
head(data_atlantic) # Permite observar las seis primeras filas de la base de datos
tail(data_atlantic) # Permite observar las seis ultimas filas de la base de datos
Observamos el resumen estadistico de la base de datos
summary(data_atlantic) # Resumen de datos estadisticos
## FIPS x y Rate
## Min. :10001 Min. : 949877 Min. : 932261 Min. : 39.00
## 1st Qu.:24010 1st Qu.:1248238 1st Qu.:1332531 1st Qu.: 62.00
## Median :37126 Median :1425055 Median :1675880 Median : 69.00
## Mean :35728 Mean :1426696 Mean :1669654 Mean : 69.97
## 3rd Qu.:51033 3rd Qu.:1617599 3rd Qu.:1970044 3rd Qu.: 77.00
## Max. :54109 Max. :1916559 Max. :2623675 Max. :120.00
## POV SMOK PM25 NO2
## Min. : 3.26 Min. :13.51 Min. : 7.323 Min. : 0.5954
## 1st Qu.:11.35 1st Qu.:24.57 1st Qu.:10.508 1st Qu.: 1.2608
## Median :14.54 Median :26.52 Median :11.652 Median : 1.7766
## Mean :15.21 Mean :26.06 Mean :11.394 Mean : 2.2015
## 3rd Qu.:18.91 3rd Qu.:28.18 3rd Qu.:12.384 3rd Qu.: 2.6475
## Max. :35.26 Max. :34.96 Max. :14.133 Max. :14.6895
## SO2
## Min. :0.001209
## 1st Qu.:0.027947
## Median :0.049337
## Mean :0.072520
## 3rd Qu.:0.091953
## Max. :0.420847
Conversión a datos espaciales
Convertimos las columnas (x,y) de la base de datos data_atlantic en un objeto espacial
coordinates(data_atlantic) = c("x", "y") # Convierte el data.frame a un objeto espacial de tipo SPACIAPOINTDATAFRAME
class(data_atlantic)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
Diagramamos la base de datos data_atlantic
plot(data_atlantic) # Diagrama la base de datos
plot(data_atlantic, pch=1, main="Dato espacial: Punto") # Diagrama la base de datos con forma de puntos
Tratamiento de los datos espaciales
Diagramamos la distrbución espacial de la variable POV, por medio de la libreria “laticce”
DISTRIBUCIÓN ESPACIAL DE LA VARIABLE POV
library(lattice) #cargar libreria lattice
spplot(data_atlantic, c("POV"))
**Diagramamos la distrbución espacial conlas coordenadas (x,y)
#Distribución espacial con coordenadas (x,y)
library(ggplot2)
methods(fortify)
## [1] fortify.cld* fortify.confint.glht*
## [3] fortify.data.frame* fortify.default*
## [5] fortify.formula* fortify.function*
## [7] fortify.glht* fortify.grouped_df*
## [9] fortify.Line* fortify.Lines*
## [11] fortify.lm* fortify.map*
## [13] fortify.NULL* fortify.Polygon*
## [15] fortify.Polygons* fortify.sfc*
## [17] fortify.sfg* fortify.SpatialLinesDataFrame*
## [19] fortify.SpatialPolygons* fortify.SpatialPolygonsDataFrame*
## [21] fortify.summary.glht* fortify.tbl*
## [23] fortify.tbl_df*
## see '?methods' for accessing help and source code
m = as(data_atlantic, "data.frame")
ggplot(m, aes(x, y)) + geom_point() + coord_equal()
#Diagramamos la caja de bloques de la variable POV
boxplot(data_atlantic$POV)# Caja de bigotes de la variable POV
Ejercicio implementar el codigo visto
implementamos lo visto la clase anterior
library(GWmodel) #Instalar GWmodel (modelos geograficamentes ponderados)
## Loading required package: maptools
## Checking rgeos availability: TRUE
## Loading required package: robustbase
## Loading required package: Rcpp
## Loading required package: spatialreg
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: Matrix
## Registered S3 methods overwritten by 'spatialreg':
## method from
## residuals.stsls spdep
## deviance.stsls spdep
## coef.stsls spdep
## print.stsls spdep
## summary.stsls spdep
## print.summary.stsls spdep
## residuals.gmsar spdep
## deviance.gmsar spdep
## coef.gmsar spdep
## fitted.gmsar spdep
## print.gmsar spdep
## summary.gmsar spdep
## print.summary.gmsar spdep
## print.lagmess spdep
## summary.lagmess spdep
## print.summary.lagmess spdep
## residuals.lagmess spdep
## deviance.lagmess spdep
## coef.lagmess spdep
## fitted.lagmess spdep
## logLik.lagmess spdep
## fitted.SFResult spdep
## print.SFResult spdep
## fitted.ME_res spdep
## print.ME_res spdep
## print.lagImpact spdep
## plot.lagImpact spdep
## summary.lagImpact spdep
## HPDinterval.lagImpact spdep
## print.summary.lagImpact spdep
## print.sarlm spdep
## summary.sarlm spdep
## residuals.sarlm spdep
## deviance.sarlm spdep
## coef.sarlm spdep
## vcov.sarlm spdep
## fitted.sarlm spdep
## logLik.sarlm spdep
## anova.sarlm spdep
## predict.sarlm spdep
## print.summary.sarlm spdep
## print.sarlm.pred spdep
## as.data.frame.sarlm.pred spdep
## residuals.spautolm spdep
## deviance.spautolm spdep
## coef.spautolm spdep
## fitted.spautolm spdep
## print.spautolm spdep
## summary.spautolm spdep
## logLik.spautolm spdep
## print.summary.spautolm spdep
## print.WXImpact spdep
## summary.WXImpact spdep
## print.summary.WXImpact spdep
## predict.SLX spdep
## Welcome to GWmodel version 2.2-3.
## The new version of GWmodel 2.2-4 now is ready
##
## Attaching package: 'GWmodel'
## The following objects are masked from 'package:stats':
##
## BIC, fitted
data(data_atlantic) #cargamos la base de datos
## Warning in data(data_atlantic): data set 'data_atlantic' not found
class(data_atlantic) #clase de la base de datos
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
Vemos que esta reconocido como un objeto espacial
Se implementa la base de datos “data_atlantic”, la observamos y vemos sus dimensiones:
head(data_atlantic) #muestra las primeras 5 observaciones
Data_atlantic es una base de datos que contiene 666 observaciones y 9 variables:
dim(data_atlantic) #muestra la dimensiòn
## [1] 666 7
?data_atlantic
## No documentation for 'data_atlantic' in specified packages and libraries:
## you could try '??data_atlantic'
De las cuales reconocemos las coordenadas (x, y) con la siguiente información:
Rate: Tasa de cáncer de pulmón. POV: Pobreza PM25: Indicador de la contaminación urbana. NO2: El dióxido de nitrógeno como consecuencia del tráfico rodado. SO2: El dióxido de azufre como gas contaminante del aire.
Realizar el codigo enviado y comentar su significado:
Primero nos indica que tenemos un archio shp con la función readShapeSpatial muestra que es un archivo de datos espaciales
leer el archivo COUNTY_ATLANTIC.shp, conociendo este formato (shp), se deduce que es un archivo de datos espaciales..
library(GWmodel)
library(rgdal)
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/asus/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/asus/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/asus/Documents/R/win-library/4.0/rgdal/proj
Condado = readShapeSpatial(file.choose())
## Warning: readShapeSpatial is deprecated; use rgdal::readOGR or sf::st_read
## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read
plot(Condado)
Estado = readOGR(file.choose()) #Cargar Estados.shp
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\asus\Documents\doctorado\DOCTORADO\Clase 3 - Docente Fabio\ejercicios\Data_GWR\Data_GWR\STATE_ATLANTIC.shp", layer: "STATE_ATLANTIC"
## with 11 features
## It has 10 fields
## Integer64 fields read as strings: ID ID2
plot(Estado)
Datos = read.csv("C:/Users/asus/Documents/data_atlantic_1998_2012.csv", header=T)
SPDF = merge(Condado, Datos, by = "FIPS") # unir base de datos
names (SPDF)
## [1] "FIPS" "ID" "x.x" "y.x" "REGION_ID"
## [6] "DIVISION_I" "STATE_ID" "COUNTY_ID" "REGION" "DIVISION"
## [11] "STATE" "COUNTY" "x.y" "y.y" "Rate"
## [16] "POV" "SMOK" "PM25" "NO2" "SO2"
Polys= Define caracteristicas para generar una función que grafique cada estadístico para cada variable.
Atributo que dibuja las delimitaciones de los estado
col.palette=Crea una paleta de colores
polys = list("sp.lines", as(Estado, "SpatialLines"), col = "transparent", lwd=0.8, lty=1)
col.palette = colorRampPalette(c("blue", "sky blue", "green", "yellow", "red"), space = "rgb", interpolate = "linear")
Ahora, utilizamos la función gwss
la función gwss Calcula medidas descriptivas del conjunto de datos geográficamente ponderadas incluye variables: Rate, PM2.5, el kernel utilizado es bisquare y ancho de banda =48
gwss.rate.pm25 = gwss(SPDF, vars = c("Rate","PM25", "POV"), kernel = "bisquare", adaptive = TRUE, bw = 48)
gwss.rate.pm25
## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
##
## ***********************Calibration information*************************
##
## Local summary statistics calculated for variables:
## Rate PM25 POV
## Number of summary points: 666
## Kernel function: bisquare
## Summary points: the same locations as observations are used.
## Adaptive bandwidth: 48 (number of nearest neighbours)
## Distance metric: Euclidean distance metric is used.
##
## ************************Local Summary Statistics:**********************
## Summary information for Local means:
## Min. 1st Qu. Median 3rd Qu. Max.
## Rate_LM 60.719 66.421 70.614 72.207 76.842
## PM25_LM 10.165 11.186 11.489 11.864 12.427
## POV_LM 10.779 12.453 15.266 16.509 20.584
## Summary information for local standard deviation :
## Min. 1st Qu. Median 3rd Qu. Max.
## Rate_LSD 9.05700 10.64575 11.07438 11.73020 15.2284
## PM25_LSD 0.54993 0.94388 1.12628 1.23341 1.6819
## POV_LSD 3.73408 4.73953 5.12185 5.39749 5.8609
## Summary information for local variance :
## Min. 1st Qu. Median 3rd Qu. Max.
## Rate_LVar 82.02917 113.33210 122.64193 137.59767 231.9056
## PM25_LVar 0.30243 0.89091 1.26851 1.52131 2.8288
## POV_LVar 13.94337 22.46319 26.23332 29.13289 34.3497
## Summary information for Local skewness:
## Min. 1st Qu. Median 3rd Qu. Max.
## Rate_LSKe -0.071523 0.216442 0.386823 0.598824 0.8325
## PM25_LSKe -1.456352 -0.890750 -0.410189 -0.102092 0.4324
## POV_LSKe -0.434795 0.205608 0.400079 0.584543 0.9860
## Summary information for localized coefficient of variation:
## Min. 1st Qu. Median 3rd Qu. Max.
## Rate_LCV 0.135837 0.149828 0.159503 0.173521 0.2031
## PM25_LCV 0.044980 0.080001 0.098187 0.110133 0.1637
## POV_LCV 0.255187 0.304386 0.331298 0.379579 0.4472
## Summary information for localized Covariance and Correlation between these variables:
## Min. 1st Qu. Median 3rd Qu. Max.
## Cov_Rate.PM25 -3.1528683 -1.1347505 -0.6292969 -0.0801118 0.8349
## Cov_Rate.POV 8.6548252 18.7118199 22.1321511 26.0467903 47.7038
## Cov_PM25.POV -2.0831025 -1.2793668 -0.8609368 -0.4255488 0.2914
## Corr_Rate.PM25 -0.2213778 -0.1009946 -0.0569110 -0.0064755 0.0680
## Corr_Rate.POV 0.1752022 0.3354824 0.3999972 0.4623835 0.6018
## Corr_PM25.POV -0.3563237 -0.2314432 -0.1609050 -0.0810170 0.0495
## Spearman_rho_Rate.PM25 -0.2001301 -0.0743685 -0.0125142 0.0456275 0.1341
## Spearman_rho_Rate.POV 0.2501326 0.3686766 0.4001836 0.4476783 0.5721
## Spearman_rho_PM25.POV -0.3639950 -0.2338469 -0.1655752 -0.0897677 0.0403
##
## ************************************************************************
Ahora, visualizamos las diferentes medidas, por ejemplo:
Media local
Para la media geográficamente ponderada de la variable Lung rate sería:
Mean.Rate = spplot(gwss.rate.pm25$SDF, "Rate_LM", main = "Media Local de tasa cancer", sp.layout = list(polys), col = "transparent", col.regions = col.palette(100))
Mean.Rate
#Ahora Para la variable PM25 material particulado, tamaño 2.5 micrómetros la media local seria:
Mean.Rate = spplot(gwss.rate.pm25$SDF, "PM25_LM", main = "material particulado, tamaño 2.5 micrómetros ", sp.layout = list(polys), col = "transparent", col.regions = col.palette(100))
Mean.Rate #media local para PM2.5
Si deseamos analizar y comparar las medias podemos realizar la siguiente visualización que permite observar mejor
Media Local Ahora analizamos y comparamos la media local para las variables Rate, PM25 y POV en un solo grafico
datos.grafico <- function(vble, title) {spplot(gwss.rate.pm25$SDF, vble, col = "transparent", col.regions = col.palette(100), cuts = 7, main = title, sp.layout = polys)}
datos.grafico
## function(vble, title) {spplot(gwss.rate.pm25$SDF, vble, col = "transparent", col.regions = col.palette(100), cuts = 7, main = title, sp.layout = polys)}
library(gridExtra) #cargar libreria
x=datos.grafico("Rate_LM","Tasa de cancer ML")
y=datos.grafico("PM25_LM","PM2.5 ML")
z=datos.grafico("POV_LM","Pobreza ML")
grid.arrange(x,y,z, nrow=3, ncol=2) #ensamblar varios gráficos en uno, especificando por numero de filas y columnas a ordenar
Análisis:
De los anteriores graficos se puede evidenciar que las tres variables tienen comportamientos similares, en el sur del mapa se observa que las tres variables es donde presentan mayor concentración, aunque la variable PM25 tiende un poco más al suroeste. Por otro lado se observa que en el norte del mapa presentan la menor concentración.
#Desviacion
Ahora visualizamos y analizamos el comportamiento de la desviación:
#Desviacion para Rate
Desviación de la variable rate (tasa de cáncer de pulmón)
spplot(gwss.rate.pm25$SDF, "Rate_LSD", key.space = "right", col = "transparent", col.regions = col.palette(100), cuts = 7, main = "Desviación GP para RATE", sp.layout = polys)
De igual forma que la anterior medida, la podemos unificar en un solo grafico que permita observar el comportamiento de la desviación de las tres variables y escribir conclusiones a partir de ella:
Desviación estándar local Para las variables (tasa de cáncer de pulmón), PM2.5 (material particulado, tamaño 2.5 micrómetros), POV (pobreza)
x=datos.grafico("Rate_LSD","Tasa de cancer DeL")
y=datos.grafico("PM25_LSD","PM2.5 DeL")
z=datos.grafico("POV_LSD","Pobreza DeL")
grid.arrange(x,y,z, nrow=3, ncol=1)
Análisis: Las anteriores graficas permite observar que el primer grafico infica que en el centro y al oeste se encuentra mayor dispersion en la tasa de cancer (Rate), es decir, la variable tasa de cáncer en esa zonas varió bastante. Por otro lado, donde se observa azul, la tasa de cáncer fue similar por lo tanto hubo menor dispersión.
el grafico 2 de la variable PM2.5 (material particulado, tamaño 2.5 micrómetros) permite observar que al norte se encuentra mayor dispersión, evidenciando que la variable PM25 varia bastante, pero en el sur el comportamiento es más similar.
Finalmente el grafico de la pobreza, muestra que en la zona del sur hubo mayor dispersión, indicando que la probabilidad que existe mayor desigualdad económica que en el norte del mapa, donde se visualiza que se presenta el menor valor de dispersión y por lo tanto, mayor similaridad.
Coeficiente de variación local
Por otro lado tambien podemos visualizar el coeficiende de variación local entre la tasa de cáncer y el PM25
spplot(gwss.rate.pm25$SDF, "Cov_Rate.PM25", key.space = "right", col = "transparent", col.regions = col.palette(100), cuts = 7, main = "Cov_Rate.PM25", sp.layout = polys)
Por otro lado tambien podemos visualizar el coeficiende de variación local entre la tasa de cáncer y la pobreza
spplot(gwss.rate.pm25$SDF, "Cov_Rate.POV", key.space = "right", col = "transparent", col.regions = col.palette(100), cuts = 7, main = "Cov_Rate.PM25", sp.layout = polys)
Correlación geográficamente ponderada
Ahora visualizamos la correlación geograficamente ponderada entre la tasa de cáncer y el PM2.5
spplot(gwss.rate.pm25$SDF, "Corr_Rate.PM25", key.space = "right", col = "transparent", col.regions = col.palette(100), cuts = 7, main = "Cov_Rate.PM25", sp.layout = polys)
Ahora visualizamos la correlación geograficamente ponderada entre la tasa de cáncer y la pobreza
spplot(gwss.rate.pm25$SDF, "Corr_Rate.POV", key.space = "right", col = "transparent", col.regions = col.palette(100), cuts = 7, main = "Cov_Rate.PM25", sp.layout = polys)