Global and Local Measures of Spatial Autocorrelation

This post aims at being a summary of the available techniques to investigate spatial autocorrelation for the social sciences, rather than presenting the theory behind spatial autocorrelation. For that there are great books available on line, like Anselin’s, Le Page, and Bivand just to cite a few.

The techniques presented here work for a spatial polygons data frame. The difference between spatial points and polygons data frames is not that big, the idea is the same and most of what I am doing here can be applied to data points.

Why do we look at spatial autocorrelation at all? Spatial autocorrelation leads to biased results in regressions, this is the reason why we want to compute Moran’s I and why we include spatial autocorrelation if its measurement proves to be significative and non-random.

Spatial autocorrelation can be investigated globally or locally. “Globally”, implies that the measure you’re going to obtain refers to the dataset as a whole, whether it is a whole country, continent or region. “Locally”, means that you are taking into consideration each and every polygon and getting a measure for each one of them.

We start by uploading the data, projecting them (important when considering the distance based measures -earth is not flat, whatever they may say!), and construct neighborhood relations (in this case Queen and Rook, but could be any other). For more detail see this post on how to construct neighbor relations.

library(maptools)
library(spdep)
NC= readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
nb.FOQ = poly2nb(NC, queen=TRUE, row.names=NC$FIPSNO)
nb.RK = poly2nb(NC, queen=FALSE, row.names=NC$FIPSNO)

Global measures of spatial autocorrelation

There are two main measures of global spatial autocorrelation: Moran’s I and Geary’s C. Moran’s I is the most used in my experience, but both work are perfectly acceptable.
Moran’s I ranges between -1 (strong negative spatial autocorrelation with a dispersed pattern) and 1 (strong positive spatial autocorrelation with a clustered pattern) with 0 being the absence of spatial autocorrelation.
Geary’s C ranges between 0 and 2, with positive spatial autocorrelation ranging from 0 to 1 and negative spatial autocorrelation between 1 and 2.
Of course being these inferential measures, if the p-value is non significant we cannot exclude that the patterns could be random(!)

In this case Moran’s I is positive and significant, the z-score (not provided by moran.test) is positive implying spatial clusters, so we can reject the null hypothesis.

library(spdep)
nwb <- NC$NWBIR74
moran.test(nwb, listw = nb2listw(nb.RK))
##
##  Moran I test under randomisation
##
## data:  nwb
## weights: nb2listw(nb.RK)
##
## Moran I statistic standard deviate = 3.0787, p-value = 0.001039
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance
##       0.185965551      -0.010101010       0.004055701
geary.test(nwb, listw = nb2listw(nb.RK))
##
##  Geary C test under randomisation
##
## data:  nwb
## weights: nb2listw(nb.RK)
##
## Geary C statistic standard deviate = 2.0324, p-value = 0.02106
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance
##        0.83274888        1.00000000        0.00677185

if you have polygons with no neighbors remember to specify zero.policy=NULL

moran.plot(nwb, listw = nb2listw(nb.RK))

moran plot

moran.mc(nwb, listw = nb2listw(nb.RK), nsim=100)
##
##  Monte-Carlo simulation of Moran I
##
## data:  nwb
## weights: nb2listw(nb.RK)
## number of simulations + 1: 101
##
## statistic = 0.18597, observed rank = 100, p-value = 0.009901
## alternative hypothesis: greater

permutation Moran_s I

plot(moran.mc(nwb, listw = nb2listw(nb.RK), nsim=100))


Same thing can be done for Geary’s C:

geary.mc(nwb, listw = nb2listw(nb.RK), nsim=100)
##
##  Monte-Carlo simulation of Geary C
##
## data:  nwb
## weights: nb2listw(nb.RK)
## number of simulations + 1: 101
##
## statistic = 0.83275, observed rank = 5, p-value = 0.0495
## alternative hypothesis: greater
plot(geary.mc(nwb, listw = nb2listw(nb.RK), nsim=100))

Local measures of spatial autocorrelation

Local Moran and Local G

locm <- localmoran(nwb, listw = nb2listw(nb.RK))
locG <- localG(nwb, listw = nb2listw(nb.RK))

