Ejercio tarea práctico con datos atlanta`

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)

1. Cargamos la base de datos “data_atlantic”, la observamos y vemos sus dimensiones

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.

2. Tarea ejercicio 24-03-2021 -Realizar y comentar el script enviado por el profesor en las notas de clase:

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)