Analyzing Spatial Autocorrelation Using CRD Census Tract Data

I. Introduction

Two different variables from the Victoria Capital Regional District (CRD) 2016 Census Data will be analyzed to determine if either of them has any spatial autocorrelation with themselves, over space.

The following method is for the first variable Median Income (MdInc), and the entire will be duplicated for the second variable Renter, after changing the new variable name, and renaming the output files.

II. Set Working Directory

Set the working directory to point to the data

## Set working directory
dir <- "/Users/geographydepartment/Documents/Geog418/Assign3-Geog418/Data/Working/Data"
setwd(dir)
getwd()

III. Load Libraries

Load the Libraries that are needed to run the code. If these packages have not previously been installed, use the function install.packages(“”) for each package.

plyr: A tool to split, apply and combine data https://www.rdocumentation.org/packages/plyr/versions/1.8.4
dplyr: A tool to transform and manipulate data https://dplyr.tidyverse.org/, https://github.com/rstudio/cheatsheets/blob/master/data-transformation.pdf
spdep: Spatial Dependence to create weighted matrix from contiguous polygons https://cran.r-project.org/web/packages/spdep/index.html, https://cran.r-project.org/web/packages/spdep/spdep.pdf
GISTools: mapping and spatial data tool, especially for choloropleth maps and legends https://cran.r-project.org/web/packages/GISTools/GISTools.pdf
raster: read, write, analyze and model spatial gridded data https://cran.r-project.org/web/packages/raster/raster.pdf
maptools: tools to manipulate geographic data https://cran.r-project.org/web/packages/maptools/maptools.pdf
rgdal: Projection transformations for PROJ.4 library https://cran.r-project.org/web/packages/rgdal/rgdal.pdf
tmap: Thematic map tool to visualize spatial data https://cran.r-project.org/web/packages/tmap/tmap.pdf
BAMMtools: ??? Not sure WHY THIS IS NEEDED ??? https://cran.r-project.org/web/packages/BAMMtools/BAMMtools.pdf
shinyjs: Allows javascript to be used https://deanattali.com/shinyjs/overview

## Load the following Libraries which will be needed to run this census analysis
# if the packages are not already installed, first run function for each: 
# install.packages("plyr")
library("plyr")
library("dplyr")
library("spdep")
library("GISTools")
library("raster")
library("maptools")
library("rgdal")
library("tmap")
library("BAMMtools")
library("shinyjs")

IV. Clean Data

1. Read in the data from the census tract boundaries shapefile, using the rgdal package (don’t include the .shp file extension), and the 2016 census data .csv file (include the .csv file extension), assigning them each a variable name.
2. Merge the census tract boundaries with the census data, assigning this merged data a new variable name.
3. Remove any missing ‘na’ values from the merged data, reassigning these

## Clean data
# Read in shapefile from working directory "."; do not include the .shp extension
# uses rgdal package
tracts <- readOGR(dsn = ".", layer = "Vic_Census")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/geographydepartment/Documents/Geog418/Assign3-Geog418/Data/Working/Data", layer: "Vic_Census"
## with 78 features
## It has 2 fields

tracts

## class       : SpatialPolygonsDataFrame 
## features    : 78 
## extent      : 3903547, 3974156, 1910542, 1962451  (xmin, xmax, ymin, ymax)
## crs         : +proj=lcc +lat_1=49 +lat_2=77 +lat_0=63.390675 +lon_0=-91.86666666666666 +x_0=6200000 +y_0=3000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## variables   : 2
## names       :      CTUID,       GUID 
## min values  : 9350001.00,    9350001 
## max values  : 9350180.06, 9350180.06

# Read in csv file from working directory; include .csv extension
census.16 <- read.csv("CensusTractData.csv")

# Summary of the data attributes
summary(census.16)

