library(sf)
## Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1
library(spdep)
## Loading required package: sp
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(tmap)
setwd("/Users/asiyavalidova/Dropbox/Demography/GIS/R&D/Cluster analysis")

Load COI data and inpute missing

setwd("/Users/asiyavalidova/Dropbox/Demography/GIS/R&D/Cluster analysis")
sa15<-get(load("sa15.Rdata"))
sa15$county=48029
title="COI"
legend_title="San Antonio Child Opportunity Index 2015"
tm_shape(sa15)+
          tm_polygons(c("r_COI_stt"), title=c(title), palette="Blues")+
          tm_layout(title=legend_title, title.size =1.5, title.color="black",legend.frame = TRUE, title.position = c('center', 'bottom'))+
          tm_scale_bar()+
          tm_compass()

mn<-mean(sa15$r_COI_stt,na.rm=TRUE)  
sa15$r_COI_stt<-ifelse(is.na(sa15$r_COI_stt),mn,sa15$r_COI_stt)
#Make a rook style weight matrix
sa15nb<-poly2nb(as(sa15, "Spatial"), queen=F)  

sa15lw<-nb2listw(sa15nb, style="W")
sa15lw$weights[1]
## [[1]]
## [1] 0.2 0.2 0.2 0.2 0.2
sa15lag <- lag.listw(sa15lw, sa15$r_COI_stt)
fit<-lm(sa15lag~sa15$r_COI_stt)

#the Moran's I coefficient
coef(fit)[2]
## sa15$r_COI_stt 
##      0.7447655
#png("Moranplot_COI.png")
#Moran Scatter plot
moran.plot(sa15$r_COI_stt, listw=sa15lw, xlab="Child opportunity index", ylab="Laggged child opportunity index", main=c("Moran Scatterplot for Child Opportunity Index", "in San Antonio 2015. I= 0.744"))

#dev.off()
#The Moran's I coefficient is 0.744. The positive (upward) slope suggests that as the percent Hispanic value of a polygon increases, so does those of its neighboring polygons.

Local Moran Cluster Map

sa15.nbs<-poly2nb(sa15, queen = T)
sa15.wts<-nb2listw(neighbours = sa15.nbs, style = "W")

locali<-localmoran(sa15$r_COI_stt,
                   listw = sa15.wts )

sa15$locali<-locali[,1]
sa15$localp<-locali[,5] 
sa15$index<-scale(sa15$r_COI_stt)
sa15$lag_index<-lag.listw(var=sa15$index, x = sa15.wts)
sa15$quad_sig <- NA
sa15$quad_sig[(sa15$index >= 0 & sa15$lag_index >= 0) & (sa15$localp <= 0.05)] <- 1 #high high
sa15$quad_sig[(sa15$index <= 0 & sa15$lag_index <= 0) & (sa15$localp <= 0.05)] <- 2 #low low
sa15$quad_sig[(sa15$index >= 0 & sa15$lag_index <= 0) & (sa15$localp <= 0.05)] <- 3 #high low
sa15$quad_sig[(sa15$index <= 0 & sa15$lag_index >= 0) & (sa15$localp <= 0.05)] <- 4 #low high

breaks <- seq(1, 5, 1)

# Set the corresponding labels for the thematic map classes
labels <- c("High-High", "Low-Low", "High-Low", "Low-High", "Not Clustered")
np <- findInterval(sa15$quad_sig, breaks)

# Assign colors to each map class
colors <- c("red", "blue", "lightpink", "skyblue2", "white")

#png("Moran_Cluster_Map_COI.png")
plot(sa15["county"],col = colors[np] , main="Moran LISA Cluster Map of San Antonio - Child Opportunity Index",sub=" San Antonio")
legend("topright", legend = labels, fill = colors, bty = "n")

#dev.off()

Load ACS data, impute missing

2015

sa15<-get(load(file="sa_acs15.Rdata"))
title="Percent Hispanic"
legend_title="San Antonio Percent Hispanic"
tm_shape(sa15)+
          tm_polygons(c("phisp"), title=c(title), palette="Blues", style="quantile", n=5)+
          tm_layout(title=legend_title, title.size =1.5, title.color="black",legend.frame = TRUE, title.position = c('center', 'bottom'))+
          tm_scale_bar()+
          tm_compass() 

mn<-mean(sa15$phisp,na.rm=TRUE)  
sa15$phisp<-ifelse(is.na(sa15$phisp),mn,sa15$phisp)

Make a rook style weight matrix

sa15nb<-poly2nb(as(sa15, "Spatial"), queen=F)       
sa15lw<-nb2listw(sa15nb, style="W")
sa15lw$weights[1]
## [[1]]
## [1] 0.2 0.2 0.2 0.2 0.2
sa15lag <- lag.listw(sa15lw, sa15$phisp)
fit<-lm(sa15lag~sa15$phisp)

#the Moran's I coefficient
coef(fit)[2]
## sa15$phisp 
##  0.7838428
#png("Moranplot_phisp15.png")
#Moran Scatter plot
moran.plot(sa15$phisp, listw=sa15lw, xlab="Percent Hispanic", ylab="Lagged percent Hispanic", main=c("Moran Scatterplot for percent Hispanic", "in San Antonio 2015. I= 0.784"))

#dev.off()
#The Moran's I coefficient is 0.784. The positive (upward) slope suggests that as the percent Hispanic value of a polygon increases, so does those of its neighboring polygons.

Local Moran Cluster Map

sa15.nbs<-poly2nb(sa15, queen = T)
sa15.wts<-nb2listw(neighbours = sa15.nbs, style = "W")

locali<-localmoran(sa15$phisp,
                   listw = sa15.wts )

sa15$locali<-locali[,1]
sa15$localp<-locali[,5] 
sa15$hisp<-scale(sa15$phisp)
sa15$lag_hisp<-lag.listw(var=sa15$hisp, x = sa15.wts)
sa15$quad_sig <- NA
sa15$quad_sig[(sa15$hisp >= 0 & sa15$lag_hisp >= 0) & (sa15$localp <= 0.05)] <- 1 #high high
sa15$quad_sig[(sa15$hisp <= 0 & sa15$lag_hisp <= 0) & (sa15$localp <= 0.05)] <- 2 #low low
sa15$quad_sig[(sa15$hisp >= 0 & sa15$lag_hisp <= 0) & (sa15$localp <= 0.05)] <- 3 #high low
sa15$quad_sig[(sa15$hisp <= 0 & sa15$lag_hisp >= 0) & (sa15$localp <= 0.05)] <- 4 #low high

breaks <- seq(1, 5, 1)

# Set the corresponding labels for the thematic map classes
labels <- c("High-High", "Low-Low", "High-Low", "Low-High", "Not Clustered")
np <- findInterval(sa15$quad_sig, breaks)

# Assign colors to each map class
colors <- c("red", "blue", "lightpink", "skyblue2", "white")

#png("Moran_Cluster_Map_phisp15.png")
plot(sa15["county"],col = colors[np] , main="Moran LISA Cluster Map of San Antonio - Percent Hispanic",sub=" San Antonio")
legend("topright", legend = labels, fill = colors, bty = "n")

#dev.off()