# Moran's I examples # Michael Minn # 15 December 2020 library(sf) library(spdep) # Create sample polygon grid polygons = list() for (x in seq(-5, 4)) { for (y in seq(-5, 4)) { polygons = c(polygons, list(matrix(c(x, y, x, y+1, x+1, y+1, x+1, y, x, y), ncol=2, byrow=T))) } # for y } # for x polygons = lapply(polygons, function(x) { st_polygon(list(x)) }) polygons = st_sfc(polygons, crs=4326) polygons = st_sf(polygons) polygons = st_transform(polygons, 3857) # https://mgimond.github.io/Spatial/spatial-autocorrelation-in-r.html png("2020-morans-i-theoretical.png", width=700, height=500, pointsize=16) layout(matrix(c(1,1,3,3,5,5,0,2,2,4,4,0), nrow=2, byrow=T)) par(mar=c(1,1,1,1)) polygons$value = rep(c(rep(c(0,1), 5), rep(c(1,0), 5)), 5) neighbors = poly2nb(as_Spatial(polygons), queen=F) weights = nb2listw(neighbors, style="W", zero.policy=T) m = moran.test(polygons$value, weights) title = paste0("Perfect Dispersion (I = ", round(m$estimate[1], 2), ")") plot(polygons, col=ifelse(polygons$value == 0, "#f0f0f0", "#a10022"), border="#a0a0a0", main=title, reset=F) polygons$value = rep(c(0,0,1,1,1,0,0), 15)[1:100] neighbors = poly2nb(as_Spatial(polygons), queen=T) weights = nb2listw(neighbors, style="W", zero.policy=T) m = moran.test(polygons$value, weights) title = paste0("Dispersed (I = ", round(m$estimate[1], 2), ")") plot(polygons, col=ifelse(polygons$value == 0, "#f0f0f0", "#a10022"), border="#a0a0a0", main=title, reset=F) set.seed(100) polygons$value = round(runif(100)) neighbors = poly2nb(as_Spatial(polygons), queen=T) weights = nb2listw(neighbors, style="W", zero.policy=T) m = moran.test(polygons$value, weights) title = paste0("Spatial Independence (I = ", round(m$estimate[1], 2), ")") plot(polygons, col=ifelse(polygons$value == 0, "#f0f0f0", "#a10022"), border="#a0a0a0", main=title, reset=F) polygons$value = c(rep(c(1,1,1,1,0,0,0,0,0,0), 3), rep(c(0,0,0,0,0,1,1,1,1,0), 3), rep(c(1,1,1,1,0,0,0,0,0,0), 4)) neighbors = poly2nb(as_Spatial(polygons), queen=T) weights = nb2listw(neighbors, style="W", zero.policy=T) m = moran.test(polygons$value, weights) title = paste0("Clustering (I = ", round(m$estimate[1], 2), ")") plot(polygons, col=ifelse(polygons$value == 0, "#f0f0f0", "#a10022"), border="#a0a0a0", main=title, reset=F) polygons$value = c(rep(1, 50), rep(0,50)) neighbors = poly2nb(as_Spatial(polygons), queen=F) weights = nb2listw(neighbors, style="W", zero.policy=T) m = moran.test(polygons$value, weights) title = paste0("Perfect Bifurcation (I = ", round(m$estimate[1], 2), ")") plot(polygons, col=ifelse(polygons$value == 0, "#f0f0f0", "#a10022"), border="#a0a0a0", main=title, reset=F) dev.off() counties = st_read("2014-2018-acs-counties.geojson", stringsAsFactors = F) counties = counties[(counties$FIPS >= "17000") & (counties$FIPS <= "17999"),] neighbors = poly2nb(as_Spatial(counties), queen=T) weights = nb2listw(neighbors, style="W", zero.policy=T) #m = moran.test(counties$Median.Household.Income, weights) #m png("2020-morans-i-mortgage.png", width=400, height=600, pointsize=16) plot(counties[,c("Median.Monthly.Mortgage", "geometry")], pal=colorRampPalette(c("#004d3c", "#f0f0f0", "#a10022"))(5), nbreaks=5, breaks="jenks", main="", key.pos=1) m = moran.test(counties$Median.Monthly.Mortgage, weights) print(m) dev.off() png("2020-morans-i-work-home.png", width=400, height=600, pointsize=16) plot(counties[,c("X..Work.at.Home", "geometry")], pal=colorRampPalette(c("#004d3c", "#f0f0f0", "#a10022"))(5), nbreaks=5, breaks="jenks", main="", key.pos=1) m = moran.test(counties$X..Work.at.Home, weights) print(m) dev.off() png("2020-morans-i-income.png", width=400, height=600, pointsize=16) plot(counties[,c("Median.Household.Income", "geometry")], pal=colorRampPalette(c("#004d3c", "#f0f0f0", "#a10022"))(5), nbreaks=5, breaks="jenks", main="", key.pos=1) m = moran.test(counties$Median.Household.Income, weights) print(m) dev.off() stop("Done") # https://stats.idre.ucla.edu/r/faq/how-can-i-calculate-morans-i-in-r/ points = list() for (x in seq(-5, 4)) { for (y in seq(-5, 4)) { polygons = c(polygons, list(matrix(c(x, y, x, y+1, x+1, y+1, x+1, y, x, y), ncol=2, byrow=T))) } # for y } # for x polygons = lapply(polygons, function(x) { st_polygon(list(x)) }) polygons = st_sfc(polygons, crs=4326) polygons = st_sf(polygons) polygons = st_transform(polygons, 3857) library(ape) #distances = st_distance(polygons, by_element=F) #distances = 1 / distances #diag(distances) = 0 #distances = drop_units(distances) #Moran.I(polygons$value, distances)