Tutorial: Spatial econometrics with R 2

Load package

require(spdep) #require("spdep)
require(spatialreg)
library(raster) 
library(GeoXp)
library(geodata)
library(cartography)
library(sp)
library(rgdal)
library(rgdal)
require(spData) 
library(GeoXp)

The class listw : small example

t=c(1,2,3,4) 
u=c(3,2,5,1) 
plot(t,u) 
co=cbind(t,u) 
W.knn=knearneigh(co,k=2,longlat=TRUE) 
plot(knn2nb(W.knn), co, add=TRUE) 

W.nb=knn2nb(W.knn) 
W.listw1=nb2listw(W.nb,style="B") 
str(W.listw1) 

## List of 3
##  $ style     : chr "B"
##  $ neighbours:List of 4
##   ..$ : int [1:2] 2 3
##   ..$ : int [1:2] 1 4
##   ..$ : int [1:2] 1 2
##   ..$ : int [1:2] 1 2
##   ..- attr(*, "region.id")= chr [1:4] "1" "2" "3" "4"
##   ..- attr(*, "call")= language knearneigh(x = co, k = 2, longlat = TRUE)
##   ..- attr(*, "sym")= logi FALSE
##   ..- attr(*, "type")= chr "knn"
##   ..- attr(*, "knn-k")= num 2
##   ..- attr(*, "class")= chr "nb"
##  $ weights   :List of 4
##   ..$ : num [1:2] 1 1
##   ..$ : num [1:2] 1 1
##   ..$ : num [1:2] 1 1
##   ..$ : num [1:2] 1 1
##   ..- attr(*, "mode")= chr "binary"
##   ..- attr(*, "B")= logi TRUE
##  - attr(*, "class")= chr [1:2] "listw" "nb"
##  - attr(*, "region.id")= chr [1:4] "1" "2" "3" "4"
##  - attr(*, "call")= language nb2listw(neighbours = W.nb, style = "B")

W.listw1$neighbours[] 

## [[1]]
## [1] 2 3
## 
## [[2]]
## [1] 1 4
## 
## [[3]]
## [1] 1 2
## 
## [[4]]
## [1] 1 2
## 
## attr(,"region.id")
## [1] "1" "2" "3" "4"
## attr(,"call")
## knearneigh(x = co, k = 2, longlat = TRUE)
## attr(,"sym")
## [1] FALSE
## attr(,"type")
## [1] "knn"
## attr(,"knn-k")
## [1] 2

W.listw1$weights[] 

## [[1]]
## [1] 1 1
## 
## [[2]]
## [1] 1 1
## 
## [[3]]
## [1] 1 1
## 
## [[4]]
## [1] 1 1
## 
## attr(,"mode")
## [1] "binary"
## attr(,"B")
## [1] TRUE

W.listw2=nb2listw(W.nb,style="W") 
W.listw2$neighbours[] 

## [[1]]
## [1] 2 3
## 
## [[2]]
## [1] 1 4
## 
## [[3]]
## [1] 1 2
## 
## [[4]]
## [1] 1 2
## 
## attr(,"region.id")
## [1] "1" "2" "3" "4"
## attr(,"call")
## knearneigh(x = co, k = 2, longlat = TRUE)
## attr(,"sym")
## [1] FALSE
## attr(,"type")
## [1] "knn"
## attr(,"knn-k")
## [1] 2

W.listw2$weights[]

## [[1]]
## [1] 0.5 0.5
## 
## [[2]]
## [1] 0.5 0.5
## 
## [[3]]
## [1] 0.5 0.5
## 
## [[4]]
## [1] 0.5 0.5
## 
## attr(,"mode")
## [1] "binary"
## attr(,"W")
## [1] TRUE
## attr(,"comp")
## attr(,"comp")$d
## [1] 2 2 2 2

Auckland example

#loading Auckland 
require("spData") 
library(rgdal) 
library(spdep) 
auckland <- readOGR(system.file("shapes/auckland.shp",package="spData")[1]) 

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\TRINHHUONG\Documents\R\win-library\4.0\spData\shapes\auckland.shp", layer: "auckland"
## with 167 features
## It has 4 fields

auckland.nb <- poly2nb(auckland) 
summary(auckland.nb) 

## Neighbour list object:
## Number of regions: 167 
## Number of nonzero links: 772 
## Percentage nonzero weights: 2.768116 
## Average number of links: 4.622754 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 
##  1  9 39 42 26 27 13  7  1  2 
## 1 least connected region:
## 57 with 1 link
## 2 most connected regions:
## 141 156 with 10 links

plot(auckland) 
plot(auckland.nb,coord=cbind(auckland$Easting,auckland$Northing), add=TRUE)

Loading columbus

Loading columbus with sf library(rgdal)

####Loading columbus with sf 
library(rgdal) 
columbus <- readOGR(system.file("shapes/columbus.shp", package = "spData")[1]) 

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\TRINHHUONG\Documents\R\win-library\4.0\spData\shapes\columbus.shp", layer: "columbus"
## with 49 features
## It has 20 fields
## Integer64 fields read as strings:  COLUMBUS_ COLUMBUS_I POLYID

library(sf) 
columbus_sf <- read_sf(system.file("shapes/columbus.shp", package = "spData")[1]) 
plot(columbus_sf[, c("HOVAL", "CRIME", "NSA", "CP")]) 

class(columbus) 

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

coord <- coordinates(columbus)[, 1:2] 
dim(columbus) 

## [1] 49 20

head(columbus@data)

##       AREA PERIMETER COLUMBUS_ COLUMBUS_I POLYID NEIG  HOVAL    INC    CRIME
## 0 0.309441  2.440629         2          5      1    5 80.467 19.531 15.72598
## 1 0.259329  2.236939         3          1      2    1 44.567 21.232 18.80175
## 2 0.192468  2.187547         4          6      3    6 26.350 15.956 30.62678
## 3 0.083841  1.427635         5          2      4    2 33.200  4.477 32.38776
## 4 0.488888  2.997133         6          7      5    7 23.225 11.252 50.73151
## 5 0.283079  2.335634         7          8      6    8 28.750 16.029 26.06666
##       OPEN    PLUMB DISCBD     X     Y NSA NSB EW CP THOUS NEIGNO
## 0 2.850747 0.217155   5.03 38.80 44.07   1   1  1  0  1000   1005
## 1 5.296720 0.320581   4.27 35.62 42.38   1   1  0  0  1000   1001
## 2 4.534649 0.374404   3.89 39.82 41.18   1   1  1  0  1000   1006
## 3 0.394427 1.186944   3.70 36.50 40.52   1   1  0  0  1000   1002
## 4 0.405664 0.624596   2.83 40.01 38.00   1   1  1  0  1000   1007
## 5 0.563075 0.254130   3.78 43.75 39.28   1   1  1  0  1000   1008

str(columbus@data)

## 'data.frame':    49 obs. of  20 variables:
##  $ AREA      : num  0.3094 0.2593 0.1925 0.0838 0.4889 ...
##  $ PERIMETER : num  2.44 2.24 2.19 1.43 3 ...
##  $ COLUMBUS_ : chr  "2" "3" "4" "5" ...
##  $ COLUMBUS_I: chr  "5" "1" "6" "2" ...
##  $ POLYID    : chr  "1" "2" "3" "4" ...
##  $ NEIG      : int  5 1 6 2 7 8 4 3 18 10 ...
##  $ HOVAL     : num  80.5 44.6 26.4 33.2 23.2 ...
##  $ INC       : num  19.53 21.23 15.96 4.48 11.25 ...
##  $ CRIME     : num  15.7 18.8 30.6 32.4 50.7 ...
##  $ OPEN      : num  2.851 5.297 4.535 0.394 0.406 ...
##  $ PLUMB     : num  0.217 0.321 0.374 1.187 0.625 ...
##  $ DISCBD    : num  5.03 4.27 3.89 3.7 2.83 3.78 2.74 2.89 3.17 4.33 ...
##  $ X         : num  38.8 35.6 39.8 36.5 40 ...
##  $ Y         : num  44.1 42.4 41.2 40.5 38 ...
##  $ NSA       : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ NSB       : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ EW        : num  1 0 1 0 1 1 0 0 1 1 ...
##  $ CP        : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ THOUS     : num  1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 ...
##  $ NEIGNO    : num  1005 1001 1006 1002 1007 ...

