Suggested citation:

Mendez C. (2020). Spatial regression analysis in R. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/tutorial-spatial-regression Acknowledgment:

Material adapted from multiple sources, in particular BurkeyAcademy’s GIS & Spatial Econometrics Project

# 1 Libraries

``````knitr::opts_chunk\$set(echo = TRUE)

library(tidyverse)  # Modern data science workflow
library(spdep)
library(spatialreg)
library(rgdal)
library(rgeos)

# Change the presentation of decimal numbers to 4 and avoid scientific notation
options(prompt="R> ", digits=4, scipen=7)``````

# 2 Tutorial objectives

• Import shapefiles into R

• Import neighbor relationship from `.gal` files

• Create neighbor relationships in R from shape files

• Create neighbor relationships in R from shape latitude and longitude

• Understand the difference between Great Circle and Euclidean distances

• Export neighbor relationships as weight matrices to plain text files

• Test for spatial dependence via the Moran’s I test

• Evaluate the four simplest models of spatial regression

# 5 Import spatial data

Let us use the `readOGR` function from the `rgdal` library to import the `.shp` file

``NCVACO <- readOGR(dsn = ".", layer = "NCVACO")``
``````## OGR data source with driver: ESRI Shapefile
## Source: "/Users/carlos/Github/QuaRCS-lab/tutorial-spatial-regression", layer: "NCVACO"
## with 234 features
## It has 49 fields
## Integer64 fields read as strings:  FIPS qtystores PCI COUNTMXBV DC GA KY MD SC TN WV VA COUNTBKGR TOTALPOP POP18OV LABFORCE HHOLDS POP25OV POP16OV``````

Note that the file is imported as a `SpatialPolygonsDataFrame` object

# 6 Import neighbor relationship: `.gal` file

Let us use the `read.gal` function from the `rgdal` library to import the `.gal` weights matrix created in GeoDa

``queen.nb <- read.gal("queen.gal", region.id=NCVACO\$FIPS)``

## 6.1 Summarize neighbor relationships

``summary(queen.nb)``
``````## Neighbour list object:
## Number of regions: 234
## Number of nonzero links: 1132
## Percentage nonzero weights: 2.067
## Average number of links: 4.838
##
##  1  2  3  4  5  6  7  8  9 11
## 16 26 22 32 36 49 37  9  6  1
## 16 least connected regions:
## 51515 51530 51595 51660 51690 51720 51820 51131 51520 51540 51580 51600 51678 51790 51840 51001 with 1 link
## 1 most connected region:
• Is the is the neighbor relationship symmetric?
``is.symmetric.nb(queen.nb)``
``##  TRUE``

# 7 Create is the neighbor relationship in R

## 7.1 From a shapefile

• For queen contiguity
``queen.R.nb <- poly2nb(NCVACO, row.names = NCVACO\$FIPS)``
• Alternatively, you can create a Rook contiguity relationship as
``rook.R.nb <- poly2nb(NCVACO, row.names = NCVACO\$FIPS, queen=FALSE)``

### 7.1.1 Summarize neighbor relationships

``summary(queen.R.nb)``
``````## Neighbour list object:
## Number of regions: 234
## Number of nonzero links: 1132
## Percentage nonzero weights: 2.067
## Average number of links: 4.838
##
##  1  2  3  4  5  6  7  8  9 11
## 16 26 22 32 36 49 37  9  6  1
## 16 least connected regions:
## 51515 51530 51595 51660 51690 51720 51820 51131 51520 51540 51580 51600 51678 51790 51840 51001 with 1 link
## 1 most connected region:
• Are the relationships symmetric?
``is.symmetric.nb(queen.R.nb)``
``##  TRUE``

## 7.2 From latitude and longitude

• Import table
``nc.cent <- read.csv(file="CenPop2010_Mean_CO37.txt")``
• Identify coordinates
``nc.coords <- cbind(nc.cent\$LONGITUDE, nc.cent\$LATITUDE)``

### 7.2.1 Identify 5 nearest neighbors

#### 7.2.1.1 The right way

Recognize that latitude and longitude are handled using great circle distances

``nc.5nn <- knearneigh(nc.coords, k=5, longlat = TRUE)``

#### 7.2.1.2 The wrong way

Fail to recognize that latitude and longitude are handled using great circle distances. Latitude and longitude should not be used to compute Euclidean distances

``nc.5nn.wrong <- knearneigh(nc.coords, k=5, longlat = FALSE)``

### 7.2.2 Create 5 nearest neighbors relationship

``nc.5nn.nb <- knn2nb(nc.5nn)``

Plot the right neighbor relationship

``plot(nc.5nn.nb, nc.coords)``