## some summary stats have been ommitted ... 
##       GUID              CUID            Pop             PopCh        
##  Min.   :9350000   Min.   :  0.0   Min.   :    94   Min.   :-60.200  
##  1st Qu.:9350102   1st Qu.:101.5   1st Qu.:  3479   1st Qu.:  1.725  
##  Median :9350128   Median :128.0   Median :  4689   Median :  4.100  
##  Mean   :9350110   Mean   :109.6   Mean   :  9430   Mean   :  4.599  
##  3rd Qu.:9350154   3rd Qu.:154.0   3rd Qu.:  6156   3rd Qu.:  7.275  
##  Max.   :9350180   Max.   :180.1   Max.   :367770   Max.   : 44.200  
##                                    NA's   :1        NA's   :1        
##       PopD             Area            Age0yrs           Age15yrs     
##  Min.   :  21.1   Min.   :  0.220   Min.   :   25.0   Min.   :    60  
##  1st Qu.: 574.3   1st Qu.:  1.350   1st Qu.:  416.2   1st Qu.:  2192  
##  Median :1974.3   Median :  2.490   Median :  552.5   Median :  3108  
##  Mean   :2137.9   Mean   : 17.624   Mean   : 1237.2   Mean   :  6198  
##  3rd Qu.:2877.9   3rd Qu.:  4.975   3rd Qu.:  735.0   3rd Qu.:  4121  
##  Max.   :7582.6   Max.   :696.150   Max.   :48255.0   Max.   :241740  
##  NA's   :1                          NA's   :1         NA's   :1       
##     Age65yrs           AvgAge          Detach            Apart        
##  Min.   :   10.0   Min.   :29.20   Min.   :   25.0   Min.   :    0.0  
##  1st Qu.:  622.5   1st Qu.:40.60   1st Qu.:  482.5   1st Qu.:    0.0  
##  Median :  930.0   Median :43.85   Median :  835.0   Median :    5.0  
##  Mean   : 1993.8   Mean   :43.65   Mean   : 1647.0   Mean   :  266.9  
##  3rd Qu.: 1418.8   3rd Qu.:46.58   3rd Qu.: 1130.0   3rd Qu.:   85.0  
##  Max.   :77770.0   Max.   :56.50   Max.   :64235.0   Max.   :10400.0  
##  NA's   :1         NA's   :1       NA's   :1         NA's   :1        
##     Divorced          Widowed            MdInc        MDAftTaxInc   
##  Min.   :    0.0   Min.   :    0.0   Min.   :12384   Min.   :12192  
##  1st Qu.:  196.2   1st Qu.:  135.0   1st Qu.:34056   1st Qu.:30667  
##  Median :  302.5   Median :  215.0   Median :37481   Median :33380  
##  Mean   :  671.4   Mean   :  476.2   Mean   :37557   Mean   :33388  
##  3rd Qu.:  473.8   3rd Qu.:  328.8   3rd Qu.:42240   3rd Qu.:37328  
##  Max.   :26195.0   Max.   :18565.0   Max.   :57557   Max.   :48939  
##  NA's   :1         NA's   :1         NA's   :4       NA's   :4      
##     MVisMin            Owner              Renter             Band        
##  Min.   :   0.00   Min.   :    20.0   Min.   :    0.0   Min.   :  0.000  
##  1st Qu.:  10.00   1st Qu.:   937.5   1st Qu.:  308.8   1st Qu.:  0.000  
##  Median :  20.00   Median :  1277.5   Median :  602.5   Median :  0.000  
##  Mean   :  38.78   Mean   :  2612.4   Mean   : 1557.2   Mean   :  2.821  
##  3rd Qu.:  28.75   3rd Qu.:  1690.0   3rd Qu.: 1071.2   3rd Qu.:  0.000  
##  Max.   :1505.00   Max.   :101885.0   Max.   :60745.0   Max.   :105.000  
##  NA's   :1         NA's   :1          NA's   :1         NA's   :1        

## merge csv join to shape file by column GUID
crd.data <- merge(tracts, census.16, by = "GUID")

## remove 'na' values in data column: MInc
# '!' means 'not'; !is.na will keep everything that is not na
crd.data <- crd.data[!is.na(crd.data$MInc),]

To understand what type of data, use class() to see this a “SpatialPolygonsDataFrame”, and can be used by the sp package. Summary() will produce a table of the dataframe summary statistics

class(crd.data)

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

summary(crd.data)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##       min     max
## x 3903547 3974156
## y 1910542 1962451
## Is projected: TRUE 
## proj4string :
## [+proj=lcc +lat_1=49 +lat_2=77 +lat_0=63.390675
## +lon_0=-91.86666666666666 +x_0=6200000 +y_0=3000000 +ellps=GRS80
## +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Data attributes:
## some summary stats have been ommitted ... 
##       GUID                CTUID         CUID            Pop       
##  Min.   :9350001   9350001.00: 1   Min.   :  1.0   Min.   :  332  
##  1st Qu.:9350101   9350002.00: 1   1st Qu.:101.2   1st Qu.: 3538  
##  Median :9350126   9350003.01: 1   Median :126.5   Median : 4744  
##  Mean   :9350108   9350003.02: 1   Mean   :108.4   Mean   : 4964  
##  3rd Qu.:9350153   9350004.00: 1   3rd Qu.:152.8   3rd Qu.: 6156  
##  Max.   :9350180   9350005.00: 1   Max.   :180.1   Max.   :10393  
##                    (Other)   :68                                  
##      PopCh             PopD             Area            Age0yrs      
##  Min.   :-2.000   Min.   :  21.1   Min.   :  0.500   Min.   : 105.0  
##  1st Qu.: 1.800   1st Qu.: 641.2   1st Qu.:  1.383   1st Qu.: 422.5  
##  Median : 4.250   Median :2025.1   Median :  2.500   Median : 560.0  
##  Mean   : 5.791   Mean   :2237.2   Mean   :  9.355   Mean   : 650.7  
##  3rd Qu.: 7.300   3rd Qu.:2899.8   3rd Qu.:  4.992   3rd Qu.: 735.0  
##  Max.   :44.200   Max.   :7582.6   Max.   :221.580   Max.   :2000.0  
##                                                                      
##     Age15yrs       Age65yrs          AvgAge          Detach      
##  Min.   : 210   Min.   :  20.0   Min.   :29.20   Min.   :  70.0  
##  1st Qu.:2232   1st Qu.: 646.2   1st Qu.:41.30   1st Qu.: 503.8  
##  Median :3115   Median : 965.0   Median :43.90   Median : 865.0  
##  Mean   :3263   Mean   :1050.1   Mean   :43.99   Mean   : 866.1  
##  3rd Qu.:4121   3rd Qu.:1418.8   3rd Qu.:46.67   3rd Qu.:1130.0  
##  Max.   :7300   Max.   :2410.0   Max.   :56.50   Max.   :2715.0  
##    Separated        Divorced         Widowed          MdInc      
##  Min.   :  5.0   Min.   :  10.0   Min.   :  5.0   Min.   :12384  
##  1st Qu.: 60.0   1st Qu.: 200.0   1st Qu.:138.8   1st Qu.:33878  
##  Median : 95.0   Median : 310.0   Median :215.0   Median :37563  
##  Mean   :117.0   Mean   : 353.6   Mean   :250.9   Mean   :37558  
##  3rd Qu.:158.8   3rd Qu.: 473.8   3rd Qu.:328.8   3rd Qu.:42250  
##  Max.   :380.0   Max.   :1065.0   Max.   :670.0   Max.   :57557  
##                                                                  
##    Vminority         MVisMin          Owner            Renter      
##  Min.   : 0.000   Min.   : 0.00   Min.   :  65.0   Min.   :  15.0  
##  1st Qu.: 0.000   1st Qu.:10.00   1st Qu.: 966.2   1st Qu.: 325.0  
##  Median : 0.000   Median :20.00   Median :1287.5   Median : 620.0  
##  Mean   : 5.608   Mean   :20.54   Mean   :1375.5   Mean   : 820.3  
##  3rd Qu.:10.000   3rd Qu.:28.75   3rd Qu.:1690.0   3rd Qu.:1071.2  
##  Max.   :35.000   Max.   :80.00   Max.   :3255.0   Max.   :3945.0  

