Tutorial: Spatial econometrics with R 3

Load package

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

Load columbus data

#------------

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 ...

Wcontig.nb =poly2nb(columbus, queen=FALSE)
Wdel.nb = tri2nb(coord) #Matrix based on Delaunay triangulation
Wdist.nb = dnearneigh(coord,0,0.7,longlat=FALSE) #Distance threshold matrix for columbus 

Wdist.listw = nb2listw(Wdist.nb) #6. comparing the weights  

neighbourmap(columbus,"X", Wcontig.nb) 
neighbourmap(columbus,"X", Wdel.nb) 
barnbmap(columbus,Wdist.nb) 
#histnbmap(columbus, Wnear.nb)


moran.plot(columbus@data$CRIME, Wdist.listw)


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

For use in interactive below function c(“orange”,“blue”,“pink”,“green”)

#
moranplotmap(columbus,"CRIME", Wcontig.listw) 
#LISA"https://fr.overleaf.com/project/5e15e2ec2ced2c0001aa9e88
columbus@data$localI = localmoran(columbus@data$CRIME,Wcontig.listw, alternative="two.sided")

columbus@data$pval = columbus@data$localI[,5]
columbus@data$pre_select = (columbus@data$pval<0.05)
columbus@data$localI[columbus@data$pre_select,1]

##          6         10         11         15         24         27         28 
## -1.5641932  1.4893917  1.3166926  1.2555476  1.2939326  0.8759771  2.0796486 
##         29         31         35         39 
##  1.2634974  1.2398550  1.2432425  1.2812882

moranplotmap(columbus,"CRIME",Wcontig.listw,identify=TRUE, 
             flower=TRUE,criteria=columbus@data$pre_select)
columbus@data$bub=abs(columbus@data$localI[,4])
columbus@data[(columbus@data$pval)<0.05,]$localI

##            Ii        E.Ii    Var.Ii      Z.Ii   Pr(z != 0)
## 6  -1.5641932 -0.02083333 0.3112215 -2.766511 5.665965e-03
## 10  1.4893917 -0.02083333 0.2283710  3.160248 1.576347e-03
## 11  1.3166926 -0.02083333 0.1786607  3.164374 1.554172e-03
## 15  1.2555476 -0.02083333 0.1218489  3.656533 2.556495e-04
## 24  1.2939326 -0.02083333 0.1786607  3.110527 1.867538e-03
## 27  0.8759771 -0.02083333 0.1218489  2.569152 1.019476e-02
## 28  2.0796486 -0.02083333 0.2283710  4.395401 1.105684e-05
## 29  1.2634974 -0.02083333 0.2283710  2.687549 7.197848e-03
## 31  1.2398550 -0.02083333 0.2283710  2.638076 8.337789e-03
## 35  1.2432425 -0.02083333 0.1786607  2.990602 2.784278e-03
## 39  1.2812882 -0.02083333 0.2283710  2.724778 6.434480e-03

#---------------

###################################### February 22 ####################
#discrete variables 
med <- median(columbus$CRIME)
HICRIME=as.factor(columbus$CRIME>med)

joincount.multi(HICRIME,Wcontig.listw) 

##             Joincount Expected Variance z-value
## FALSE:FALSE   9.51944  6.25000  0.54227  4.4398
## TRUE:TRUE     9.17619  5.75000  0.52929  4.7094
## TRUE:FALSE    5.80437 12.50000  1.53304 -5.4077
## Jtot          5.80437 12.50000  1.53304 -5.4077

HICRIME <- c("High", "Low")[(columbus$CRIME < med) + 1]

HICRIME=as.factor(HICRIME) 
s.HICRIME<-sample(HICRIME)

#Moran tests #Free sampling model - Gaussian test #one sided test 
moran.test(columbus@data$CRIME, Wcontig.listw, randomisation=FALSE)

## 
##  Moran I test under normality
## 
## data:  columbus@data$CRIME  
## weights: Wcontig.listw    
## 
## Moran I statistic standard deviate = 5.4978, p-value = 1.923e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.52367021       -0.02083333        0.00980890

#two sided test 
moran.test(columbus$CRIME, Wcontig.listw, randomisation=TRUE, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  columbus$CRIME  
## weights: Wcontig.listw    
## 
## Moran I statistic standard deviate = 5.4579, p-value = 4.819e-08
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.523670213      -0.020833333       0.009952987

#Non-free sampling model - gaussian test #one sided test 
moran.test(columbus$CRIME, Wcontig.listw, randomisation=TRUE) #two sided test 

## 
##  Moran I test under randomisation
## 
## data:  columbus$CRIME  
## weights: Wcontig.listw    
## 
## Moran I statistic standard deviate = 5.4579, p-value = 2.409e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.523670213      -0.020833333       0.009952987

moran.test(columbus$CRIME, Wcontig.listw, randomisation=TRUE, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  columbus$CRIME  
## weights: Wcontig.listw    
## 
## Moran I statistic standard deviate = 5.4579, p-value = 4.819e-08
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.523670213      -0.020833333       0.009952987

#Non-free sampling model - Permutation test #one sided test 
moran.mc(columbus$CRIME, Wcontig.listw, 9999) 

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  columbus$CRIME 
## weights: Wcontig.listw  
## number of simulations + 1: 10000 
## 
## statistic = 0.52367, observed rank = 10000, p-value = 1e-04
## alternative hypothesis: greater

#two sided test

#moran.mc(columbus$CRIME, Wcontig.listw, 9999, 
 #        alternative = "two.sided") # Huong's note: Does not work

#small exercise #permute the values of CRIME and re-run the moran.test 
S=sample(columbus$CRIME, length(columbus$CRIME)) 
moran.test(S, Wcontig.listw, randomisation=TRUE, alternative="two.sided") 

## 
##  Moran I test under randomisation
## 
## data:  S  
## weights: Wcontig.listw    
## 
## Moran I statistic standard deviate = -0.54151, p-value = 0.5882
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.074856867      -0.020833333       0.009952987

#alternative with Geary’s test 
geary.test(columbus@data$CRIME, Wcontig.listw, randomisation=FALSE)

## 
##  Geary C test under normality
## 
## data:  columbus@data$CRIME 
## weights: Wcontig.listw 
## 
## Geary C statistic standard deviate = 4.5377, p-value = 2.844e-06
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.51544081        1.00000000        0.01140311

lmmod=lm(HOVAL~CRIME+INC,data=columbus)
lm.morantest(lmmod,Wcontig.listw, alternative="two.sided")

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = HOVAL ~ CRIME + INC, data = columbus)
## weights: Wcontig.listw
## 
## Moran I statistic standard deviate = 2.5033, p-value = 0.0123
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.207051059     -0.035519797      0.009389584

#end