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!

Setup

library(expm)
library(extrafont)
library(httr)
library(jsonlite)
library(lubridate)
library(magick)
library(rgdal)
library(spdep)
library(spdplyr)
library(tidyverse)

loadfonts(quiet=TRUE)

Data

Data from St. Louis FRED

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 from USDA NASS

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')

Shapefile from US Census Bureau

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/')

Manipulation

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')

Choropleth Maps

Farmland Values

Legend Preparation

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
)

Map Generation

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)
  })

Animated Map Compilation

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')

Appreciation Rates

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.

Legend Preparation

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

Map Generation

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)
})

Animated Map Compilation

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')

Spatial Correlation

Global: Moran’s I

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)

Local: LISA Clusters

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:

  • High-High: Above yearly average appreciation and positive Local Moran
  • High-Low: Above yearly average appreciation and negative Local Moran
  • Low-High: Below yearly average appreciation and negative Local Moran
  • Low-Low: Below yearly average appreciation and positive Local Moran

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')

Spatial Autoregression

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')

Moran’s I

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.