#Evaluamos la existencia de autocorrelaciĆ³n espacial
1. Carga de libreria y datos
library(sp) #Libraria para el tratamiento de datos espaciales
library(readr)
library(spdep)
## 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: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(maptools)
## Checking rgeos availability: TRUE
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
data_atlantic<- read_csv("data_atlantic_1998_2012.csv")
##
## -- Column specification --------------------------------------------------------
## cols(
## FIPS = col_double(),
## x = col_double(),
## y = col_double(),
## Rate = col_double(),
## POV = col_double(),
## SMOK = col_double(),
## PM25 = col_double(),
## NO2 = col_double(),
## SO2 = col_double()
## )
#Cargamos el archivo County Atlantic
Condado = readOGR(file.choose())
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\asus\Documents\DOCTORADO\MATERIAS\PROFUNDIZACION_I\R\AUTOCORRELACION\COUNTY_ATLANTIC.shp", layer: "COUNTY_ATLANTIC"
## with 666 features
## It has 12 fields
## Integer64 fields read as strings: ID
#Ploteamos el archivo condado
plot(Condado)
class(Condado) #Observamos el Tipo de datos de condado ("SpatialPolygonsDataFrame")
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
head(Condado@data)
Datos = read.csv(file.choose(), header=T) # Cargamos la base de datos
head(Datos)
2. Union de datos comando āmergeā
SPDF = merge(Condado, Datos, by = "FIPS") # unir base de datos con el SpatialPolygonsDataFrame "Condado"
head(SPDF@data)
class(SPDF)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
dim(SPDF)
## [1] 666 20
names (SPDF) # Observamos los nombres de la base de datos creada
## [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"
spplot(SPDF,"PM25")#Grafica de distribuciĆ³n espacial de la variable PM25
3. ConstrucciĆ³n de una lista de vecindades
X_nbq = poly2nb(SPDF)
coords = coordinates(SPDF)
class(coords)
## [1] "matrix" "array"
head(coords)
## [,1] [,2]
## 0 1056524 1376613
## 1 1653442 2267301
## 2 1633708 2096683
## 3 1584049 1901443
## 4 1735811 2409536
## 5 1003647 1193902
dim(coords)
## [1] 666 2
plot(SPDF, main = "Mapa de contactos usando criterio de contigĆ¼idad tipo reina")
plot(X_nbq, coords, add=T)
######### K- vecinos
IDs = row.names(as(SPDF, "data.frame"))
class(IDs)
## [1] "character"
head(IDs)
## [1] "59" "428" "408" "594" "237" "78"
length(IDs)
## [1] 666
X_kn4 = knn2nb(knearneigh(coords, k = 4), row.names=IDs)
?knearneigh
## starting httpd help server ... done
plot(SPDF, main = "Mapa de contacto usando criterio de contigĆ¼idad 4 vecinos mĆ”s cercanos")
plot(X_kn4, coords, add=T)
4. ConstrucciĆ³n de la matriz w para el criterio QUEEN(w)
#ConstrucciĆ³n de la matriz w
# Matriz de pesos normalizados en fila usando QUEEN(W) (Tipo reina):
X_nbq_w = nb2listw(X_nbq)
#?nb2listw
class(X_nbq_w)
## [1] "listw" "nb"
#head(X_nbq_w)
summary(unlist(X_nbq_w$weights))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.09091 0.14286 0.16667 0.18439 0.20000 1.00000
X_nbq_w$style
## [1] "W"
X_nbq_w$neighbours
## Neighbour list object:
## Number of regions: 666
## Number of nonzero links: 3612
## Percentage nonzero weights: 0.8143278
## Average number of links: 5.423423
#X_nbq_w$weights
5. ConstrucciĆ³n de la matriz w para k_vecinos = 4
# Matriz de pesos normalizados en fila usando k = 4(W):
X_kn4_w = nb2listw(X_kn4)
class(X_kn4_w)
## [1] "listw" "nb"
#head(X_kn4_w)
summary(unlist(X_kn4_w$weights))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.25 0.25 0.25 0.25 0.25 0.25
6. Indice de Moran para el criterio tipo Queen
# Indice de Moran (TIPO QUEEN)
Ind_Moran_reina = moran.test(SPDF$PM25, listw = X_nbq_w, alternative = "two.side")
#?moran.test
Ind_Moran_reina
##
## Moran I test under randomisation
##
## data: SPDF$PM25
## weights: X_nbq_w
##
## Moran I statistic standard deviate = 37.413, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.9185498483 -0.0015037594 0.0006047565
scatterplot_moran = moran.plot(SPDF$PM25, listw = X_nbq_w, main = "Scatterplot de Moran para EL PM25 DE ATLANTIC (QUEEN)")
#Prueba de Monte Carlo para Queen
moran.mc(SPDF$PM25, listw = X_nbq_w, nsim=999)
##
## Monte-Carlo simulation of Moran I
##
## data: SPDF$PM25
## weights: X_nbq_w
## number of simulations + 1: 1000
##
## statistic = 0.91855, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
#?moran.mc
7.Indice de Moran para k_vecinos = 4
# Indice de Moran (k_vecions=4)
Ind_Moran_kn4 = moran.test(SPDF$PM25, listw = X_kn4_w, alternative = "two.side")
#?moran.test
Ind_Moran_kn4
##
## Moran I test under randomisation
##
## data: SPDF$PM25
## weights: X_kn4_w
##
## Moran I statistic standard deviate = 35.484, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.9261589438 -0.0015037594 0.0006834622
scatterplot_moran = moran.plot(SPDF$PM25, listw = X_kn4_w, main = "Scatterplot de Moran para EL PM25 DE ATLANTIC(K=4)")
#Prueba de Monte Carlo para k=4
moran.mc(SPDF$PM25, listw = X_kn4_w, nsim=999)
##
## Monte-Carlo simulation of Moran I
##
## data: SPDF$PM25
## weights: X_kn4_w
## number of simulations + 1: 1000
##
## statistic = 0.92616, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater