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/

Advertisements

The marriage market in XIX century Spain. La soltería en España en el siglo XIX

Using data on single men (21-35 years old) and women (16-30 years old) you can map where the unbalances in the marriage market are caused by “excess” of male or female population

Using data on single men (21-35 years old) and women (16-30 years old) you can map where the unbalances in the marriage market are caused by “excess” of male or female population.

Untitled.png

A match made in R: checking the order of geographical areas in shape files and in your data frames

Not every shape file is as nice as those provided in libraries. Sometimes we have to deal with historical maps, which have been hand-drawn, re-touched and what not. To work with geo-referenced data it is essential to have a variable in both shape file and dataframe with unique coding that has exactly the same number of areas and the same ordering in both files.

A quick way to check if shapefile and dataframe have the same number of areas:

nrow(df) == length(shape.file$Code)

In the shapefile, one can also select a couple of areas big enough so that they can easily be located, and plot them as “control” areas.
For instance, I want to select the area with code “15078” in the shapefile:
>which(shape.file$Code=="15078",arr.Ind=T)
[1] 271

which is the area in the 271-th position (same way shape.file$Code[271] gives the code of area 271).
plot(shape.file)
plot(shape.file[c(271,898),],col="red",border="red",add=T)

this is an easy way to locate your “control” area(s).
Rplot
Ideally, you should have some variable that is identical to the one in the shapefile, a codification of some sort, providing a unique Code, the name of the area or some factors that allow you to locate the area in space.

An easy way to check if both shape file and data frame have the same ordering of geographical areas is to test it:
>code.sh <- cbind(c(1:length(shape.file$Code)),as.vector(shape.file$Code))
>code.df <- cbind(c(1:nrow(df)),df$Code)
>code.df==code
.sh
[,1]  [,2]
[1,] TRUE  TRUE
[2,] TRUE  TRUE
[3,] TRUE  TRUE

What if it’s not?
First option: the inelegant solution
Manually change the order of the areas in a csv file according to the exact order they have in the shape file. It’s easy as you can create an ordinal index for the shapefile codes, paste it in excel, and assign it with a vlookup function.
Second option: the smart R match
In R there is a function called match that returns a vector of the positions of first matches of the first argument in its second:
>my.match <- match(df$Code, shape.file$Code)
NB: to use match the two variables providing the code for the areas have to have the very same unique and identical codes, or else funny stuff happens. To check that everything is in its right place, you can plot the two “control” spatial polygons we chose in the beginning, using their position in the dataframe rather than in the shapefile:
>plot(shape.file)
>plot(shape.file[c(which(df$Code=="305"),which(df$Code=="15078")),],col="orange",add=T)

Creating neighborhood matrices for Spatial Polygons in R (updated)

One of the first steps in spatial analysis is to create a neighborhood matrix, that is to say create a relationship/connection between each and (ideally!) every polygon. Why? Well, given that the premise for spatial analysis is that neighboring locations are more similar than far away locations, we need to define what is “near”, a set of neighbors for each location capturing such dependence.

There are many ways to define neighbors, and usually, they are not interchangeable, meaning that one neighborhood definition will capture spatial autocorrelation differently from another.

In R the package spdep allows to create a neighbor matrix according to a wide range of definitions: contiguity, radial distance, graph based, and triangulation (and more). There are 3 main and most used neighbors:

A) Contiguity based of order 1 or higher (most used in social sciences)

B) Distance based

C) Graph based

Install and load the maptools and spdep libraries shapefile from North Carolina counties:

library(maptools)
library(spdep)
NC= readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))

A. Contiguity based relations

are the most used in the presence of irregular polygons with varying shape and surface, since contiguity ignores distance and focuses instead on the location of an area. The function poly2nb allows to create 2 types of contiguity based relations:

1. First Order Queen Contiguity

FOQ contiguity defines a neighbor when at least one point on the boundary of one polygon is shared with at least one point of its neighbor (common border or corner);