Get the neighbor matrix into a listwise format with listw: there’s two options here, row-standardized weights matrix style = "W" creates proportional weights when polygons have an unequal number of neighbors, balancing out observations with few neighbors. Binary weights style = "B" upweight observations with many neighbors.

library(classInt)
library(dplyr)
myvar <- NC$NWBIR74
nb <- nb.RK
# Define weight style
ws <- c("W")

# Define significance for the maps
significance <- 0.05
plot.only.significant <- TRUE

# Transform the neigh mtx into a listwise object
listw <- nb2listw(nb, style=ws)

# Create the lagged variable
lagvar <- lag.listw(listw, myvar)

# get the mean of each
m.myvar <- mean(myvar)
m.lagvar <- mean(lagvar)

The next step is to derive the quadrants and set the coloring scheme. I like to color the border of each polygon with the color of their local moran score, regardless of their pvalue, and then fill only the significant ones.

n <- length(NC)
#
vec <- c(1:n)
vec <- ifelse(locm[,5] < significance, 1,0)

# Derive quadrants
q <- c(1:n) for (i in 1:n) {   if (myvar[[i]]>=m.myvar & lagvar[[i]]>=m.lagvar) q[i] <- 1
  if (myvar[[i]]<m.myvar & lagvar[[i]]<m.lagvar) q[i] <- 2
  if (myvar[[i]]<m.myvar & lagvar[[i]]>=m.lagvar) q[i] <- 3   if (myvar[[i]]>=m.myvar & lagvar[[i]]<m.lagvar) q[i] <- 4
}

# set coloring scheme
q.all <- q
colors <- c(1:n)
for (i in 1:n) {
  if (q.all[i]==1) colors[i] <- "red"
  if (q.all[i]==2) colors[i] <- "blue"
  if (q.all[i]==3) colors[i] <- "lightblue"
  if (q.all[i]==4) colors[i] <- "pink"
  if (q.all[i]==0) colors[i] <- "white"   if (q.all[i]>4) colors[i] <- "white"
}

# Mark all non-significant regions white
locm.dt <- q*vec
colors1 <- colors
for (i in 1:n)
{
  if ( !(is.na (locm.dt[i])) )  {
  if (locm.dt[i]==0) colors1[i] <- "white"
}
}
colors2 <- colors
colors2 <- paste(colors2,vec)
pos = list()
for (i in 1:n) {
  pos[[i]] <- c(which(NC$NWBIR74==colors2["blue 0"]))
}

blue0 <- which(colors2=="blue 0")
red0 <- which(colors2=="red 0")
lightblue0 <- which(colors2=="lightblue 0")
pink0 <- which(colors2=="pink 0")
lb <- 6
labels=c("High-High", "High-Low", "Low-High", "Low-Low")

# plot the map
if (plot.only.significant==TRUE) plot(NC, col=colors1,border=F) else
  plot(NC, col=colors,border=F)
plot(NC[blue0,],border="blue",lwd=0.5,add=T)
plot(NC[lightblue0,],border="lightblue",add=T,lwd=0.5)
plot(NC[red0,],border="red",add=T,lwd=0.5)
plot(NC[pink0,],border="pink",add=T,lwd=0.5)
legend("bottomleft", legend = labels, fill = c("red", "pink", "lightblue", "blue"), bty = "n")

lisa map
Local G gives back z-scores values and indicate the posibility of a local cluster of high values of the variable being analysed, very low values indicate a similar cluster of low values.

library(RColorBrewer)

nclassint <- 3
colpal <- brewer.pal(nclassint,"PiYG")
cat <- classIntervals(locG, nclassint, style = "jenks", na.ignore=T)
color.z <- findColours(cat, colpal)

plot(NC, col= color.z, border=T)

Rplot

# color only significant polygons
plot(NC, border=T)
plot(NC[which(vec==1),], col=color.z[which(vec==1)], border=T, add=T)

lisa map2

Advertisements

A nice example of hafen/geofacet from Washington Post to ggplot2

I’ve recently came across the hafen/geofacet function and was pondering to blog an example. Then, I came across a perfect example, thanks to  kanishkamisra for working on the dataset & code and making it available via github here!

 

usa_vs_state1

High Resolution Mapping of Fertility and Mortality from Household Survey Data in Low Income Settings – PAA presentation

