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"
# dir <- "Z:/Geog418-Assign-3/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, using the common column: GUID, assigning this merged data a new variable name.
   
3. Remove any missing ‘na’ values from the merged data, reassigning these
   
4. Summary(census.16) provides some summary statistics for each of the variable columns of the Census Tract Data. Note the summary also records the number of missing NA values in each column of data.

## 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        
##     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  
##    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 …

1. To select a colour palette to create the chloroplet map, use tmaptools::palette_explorer(). To generate the colour code needed for the map palette, check the Code Generator button for tmap layer function code - the code will be found under the colour palette. Note, the palette explorer window must be closed in order for R to continue processing the code.
2. The shape of the polygons comes from the merged crd.data file. The classified polygons are coloured by the values in the MdInc column. The title is for the legend. The classification style is Jenks uses natural breaks to set up the best classification of variable values. The layout includes the title of the plot, and the position of the title (note the positions need to be in capital letter format).

1. To place scale bar position https://www.rdocumentation.org/packages/tmap/versions/2.3/topics/tm_scale_bar 2) To make room for north arrow and scale bar, place margins around map and image boundary (bottom, left, top, right) https://www.rdocumentation.org/packages/tmap/versions/1.2/topics/tm_layout

3. map_PopDens prints the plot in R window. In order to save the map, the file object png is called using a filename to save to, call the plot, then close the graphic device connection with dev.off() in order to return to R.

1. dev.off() must be used when saving a plot, in order to close the plot window, and to save the file. For Mac OSX users, ‘quartz_off_screen’ is a graphic device driver. Calling dev.off() again should return user to the RStudio graphic device or ‘null device’. https://support.rstudio.com/hc/en-us/community/posts/200663616-plots-not-showing-up, https://stackoverflow.com/questions/23621012/display-and-save-the-plot-simultaneously-in-r-rstudio, https://www.rdocumentation.org/packages/grDevices/versions/3.6.1/topics/quartz

# tmap tool for selecting a colour palette
# the code for palette explorer is commented-out from the general code, as the code will not continue until the palette explorer window is closed
#
# tmaptools::palette_explorer()

# make Chloropleth map
map_PopDens <- tm_shape(crd.data) + 
  tm_polygons(col = "MdInc", 
              title = "Median $ Income", 
              style = "jenks", 
              palette = "viridis", n = 6) + 
  # add scale bar - first value is TOP, second value is BOTTOM
  tm_scale_bar(width = 0.22, position = c("RIGHT", "BOTTOM")) + 
  # add compass
  tm_compass(position = c("RIGHT", "TOP")) + 
  tm_layout(title = "CRD Census Tracts Median $ Income 2016", title.position = c("LEFT", "TOP"), inner.margins = c(.08, .03, .08, .03))

map_PopDens

# create png of Median Income 
png("map_PopDens_MdInc_2.png")  
map_PopDens
dev.off()

## quartz_off_screen 
##                 2

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') +
  # add scale bar - first value is TOP, second value is BOTTOM
  tm_scale_bar(width = 0.22, position = c("RIGHT", "BOTTOM")) + 
  # add compass
  tm_compass(position = c("RIGHT", "TOP")) + 
  tm_layout(title = "CRD Census Tracts Median Income\nQueen Neighbours 2016", title.position = c("LEFT", "TOP"), inner.margins = c(.08, .03, .08, .03))

##
# 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='red', lwd = 2) +
  tm_shape(crd.net2) + tm_lines(col='blue', lwd = 2) + # add scale bar - first value is TOP, second value is BOTTOM
  tm_scale_bar(width = 0.22, position = c("RIGHT", "BOTTOM")) + 
  # add compass
  tm_compass(position = c("RIGHT", "TOP")) + 
  tm_layout(title = "CRD Census Tracts Median Income\nRook Neighbours 2016", title.position = c("LEFT", "TOP"), inner.margins = c(.08, .03, .08, .03))

map_crd.net

## Warning: Current projection of shape crd.net unknown and cannot be
## determined.

map_crd.net2

## Warning: Current projection of shape crd.net unknown and cannot be
## determined.

# create png of Median Income CRD net
png("map_map_crd.net_MdInc_2.png")  
map_crd.net

## Warning: Current projection of shape crd.net unknown and cannot be
## determined.