Loading columbus the old way

library("spData") 
library("rgdal") 
library("spdep")
columbus <- readOGR(system.file("shapes/columbus.shp",package="spData")[1]) 

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\TRINHHUONG\Documents\R\win-library\4.0\spData\shapes\columbus.shp", layer: "columbus"
## with 49 features
## It has 20 fields
## Integer64 fields read as strings:  COLUMBUS_ COLUMBUS_I POLYID

coord=coordinates(columbus)
if (requireNamespace("spdep", quietly = TRUE)) { 
  library(spdep) 
  col.gal.nb <- read.gal(system.file("weights/columbus.gal",package="spData")[1]) }

str(columbus@data)

## 'data.frame':    49 obs. of  20 variables:
##  $ AREA      : num  0.3094 0.2593 0.1925 0.0838 0.4889 ...
##  $ PERIMETER : num  2.44 2.24 2.19 1.43 3 ...
##  $ COLUMBUS_ : chr  "2" "3" "4" "5" ...
##  $ COLUMBUS_I: chr  "5" "1" "6" "2" ...
##  $ POLYID    : chr  "1" "2" "3" "4" ...
##  $ NEIG      : int  5 1 6 2 7 8 4 3 18 10 ...
##  $ HOVAL     : num  80.5 44.6 26.4 33.2 23.2 ...
##  $ INC       : num  19.53 21.23 15.96 4.48 11.25 ...
##  $ CRIME     : num  15.7 18.8 30.6 32.4 50.7 ...
##  $ OPEN      : num  2.851 5.297 4.535 0.394 0.406 ...
##  $ PLUMB     : num  0.217 0.321 0.374 1.187 0.625 ...
##  $ DISCBD    : num  5.03 4.27 3.89 3.7 2.83 3.78 2.74 2.89 3.17 4.33 ...
##  $ X         : num  38.8 35.6 39.8 36.5 40 ...
##  $ Y         : num  44.1 42.4 41.2 40.5 38 ...
##  $ NSA       : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ NSB       : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ EW        : num  1 0 1 0 1 1 0 0 1 1 ...
##  $ CP        : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ THOUS     : num  1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 ...
##  $ NEIGNO    : num  1005 1001 1006 1002 1007 ...

Some exploratory analysis

str(columbus@data)

## 'data.frame':    49 obs. of  20 variables:
##  $ AREA      : num  0.3094 0.2593 0.1925 0.0838 0.4889 ...
##  $ PERIMETER : num  2.44 2.24 2.19 1.43 3 ...
##  $ COLUMBUS_ : chr  "2" "3" "4" "5" ...
##  $ COLUMBUS_I: chr  "5" "1" "6" "2" ...
##  $ POLYID    : chr  "1" "2" "3" "4" ...
##  $ NEIG      : int  5 1 6 2 7 8 4 3 18 10 ...
##  $ HOVAL     : num  80.5 44.6 26.4 33.2 23.2 ...
##  $ INC       : num  19.53 21.23 15.96 4.48 11.25 ...
##  $ CRIME     : num  15.7 18.8 30.6 32.4 50.7 ...
##  $ OPEN      : num  2.851 5.297 4.535 0.394 0.406 ...
##  $ PLUMB     : num  0.217 0.321 0.374 1.187 0.625 ...
##  $ DISCBD    : num  5.03 4.27 3.89 3.7 2.83 3.78 2.74 2.89 3.17 4.33 ...
##  $ X         : num  38.8 35.6 39.8 36.5 40 ...
##  $ Y         : num  44.1 42.4 41.2 40.5 38 ...
##  $ NSA       : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ NSB       : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ EW        : num  1 0 1 0 1 1 0 0 1 1 ...
##  $ CP        : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ THOUS     : num  1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 ...
##  $ NEIGNO    : num  1005 1001 1006 1002 1007 ...