V.i. Make Chloropleth Map: Median Income

Using the tmap package …

# make Chloropleth map
map_PopDens <- tm_shape(crd.data) + 
  tm_polygons(col = "MdInc", 
              title = "Median Income", 
              style = "jenks", 
              palette = "viridis", n = 6) + 
  tm_layout(title = "CRD Census Tracts Median Income", title.position = c("LEFT", "TOP"))

map_PopDens

# create png of Median Income Population Density
png("map_PopDens_MdInc-2.png")  
map_PopDens
dev.off()

VI.i. Defining Neighbourhoods: Median Income

1. Using the spdep package,
   

1. poly2nb builds a list of neighbours from polygon regions whose boundaries are continguous: https://www.rdocumentation.org/packages/spdep/versions/1.1-3/topics/poly2nb
2. nb2lines stores spatial coordinates as vectors (needed for weights): https://www.rdocumentation.org/packages/spdep/versions/1.1-3/topics/nb2lines
3. queen - FALSE defines the neighbourhood as Rook https://cran.r-project.org/web/packages/spdep/vignettes/nb.pdf

## Defining Neighbourhoods
# Queen Neighbours
crd.nb <- poly2nb(crd.data)
crd.net <- nb2lines(crd.nb,coords=coordinates(crd.data))

map_crd.net <- tm_shape(crd.data) + tm_borders(col='lightgrey') + 
  tm_shape(crd.net) + tm_lines(col='red') +
  tm_layout(title = "CRD Census Tracts Median Income Queen Neighbours", title.position = c("LEFT", "TOP"))

##
# Rook Neighbours
crd.nb2 <- poly2nb(crd.data, queen = FALSE)
crd.net2 <- nb2lines(crd.nb2,coords=coordinates(crd.data))

map_crd.net2 <- tm_shape(crd.data) + tm_borders(col='lightgrey') + 
  tm_shape(crd.net) + tm_lines(col='blue', lwd = 2) +
  tm_shape(crd.net2) + tm_lines(col='yellow', lwd = 2) +
    tm_layout(title = "CRD Census Tracts Median Income Rook Neighbours", title.position = c("LEFT", "TOP"))

map_crd.net

map_crd.net2

# create png of Median Income CRD net
png("map_map_crd.net_MdInc-2.png")  
map_crd.net

dev.off()

# create png of Median Income CRD net 2
png("map_map_crd.net2_MdInc-2.png")  
map_crd.net2

dev.off()

VII.i. Weight Matrix: Median Income

1. Using the spdep package, nb2listw uses the neighbours list previously created and adds weight: https://www.rdocumentation.org/packages/spdep/versions/1.1-3/topics/nb2listw

# Weight Matrix
crd.lw <- nb2listw(crd.nb, zero.policy = TRUE, style = "W")
print.listw(crd.lw, zero.policy = TRUE)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 74 
## Number of nonzero links: 424 
## Percentage nonzero weights: 7.742878 
## Average number of links: 5.72973 
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 74 5476 74 29.46931 322.5924

VIII.i. Calculate Lag Means: Median Income

1. Using the spdep package, lag.listw : https://www.rdocumentation.org/packages/spdep/versions/1.1-3/topics/lag.listw
   
2. Creates a map using tmap

# Calculate Lag Means
crd.data$IncLagMeans = lag.listw(crd.lw, crd.data$MdInc, zero.policy = TRUE)

map_LagMean <- tm_shape(crd.data) + 
  tm_polygons(col = "IncLagMeans", 
              title = "Median Income\nLagged Means", 
              style = "jenks", 
              palette = "viridis", n = 6) +
  tm_layout(title = "CRD Census Tracts Median Income Lagged Means", title.position = c("LEFT", "TOP"))