I will present at PAA my WorldPop mapping of Demographic indicators in low-income settings at PAA in Chicago.  “Advances in Mathematical, Spatial, and Small-Area Demography”, Thursday, April 27, 2017: 10:15 AM – 11:45 AM, Hilton, Joliet Room.

Find color breaks for mapping (fast)

I’ve stumbled upon a little trick to compute jenks breaks faster than with the classInt package, just be sure to use n+1 instead of n as the breaks are computed a little bit differently. That is to say, if you want 5 breaks, n=6, no biggie there.

For more on the Bayesian Analysis of Macroevolutionary Mixtures see BAMMtools library

install.packages("BAMMtools")
library(BAMMtools)
system.time(getJenksBreaks(mydata$myvar, 6))
> user system elapsed
> 0.970 0.001 0.971

On the other hand this takes way more time with large datasets
library(classInt)
system.time(classIntervals(mydata$myvar, n=5, style="jenks"))
> Timing stopped at: 1081.894 1.345 1083.511

A map of the US election results

  1. Upload libraries:
rm(list = ls(all=T)) #clear workspace
library(dplyr)
library(readr)
library(stringr)
library(tidyr)
library(readxl)
library(classInt)
library(RColorBrewer)
library(maptools) #to read shapefiles

2. Download the data files (note they are not ready for use but need some cleaning as there are more areas in the excel files than polygons in the shape file). I copy here the code as I have used it in my script but it’s available at RPubs thanks to David Robinson.

download.file("http://www2.census.gov/prod2/statcomp/usac/excel/LND01.xls", "LND01.xls")
download.file("http://www2.census.gov/prod2/statcomp/usac/excel/POP01.xls", "POP01.xls")

according to metadata, this is Land Area in 2010 and resident population in 2010:

us_county_area <- read_excel("LND01.xls") 
transmute(CountyCode = as.character(as.integer(STCOU)), Area = LND110210D)

us_county_population <- read_excel("POP01.xls") 
transmute(CountyCode = as.character(as.integer(STCOU)),Population = POP010210D)

3. Adjust data

election_url <- "https://raw.githubusercontent.com/Prooffreader/election_2016_data/master/data/presidential_general_election_2016_by_county.csv"
county_data <- read_csv(election_url) 
group_by(CountyCode = as.character(fips)) 
mutate(TotalVotes = sum(votes)) 
ungroup() 
mutate(name = str_replace(name, ".\\. ", "")) 
filter(name %in% c("Trump", "Clinton", "Johnson", "Stein")) 
transmute(County = str_replace(geo_name, " County", ""),
State = state,
CountyCode = as.character(fips),
Candidate = name,
Percent = vote_pct / 100,
TotalVotes) 
spread(Candidate, Percent, fill = 0) 
inner_join(us_county_population, by = "CountyCode") 
inner_join(us_county_area, by = "CountyCode")

you can save the data into a csv file:

# write_csv(county_data, "county_election_2016.csv")

You can download the cleaned datafile here: data_election_2016_by_county

4. Upload data and shape files

setwd("/Users/...")
dt <- read.csv("new_county_election_2016.csv", header=T)
us <- readShapePoly("./USA_adm/USA_adm2.shp")
us0 <- readShapePoly("./USA_adm/USA_adm0.shp")
us.m <- us[-c(which(us$NAME_1=="Alaska")),] #get rid of Alaska
us.d <- us.m[-c(67:71),]

5. Prepare the color palette(s)

nclassint <- 5 #number of colors to be used in the palette
cat.T <- classIntervals(dt$Trump[-c(67:71)], nclassint,style = "jenks") #style refers to how the breaks are created
colpal.T <- brewer.pal(nclassint,"Reds")
color.T <- findColours(cat.T,colpal.T) #sequential
bins.T <- cat.T$brks
lb.T <- length(bins.T)

5. Plot the maps with map basic

# pdf("Where are the trump voters.pdf")
# plot(us.d, col=color.T, border=F)
# plot(us0,add=T, lwd=0.1)
# legend("bottomleft",fill=colpal.T,legend=paste(round(bins[-length(bins.T)],1),":",round(bins.T[-1],1)),cex=1, bg="white")
# dev.off()
clinton-voters
% Votes for Clinton
where-are-the-trump-voters
% Votes for Trump