dev.off()

## quartz_off_screen 
##                 2

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

## Warning: Current projection of shape crd.net unknown and cannot be
## determined.

dev.off()

## quartz_off_screen 
##                 2

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) +
  # add scale bar - first value is TOP, second value is BOTTOM
  tm_scale_bar(width = 0.22, position = c("RIGHT", "BOTTOM")) + 
  # add compass
  tm_compass(position = c("RIGHT", "TOP")) + 
  tm_layout(title = "CRD Census Tracts Median Income\nLagged Means 2016", title.position = c("LEFT", "TOP"), inner.margins = c(.08, .03, .08, .03))

map_LagMean

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

## quartz_off_screen 
##                 2

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)  + 
  # add scale bar - first value is TOP, second value is BOTTOM
  tm_scale_bar(width = 0.22, position = c("RIGHT", "BOTTOM")) + 
  # add compass
  tm_compass(position = c("RIGHT", "TOP")) + 
  tm_layout(title = "CRD Census Tracts Median Income\nLISA Local Moran's i 2016", title.position = c("LEFT", "TOP"), inner.margins = c(.08, .03, .08, .03))

map_LISA

# create png of Map LISA MdInc
png("map_LISA_MdInc_2.png")  
map_LISA
dev.off()

## quartz_off_screen 
##                 2

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="Median Income", 
           ylab="Spatially Lagged Mean Median Income", quiet=NULL, main ="Local Moran's i 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="Median Income", 
           ylab="Spatially Lagged Median Income", quiet=NULL, main ="Local Moran's i Plot for CRD Median Income 2016")
dev.off()

## quartz_off_screen 
##                 2

---

Variable II: Renter

IV.ii. Clean Data

1. Remove the na values in the new column: Renter

## 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$Renter),]
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  
##                                                                  
##   MDAftTaxInc       GovTrans          EInc            MInc      
##  Min.   :12192   Min.   :  720   Min.   :11616   Min.   :11200  
##  1st Qu.:30491   1st Qu.: 3406   1st Qu.:29105   1st Qu.:31412  
##  Median :33420   Median : 5824   Median :31669   Median :34421  
##  Mean   :33389   Mean   : 5454   Mean   :31549   Mean   :34637  
##  3rd Qu.:37373   3rd Qu.: 7254   3rd Qu.:34292   3rd Qu.:38460  
##  Max.   :48939   Max.   :12079   Max.   :52395   Max.   :53632  
##                                                                   
##    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  
##                                                                  
##        Ii                E.Ii             Var.Ii             Z.Ii         
##  Min.   :-3.72431   Min.   :-0.0137   Min.   :0.02923   Min.   :-3.82072  
##  1st Qu.:-0.01818   1st Qu.:-0.0137   1st Qu.:0.12419   1st Qu.:-0.01247  
##  Median : 0.10445   Median :-0.0137   Median :0.16286   Median : 0.28519  
##  Mean   : 0.14317   Mean   :-0.0137   Mean   :0.20140   Mean   : 0.45976  
##  3rd Qu.: 0.32975   3rd Qu.:-0.0137   3rd Qu.:0.22656   3rd Qu.: 0.98397  
##  Max.   : 2.10785   Max.   :-0.0137   Max.   :0.94320   Max.   : 4.45715  
##                                                                           
##        P            
##  Min.   :0.0000042  
##  1st Qu.:0.1626611  
##  Median :0.3878925  
##  Mean   :0.3611524  
##  3rd Qu.:0.5049726  
##  Max.   :0.9999335  
## 

2. 1. In order to compare the Renter variable between different polygons, the variable needs to be ‘normalized’, by dividing each count of Renters in a polygon by the population of that polygon, and creating a new column with the % of population that is a renter in any polygon … these values can now be compared across different polygons.2) A new column is created $RenterNorm, and this new variable will now be used for the analysis.

## 
crd.data$RenterNorm <- crd.data$Renter / crd.data$Pop
class(crd.data$RenterNorm)

## [1] "numeric"

# new column

V.ii. Make Chloropleth Map: Renter

Using the tmap package …