map_LagMean

# create png of map LagMean Median Income CRD net 2
png("map_LagMean_MdInc-2.png")  
map_LagMean
dev.off()

IX.i. Morans I test: Median Income

Morans Global I Test creates one i value for every polygon as if they were all one polygon

# Morans I test
mi <- moran.test(crd.data$MdInc, crd.lw, zero.policy = TRUE)
mi

## 
##  Moran I test under randomisation
## 
## data:  crd.data$MdInc  
## weights: crd.lw    
## 
## Moran I statistic standard deviate = 2.2665, p-value = 0.01171
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.143170949      -0.013698630       0.004790316

moran.range <- function(lw) {
  wmat <- listw2mat(lw)
  return(range(eigen((wmat + t(wmat))/2)$values))
}
moran.range(crd.lw)

## [1] -0.8193252  1.0821054

## calculate z value
mI <- mi$estimate[[1]]
eI <- mi$estimate[[2]]
var <- mi$estimate[[3]]

# calculate z value 
z <- (mI - eI) / var^(1/2)
z

## [1] 2.266505

# Add Stats data objects to data.frame ########
data.for.table1 = data.frame(mI, eI, var, z)

data.for.table1

##          mI          eI         var        z
## 1 0.1431709 -0.01369863 0.004790316 2.266505

table1 <- data.for.table1 

# write csv file with Morans Global results for I, E, Var, Z, P
write.csv(table1, "MoransGlobal-MdInc.csv", row.names = FALSE)

X.i. Local Moran’s i Test (Lisa Test): Median Income

Separate i values are created for every single polygon

# Lisa Test
lisa.test <- localmoran(crd.data$MdInc, crd.lw)

crd.data$Ii <- lisa.test[,1]
crd.data$E.Ii<- lisa.test[,2]
crd.data$Var.Ii<- lisa.test[,3]
crd.data$Z.Ii<- lisa.test[,4]
crd.data$P<- lisa.test[,5]

# make a data frome from the Lisa Test Data
data.for.table2 = data.frame(crd.data$Ii, crd.data$E.Ii, crd.data$Var.Ii, crd.data$Z.Ii, crd.data$P)
table2 <- data.for.table2
table2