histomap(columbus,"HOVAL", col="green",pch=5)

histomap(columbus,"HOVAL",col="green",pch=5,identify="TRUE")

histomap(columbus,"HOVAL",type="percent",xlab="House Price")

histomap(columbus,"CRIME",criteria=(columbus@data$HOVAL>38))

densitymap(columbus,"HOVAL", col="green",pch=5,xlab="House Price")

Contiguity matrix for columbus

#Contiguity matrix for columbus

Wcontig.nb=poly2nb(columbus,queen=FALSE) 

plot(columbus, border= "grey", axes=TRUE, main="Contiguity") 

plot(Wcontig.nb,coords=coord,add=TRUE)

#Distance threshold matrix for columbus 
max(dist(coord)) 

## [1] 4.732267

Wdist.nb=dnearneigh(coord,0,0.7,longlat=FALSE) 
plot(columbus, border="grey",axes=TRUE,main="Distance threshold") 
plot(Wdist.nb,coords=coord,add=TRUE)

Wdist_inf.nb=dnearneigh(coord,0,4,longlat=FALSE) 
dlist_inf <- nbdists(Wdist_inf.nb, coord) 
Wdist_inv <- lapply(dlist_inf, function(x) 1/x) 
Wdist_inv.listw = nb2listw(Wdist_inf.nb, glist=Wdist_inv)

plot(columbus, border="grey",axes=TRUE,main="Inverse Distance ") 
plot(Wdist_inv.listw,coords=coord,add=TRUE)

Comparing the weights

#comparing the weights 
Wdist.listw=nb2listw(Wdist.nb) 
head(Wdist_inv.listw$weights[])