# make Chloropleth map
map_PopDens <- tm_shape(crd.data) + 
  tm_polygons(col = "RenterNorm", 
              title = "Renter (% of Population)", 
              style = "jenks", 
              palette = "viridis", n = 6) + 
  # add scale bar - first value is TOP, second value is BOTTOM
  tm_scale_bar(width = 0.22, position = c("RIGHT", "BOTTOM")) + 
  # add compass
  tm_compass(position = c("RIGHT", "TOP")) + 
  tm_layout(title = "CRD Census Tracts Renter % of Population 2016", title.position = c("LEFT", "TOP"), inner.margins = c(.08, .03, .08, .03))

map_PopDens

# create png of RenterNormalized % of Population
png("map_PopDens_RenterNorm_2.png")  
map_PopDens
dev.off()

## quartz_off_screen 
##                 2

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') + 
  # add scale bar - first value is TOP, second value is BOTTOM
  tm_scale_bar(width = 0.22, position = c("RIGHT", "BOTTOM")) + 
  # add compass
  tm_compass(position = c("RIGHT", "TOP")) + 
  tm_layout(title = "CRD Census Tracts Renter % of Population\nQueen Neighbours 2016", title.position = c("LEFT", "TOP"), inner.margins = c(.08, .03, .08, .03))

##

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='red', lwd = 2) +
  tm_shape(crd.net2) + tm_lines(col='blue', lwd = 2) + 
  # add scale bar - first value is TOP, second value is BOTTOM
  tm_scale_bar(width = 0.22, position = c("RIGHT", "BOTTOM")) + 
  # add compass
  tm_compass(position = c("RIGHT", "TOP")) + 
  tm_layout(title = "CRD Census Tracts Renter % of Population\nRook Neighbours 2016", title.position = c("LEFT", "TOP"), inner.margins = c(.08, .03, .08, .03))

map_crd.net

## Warning: Current projection of shape crd.net unknown and cannot be
## determined.

map_crd.net2

## Warning: Current projection of shape crd.net unknown and cannot be
## determined.

# create png of Renter CRD net
png("map_map_crd.net_RenterNorm_2.png")  
map_crd.net

## Warning: Current projection of shape crd.net unknown and cannot be
## determined.

dev.off()

## quartz_off_screen 
##                 2

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

## Warning: Current projection of shape crd.net unknown and cannot be
## determined.

## Warning: Current projection of shape crd.net2 unknown and cannot be
## determined.

dev.off()

## quartz_off_screen 
##                 2

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.listw2. Creates a map using tmap.

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

map_LagMean <- tm_shape(crd.data) + 
  tm_polygons(col = "RenterNormLagMeans", 
              title = "Renter % of Population\nLagged Means", 
              style = "jenks", 
              palette = "viridis", n = 6) + 
  # add scale bar - first value is TOP, second value is BOTTOM
  tm_scale_bar(width = 0.22, position = c("RIGHT", "BOTTOM")) + 
  # add compass
  tm_compass(position = c("RIGHT", "TOP")) + 
  tm_layout(title = "CRD Census Tracts Renter % of Population\nLagged Means 2016", title.position = c("LEFT", "TOP"), inner.margins = c(.08, .03, .08, .03))

map_LagMean

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

## quartz_off_screen 
##                 2

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$RenterNorm, crd.lw, zero.policy = TRUE)
mi

## 
##  Moran I test under randomisation
## 
## data:  crd.data$RenterNorm  
## weights: crd.lw    
## 
## Moran I statistic standard deviate = 8.6868, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.595415952      -0.013698630       0.004916714

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] 8.686834

# 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.595416 -0.01369863 0.004916714 8.686834

table3 <- data.for.table3 