##      crd.data.Ii crd.data.E.Ii crd.data.Var.Ii crd.data.Z.Ii   crd.data.P
## 1   0.2744151305   -0.01369863      0.14693811   0.751617311 2.261406e-01
## 2   0.0557309861   -0.01369863      0.22656398   0.145864407 4.420142e-01
## 3  -0.2349448907   -0.01369863      0.14693811  -0.577176596 7.180899e-01
## 4   0.3326591723   -0.01369863      0.22656398   0.727661741 2.334103e-01
## 5   0.1990156088   -0.01369863      0.17878846   0.503068148 3.074582e-01
## 6  -0.0338955771   -0.01369863      0.14693811  -0.052688823 5.210101e-01
## 7   0.0128136506   -0.01369863      0.17878846   0.062701416 4.750021e-01
## 8   0.7919282302   -0.01369863      0.17878846   1.905303636 2.837031e-02
## 9  -0.0464760893   -0.01369863      0.17878846  -0.077518533 5.308945e-01
## 10  0.4657065327   -0.01369863      0.12418786   1.360389058 8.685342e-02
## 11  0.0067879072   -0.01369863      0.10712518   0.062592570 4.750455e-01
## 12  0.2947025748   -0.01369863      0.09385420   1.006674841 1.570455e-01
## 13 -0.0923032394   -0.01369863      0.46544159  -0.115216696 5.458633e-01
## 14 -0.0285935895   -0.01369863      0.30618985  -0.026918070 5.107375e-01
## 15  0.4126788509   -0.01369863      0.09385420   1.391769798 8.199606e-02
## 16 -0.2187477999   -0.01369863      0.10712518  -0.626487255 7.345023e-01
## 17  0.4433240328   -0.01369863      0.12418786   1.296875124 9.733708e-02
## 18  1.2205874502   -0.01369863      0.30618985   2.230593514 1.285403e-02
## 19  0.4616826274   -0.01369863      0.17878846   1.124274380 1.304483e-01
## 20  0.2638695873   -0.01369863      0.12418786   0.787644346 2.154524e-01
## 21  2.1078481863   -0.01369863      0.22656398   4.457149340 4.152837e-06
## 22 -0.0124306146   -0.01369863      0.17878846   0.002998851 4.988036e-01
## 23  0.0980896658   -0.01369863      0.08323742   0.387469072 3.492045e-01
## 24 -0.0231942136   -0.01369863      0.22656398  -0.019949234 5.079581e-01
## 25  0.1047856106   -0.01369863      0.12418786   0.336218041 3.683532e-01
## 26  0.0246750827   -0.01369863      0.10712518   0.117243303 4.533336e-01
## 27  0.1900350829   -0.01369863      0.14693811   0.531490704 2.975394e-01
## 28  0.0019320024   -0.01369863      0.08323742   0.054177288 4.783970e-01
## 29  0.9284834191   -0.01369863      0.17878846   2.228255999 1.293172e-02
## 30  0.2032721811   -0.01369863      0.17878846   0.513134921 3.039285e-01
## 31  0.1259530137   -0.01369863      0.12418786   0.396283943 3.459478e-01
## 32  0.4535832667   -0.01369863      0.17878846   1.105119431 1.345539e-01
## 33  0.0760594107   -0.01369863      0.14693811   0.234156457 4.074318e-01
## 34  0.6622247407   -0.01369863      0.30618985   1.221524176 1.109438e-01
## 35  0.2429889826   -0.01369863      0.17878846   0.607064965 2.719039e-01
## 36  0.0356433844   -0.01369863      0.17878846   0.116693626 4.535514e-01
## 37 -0.0004154355   -0.01369863      0.17878846   0.031414691 4.874694e-01
## 38 -0.0193496543   -0.01369863      0.12418786  -0.016035688 5.063970e-01
## 39  0.0074249521   -0.01369863      0.17878846   0.049957170 4.800783e-01
## 40  0.1856654215   -0.01369863      0.22656398   0.418843150 3.376654e-01
## 41  0.0595949852   -0.01369863      0.12418786   0.207982392 4.176214e-01
## 42 -0.2359649540   -0.01369863      0.12418786  -0.630716351 7.358870e-01
## 43 -0.0579399062   -0.01369863      0.14693811  -0.115414512 5.459417e-01
## 44  0.1041228751   -0.01369863      0.30618985   0.212926233 4.156923e-01
## 45  0.1458425508   -0.01369863      0.14693811   0.416203354 3.386306e-01
## 46 -0.1574332175   -0.01369863      0.17878846  -0.339931605 6.330460e-01
## 47  0.0536021511   -0.01369863      0.22656398   0.141391946 4.437802e-01
## 48  0.6232149028   -0.01369863      0.14693811   1.661549369 4.830158e-02
## 49  0.5484882380   -0.01369863      0.10712518   1.717650982 4.293015e-02
## 50 -0.0440773555   -0.01369863      0.02923030  -0.177685867 5.705152e-01
## 51 -1.9458534186   -0.01369863      0.30618985  -3.491777155 9.997601e-01
## 52 -0.3484095342   -0.01369863      0.14693811  -0.873177696 8.087169e-01
## 53 -0.0252503561   -0.01369863      0.12418786  -0.032779876 5.130749e-01
## 54  1.0843288469   -0.01369863      0.14693811   2.864481232 2.088465e-03
## 55  0.2836125755   -0.01369863      0.10712518   0.908375690 1.818399e-01
## 56  0.0384246510   -0.01369863      0.22656398   0.109505596 4.564007e-01
## 57  0.0025467591   -0.01369863      0.14693811   0.042380189 4.830978e-01
## 58 -0.0346348845   -0.01369863      0.22656398  -0.043984894 5.175418e-01
## 59 -3.7243153785   -0.01369863      0.94319681  -3.820717593 9.999335e-01
## 60  0.0117550365   -0.01369863      0.09385420   0.083085168 4.668919e-01
## 61  0.2614551227   -0.01369863      0.22656398   0.578069434 2.816086e-01
## 62  0.1399328018   -0.01369863      0.46544159   0.225189161 4.109161e-01
## 63  0.4384667710   -0.01369863      0.46544159   0.662772886 2.537380e-01
## 64 -0.3612599654   -0.01369863      0.94319681  -0.357874121 6.397812e-01
## 65  0.2206921127   -0.01369863      0.14693811   0.611467288 2.704451e-01
## 66  0.4853407163   -0.01369863      0.14693811   1.301869827 9.648044e-02
## 67  0.3210323517   -0.01369863      0.22656398   0.703235000 2.409547e-01
## 68 -0.0802313849   -0.01369863      0.14693811  -0.173567448 5.688973e-01
## 69  0.3159098132   -0.01369863      0.12418786   0.935316835 1.748125e-01
## 70  0.6415313555   -0.01369863      0.12418786   1.859320199 3.149088e-02
## 71  0.4949972138   -0.01369863      0.17878846   1.203063216 1.144759e-01
## 72 -0.0146686830   -0.01369863      0.30618985  -0.001753073 5.006994e-01
## 73  0.2927142780   -0.01369863      0.09385420   1.000184696 1.586106e-01
## 74  1.0828675345   -0.01369863      0.30618985   1.981707007 2.375602e-02

head(table2)

##   crd.data.Ii crd.data.E.Ii crd.data.Var.Ii crd.data.Z.Ii crd.data.P
## 1  0.27441513   -0.01369863       0.1469381    0.75161731  0.2261406
## 2  0.05573099   -0.01369863       0.2265640    0.14586441  0.4420142
## 3 -0.23494489   -0.01369863       0.1469381   -0.57717660  0.7180899
## 4  0.33265917   -0.01369863       0.2265640    0.72766174  0.2334103
## 5  0.19901561   -0.01369863       0.1787885    0.50306815  0.3074582
## 6 -0.03389558   -0.01369863       0.1469381   -0.05268882  0.5210101

# write csv file with Lisa Test results for I, E, Var, Z, P
write.csv(table2, "LisaTest-MdInc.csv", row.names = FALSE)

XI.i. Map Lisa: Median Income