nb.FOQ = poly2nb(NC, queen=TRUE, row.names=NC$FIPSNO)
#row.names refers to the unique names of each polygon
nb.FOQ
## Neighbour list object:
## Number of regions: 100
## Number of nonzero links: 490
## Percentage nonzero weights: 4.9
## Average number of links: 4.9

Calling nb.FOQ you get a summary of the neighbor matrix, including the total number of areas/counties, and average number of links.

2. First Order Rook Contiguity

FOR contiguity does not include corners, only borders, thus comprising only polygons sharing more than one boundary point;

nb.RK = poly2nb(NC, queen=FALSE, row.names=NC$FIPSNO)
nb.RK
## Neighbour list object:
## Number of regions: 100
## Number of nonzero links: 462
## Percentage nonzero weights: 4.62
## Average number of links: 4.62

NB: if there is a region without any link, there will be a message like this:
Neighbour list object:
Number of regions: 910
Number of nonzero links: 4620
Percentage nonzero weights: 0.5924405
Average number of links: 5.391209
10 regions with no links:
1014 3507 3801 8245 9018 10037 22125 30005 390299 390399

where you can identify the regions with no links (1014, 3507,…) using which(…), and in R it is possible to “manually” connect them or change the neighbor matrix so that they can be included (such as graph or distance based neighbors).
Sometimes, it also happens that some polygons that have been retouched (sounds like a blasphemy but it happens a lot with historical maps) may not recognize shared borders. This is when manually setting up neighbors comes in handy (you can’t do that in Geoda).

Contiguity

Higher order neighbors are useful when looking at the effect of lags on spatial autocorrelation and in spatial autoregressive models like SAR with a more global spatial autocorrelation:

nb.SRC = nblag(nb.RK,2) #second order rook contiguity
nb.SRC
## [[1]]
## Neighbour list object:
## Number of regions: 100
## Number of nonzero links: 490
## Percentage nonzero weights: 4.9
## Average number of links: 4.9
##
## [[2]]
## Neighbour list object:
## Number of regions: 100
## Number of nonzero links: 868
## Percentage nonzero weights: 8.68
## Average number of links: 8.68
##
## attr(,"call")
## nblag(neighbours = nb.RK, maxlag = 2)

Contiguity2

B. Distance based neighbors

DBN defines a set of connections between polygons either based on a (1) defined Euclidean distance between centroids dnearneigh or a certain (2) number of neighbors knn2nb (e.g. 5 nearest neighbors);

coordNC = coordinates(NC) #get centroids coordinates
d05m = dnearneigh(coordNC, 0, 0.5, row.names=NC$FIPSNO)
nb.5NN = knn2nb(knearneigh(coordNC,k=5),row.names=NC$FIPSNO) #set the number of neighbors (here 5)
d05m
## Neighbour list object:
## Number of regions: 100
## Number of nonzero links: 430
## Percentage nonzero weights: 4.3
## Average number of links: 4.3
nb.5NN
## Neighbour list object:
## Number of regions: 100
## Number of nonzero links: 500
## Percentage nonzero weights: 5
## Average number of links: 5
## Non-symmetric neighbours list

a little trick: if you want information on neighbor distances whatever the type of neighborhood may be:

distance = unlist(nbdists(nb.5NN, coordNC))
distance
##   [1] 0.3613728 0.3693554 0.3864847 0.2766561 0.5168459 0.3709748 0.2607982
##   [8] 0.3232974 0.4376632 0.2862144 0.5773310 0.3778483 0.4463538 0.2914539
## ...
## [498] 0.3407192 0.3995114 0.1838115

Distance

C. Graph based (I’ve never used them, but it’s good to know that they exist)