# write csv file with Morans Global results for I, E, Var, Z, P
write.csv(table3, "MoransGlobal-RenterNorm.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$RenterNorm, 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.233821626   -0.01369863      0.15045712    0.63812297 2.616968e-01
## 2   0.032943097   -0.01369863      0.23218304    0.09679644 4.614440e-01
## 3   0.549070589   -0.01369863      0.15045712    1.45085485 7.341014e-02
## 4   0.503692530   -0.01369863      0.23218304    1.07375144 1.414671e-01
## 5   0.060472308   -0.01369863      0.18314749    0.17331386 4.312024e-01
## 6   1.556366078   -0.01369863      0.15045712    4.04772671 2.585874e-05
## 7   3.391433333   -0.01369863      0.18314749    7.95670883 8.833791e-16
## 8   1.553256264   -0.01369863      0.18314749    3.66147450 1.253839e-04
## 9   1.140736503   -0.01369863      0.18314749    2.69754721 3.492619e-03
## 10 -0.071622535   -0.01369863      0.12710685   -0.16247006 5.645321e-01
## 11  1.333661488   -0.01369863      0.10959416    4.06995858 2.351075e-05
## 12  0.360821928   -0.01369863      0.09597317    1.20892885 1.133451e-01
## 13  0.205149532   -0.01369863      0.47736081    0.31675212 3.757159e-01
## 14  0.316206463   -0.01369863      0.31390896    0.58882625 2.779889e-01
## 15  3.643380638   -0.01369863      0.09597317   11.80482233 1.842843e-32
## 16  1.379862981   -0.01369863      0.10959416    4.20951901 1.279575e-05
## 17  1.839610501   -0.01369863      0.12710685    5.19832443 1.005465e-07
## 18  0.359894929   -0.01369863      0.31390896    0.66680297 2.524490e-01
## 19 -0.045775803   -0.01369863      0.18314749   -0.07495414 5.298744e-01
## 20 -0.018460853   -0.01369863      0.12710685   -0.01335750 5.053287e-01
## 21  0.640176636   -0.01369863      0.23218304    1.35699942 8.739070e-02
## 22  1.897336263   -0.01369863      0.18314749    4.46547986 3.994474e-06
## 23  0.098397443   -0.01369863      0.08507638    0.38431373 3.503730e-01
## 24 -0.026062633   -0.01369863      0.23218304   -0.02565924 5.102354e-01
## 25  0.147756323   -0.01369863      0.12710685    0.45286305 3.253237e-01
## 26  0.431022202   -0.01369863      0.10959416    1.34336421 8.957703e-02
## 27  0.356600204   -0.01369863      0.15045712    0.95465395 1.698764e-01
## 28  0.304032189   -0.01369863      0.08507638    1.08931840 1.380067e-01
## 29  0.577861864   -0.01369863      0.18314749    1.38228846 8.344157e-02
## 30 -0.168322406   -0.01369863      0.18314749   -0.36130652 6.410648e-01
## 31 -0.001621923   -0.01369863      0.12710685    0.03387381 4.864889e-01
## 32  4.474925356   -0.01369863      0.18314749   10.48848459 4.879043e-26
## 33  0.064613726   -0.01369863      0.15045712    0.20189424 4.199997e-01
## 34 -0.904832954   -0.01369863      0.31390896   -1.59052798 9.441421e-01
## 35  0.011019807   -0.01369863      0.18314749    0.05775911 4.769703e-01
## 36  3.129918431   -0.01369863      0.18314749    7.34563180 1.023947e-13
## 37  0.136728114   -0.01369863      0.18314749    0.35149939 3.626069e-01
## 38  0.199580888   -0.01369863      0.12710685    0.59822515 2.748449e-01
## 39  0.097877099   -0.01369863      0.18314749    0.26071694 3.971554e-01
## 40  0.529096109   -0.01369863      0.23218304    1.12647195 1.299829e-01
## 41  0.236191920   -0.01369863      0.12710685    0.70091499 2.416780e-01
## 42  0.226356428   -0.01369863      0.12710685    0.67332754 2.503695e-01
## 43  0.604421772   -0.01369863      0.15045712    1.59355372 5.551801e-02
## 44  0.442216620   -0.01369863      0.31390896    0.81373362 2.078988e-01
## 45  0.147383189   -0.01369863      0.15045712    0.41527918 3.389688e-01
## 46  0.713133316   -0.01369863      0.18314749    1.69837476 4.471853e-02
## 47  0.387150645   -0.01369863      0.23218304    0.83188991 2.027355e-01
## 48  0.612030700   -0.01369863      0.15045712    1.61317003 5.335379e-02
## 49  0.156866528   -0.01369863      0.10959416    0.51522464 3.031980e-01
## 50  0.234052890   -0.01369863      0.02964488    1.43893589 7.508434e-02
## 51  0.662880165   -0.01369863      0.31390896    1.20758171 1.136041e-01
## 52  0.630290771   -0.01369863      0.15045712    1.66024565 4.843252e-02
## 53  0.223824443   -0.01369863      0.12710685    0.66622560 2.526335e-01
## 54  1.838131703   -0.01369863      0.15045712    4.77413642 9.024001e-07
## 55  0.305859518   -0.01369863      0.10959416    0.96528642 1.672007e-01
## 56  0.387356307   -0.01369863      0.23218304    0.83231672 2.026151e-01
## 57  0.099487976   -0.01369863      0.15045712    0.29180227 3.852189e-01
## 58  0.046748357   -0.01369863      0.23218304    0.12544675 4.500849e-01
## 59  0.998762311   -0.01369863      0.96771634    1.02921058 1.516904e-01
## 60  0.043229883   -0.01369863      0.09597317    0.18376167 4.271002e-01
## 61 -0.252114651   -0.01369863      0.23218304   -0.49478917 6.896255e-01
## 62  0.355471073   -0.01369863      0.47736081    0.53432153 2.965595e-01
## 63  0.645739004   -0.01369863      0.47736081    0.95444377 1.699295e-01
## 64  0.290854679   -0.01369863      0.96771634    0.30959168 3.784357e-01
## 65  0.299887489   -0.01369863      0.15045712    0.80844497 2.094172e-01
## 66  0.601525828   -0.01369863      0.15045712    1.58608780 5.635970e-02
## 67  0.811880946   -0.01369863      0.23218304    1.71334055 4.332495e-02
## 68 -0.037825907   -0.01369863      0.15045712   -0.06220166 5.247989e-01
## 69  0.151580204   -0.01369863      0.12710685    0.46358861 3.214713e-01
## 70  0.262426181   -0.01369863      0.12710685    0.77449915 2.193178e-01
## 71  0.584495330   -0.01369863      0.18314749    1.39778875 8.108826e-02
## 72  0.157706771   -0.01369863      0.31390896    0.30593041 3.798288e-01
## 73  0.104039684   -0.01369863      0.09597317    0.38005189 3.519534e-01
## 74  0.768114071   -0.01369863      0.31390896    1.39540689 8.144659e-02

head(table4)

##   crd.data.Ii crd.data.E.Ii crd.data.Var.Ii crd.data.Z.Ii   crd.data.P
## 1  0.23382163   -0.01369863       0.1504571    0.63812297 2.616968e-01
## 2  0.03294310   -0.01369863       0.2321830    0.09679644 4.614440e-01
## 3  0.54907059   -0.01369863       0.1504571    1.45085485 7.341014e-02
## 4  0.50369253   -0.01369863       0.2321830    1.07375144 1.414671e-01
## 5  0.06047231   -0.01369863       0.1831475    0.17331386 4.312024e-01
## 6  1.55636608   -0.01369863       0.1504571    4.04772671 2.585874e-05

# write csv file with Lisa Test results for I, E, Var, Z, P
write.csv(table4, "LisaTest-RenterNorm.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)  + 
  # add scale bar - first value is TOP, second value is BOTTOM
  tm_scale_bar(width = 0.22, position = c("RIGHT", "BOTTOM")) + 
  # add compass
  tm_compass(position = c("RIGHT", "TOP")) + 
  tm_layout(title = "CRD Census Tracts Renter % of Population\nLISA Local Moran's i 2016", title.position = c("LEFT", "TOP"), inner.margins = c(.08, .03, .08, .03))

map_LISA

# create png of Map LISA Renter Normalized
png("map_LISA_RenterNorm_2.png")  
map_LISA
dev.off()

## quartz_off_screen 
##                 2

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$RenterNorm, crd.lw, zero.policy=NULL, spChk=NULL, labels=NULL, xlab="Renter % of Population", 
           ylab="Spatially Lagged Renter % of Population", quiet=NULL, main ="Local Moran's i Plot for CRD Renter % of Population 2016")

# create png of Moran scatter plot
png("moran.plot_RenterNorm_2.png")  
moran.plot(crd.data$RenterNorm, crd.lw, zero.policy=NULL, spChk=NULL, labels=NULL, xlab="Renter % of Population", 
           ylab="Spatially Lagged Renter % of Population", quiet=NULL, main ="Local Moran's i Plot for CRD Renter % of Population 2016")
dev.off()

## quartz_off_screen 
##                 2

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.