# Map Lisa
map_LISA <- tm_shape(crd.data) + 
  tm_polygons(col = "Ii", 
              title = "Local Moran's I", 
              style = "jenks", 
              palette = "viridis", n = 6)  +
  tm_layout(title = "CRD Census Tracts Median Income LISA Local Moran's I", title.position = c("LEFT", "TOP"))

map_LISA

# create png of Map LISA MdInc
png("map_LISA_MdInc-2.png")  
map_LISA
dev.off()

XII.i. Scatter Plot: Median Income

when running the code for the second variable, the new column name will need to replace this item crd.data$MdInc

# Scatter Plot
moran.plot(crd.data$MdInc, crd.lw, zero.policy=NULL, spChk=NULL, labels=NULL, xlab="Population Density", 
           ylab="Spatially Lagged Population Density", quiet=NULL, main ="Moran Plot for CRD Median Income")

# create png of Moran scatter plot
png("moran.plot_MdInc-2.png")  
moran.plot(crd.data$MdInc, crd.lw, zero.policy=NULL, spChk=NULL, labels=NULL, xlab="Population Density", 
           ylab="Spatially Lagged Population Density", quiet=NULL, main ="Moran Plot for CRD Median Income")
dev.off()

V.ii. Make Chloropleth Map: Renter

Using the tmap package …

# make Chloropleth map
map_PopDens <- tm_shape(crd.data) + 
  tm_polygons(col = "Renter", 
              title = "Renter", 
              style = "jenks", 
              palette = "viridis", n = 6) + 
  tm_layout(title = "CRD Census Tracts Renter", title.position = c("LEFT", "TOP"))

map_PopDens

# create png of Renter Population Density
png("map_PopDens_Renter-2.png")  
map_PopDens
dev.off()

VI.ii. Defining Neighbourhoods: Renter

1. Using the spdep package,
   

1. poly2nb builds a list of neighbours from polygon regions whose boundaries are continguous: https://www.rdocumentation.org/packages/spdep/versions/1.1-3/topics/poly2nb
   
2. nb2lines stores spatial coordinates as vectors (needed for weights): https://www.rdocumentation.org/packages/spdep/versions/1.1-3/topics/nb2lines
   

2. Using the tmap package …

# Defining Neighbourhoods
crd.nb <- poly2nb(crd.data)
crd.net <- nb2lines(crd.nb,coords=coordinates(crd.data))

map_crd.net <- tm_shape(crd.data) + tm_borders(col='lightgrey') + 
  tm_shape(crd.net) + tm_lines(col='red')

##

crd.nb2 <- poly2nb(crd.data, queen = FALSE)
crd.net2 <- nb2lines(crd.nb2,coords=coordinates(crd.data))

map_crd.net2 <- tm_shape(crd.data) + tm_borders(col='lightgrey') + 
  tm_shape(crd.net) + tm_lines(col='blue', lwd = 2) +
  tm_shape(crd.net2) + tm_lines(col='yellow', lwd = 2)

map_crd.net

map_crd.net2

# create png of Renter CRD net
png("map_map_crd.net_Renter-2.png")  
map_crd.net

dev.off()

# create png of Median Income CRD net 2
png("map_map_crd.net2_Renter-2.png")  
map_crd.net2

dev.off()

VII.ii. Weight Matrix: Renter

1. Using the spdep package, nb2listw uses the neighbours list previously created and adds weight: https://www.rdocumentation.org/packages/spdep/versions/1.1-3/topics/nb2listw

# Weight Matrix
crd.lw <- nb2listw(crd.nb, zero.policy = TRUE, style = "W")
print.listw(crd.lw, zero.policy = TRUE)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 74 
## Number of nonzero links: 424 
## Percentage nonzero weights: 7.742878 
## Average number of links: 5.72973 
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 74 5476 74 29.46931 322.5924

VIII.ii. Calculate Lag Means: Renter

1. Using the spdep package, lag.listw : https://www.rdocumentation.org/packages/spdep/versions/1.1-3/topics/lag.listw
   
2. Creates a map using tmap

# Calculate Lag Means
crd.data$RentLagMeans = lag.listw(crd.lw, crd.data$Renter, zero.policy = TRUE)

map_LagMean <- tm_shape(crd.data) + 
  tm_polygons(col = "RentLagMeans", 
              title = "Renter\nLagged Means", 
              style = "jenks", 
              palette = "viridis", n = 6) +
  tm_layout(title = "CRD Census Tracts Renter Lagged Means", title.position = c("LEFT", "TOP"))

map_LagMean

# create png of map LagMean Renter CRD net 2
png("map_LagMean_Renter-2.png")  
map_LagMean
dev.off()

IX.ii. Morans I test: Renter

Morans Global I Test creates one i value for every polygon as if they were all one polygon.

# Morans I test
mi <- moran.test(crd.data$Renter, crd.lw, zero.policy = TRUE)
mi

## 
##  Moran I test under randomisation
## 
## data:  crd.data$Renter  
## weights: crd.lw    
## 
## Moran I statistic standard deviate = 7.1205, p-value = 5.378e-13
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.472963854      -0.013698630       0.004671315

moran.range <- function(lw) {
  wmat <- listw2mat(lw)
  return(range(eigen((wmat + t(wmat))/2)$values))
}
moran.range(crd.lw)

## [1] -0.8193252  1.0821054

## calculate z value
mI <- mi$estimate[[1]]
eI <- mi$estimate[[2]]
var <- mi$estimate[[3]]