… or ggplot2

library(ggplot2)
library(scales)
theme_set(theme_bw())

ggplot(county_data, aes(Population / Area, Trump)) +
  geom_point() +
  geom_point(data=county_data[which(county_data$State=="Texas"),], aes(x=Population/Area, y=Trump), colour="red")+
  scale_x_log10() +
  scale_y_continuous(labels = percent_format()) +
  xlab("Population density (ppl / square mile)") +
  ylab("% of votes going to Trump") +
  geom_text(aes(label = County), vjust = 1, hjust = 1, check_overlap = TRUE) +
  geom_smooth(method = "lm") +
  ggtitle("Population density vs Trump voters by county (Texas Counties in red)")

This is the code to plot in red points according to State (in red) and to add red labels to those points. The check_overlap=T avoids overlapping labels.

# ggplot(county_data, aes(Population / Area, Trump)) +
#   geom_point() +
#   geom_point(data=county_data[which(county_data$State=="California"),], aes(x=Population/Area, y=Trump), colour="red")+
#   scale_x_log10() +
#   scale_y_continuous(labels = percent_format()) +
#   xlab("Population density (ppl / square mile)") +
#   ylab("% of votes going to Trump") +
#   geom_text(data=county_data[which(county_data$State=="California"),], aes(label = ifelse(Trump&amp;gt;.5, as.character(dt$County), "" )), color= "red",size=5,vjust = 1, hjust = 1, check_overlap = TRUE) +
#   geom_smooth(method = "lm") +
#   ggtitle("Population density vs Trump voters by county (California in red)")

rplot1geom_point_texas

californiaclinton

How to get good maps in R and avoid the expensive softwares

How to convey as much information as possible in a clear and simple way? Producing maps for social sciences is not difficult, there are a plethora of softwares that can help us. But there are a few issues to consider when choosing your to go program:
(1) Do I want to do all my analysis in one (or more) program(s) and then switch to another one to make those maps?
(2) Are those programs freely accessible to me?
1. Using more than one softwares usually implies spending time to learn different syntax:  why do your analysis in (insert name here ____) and then plot in R when you can do everything in R?
2. The availability of mapping softwares is no trivial issue. Not all researchers have powerful computers, not all institutes have bottomless funds to buy licences, and sometimes having the possibility to map on your laptop while bingeing on Netflix is way nicer than waiting for the one computer with the one licence.

Probably the best and most elegant mapping tool available to Geographers is ArcGIS (to my knowledge, but again, I use R and own a Mac), however it does not come for free. What to do? Well, R is a very good alternative, you can produce elegant maps, customizable to the very last detail. The only drawback I have encountered is the time you would spend to get the first map, but then you would have the syntax and any other map would be pretty quick to plot, and you can always for loop all graphics (although I do not recommend it). Moreover, R runs on your Mac (and Linux), it allows for way more control over features, and has great color palettes (see here and here).

Here are some useful libraries:
library(maps) #for creating geographical maps
library(RColorBrewer) #contains color palettes
library(classInt) #defines the class intervals for the color palettes
library(maptools) #tools for handling spatial objects
library(raster) #tools to deal with raster maps
library(ggplot2) #to create maps, quick and painless

Some stuff to keep in mind:
(1) add a scale with scale.map (or a nice  scalebar);
(2) it is sometimes required to add a north arrow, you can find many versions for that (see this document on page 4 for  examples, I use the same with no labels);
(3) locator() is a very useful tool to get the coordinates when adding labels, arrows, scales and so on.

Part 1: get a plain map.

Below is a very simple example produced using EUROSTAT shape files for world countries (world) and DIVA-GIS for Spain at NUTS3 level (spain). In this map I have removed the Canary Islands, but you can always cut it and paste it in the map using either par(fig=c(…)) or par(fin(…)), inset, or something more elaborated with layout, and framing it using box() or rectangle.

world is the shapefile for the whole world, where I select the neighboring countries I want to appear in the map, in this case Spain, France, Portugal, France, Morocco, and Algeria.
spain is the Spain NUTS3 shapefile where I remove the Canary Islands (45)


