I’ve been working on a project on state-level farmland values and appreciation rates, and colleagues have expressed interest in knowing how I achieved some of the elementary spatial data exploration for it, particularly the eye-catching animated choropleth maps that I built. You can see the code and results below!
library(expm)
library(extrafont)
library(httr)
library(jsonlite)
library(lubridate)
library(magick)
library(rgdal)
library(spdep)
library(spdplyr)
library(tidyverse)
loadfonts(quiet=TRUE)
As I’m analyzing real farmland values and appreciation rates, I’ll need data on the general price level in the US over time. I choose the GDP deflator as it is a good broad indicator of the price level independent of particular consumer goods baskets.
fred.url <- paste(
'https://api.stlouisfed.org/fred/series/observations',
'?series_id=GDPDEF',
'&api_key=', api.token.fred,
'&file_type=json',
sep=''
)
fred.return <- GET(fred.url)
fred.data <-
fromJSON(content(fred.return, 'text'))$observations %>%
type_convert() %>%
mutate(year = year(date)) %>%
group_by(year) %>%
summarize(gdpdef = mean(value))
write_csv(fred.data, '../data/fred-data.csv')
Data on farmland value is available directly from the USDA National Agricultural Statistics Service via their API:
# NASS Data
nass.url <- paste(
'http://quickstats.nass.usda.gov/api/api_GET/?key=', api.token.nass,
'&source_desc=SURVEY',
'&reference_period_desc=YEAR',
'&agg_level_desc=STATE',
'&commodity_desc=AG LAND',
'&statisticcat=ASSET VALUE',
'&short_desc=AG LAND, INCL BUILDINGS - ASSET VALUE, MEASURED IN $ / ACRE',
'&domain_desc=TOTAL',
sep=''
)
nass.return <- GET(nass.url)
nass.data <-
fromJSON(content(nass.return, 'text'))$data %>%
as_tibble() %>%
type_convert() %>%
mutate(year = parse_integer(year)) %>%
rename(landvalac = Value)
write_csv(nass.data, '../data/nass-data.csv')
A shapefile of the United States is publicly available from the US Census Bureau:
download.file(
'http://www2.census.gov/geo/tiger/GENZ2016/shp/cb_2016_us_state_20m.zip',
destfile = '../data/cb_2016_us_state_20m.zip')
unzip('../data/cb_2016_us_state_20m.zip', exdir = '../data/')
With the data downloaded from their respective sources, I can prepare it for my map generating procedures and spatial regressions. First, I can access the GDP deflator for 2010 (our base year for real values). Then generating real values and appreciation rates is fairly straightfoward.
For the shapefile, there are some districts that should be removed if they are not included in the NASS data. Then for the spatial regressions later on, I make sure that the states are alphabetized in order to ensure the spatial weights matrices I generate are in accordance with the order of the NASS data.
# Set a range of years
year.range <- 1960:2016
# Read in FRED data from local copy
fred.data <- read_csv('../data/fred-data.csv')
gdpdef2010 <-
fred.data %>%
filter(year == 2010) %>%
select(gdpdef) %>%
as_vector()
# Function to calculate appreciation rates
pct <- function(x) 100 * (x / lag(x) - 1)
# Read in NASS data from local copy, adjust land values to real (2010) values,
# and add appreciation rates
nass.data <- read_csv('../data/nass-data.csv') %>%
select(state_alpha, year, landvalac) %>%
arrange(state_alpha, year) %>%
inner_join(fred.data, by = 'year') %>%
mutate(realvalac = landvalac * gdpdef / gdpdef2010) %>%
group_by(state_alpha) %>%
mutate(realvalac.PD1 = pct(realvalac))
# Read in US state shapefile
states <-
readOGR(dsn='../data', layer='cb_2016_us_state_20m',
stringsAsFactors = FALSE) %>%
filter(!(STUSPS %in% c('DC', 'AK', 'HI', 'PR'))) %>%
arrange(STUSPS)
## OGR data source with driver: ESRI Shapefile
## Source: "../data", layer: "cb_2016_us_state_20m"
## with 52 features
## It has 9 fields
## Integer64 fields read as strings: ALAND AWATER
states.tbl <- fortify(states, region='STUSPS')
A problem with an animated map is that I’m making 57 separate maps (one for each year) and stitching them together. Because of that, I need to explicitly tell ggplot2
that I want the legend to use the data available for all years in the animation, so that the colors have constant meaning year by year. There’s no built-in color scale in ggplot2
that accounts for the legend breaks that I want to hold constant across year, so I will use the scale_fill_manual
function. This option, however, requires us to explicitly state my desired breakpoints and the colors I will use between them, so I define both below:
# Create bins and color palette for legend
# Use K-State Purple and an accent tan from branding guide
quantile(nass.data$realvalac, seq(0, 1, 0.1))
## 0% 10% 20% 30% 40%
## 1.527728 11.524702 22.362028 45.977944 118.194399
## 50% 60% 70% 80% 90%
## 300.622549 552.699499 932.764775 1698.832913 3326.439259
## 100%
## 16473.475800
realvalac.levels <- c(0, 25, 50, 125, 275, 450,
750, 1150, 2000, 3500, 17000)
realvalac.labels <- c('0 - 25', '25 - 50', '50 - 125', '125 - 275',
'275 - 450', '450 - 750', '750 - 1150', '1150 - 2000',
'2000 - 3500', '3500 - 17000')
colfunc <- colorRampPalette(c('#ffffff', '#cfa35c','#512888'))
realvalac.colors <- colfunc(10)
names(realvalac.colors) <- realvalac.labels
nass.data$realvalac.fc <- factor(
cut(nass.data$realvalac, realvalac.levels),
labels = realvalac.labels
)
Now that the legend breaks and colors have been created, I can generate my year-by-year choropleth maps using ggplot2
:
# Generate a choropleth map for each year
maplist <- lapply(year.range, function(yr) {
ggmap <- ggplot(nass.data %>% filter(year == yr)) +
geom_map(aes(map_id=state_alpha, fill=realvalac.fc),
map=states.tbl, color='#cccccc', size=0.20) +
coord_map() +
expand_limits(x = states.tbl$long, y = states.tbl$lat) +
labs(title=paste('Real Farmland Value (', yr, ')', sep=''),
subtitle='Including Buildings, Constant 2010 US Dollars/Acre',
caption='Data: USDA National Agricultural Statistics Service') +
scale_fill_manual(element_blank(),
values = realvalac.colors, drop=FALSE) +
theme_void() +
theme(text=element_text(family='Myriad Pro'),
plot.title=element_text(face='bold'),
plot.subtitle=element_text(size=8),
plot.caption=element_text(size=6),
legend.margin=margin(t=-0.6, r=0, b=-0, l=0, unit='cm'))
ggsave(paste('../png/landvalues-levels-01-choropleth-', yr, '.png', sep=''),
width = 6, height = 3.1, units = 'in', dpi=150, pointsize=14)
return(ggmap)
})
With a choropleth map generated and stored for each year in my sample, I can use ImageMagick to stitch them together into an animated GIF.
# Read the individual maps into a data structure for use with 'magick'
imglayers <- sapply(year.range, function(yr) {
image_read(paste('../png/landvalues-levels-01-choropleth-', yr, '.png', sep=''))
})
# Generate an animated GIF with the individual maps and write to a file
imganim <- image_animate(image_join(imglayers), fps = 1, dispose = "previous")
image_write(imganim, '../gif/landvalues-levels-01-choropleth-anim.gif')
Notice that there is a clear upward trend in farmland values for each state over my sample, and the legend breaks indicate that this growth is exponential in nature. As stationary variables are generally preferred for statistical analysis, one logical step toward a more informative map is to plot farmland appreciation rates instead of value levels. I follow a very similar procedure as followed for the farmland levels, only here I account for the imbalance of negative and positive appreciation rates by separately creating breakpoints and color scales, then stitching them together.
# Create bins and color palette for legend
# Use K-State Purple and an accent goldenrod from branding guide
quantile(nass.data$realvalac.PD1[nass.data$realvalac.PD1 < 0],
seq(0, 1, 0.5), na.rm = TRUE)
## 0% 50% 100%
## -34.562556217 -3.490169972 -0.001892014
quantile(nass.data$realvalac.PD1[nass.data$realvalac.PD1 > 0],
seq(0, 1, 0.125), na.rm = TRUE)
## 0% 12.5% 25% 37.5% 50%
## 0.02342932 3.17082868 5.47912534 7.17475612 9.10742226
## 62.5% 75% 87.5% 100%
## 11.34159156 14.57069859 20.79343575 161.78167924
realvalac.PD1.levels <- c(-26, -4, 0, 4, 6, 8,
10, 12, 16, 23, 163)
realvalac.PD1.labels <- c('-26 to -4', '-4 to 0', '0 to 4', '4 to 6',
'6 to 8', '8 to 10', '10 to 12', '12 to 16',
'16 to 23', '23 to 163')
colfunc <- colorRampPalette(c('#f0ad00', '#ffffff'))
realvalac.PD1.colors <- colfunc(3)
colfunc <- colorRampPalette(c('#ffffff', '#512888'))
realvalac.PD1.colors <- c(realvalac.PD1.colors, colfunc(8)[2:8])
names(realvalac.PD1.colors) <- realvalac.PD1.labels
nass.data$realvalac.PD1.fc <- factor(
cut(nass.data$realvalac.PD1, realvalac.PD1.levels),
labels = realvalac.PD1.labels
)
With the legends constructed, the maps are created in exactly the same manner as for the levels above:
# Generate a choropleth map for each year
maplist <- lapply(year.range, function(yr) {
ggmap <- ggplot(nass.data %>% filter(year == yr)) +
geom_map(aes(map_id=state_alpha, fill=realvalac.PD1.fc),
map=states.tbl, color='#cccccc', size=0.20) +
coord_map() +
expand_limits(x = states.tbl$long, y = states.tbl$lat) +
labs(title=paste('Real Farmland Appreciation in', yr),
subtitle='Including Buildings, Annual Percentage Change',
caption='Data: USDA National Agricultural Statistics Service') +
scale_fill_manual(element_blank(),
values = realvalac.PD1.colors, drop=FALSE) +
theme_void() +
theme(text=element_text(family='Myriad Pro'),
plot.title=element_text(face='bold'),
plot.subtitle=element_text(size=8),
plot.caption=element_text(size=6),
legend.margin=margin(t=-0.6, r=0, b=-0, l=0, unit='cm'))
ggsave(paste('../png/landvalues-pctchange-01-choropleth-', yr, '.png', sep=''),
width = 6, height = 3.1, units = 'in', dpi=150, pointsize=14)
return(ggmap)
})
The animated maps are likewise created in exactly the same manner as for the levels:
# Read the individual maps into a data structure for use with 'magick'
imglayers <- sapply(year.range, function(yr) {
image_read(paste('../png/landvalues-pctchange-01-choropleth-', yr, '.png', sep=''))
})
# Generate an animated GIF with the individual maps and write to a file
imganim <- image_animate(image_join(imglayers), fps = 1, dispose = "previous")
image_write(imganim, '../gif/landvalues-pctchange-01-choropleth-anim.gif')
Note how in the animated choropleth maps for both the levels and appreciation rates, there appears to be a high degree of correlation between states, most notably in the 1970s when most states appear to be appreciating very rapidly and the 1980s when most states appear to be depreciating rapidly. As the eye alone is inadequate for determining spatial correlation, however, a statistical measure would be useful. Here I calculate Moran’s I Statistic for each year using Queen contiguity matrices of order 1-3, where Order 1 means I am looking for the global correlation between states and their immediate neighbors, Order 2 means I am looking for correlation with neighbors of neighbors, and so on.
I begin by constructing a dataframe of the Global Moran’s I statistics for each year for both farmland value levels and appreciation rates:
# Generate Spatial Weights Matrices
states.q1.nb <- poly2nb(states, row.names=states$STUSPS)
states.q1.W <- nb2mat(states.q1.nb)
states.queen.W <- lapply(1:3, function(i) states.q1.W %^% i)
states.queen.listW <- lapply(states.queen.W, function(X) mat2listw(X))
# Yearly Moran's I tests on levels and
# appreciation rates of farmland value
moran.stats <- lapply(1:3, function(i) {
W <- states.queen.listW[[i]]
get.morans <- function(series.in, series.name) {
morans.out <- sapply(year.range, function(yr) {
data.yr <- nass.data %>% filter(year == yr)
series.yr <- data.yr[, series.in] %>% as_vector()
moran.yr <- moran.test(series.yr, W)
moran.out <- c(yr, series.name, moran.yr$estimate)
names(moran.out) <- c('year', 'series', 'moran', 'expected', 'variance')
return(moran.out)
}) %>% t() %>% as_tibble() %>% type_convert()
}
moran.out <-
bind_rows(get.morans('realvalac', 'Val/Acre'),
get.morans('realvalac.PD1', 'Appr Rate')) %>%
mutate(series = factor(series, levels = c('Val/Acre', 'Appr Rate')),
moran.lower = moran - 1.96 * sqrt(variance),
moran.upper = moran + 1.96 * sqrt(variance),
queen.order = i)
return(moran.out)
}) %>%
bind_rows() %>%
mutate(queen.order = factor(queen.order,
levels = 1:3,
labels = c('First-Order', 'Second-Order',
'Third-Order')))
With the Moran’s I statistics in hand, I can plot them as time series. Here I plot them using two methods: The first keeps all the time series plots separate but includes confidence bands to demonstrate that the spatial correlation detected is in fact statistically significant. The second puts the Moran’s I plots for each contiguity definition on the same set of axes for easier comparison.
# Moran's I Plots (with confidence bands)
moran.plot.noband <- ggplot(moran.stats) +
geom_ribbon(aes(x=year, ymin=moran.lower, ymax=moran.upper), alpha=0.25) +
geom_line(aes(x=year, y=moran)) +
geom_hline(yintercept = -1/47, linetype='dashed') +
facet_grid(series ~ queen.order, scales='free_y') +
labs(title = "Moran's I Statistics: 1960 - 2016",
subtitle = paste('Farmland Values and Appreciation Rates\n',
'Queen Contiguity Spatial Weights', sep=''),
caption = 'Data: USDA National Agriculutural Statistics Service') +
theme(text = element_text(family='Myriad Pro'),
plot.title = element_text(face='bold'),
plot.subtitle = element_text(size=8),
plot.caption = element_text(size=6),
axis.title.x = element_blank(),
axis.title.y = element_blank())
ggsave(paste('../png/moranstats-bands.png', sep=''),
width = 6, height = 3.1, units = 'in', dpi=150, pointsize=14)
# Moran's I Plots (no confidence bands)
moran.plot.bands <- ggplot(moran.stats) +
geom_line(aes(x=year, y=moran, color=queen.order)) +
geom_hline(yintercept = -1/47, linetype='dashed') +
facet_grid(series ~ ., scales='free_y') +
labs(title = "Moran's I Statistics: 1960 - 2016",
subtitle = paste('Farmland Values and Appreciation Rates\n',
'Queen Contiguity Spatial Weights', sep=''),
caption = 'Data: USDA National Agriculutural Statistics Service') +
theme(text = element_text(family='Myriad Pro'),
plot.title = element_text(face='bold'),
plot.subtitle=element_text(size=8),
plot.caption=element_text(size=6),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank(),
legend.position = 'bottom')
ggsave(paste('../png/moranstats-nobands.png', sep=''),
width = 6, height = 3.1, units = 'in', dpi=150, pointsize=14)
One fact made clear by the Global Moran plots is that there is in fact a significant statistical correlation between the farmland values and appreciation rates seen in one state and those of their neighbors. But what drives this global correlation? That is to ask, where is the correlation strongest?
Here I will focus on appreciation rates. As there is almost always a significant positive global correlation, I would expect to see clusters of slowly-appreciating (or depreciating) states and clusters of rapidly appreciating states. Anselin’s idea of a LISA, or Local Indicator of Spatial Association, allows us to find and map these clusters.
This task requires the calculation of the Local Moran’s I for each state in each year and note its significance. This involves breaking the Local Moran statistics with p-values less than 0.05 into four groups:
The last two categories may seem counterintuitively named. The reason a below-average appreciation rate with a negative Local Moran is called “Low-High” is that the negative Local Moran indicates that neighbors’ values are the opposite of low, i.e. they are high. A similar argument explains the “High-Low” category.
The code below computes the Local Moran’s I for each state, for each year, and for each of the three orders of queen contiguity and stores them as static images using a similar procedure as for the choropleth maps above.
# Yearly LISA: Local Indicators of Spatial Association
lisa.tbl <- lapply(1:3, function(i) {
local.moran.tbl <- lapply(year.range, function(yr) {
W <- states.queen.listW[[i]]
states.yr <- rownames(listw2mat(W))
data.yr <- nass.data %>% filter(year == yr)
series.yr <- data.yr$realvalac.PD1
dev.yr <- (series.yr - mean(series.yr)) / sd(series.yr)
local.moran.yr <-
localmoran(series.yr, W) %>%
as_tibble()
local.moran.yr$queen <- i
local.moran.yr$year <- rep(yr, length(states.yr))
local.moran.yr$state <- states.yr
local.moran.yr$dev <- dev.yr
local.moran.yr <- select(local.moran.yr, queen, year, state,
dev, Ii, `Pr(z > 0)`)
colnames(local.moran.yr) <- c('queen', 'year', 'state_alpha',
'dev', 'local.moran', 'pvalue')
return(local.moran.yr)
}) %>% bind_rows()
}) %>% bind_rows()
lisa.tbl$queen <- factor(lisa.tbl$queen, levels = c(1, 2, 3),
labels = c('First-Order', 'Second-Order', 'Third-Order'))
lisa.tbl$lisa.type <- apply(lisa.tbl, 1, function(x) {
dev <- as.numeric(x[4])
local.moran <- as.numeric(x[5])
pvalue <- as.numeric(x[6])
if (dev > 0 && local.moran > 0 && pvalue < 0.05) return('High-High')
else if (dev > 0 && local.moran < 0 && pvalue < 0.05) return('High-Low')
else if (dev < 0 && local.moran < 0 && pvalue < 0.05) return('Low-High')
else if (dev < 0 && local.moran > 0 && pvalue < 0.05) return('Low-Low')
else return('Unclustered')
}) %>% factor(., levels=c('High-High', 'High-Low', 'Unclustered', 'Low-High', 'Low-Low'))
# Generate a choropleth map for each year
colfunc <- colorRampPalette(c('#aa0000', '#ffffff','#0000aa'))
lisa.colors <- colfunc(5)
maplist <- lapply(year.range, function(yr) {
ggmap <- ggplot(lisa.tbl %>% filter(year == yr)) +
geom_map(aes(map_id=state_alpha, fill=lisa.type),
map=states.tbl, color='#cccccc', size=0.20) +
coord_map() +
facet_grid(queen ~ .) +
expand_limits(x = states.tbl$long, y = states.tbl$lat) +
labs(title=paste('LISA Clusters: Farmland Appreciation', yr),
subtitle=paste('Including Buildings, Year/Year Appreciation Rate\n',
'Queen-Contiguity Spatial Weights Matrices (First-Third Order)',
sep=''),
caption='Data: USDA National Agricultural Statistics Service') +
scale_fill_manual(values = lisa.colors, drop=FALSE) +
theme_void() +
theme(text=element_text(family='Myriad Pro'),
plot.title=element_text(face='bold'),
plot.subtitle=element_text(size=8),
plot.caption=element_text(size=6),
legend.title=element_blank(),
legend.position='bottom')
ggsave(paste('../png/landvalues-pctchange-02-lisamaps-', yr, '.png', sep=''),
width = 6, height = 9.3, units = 'in', dpi=150, pointsize=14)
return(ggmap)
})
# Read the individual maps into a data structure for use with 'magick'
imglayers <- sapply(year.range, function(yr) {
image_read(paste('../png/landvalues-pctchange-02-lisamaps-', yr, '.png', sep=''))
})
# Generate an animated GIF with the individual maps and write to a file
imganim <- image_animate(image_join(imglayers), fps = 1, dispose = "previous")
image_write(imganim, '../gif/landvalues-pctchange-02-lisamaps-anim.gif')
The existence of spatial correlation in farmland appreciation, significant globally across the country and in at least one locale nearly every year, indicates that the appreciation rate in each state has a component that depends on its neighbors in addition to its own endogenous component. I can use a simple spatial autoregressive model (SAR) with no exogenous regressors to separate the observed appreciation rates into these components: \[ y = y' W y + \varepsilon \]
The residuals \(\varepsilon\) represent the estimated component of (the deviation from the mean of) farmland appreciation endogenous to the state’s own properties, separated from the “contaminating” effect of neighboring state appreciation rates. I can use the lagsarlm
function in the spdep
package to perform the relevant regression (using a vector of zeroes as my only exogenous regressor \(X\)): \[
y = y' W y + X \beta + \varepsilon \quad , \quad X = 0
\]
sar.tbl <- lapply(year.range, function(yr) {
W <- states.queen.listW[[1]]
sar.yr <- lagsarlm(realvalac.PD1 ~ 1,
data = nass.data %>% filter(year == yr),
listw = W)
tbl.out <- tibble(
year = yr,
state_alpha = states$STUSPS,
sar.res = residuals(sar.yr)
)
}) %>% bind_rows()
# Create bins and color palette for legend
# Use K-State Purple and an accent goldenrod from branding guide
quantile(sar.tbl$sar.res[sar.tbl$sar.res < 0], seq(0, 1, 0.2))
## 0% 20% 40% 60% 80%
## -17.988460381 -5.160907919 -3.070779398 -1.964696624 -0.982230384
## 100%
## -0.003456114
quantile(sar.tbl$sar.res[sar.tbl$sar.res > 0], seq(0, 1, 0.2))
## 0% 20% 40% 60% 80%
## 3.372423e-03 7.829704e-01 1.868054e+00 3.176893e+00 5.496989e+00
## 100%
## 1.500732e+02
sar.res.levels <- c(-18, -5, -3, -2, -1, 0,
1, 2, 3, 5, 163)
sar.res.labels <- c('-18 to -5', '-5 to -3', '-3 to -2', '-2 to -1', '-1 to 0',
'0 to 1', '1 to 2', '2 to 3', '3 to 5', '5 to 163')
colfunc <- colorRampPalette(c('#f0ad00', '#ffffff', '#512888'))
sar.res.colors <- colfunc(10)
names(sar.res.colors) <- sar.res.labels
sar.tbl$sar.res.fc <- factor(
cut(sar.tbl$sar.res, sar.res.levels),
labels = sar.res.labels
)
maplist <- lapply(year.range, function(yr) {
ggmap <- ggplot(sar.tbl %>% filter(year == yr)) +
geom_map(aes(map_id=state_alpha, fill=sar.res.fc),
map=states.tbl, color='#cccccc', size=0.20) +
coord_map() +
expand_limits(x = states.tbl$long, y = states.tbl$lat) +
labs(title=paste('Farmland Appreciation: SAR Residuals (', yr, ')', sep=''),
subtitle=paste('Including Buildings, Year/Year Appreciation Rate\n',
'First-Order Queen-Contiguity Spatial Weights Matrix',
sep=''),
caption='Data: USDA National Agricultural Statistics Service') +
scale_fill_manual(element_blank(),
values = sar.res.colors, drop=FALSE) +
theme_void() +
theme(text=element_text(family='Myriad Pro'),
plot.title=element_text(face='bold'),
plot.subtitle=element_text(size=8),
plot.caption=element_text(size=6),
legend.title=element_blank(),
legend.position='right',
legend.margin=margin(t=-0.6, r=0, b=-0, l=0, unit='cm'))
ggsave(paste('../png/landvalues-pctchange-03-sarresids-', yr, '.png', sep=''),
width = 6, height = 3.1, units = 'in', dpi=150, pointsize=14)
return(ggmap)
})
# Read the individual maps into a data structure for use with 'magick'
imglayers <- sapply(year.range, function(yr) {
image_read(paste('../png/landvalues-pctchange-03-sarresids-', yr, '.png', sep=''))
})
# Generate an animated GIF with the individual maps and write to a file
imganim <- image_animate(image_join(imglayers), fps = 1, dispose = "previous")
image_write(imganim, '../gif/landvalues-pctchange-03-sarresids-anim.gif')
Given the goal was to extract the component of farmland appreciation rates independent of neighbors’ rates, it is worthwhile to check the residuals themselves for spatial correlation. I do this in exactly the same manner as before (using only first-order Queen contiguity):
# Yearly Moran's I tests on levels and
# appreciation rates of farmland value
W <- states.queen.listW[[1]]
get.morans <- function(data.in, series.in) {
morans.out <- sapply(year.range, function(yr) {
data.yr <- data.in %>% filter(year == yr)
series.yr <- data.yr[, series.in] %>% as_vector()
moran.yr <- moran.test(series.yr, W)
moran.out <- c(yr, moran.yr$estimate)
names(moran.out) <- c('year', 'moran', 'expected', 'variance')
return(moran.out)
}) %>% t() %>% as_tibble() %>% type_convert()
}
sar.res.morans <-
get.morans(sar.tbl, 'sar.res') %>%
mutate(moran.lower = moran - 1.96 * sqrt(variance),
moran.upper = moran + 1.96 * sqrt(variance))
# Moran's I Plots (with confidence bands)
moran.plot.noband <- ggplot(sar.res.morans) +
geom_ribbon(aes(x=year, ymin=moran.lower, ymax=moran.upper), alpha=0.25) +
geom_line(aes(x=year, y=moran)) +
geom_hline(yintercept = -1/47, linetype='dashed') +
labs(title = "Moran's I Statistics: 1960 - 2016",
subtitle = paste('SAR Residuals: Farmland Appreciation Rates\n',
'Queen Contiguity Spatial Weights', sep=''),
caption = 'Data: USDA National Agriculutural Statistics Service') +
theme(text = element_text(family='Myriad Pro'),
plot.title = element_text(face='bold'),
plot.subtitle=element_text(size=8),
plot.caption=element_text(size=6),
axis.title.x = element_blank(),
axis.title.y = element_blank())
ggsave(paste('../png/landvalues-pctchange-04-sarresids-moranstats.png', sep=''),
width = 6, height = 3.1, units = 'in', dpi=150, pointsize=14)
From the plot, I find the Global Moran’s I statistics to be small and insignificant. I therefore conclude that I have successfully extracted the neighbor-independent farmland appreciation rates. This conclusion, however, may be dependent upon my choice of spatial weights matrix! I intend to explore this idea further as my research progresses.