# calculate z value 
z <- (mI - eI) / var^(1/2)
z

## [1] 7.120464

# Add Stats data objects to data.frame ########
data.for.table3 = data.frame(mI, eI, var, z)

data.for.table3

##          mI          eI         var        z
## 1 0.4729639 -0.01369863 0.004671315 7.120464

table3 <- data.for.table3 

# write csv file with Morans Global results for I, E, Var, Z, P
write.csv(table3, "MoransGlobal-Renter.csv", row.names = FALSE)

X.ii. Local Moran’s i Test (Lisa Test): Renter

Separate i values are created for every single polygon.

# Lisa Test
lisa.test <- localmoran(crd.data$Renter, crd.lw)

crd.data$Ii <- lisa.test[,1]
crd.data$E.Ii<- lisa.test[,2]
crd.data$Var.Ii<- lisa.test[,3]
crd.data$Z.Ii<- lisa.test[,4]
crd.data$P<- lisa.test[,5]

# make a data frome from the Lisa Test Data
data.for.table4 = data.frame(crd.data$Ii, crd.data$E.Ii, crd.data$Var.Ii, crd.data$Z.Ii, crd.data$P)
table4 <- data.for.table4
table4

##     crd.data.Ii crd.data.E.Ii crd.data.Var.Ii crd.data.Z.Ii   crd.data.P
## 1  -0.051377674   -0.01369863      0.14362507   -0.09942251 5.395986e-01
## 2   0.049185253   -0.01369863      0.22127380    0.13368244 4.468268e-01
## 3   0.386073232   -0.01369863      0.14362507    1.05486551 1.457434e-01
## 4   0.348725220   -0.01369863      0.22127380    0.77046297 2.205127e-01
## 5  -0.140094336   -0.01369863      0.17468456   -0.30241620 6.188326e-01
## 6   0.608795762   -0.01369863      0.14362507    1.64255649 5.023736e-02
## 7   2.874422715   -0.01369863      0.17468456    6.91016106 2.420519e-12
## 8   0.311890329   -0.01369863      0.17468456    0.77900887 2.179872e-01
## 9   0.208164120   -0.01369863      0.17468456    0.53083204 2.977676e-01
## 10 -0.038749345   -0.01369863      0.12143972   -0.07188524 5.286534e-01
## 11 -0.037441635   -0.01369863      0.10480070   -0.07334214 5.292331e-01
## 12 -0.024326974   -0.01369863      0.09185925   -0.03506745 5.139870e-01
## 13  0.097416695   -0.01369863      0.45422001    0.16486969 4.345233e-01
## 14  0.107618104   -0.01369863      0.29892254    0.22189186 4.121990e-01
## 15  5.129784538   -0.01369863      0.09185925   16.97054842 6.782961e-65
## 16  2.957535302   -0.01369863      0.10480070    9.17814141 2.192979e-20
## 17  2.348206689   -0.01369863      0.12143972    6.77769634 6.105358e-12
## 18  0.499077361   -0.01369863      0.29892254    0.93788233 1.741524e-01
## 19 -0.051886856   -0.01369863      0.17468456   -0.09136970 5.364006e-01
## 20 -0.068072743   -0.01369863      0.12143972   -0.15603133 5.619958e-01
## 21  0.520825597   -0.01369863      0.22127380    1.13632456 1.279104e-01
## 22  1.950423301   -0.01369863      0.17468456    4.69938665 1.304720e-06
## 23  0.371016761   -0.01369863      0.08150608    1.34754893 8.890175e-02
## 24 -0.123331272   -0.01369863      0.22127380   -0.23306383 5.921441e-01
## 25  0.355195058   -0.01369863      0.12143972    1.05857308 1.448971e-01
## 26  0.579965194   -0.01369863      0.10480070    1.83382751 3.333980e-02
## 27  0.221349851   -0.01369863      0.14362507    0.62021508 2.675581e-01
## 28  0.230017806   -0.01369863      0.08150608    0.85366957 1.966440e-01
## 29  0.463223116   -0.01369863      0.17468456    1.14108989 1.269163e-01
## 30  0.048185213   -0.01369863      0.17468456    0.14806418 4.411461e-01
## 31  0.004901072   -0.01369863      0.12143972    0.05337349 4.787172e-01
## 32  1.831026031   -0.01369863      0.17468456    4.41371501 5.080585e-06
## 33 -0.006214897   -0.01369863      0.14362507    0.01974709 4.921226e-01
## 34  0.111073161   -0.01369863      0.29892254    0.22821127 4.097410e-01
## 35  0.081964112   -0.01369863      0.17468456    0.22888407 4.094795e-01
## 36  4.103911859   -0.01369863      0.17468456    9.85185464 3.364532e-23
## 37  0.026249068   -0.01369863      0.17468456    0.09557944 4.619273e-01
## 38  0.025537837   -0.01369863      0.12143972    0.11259251 4.551768e-01
## 39  0.021262740   -0.01369863      0.17468456    0.08364908 4.666677e-01
## 40  0.560730223   -0.01369863      0.22127380    1.22115627 1.110134e-01
## 41  0.386694243   -0.01369863      0.12143972    1.14896278 1.252857e-01
## 42  0.116277949   -0.01369863      0.12143972    0.37297930 3.545819e-01
## 43  0.499152713   -0.01369863      0.14362507    1.35324480 8.798872e-02
## 44  0.424776079   -0.01369863      0.29892254    0.80198310 2.112814e-01
## 45 -0.113705920   -0.01369863      0.14362507   -0.26388611 6.040662e-01
## 46  0.590842908   -0.01369863      0.17468456    1.44643486 7.402763e-02
## 47  0.254866474   -0.01369863      0.22127380    0.57093226 2.840228e-01
## 48  0.416745105   -0.01369863      0.14362507    1.13579843 1.280205e-01
## 49  0.071323146   -0.01369863      0.10480070    0.26263226 3.964170e-01
## 50  0.172584570   -0.01369863      0.02883999    1.09692304 1.363375e-01
## 51  0.614914252   -0.01369863      0.29892254    1.14975140 1.251231e-01
## 52  0.424731981   -0.01369863      0.14362507    1.15687315 1.236621e-01
## 53  0.175231168   -0.01369863      0.12143972    0.54215078 2.938573e-01
## 54  1.193773325   -0.01369863      0.14362507    3.18611849 7.209776e-04
## 55 -0.021993901   -0.01369863      0.10480070   -0.02562409 5.102214e-01
## 56  0.371705203   -0.01369863      0.22127380    0.81931523 2.063033e-01
## 57 -0.081655662   -0.01369863      0.14362507   -0.17931610 5.711552e-01
## 58 -0.251859397   -0.01369863      0.22127380   -0.50629684 6.936759e-01
## 59  0.865291582   -0.01369863      0.92011241    0.91635467 1.797405e-01
## 60  0.072817501   -0.01369863      0.09185925    0.28545368 3.876483e-01
## 61 -0.475737466   -0.01369863      0.22127380   -0.98223064 8.370069e-01
## 62  0.009473931   -0.01369863      0.45422001    0.03438277 4.862860e-01
## 63  0.335475911   -0.01369863      0.45422001    0.51809505 3.021960e-01
## 64  0.218231309   -0.01369863      0.92011241    0.24178891 4.044719e-01
## 65  0.013301895   -0.01369863      0.14362507    0.07124544 4.716012e-01
## 66  0.262752096   -0.01369863      0.14362507    0.72946189 2.328596e-01
## 67  0.675204036   -0.01369863      0.22127380    1.46451176 7.152709e-02
## 68 -0.063642698   -0.01369863      0.14362507   -0.13178585 5.524232e-01
## 69 -0.008119337   -0.01369863      0.12143972    0.01601028 4.936131e-01
## 70  0.003457398   -0.01369863      0.12143972    0.04923074 4.803677e-01
## 71  0.239976334   -0.01369863      0.17468456    0.60694640 2.719432e-01
## 72  0.280314447   -0.01369863      0.29892254    0.53775854 2.953719e-01
## 73  0.038540937   -0.01369863      0.09185925    0.17236065 4.315770e-01
## 74  0.395325464   -0.01369863      0.29892254    0.74811706 2.271948e-01