## [[1]]
##  [1] 0.065765426 0.067924681 0.052639837 0.036195766 0.030811909 0.030615766
##  [7] 0.038546729 0.024051143 0.019126173 0.028168153 0.027271712 0.027084012
## [13] 0.024648655 0.024096595 0.022497810 0.016267480 0.021226074 0.022139343
## [19] 0.015662457 0.015842761 0.018945065 0.014174148 0.018831119 0.019769608
## [25] 0.019680857 0.016415682 0.016155917 0.016545237 0.016175060 0.012597256
## [31] 0.012374896 0.014989890 0.012626912 0.014034614 0.011422947 0.014369959
## [37] 0.014390284 0.009972995 0.011891750 0.010974573 0.011183582 0.013294662
## [43] 0.012753172 0.012222519 0.009937377 0.009928746 0.012243483 0.011515208
## 
## [[2]]
##  [1] 0.055668142 0.046819204 0.098099102 0.033385235 0.022021694 0.043489680
##  [7] 0.051276456 0.019670868 0.014621046 0.030184070 0.030321564 0.032181852
## [13] 0.028520917 0.023138671 0.023129091 0.012839848 0.022303417 0.024416274
## [19] 0.012972003 0.016613328 0.016242760 0.011433679 0.019424242 0.018949696
## [25] 0.017916477 0.014193898 0.014686632 0.015594056 0.015699428 0.013015180
## [31] 0.010191907 0.013019390 0.012817829 0.012373917 0.011555567 0.013492120
## [37] 0.013178754 0.009885526 0.010032095 0.009110408 0.011121294 0.012144846
## [43] 0.011308333 0.011335449 0.009804993 0.008357551 0.011128580 0.010312930
## 
## [[3]]
##  [1] 0.052241493 0.042540531 0.053985732 0.057862759 0.037603499 0.029235870
##  [7] 0.045654092 0.028218846 0.019143568 0.033396885 0.030759081 0.028641159
## [13] 0.025309150 0.028070809 0.024333468 0.015800423 0.021975987 0.022452283
## [19] 0.015472532 0.014493073 0.020167845 0.013425758 0.018763023 0.021014466
## [25] 0.021236069 0.016634760 0.016246952 0.016584406 0.015866741 0.010624420
## [31] 0.011471094 0.014787412 0.010985307 0.013604341 0.009581409 0.013817506
## [37] 0.013962283 0.008203130 0.011029508 0.009959060 0.009581124 0.012641968
## [43] 0.012070135 0.011329357 0.008252929 0.008894490 0.011430678 0.010642591
## 
## [[4]]
##  [1] 0.036169462 0.079631236 0.048230142 0.039201590 0.020781536 0.049858185
##  [7] 0.085517234 0.019284090 0.013292683 0.035324919 0.035206663 0.037089496
## [13] 0.031110030 0.024553099 0.024440625 0.011621461 0.023085788 0.025507260
## [19] 0.011952156 0.015518282 0.015620987 0.010285957 0.019331604 0.019061606
## [25] 0.017752673 0.013385979 0.014016842 0.015040058 0.015091003 0.011322415
## [31] 0.009115011 0.012135982 0.011413525 0.011477790 0.009988289 0.012672415
## [37] 0.012356788 0.008401869 0.009012914 0.008076532 0.009757363 0.011248212
## [43] 0.010371357 0.010376366 0.008393563 0.007376940 0.010185356 0.009354665
## 
## [[5]]
##  [1] 0.022226309 0.024218928 0.046197740 0.035033678 0.032095115 0.027190329
##  [7] 0.046282476 0.033758998 0.017395209 0.052215012 0.041951647 0.033397325
## [13] 0.028923248 0.043498076 0.031979979 0.014624395 0.026467664 0.025769513
## [19] 0.015356368 0.014412644 0.022909659 0.012441834 0.021189483 0.026324794
## [25] 0.026412831 0.017794043 0.017908808 0.018560275 0.017336228 0.009567975
## [31] 0.010612216 0.015325751 0.010248732 0.013940418 0.008571600 0.014488588
## [37] 0.014662410 0.007177229 0.010400063 0.009082158 0.008767095 0.012895321
## [43] 0.012040981 0.011246107 0.007291887 0.008085650 0.011356018 0.010367196
## 
## [[6]]
##  [1] 0.025241620 0.021312787 0.040053355 0.024777001 0.042818159 0.019776956
##  [7] 0.025898075 0.059405773 0.041100914 0.027361189 0.024740945 0.021932354
## [13] 0.020534005 0.030414219 0.023443822 0.028223381 0.020687929 0.019664850
## [19] 0.024842448 0.013738474 0.031767265 0.021422477 0.018385855 0.023786453
## [25] 0.027613899 0.024345658 0.020667703 0.019370705 0.017331078 0.009934268
## [31] 0.016756191 0.020776580 0.010651830 0.018085801 0.009140661 0.015956922
## [37] 0.017063075 0.015313211 0.013811459 0.009444392 0.015298930 0.015464115
## [43] 0.013061971 0.011775947 0.013757143 0.013048154

#Matrix based on nearest neighbors 
Wnear.knn=knearneigh(coord, k=5)
#Matrix based on Delaunay triangulation
Wdel.nb=tri2nb(coord)

plot(columbus, border="grey",axes=TRUE, main="5 Nearest neighbors") 

plot(Wdel.nb,coord,add=TRUE)

gabriel.nb=graph2nb(gabrielneigh(coord), sym=TRUE) 
plot(columbus) 
plot(gabriel.nb,coords=coord, add=TRUE)

columbus@data$CRIME_lag=lag.listw(nb2listw(knn2nb(Wnear.knn)),style="W",columbus@data$CRIME)

columbus_sf$WCRIME=columbus@data$CRIME_lag 

plot(columbus_sf[,c("CRIME","WCRIME")]) 
columbus_sf$WCRIME=columbus@data$CRIME_lag

plot(columbus_sf[,c("CRIME","WCRIME")])

is.symmetric.nb(Wcontig.nb) 

## [1] TRUE

is.symmetric.nb(Wdist.nb) 

## [1] TRUE

#is.symmetric.nb(Wnear.nb) 
is.symmetric.nb(Wdel.nb)

## [1] TRUE