Delauney triangulation tri2nb constructs neighbors through Voronoi triangles such that each centroid is a triangle node. As a consequence, DT ensures that every polygon has a neighbor, even in presence of islands. The “problem” with this specification is that it treats our area of study as if it were an island itself, without any neighbors (as if North Carolina were an island with no Virginia or South Carolina)… Therefore, distant points that would not be neighbors (such as Cherokee and Brunswick counties) become such;
Gabriel Graph gabrielneigh is a particular case of the DT, where a and b are two neighboring points/centroids if in the circles passing by a and b with diameter ab does not lie any other point/centroid;
Sphere of Influence soi.graph: twopoints a and b are SOI neighbors if the circles centered on a and b, of radius equal to the a and b nearest neighbour distances, intersect twice. It is a sort of Delauney triangulation without the longest connections;
Relative Neighbors relativeneigh is a particular case of GG. A border belongs to RN if the intersection formed by the two circles centered in a and b with radius ab does not contain any other point.

delTrinb = tri2nb(coordNC, row.names=NC$FIPSNO) #delauney triangulation
summary(delTrinb)
## Neighbour list object:
## Number of regions: 100
## Number of nonzero links: 574
## Percentage nonzero weights: 5.74
## Average number of links: 5.74
## Link number distribution:
##
##  2  3  4  5  6  7  8  9 10
##  1  2 13 29 27 22  3  1  2
## 1 least connected region:
## 37039 with 2 links
## 2 most connected regions:
## 37005 37179 with 10 links
GGnb = graph2nb(gabrielneigh(coordNC), row.names=NC$FIPSNO) #gabriel graph
summary(GGnb)
## Neighbour list object:
## Number of regions: 100
## Number of nonzero links: 204
## Percentage nonzero weights: 2.04
## Average number of links: 2.04
## 20 regions with no links:
## 37109 37131 37137 37141 37145 37147 37151 37159 37161 37165 37173 37175 37179 37183 37185 37187 37189 37195 37197 37199
## Non-symmetric neighbours list
## Link number distribution:
##
##  0  1  2  3  4  5  6  7
## 20 27 16 15 13  7  1  1
## 27 least connected regions:
## 37047 37053 37055 37075 37091 37105 37107 37113 37115 37117 37119 37121 37129 37133 37135 37139 37143 37149 37153 37155 37157 37163 37167 37177 37181 37191 37193 with 1 link
## 1 most connected region:
## 37057 with 7 links
SOInb = graph2nb(soi.graph(delTrinb, coordNC), row.names=NC$FIPSNO) #sphere of influence
summary(SOInb)
## Neighbour list object:
## Number of regions: 100
## Number of nonzero links: 470
## Percentage nonzero weights: 4.7
## Average number of links: 4.7
## Link number distribution:
##
##  1  2  3  4  5  6  7  9
##  1  5 12 26 30 15 10  1
## 1 least connected region:
## 37031 with 1 link
## 1 most connected region:
## 37097 with 9 links
RNnb = graph2nb(relativeneigh(coordNC), row.names=NC$FIPSNO) #relative graph
summary(RNnb)
## Neighbour list object:
## Number of regions: 100
## Number of nonzero links: 133
## Percentage nonzero weights: 1.33
## Average number of links: 1.33
## 31 regions with no links:
## 37047 37053 37097 37107 37109 37115 37131 37137 37141 37143 37145 37147 37151 37155 37159 37161 37163 37165 37167 37173 37175 37179 37183 37185 37187 37189 37191 37193 37195 37197 37199
## Non-symmetric neighbours list
## Link number distribution:
##
##  0  1  2  3  4
## 31 30 18 17  4
## 30 least connected regions:
## 37009 37027 37031 37035 37037 37039 37055 37073 37075 37083 37091 37095 37105 37113 37117 37119 37121 37125 37127 37129 37133 37135 37139 37149 37153 37157 37169 37171 37177 37181 with 1 link
## 4 most connected regions:
## 37001 37003 37059 37079 with 4 links

GraphBased

What to do with all this stuff? …

compute and compare global Moran’s I
LISA maps
Variograms and correlograms
…?

Mean Age at Childbearing in Spain 2011

TFR 2011 fixed

A space-time box plot of Spain’s TFR for 910 comarcas.