head(table4)

##   crd.data.Ii crd.data.E.Ii crd.data.Var.Ii crd.data.Z.Ii crd.data.P
## 1 -0.05137767   -0.01369863       0.1436251   -0.09942251 0.53959860
## 2  0.04918525   -0.01369863       0.2212738    0.13368244 0.44682685
## 3  0.38607323   -0.01369863       0.1436251    1.05486551 0.14574342
## 4  0.34872522   -0.01369863       0.2212738    0.77046297 0.22051266
## 5 -0.14009434   -0.01369863       0.1746846   -0.30241620 0.61883260
## 6  0.60879576   -0.01369863       0.1436251    1.64255649 0.05023736

# write csv file with Lisa Test results for I, E, Var, Z, P
write.csv(table4, "LisaTest-Renter.csv", row.names = FALSE)

XI.ii. Map Lisa: Renter

Using the tmap package …

# Map Lisa
map_LISA <- tm_shape(crd.data) + 
  tm_polygons(col = "Ii", 
              title = "Local Moran's I", 
              style = "jenks", 
              palette = "viridis", n = 6)  +
  tm_layout(title = "CRD Census Tracts Renter LISA Local Moran's I", title.position = c("LEFT", "TOP"))

map_LISA

# create png of Map LISA Renter
png("map_LISA_Renter-2.png")  
map_LISA
dev.off()

XII.ii. Scatter Plot: Renter

when running the code for the second variable, the new column name will need to replace this item crd.data$Renter

# Scatter Plot
moran.plot(crd.data$Renter, crd.lw, zero.policy=NULL, spChk=NULL, labels=NULL, xlab="Population Density", 
           ylab="Spatially Lagged Population Density", quiet=NULL, main ="Moran Plot for CRD Renter")

# create png of Moran scatter plot
png("moran.plot_Renter-2.png")  
moran.plot(crd.data$Renter, crd.lw, zero.policy=NULL, spChk=NULL, labels=NULL, xlab="Population Density", 
           ylab="Spatially Lagged Population Density", quiet=NULL, main ="Moran Plot for CRD Renter")
dev.off()

References

3. The R script and this document were created with R Markdown, using RStudio v.1.1.463 and R v.3.6.1.