plot(spain[-c(45),], border=F) #this first line does not plot anything, it just centers my graph on Spain, the -c(45) removes the Canary Islands
plot(world[c(6, 67,74, 132, 177),], border="lightblue",add=T, col="beige") #plotting the countries appearing in the map
plot(spain[-c(45),], border="brown", lwd=0.2, add=T, col="lightblue") #plot spain, removing the Canary Islands
map.scale(3,35.81324, ratio=F, cex=0.7, relwidth=0.1) # scale map
map.axes(cex.axis=0.9)
northarrow(c(4.8,42.9),0.7, cex=0.8)

Spain

Part 2: Add labels

Using the function shadow text to avoid labels overlapping.

coords<- coordinates(spain) # get goordinates of the centroids, it's where you center your labels
# p.names is a data frame containing the coordinates and all the names of the provinces (remember to get rid of those you don't want to use if using only a selection). Usually you can find the names in the shapefile, but I didn't have them.
shadowtext(p.names[,1],p.names[,2], label=paste(p.names[,3]), cex=0.7,col="black", bg="white",r=0.1)

Spain3

 

Ubicación, ubicación, ubicación! ¿Por qué asuntos espaciales en la demografía y por qué debemos cuidar.

Me he dado cuenta solo ahora que mi post en Demotrends sobre la dimension espacial de los fenomenos demograficos ha sido traducido en español por el grupo “Población y Desarrollo en Honduras”, muchas gracias! Aquí esta:

Los fenómenos demográficos son inherentemente espaciales, así como las poblaciones humanas no se encuentran al azar en los patrones espaciales y liquidación dependen de atributos geográficos estructurales. En este contexto, el análisis espacial se centra en el papel del espacio en la explicación del fenómeno que se investiga, ejemplificada por la Primera Ley de la Geografía de Tobler : “todo está relacionado con todo, y los lugares más que cerca están más relacionados de lugares lejanos” (Tobler, 1970). La dimensión espacial de los fenómenos demográficos ha demostrado ser de gran importancia en la comprensión del papel de las características personales y el impacto del medio ambiente en este tipo de atributos. Sin embargo, la mayoría de los estudios tienden a ignorar esta dependencia espacial. Por ejemplo, si tenemos en cuenta el nivel de la tasa global de fecundidad (TGF), podemos decir que la TGF se autocorrelaciona espacialmente, es decir grupos de áreas muestran algún grado de dependencia, con valores similares para las zonas vecinas. Este es un punto importante, ya que la presencia de autocorrelación espacial puede sugerir la existencia de variables no observadas o no incluidas en el modelo.

Recordando la Primera Ley de la Geografía de Tobler, relaciones de distancia y vecinos entre diferentes áreas pueden ser particularmente importantes para comprender hasta qué punto es la dependencia espacial que existe y para entender “cómo establecer relaciones de vecindad” con el fin de estar relacionado, o espacialmente autocorrelacionados. De los diversos instrumentos utilizados en econometría espacial para comprender la dependencia espacial, índice I de Moran (Moran, 1950) es una de las estadísticas más utilizadas, ya que ayuda a cuantificar el nivel global de autocorrelación y discernir si se trata de un fenómeno aleatorio. (Gráfico 1) Sin embargo, el I de Moran no nos dice la “historia total”, y tenemos que complementarlo con otras herramientas como (semi) variograma, correlograma o análisis de variograma, que se refieren a la dependencia espacial a distancia por medio de covarianza, correlación y semivarianza a través de valores diferenciales observados entre vecinos ( Griffith y Paelinck, 2011: capítulo 3 ) y las medidas locales de asociación espacial, tales como I de Moran a nivel local para evaluar la agrupación y el significado de cada unidad espacial.