The idea behind spatial analysis is that space matters and near things are more similar: a variable measured in city A is (ideally) different from the same variable measured in city B. A simple way to get a feeling and to represent this hypothesis is through graphical visualization, usually a map(s).

TFRG_all_4years_Spain

However, when dealing with time series maps are cumbersome and  with sometimes some information is lost, such as the national average or path convergence. Box plots are a simple yet very effective way to synthesize a lot of information in one graph. The following plot depicts TFR over a 30 years period for 910 Spanish areas with respect to the national average value (thick black line in the middle of the boxes).

p <- ggplot(dat, aes(x=factor(YEAR), y=dat$TFR))
p <- p + geom_boxplot()
p <- p + scale_y_continuous(limits=c(0,2.5)) + scale_x_discrete("YEAR", breaks=seq(1981,2011,by=5))

TFRG

Moran plots in ggplot2

Moran plots are one of the many way to depict spatial autocorrelation:
moran.test(varofint,listw)
where “varofint” is the variable we are studying, “listw” a listwise neighbourhood matrix, and the function “moran.test” performs the Moran’s test (duh!) for spatial autocorrelation and is included in the spdep funtionality. The same plot can be done using ggplo2 library. Provided that we already have our listwise matrix of neighborhood relationships listw, we first define the variable and the lagged variable under study, computing their mean and saving them into a data frame (there are a lot of datasets you can find implemented in R: afcon, columbus, syracuse, just to cite a few). The purpose is to obtain something that looks like this (I have used my own *large* set of Spanish data to obtain it):

ggplot2.moranplot1

Upload your data. Here is Anselin (1995) data on African conflicts, afcon:

data(afcon)
varofint listw varlag var.name <- "Total Conflicts"
m.varofint m.varlag
and compute the local Moran's statistic using localmoran:

lisa
and save everything into a dataframe:
df

use these variables to derive the four sectors "High-High"(red), "Low-Low"(blue), "Low-High"(lightblue), "High-Low"(pink):
df$sector significance vec =df$m.varofint & df$varlag>=df$m.varlag]  df$sector[df$varofint<df$m.varofint & df$varlag<df$m.varlag]  df$sector[df$varofint<df$m.varofint & df$varlag>=df$m.varlag]  =df$m.varofint & df$varlag<df$m.varlag]

df$sec.data

df$sector.col[df$sec.data==1] <- "red"
df$sector.col[df$sec.data==2] <- "blue"
df$sector.col[df$sec.data==3] <- "lightblue"
df$sector.col[df$sec.data==4] <- "pink"
df$sector.col[df$sec.data==0] <- "white"

df$sizevar df$sizevar 0.1)
df$FILL df$BORDER
to get the ggplot graph:
p 0.05", "High-High", "Low-Low","Low-High","High-Low"))+
scale_x_continuous(name=var.name)+
scale_y_continuous(name=paste("Lagged",var.name))+
theme(axis.line=element_line(color="black"),
axis.title.x=element_text(size=20,face="bold",vjust=0.1),
axis.title.y=element_text(size=20,face="bold",vjust=0.1),
axis.text= element_text(colour="black", size=20, angle=0,face = "plain"),
plot.margin=unit(c(0,1.5,0.5,2),"lines"),
panel.background=element_rect(fill="white",colour="black"),
panel.grid=element_line(colour="grey"),
axis.text.x  = element_text(hjust=.5, vjust=.5),
axis.text.y  = element_text(hjust=1, vjust=1),
strip.text.x  = element_text(size = 20, colour ="black", angle = 0),
plot.title= element_text(size=20))+
stat_smooth(method="lm",se=F,colour="black", size=1)+
geom_vline(xintercept=m.varofint,colour="black",linetype="longdash")+
geom_hline(yintercept=m.varlag,colour="black",linetype="longdash")+
theme(legend.background =element_rect("white"))+
theme(legend.key=element_rect("white",colour="white"),
legend.text =element_text(size=20))

Check out the interactive shiny version on pracademic