Setup
#devtools::install_github("amvallone/estdaR")
#install.packages("ggrepel")
suppressWarnings(suppressMessages({
library(estdaR)
library(sf)
library(spdep)
library(rgeos)
library(tidyverse)
library(ggrepel)
library(ExPanDaR)
library(plotly)
}))
options(repr.plot.width = 6, repr.plot.height = 4)
Import data
gdppc <- st_read("https://raw.githubusercontent.com/quarcs-lab/data-open/master/mexico/mexico.geojson")
glimpse(gdppc)
y1940 <- gdppc$PCGDP1940/mean(gdppc$PCGDP1940)
y1950 <- gdppc$PCGDP1950/mean(gdppc$PCGDP1950)
y1960 <- gdppc$PCGDP1960/mean(gdppc$PCGDP1960)
y1970 <- gdppc$PCGDP1970/mean(gdppc$PCGDP1970)
y1980 <- gdppc$PCGDP1980/mean(gdppc$PCGDP1980)
y1990 <- gdppc$PCGDP1990/mean(gdppc$PCGDP1990)
regionName <- gdppc$NAME
Weights matrix
Wqueen <- nb2listw(poly2nb(gdppc))
Spatial lags
Wy1940 <- lag.listw(Wqueen, y1940)
Wy1950 <- lag.listw(Wqueen, y1950)
Wy1960 <- lag.listw(Wqueen, y1960)
Wy1970 <- lag.listw(Wqueen, y1970)
Wy1980 <- lag.listw(Wqueen, y1980)
Wy1990 <- lag.listw(Wqueen, y1990)
RegimeMaxP <- gdppc$MAXP
Moran scatter plot
dy <- y1950 - y1940
dW <- Wy1950 - Wy1940
NofRows <- nrow(gdppc)
y0 <- rep(0, NofRows)
W0 <- rep(0, NofRows)
df <- as.data.frame(cbind(y1940, y1950, dy, Wy1940, Wy1950, dW, y0, W0))
df <- cbind(regionName, df)
head(df)
ggplot(df, aes(x = y1940, y = Wy1940)) +
geom_point(color = "dodgerblue4") +
geom_text_repel(aes(x = y1940, y = Wy1940, label = regionName), size = 2.5) +
geom_hline(yintercept=mean(Wy1940), lty=2, color = "dodgerblue4") +
geom_vline(xintercept=mean(y1940), lty=2, color = "dodgerblue4") +
labs(title="Moran scatter plot in 1940", y="Relative income of neighboring regions", x= "Relative income of focal region")
ggplot(df, aes(x = y1950, y = Wy1950)) +
geom_point(color = "dodgerblue4") +
geom_text_repel(aes(x = y1950, y = Wy1950, label = regionName), size = 2.5) +
geom_hline(yintercept=mean(Wy1950), lty=2, color = "dodgerblue4") +
geom_vline(xintercept=mean(y1950), lty=2, color = "dodgerblue4") +
labs(title="Moran scatter plot in 1950", y="Relative income of neighboring regions", x= "Relative income of focal region")
Directional Moran
ggplot(df) +
geom_segment(aes(x = y1940, y = Wy1940, xend = y1950, yend = Wy1950, color = y1940), arrow = arrow(length = unit(0.1, "cm"))) +
scale_color_gradient(low='#2c7bb6', high='#d7191c', name='Initial \nIncome') +
geom_text_repel(aes(x = y1950, y = Wy1950, label = regionName), size = 2.5) +
geom_hline(yintercept=1, lty=2, color = "dodgerblue4") +
geom_vline(xintercept=1, lty=2, color = "dodgerblue4") +
labs(title="Directional Moran scatter plot 1940-1950", y="Relative income of neighboring regions", x= "Relative income of focal region")
Standardized directional Moran
ggplot(df) +
geom_segment(aes(x = y0, y = W0, xend = dy, yend = dW, color = y1940), arrow = arrow(length = unit(0.1, "cm"))) +
scale_color_gradient(low='#2c7bb6', high='#d7191c', name='Initial \nIncome') +
geom_text_repel(aes(x = dy, y = dW, label = regionName), size = 2.5) +
geom_hline(yintercept=0, lty=2, size = 0.3) +
geom_vline(xintercept=0, lty=2, size = 0.3) +
geom_abline(intercept = 0, slope = 1, lty=2, size = 0.3) +
geom_abline(intercept = 0, slope = -1, lty=2, size = 0.3) +
labs(title="Standardized directional Moran scatter plot 1940-1950", y="Change in rel. income of neighboring regions", x= "Change in rel. income of focal region")
Using the estdaR package
directionalLISA1940_1950 <- d.LISA(y1940, y1950, Wqueen, k=4, nsim=999)
directionalLISA1940_1950$lisamap
directionalLISA1940_1950$rose
directionalLISA1940_1950$counts
directionalLISA1940_1950$p.value