Obras recientes en el campo de la demografía espacial han evidenciado que la adición de la dimensión tiempo para el análisis espacial puede proporcionar información sobre la adopción de un nuevo régimen demográfico y cómo sus variables constitutivas son impactados a través del tiempo. Esta es una cuestión importante, ya que nos enteramos del proyecto de Princeton que la dimensión espacial es crucial para entender los procesos de difusión durante la primera transición demográfica en Europa ( Coale y Watkins, 1986 ). Sin embargo en la mayoría de los estudios de la Segunda Transición Demográfica, el componente espacial es a menudo pasado por alto. Esto es en parte debido a la disponibilidad de datos y también porque las transiciones demográficas son considerados como el resultado de un país procesos específicos. Pero centrarse en el nivel nacional en vez de la local al analizar los cambios en el régimen demográfico, por lo general pierden precursores, así como los rezagados. Un ejemplo clásico en España es la región de Cataluña, que fue un precursor de la Primera y la Segunda transiciones demográficas en comparación con el resto del país y de las regiones del Sur específicamente. Mapa 1. clustrs significativas para el índice de Princeton, 1981Mapa 2. agrupaciones significativas para el índice de Princeton, 2011

La forma más sencilla y práctica de la comprensión de cómo la dependencia espacial ha evolucionado a través del tiempo es por medio de las estadísticas locales de asociación espacial, en el que probar si y donde existen grupos de áreas con características similares. Anselin (1995) sugirió que los indicadores locales de asociación espacial , LISA, una técnica similar a la I de Moran, pero computarizada y evaluado para cada unidad espacial, comparable a una regresión lineal entre la variable medida en una cierta ubicación y la misma magnitud de medida en cada ubicación.

Por lo tanto, es muy fácil de ver cómo espacial autocorrelación puede alterar el resultado de un estudio que no toma en cuenta el espacio, por lo tanto, el viejo adagio de la propiedad, “ubicación, ubicación, ubicación!” También se puede aplicar también a la demografía. En el contexto de la heterogeneidad espacial de la fertilidad, España es un país único en Europa, con una larga y bien documentada de la diversidad regional y provincial por más de dos siglos. Tener una mirada a los cambios de fertilidad municipales más de las tres últimas décadas puede ser muy indicativo de cómo 1. La fertilidad disminuye con diferentes trayectorias, 2. La reciente recuperación de la fecundidad ha interesado sólo determinadas zonas hasta el inicio de la reciente crisis económica, 3. Migración ha afectado profundamente los patrones de fecundidad en las grandes ciudades, pero dejó otras regiones afectados. En cuanto a las medidas globales de autocorrelación (ver Gráfico 1), podemos entender por qué la dependencia espacial es un fenómeno que evoluciona el tiempo que puede cambiar y revertir su camino. Por ejemplo, el gráfico 1 muestra cómo en tiempos de expansión económica -y Fertilidad, autocorrelación espacial alcanza su pico, mientras que en tiempos de recesión económica -y la fertilidad declinación- que cae en picado, estabilizándose gran parte de las diferencias de fertilidad entre las regiones. Esto se debe principalmente a la forma en que las personas tratan de hacer frente a veces en dificultades al retrasar los nacimientos hasta que vengan tiempos mejores. Los dos mapas LISA grupo de clúster de la variable de interés, en este caso de Princeton Índice, utilizando una estadística local de cuatro grupos divididos como: rojo alta altos cúmulos de áreas con -relativamente- alta fertilidad rodeadas de alta fertilidad, azul oscuro bajo-bajo clusters, la luz roja de alta bajo racimos de altas áreas de fertilidad rodeadas de baja fertilidad y de color azul claro bajo altos cúmulos. La tradicional división de España en la alta fertilidad del Sur y la baja fertilidad del Norte se ha desplazado desde mediados de los años 90 en una división Este-Oeste con grandes ciudades como puntos calientes de la alta fecundidad, como se muestra en los mapas de la LISA 1 y 2.

Aunque las técnicas espaciales en la demografía a menudo se aplican a áreas pequeñas, el enfoque a gran escala puede abordar grandes cuestiones cuando un método más heurística falla. Mapeo puede ser una poderosa herramienta para entender la dinámica geográfica, pero sin necesidad de herramientas econométricas, temas tan importantes como la aleatoriedad y la significación estadística puede sesgar sustancialmente nuestros resultados. Por otra parte, la recolección de datos SIG se está volviendo más y más común en la demografía y en el espacio definitivamente arrojar nueva luz sobre los fenómenos demográficos.

Blogpost on Demotrends: Location, location, location! Why space matters in demography and why we should care. https://demotrends.wordpress.com/2014/11/06/location-location-location/