Contar debería ser fácil. Sobre todo, cifras relativamente pequeñas. Y más para el Estado, con todo su —aunque menguante, según algunos— poder.

Por otro lado, dicen, las cifras son las que son. Y son incuestionables. Además, sirven para cuestionar. Por ejemplo, para cuestionar la labor de la Dirección General de Tráfico, cuyo rendimiento, en gran medida, se mide por el recuento anual de muertes en la carretera.

Y en este contexto, y a la luz de lo publicado recientemente en la prensa, diríase que sucede un fenómeno que me parece la mar de sorprendente: el número cobra una importancia que la metodología no puede matizar. Directores, exdirectores y asociaciones de no sé muy bien qué, se atizan —cual garroteros gollescos— con alternativas metodológicas persiguiendo intenciones que, uno sospecha, tienen poco que ver con la creación de un cuadro lo más fidedigno posible de la realidad. (A veces tengo la tentación de promover para este tipo de discusiones —así como existen las bizantinas— el término de discusión argentina).

El quid de la cuestión reside en si se cuentan como muertes relacionadas con el tráfico los fallecidos a consecuencia de accidentes hasta treinta días después y la dificultad que tiene la DGT para recopilar ese dato. También, el problema que supone obtener información sobre accidentes y fallecidos en zonas urbanas, en los que la gestión corresponde a los ayuntamientos. Y a mí me da por pensar, desde la perspectiva más alejada y neutra posible (puesto que ni coche tengo) lo siguiente:

• Igual que las cifras del paro, etc., suelen estacionalizarse, ¿no deberían igualmente tenerse en cuenta circunstancias como puentes, semanas santas, etc. para matizar las cifras de fallecidos en accidentes y permitir así la comparación entre periodos (p.e., años distintos)?
• De nuevo, para facilitar la comparación, ¿no deberían corregirse las cifras por los kilómetros recorridos? ¿Seguro que no se conducen hoy en día muchos más kilómetros que hace 10 años?
• Igualmente, ¿no deberían introducirse factores correctores por las condiciones meteorológicas?
• La incertidumbre que existe sobre el número definitivo de accidentes, ¿no podría (o debería) reflejarse en unos márgenes de confianza?
• Un camión se sale de su carril y aplasta un vehículo en el que mueren todos sus ocupantes. Pueden ser uno. O cinco. Incluso ninguno por un felicísimo azar. ¿No debería esta variabilidad corregir de alguna manera las cifras finales de fallecidos?
• Un desglose de fallecidos por tipo de accidente (urbano contra, por exceso de alcohol, por la lluvia, por exceso de velocidad, etc.) que pudiera compararse de año en año, ¿no serviría para comprender mejor el fenómeno de la siniestralidad que un agregado absoluto? ¿Por qué no se publican y discuten estas cifras desagregadas? ¿Se nos considera demasiado inmaduros para comprenderlas adecuadamente?
• Finalmente, para poder decidir sobre la bondad de un determinada política o un equipo gestor de la DGT, no debería haber un estadístico que calculase el p-valor asociado a la hipótesis nula nuestra nueva política circulatoria es igual de buena que la antigua?

Una de las principales objeciones que se le pueden hacer a mi entrada de ayer es que puede estar confundiendo la causa con efecto: puede que parte de la radialidad de la red que obtuve tenga que ver con el tamaño desproporcionado de Madrid que, a su vez, podría haber sido causado por la radialidad de la red tradicional de las comunicaciones españolas.

Así que enviemos una partida de pescado en malas condiciones a Mercamadrid, convidemos a toda la provincia, veámosla fenecer víctima de contumaces diarreas y rehagamos la simulación suponiendo que

nodos.alt <- nodos
nodos.alt$pop[nodos.alt$prov == "Madrid"] <- 0

¿Qué forma tendría ahora la red? Ejecutando

res  <- do.call(rbind, apply(aristas, 1, function(x) peso.tramos(x[1], x[2], g2, nodos.alt)))
peso <- tapply(res$pop / (res$distancia)^(1), res$tramo, sum)
 
col <- peso
col[col < median(col)] <- 0
col <- rgb(0,0,0, 255 * col/max(col), maxColorValue=255)
plot.grafo(g2, nodos, col = col)
 
g3 <- delete.edges(g2,edges=E(g2)[peso < median(peso)])
col3 <- col[peso >= median(peso)]
plot.grafo(g3, nodos, col = col3)
 
E(g3)$weight <- peso[peso >= median(peso)]
 
centralidad <- data.frame(nodo = V(g3)$name, centralidad = alpha.centrality(g3) )
centralidad <- centralidad[order(-centralidad$centralidad),]
centralidad

se obtiene

Esta vez la península se parte en dos reeditando una suerte de Hispania Tarraconensis en la que las capitales más centrales son Tarragona, Valencia, Barcelona, Castellón, Jaén y Lérida.

Si en lugar de esta versión tan extrema su ponemos que Madrid tiene una poblacion promedio, es decir,

nodos.alt <- nodos
nodos.alt$pop[nodos.alt$prov == "Madrid"] <- median( nodos$pop )

se obtiene una configuración prácticamente similar:

Y, finalmente, si toda la población está distribuida uniformemente en las provincias, es decir,

nodos.alt <- nodos
nodos.alt$pop <- mean( nodos$pop )

las cosas cambian de manera bastante sorprendente:

La vieja Castilla serviría de nexo de comunicaciones, siendo las provincias con mayor centralidad Palencia, Burgos, Segovia, Guadalajara, Toledo, Soria y Madrid. Tal vez porque la densidad de capitales de provincia (al suponer la población igual en todas las provincias) favorece esa zona en que se agrupan de manera un poco más compacta.

¿Será que consigo que alguien se anime a introducir correcciones por superficie (suponiendo una distribución uniforme de la población por kilómetro cuadrado), PIB, acceso a puertos de mar o salidas a Portugal o Francia?

Me propuse hace un tiempo combinar lo que aprendí creando rutas callejeras por Zaragoza con una entrada que escribí sobre la estructura radial de las vías de transporte de España. El problema que me planteo es si tiene sentido que la red de carreteras Española tenga estructura radial habida cuenta de la geometría peninsular bajo ciertas hipótesis, siempre discutibles y mejorables, de partida.

Así que, en primer lugar, cargué los paquetes de R necesarios, un fichero que creé que contenía las capitales de provincia, su latitud, su longitud y la población de las respectivas provincias y fabriqué una red de carreteras muy ineficiente que unía todos los nodos entre sí:

library(maptools)
library(pxR)
library(igraph)
library( geosphere)
 
## datos: provincias y población
 
nodos <- read.table( "http://www.datanalytics.com/uploads/prov_pop_lat_lon.txt",
             sep = ",", dec = ",", header = T, encoding = "latin1")
 
## aristas
 
distancia <- function(p1, p2, nodos){
  alfa  <- c(nodos$lon[nodos$prov==p1], nodos$lat[nodos$prov==p1] )
  omega <- c(nodos$lon[nodos$prov==p2], nodos$lat[nodos$prov==p2] )
  distCosine( alfa, omega ) / 1000	# kms.
}
 
aristas <- expand.grid(nodos$prov, nodos$prov)
colnames(aristas) <- c("desde", "hasta")
aristas <- aristas[ as.character(aristas$desde) < as.character(aristas$hasta), ]
 
aristas$weight <- apply(aristas,1, function(x) distancia(x[1], x[2], nodos))
 
 
## grafo
 
grafo <- graph.data.frame(aristas, directed = F)
 
plot.grafo <- function(g, nodos, col = "black", text = F){
  tmp <- get.edges(g, V(g))
  vertices <- data.frame( (V(g)[get.edges(g, E(g))[,1]])$name,
                          (V(g)[get.edges(g, E(g))[,2]])$name, col = col )
 
  plot(nodos$lon, nodos$lat, xlab ="", ylab = "", xaxt = "n", yaxt="n")
  if( text )
    text(nodos$lon, nodos$lat,nodos$prov)
 
  apply(vertices, 1, function(x){
    x1 <- nodos$lon[nodos$prov == x[1]]
    y1 <- nodos$lat[nodos$prov == x[1]]
    x2 <- nodos$lon[nodos$prov == x[2]]
    y2 <- nodos$lat[nodos$prov == x[2]]
    lines( c(x1,x2), c(y1,y2), col = x[3])
  })
}
 
plot.grafo(grafo, nodos)	# pequeño caos

El resultado es este pequeño caos:

Por simplificar, eliminé todas las autovías que unían capitales de provincia cuando pudiera encontrar una ruta alternativa cuya longitud no excediese a la original por un factor de 1.2 haciendo

exceso.edge <- function(g, e){
  a <- shortest.paths(g)
  b <- shortest.paths(delete.edges(g, e))
  max( incr <- b / a, na.rm = T )
}
 
tmp <- sapply(E(grafo), function(e) exceso.edge(grafo,e))
g2  <- delete.edges(grafo, E(grafo)[tmp < 1.2])
 
plot.grafo(g2, nodos)

para obtener

Finalmente, simulé trayectos entre provincias con este criterio: una persona viaja de A a B con una probabilidad directamente proporcional al producto de las poblaciones de dichas provincias e inversamente proporcional a la distancia (en línea recta) entre ellas. La regla del producto de la población de las provincias es compatible con una muestra aleatoria de parejas de personas sobre la población total modificada en segunda instancia por la distancia entre ellas. Así que haciendo

peso.tramos <- function(a, b, g, nodos){
  data.frame(
    tramo = as.numeric(E(g2, path = get.shortest.paths(g2, from=a, to = b)[[1]])),
    pop = nodos$pop[nodos$prov == a] / 1e6 * nodos$pop[nodos$prov == b] / 1e6,
    distancia = distancia(a,b,nodos)
  )
}
 
res  <- do.call(rbind, apply(aristas, 1, function(x) peso.tramos(x[1], x[2], g2, nodos)))
peso <- tapply(res$pop / (res$distancia)^(1), res$tramo, sum)
 
## un primer gráfico
 
col <- peso
col[col < median(col)] <- 0
col <- rgb(0,0,0, 255 * col/max(col), maxColorValue=255)
plot.grafo(g2, nodos, col = col)

obtuve

En este mapa sólo se han representado la mitad de los tramos de mayor importancia (de acuerdo con el criterio arriba especificado) y en el resto se ha modulado la intensidad en función también de ese criterio.

¿Es una estructura radial? Podemos recurrir de nuevo a la teoría de grafos para medir la centralidad de los nodos:

g3 <- delete.edges(g2,edges=E(g2)[peso < median(peso)])
col3 <- col[peso >= median(peso)]
plot.grafo(g3, nodos, col = col3)
 
E(g3)$weight <- peso[peso >= median(peso)]
 
centralidad <- data.frame(nodo = V(g3)$name, centralidad = alpha.centrality(g3) )
centralidad <- centralidad[order(-centralidad$centralidad),]
centralidad

El resultado da preeminencia a Madrid y otras capitales de su entorno:

La cuestión es: ¿está Madrid en el centro a causa de su población? ¿Es esta población de Madrid grande entre otras cosas, gracias a la estructura radial de las comunicaciones? En una nueva entrega sobre este asunto volveré a analizar el problema con hipótesis de partida distintas.

Aprovecho para recordar a los usuarios de R de Madrid que el jueves 26 de abril, a las siete de la tarde, tendrá lugar la segunda reunión del grupo de usuarios de R de Madrid en la sala Metrópolis de La Tabacalera (glorieta de Embajadores).

El programa, como siempre, puede consultarse en la página del grupo.

No sé por qué se llaman así. Ni idea. Vine a enterarme de tal nombre a través de un comentario de Carlos Ortega en la lista de usuarios de R. Parece que se usan habitualmente en diversas áreas y tienen el siguiente aspecto:

Se trata de diagramas de barras a los que se añaden unos apéndices que tratan de medir la variabilidad a la que se entiende que están sujetas aquellas. Una representación alternativa del mismo conjunto de datos ilustra el motivo por el que se desaconseja su uso:

¿Hacen falta palabras?

En Dynamite Plots: Unmitigated Evil se enumera una serie de motivos por los que es preferible no usarlos:

• su ratio tinta/datos es bajo
• esconden los datos originales, como se parecia en la comparación anterior
• suelen asumir que los intervalos de confianza son simétricos
• pueden esconder aspectos de los datos, como la bimodalidad
• cuando muestran sólo el tramo superior del intervalo de confianza dificultan la comparación entre grupos
• aparentemente, crean una ilusión visual por la que el lector tiende a agregar la longitud correspondiente al intervalo de confianza a la barra

En ese mismo enlace aparece el gráfico

en el que se comparan seis representaciones gráficas alternativas de los mismos datos. ¿Puede decirse que la de la esquina superior izquierda es la mejor?

Tras una sobremesa en la que tratamos el ya manido tema de los gráficos de tarta, me hace llegar mi tertuliano Jorge Sobrino una solapilla que le adjunta Iberdrola al recibo de la luz que parece una broma de mal gusto. Es muy parecida a la siguiente:

La gráfica compara el llamado mix de producción eléctrica de la compañía con el nacional en el periodo de referencia. O, al menos, ese parece ser su objetivo. Pero advertirán los más agudos de mis lectores cómo la elección de una gráfica de tartas dificulta la comparación, que es el presunto objetivo de esas manchas de colores con pretensiones cuantitativas. ¿No podían haber preferido, por ejemplo y sin pensar demasiado, barras adyacentes?

(Nota: la gráfica anterior está extraída de un documento de la página de Iberdrola y se refiere a la electricidad comercializada a empresas; en la factura de mi amigo, que corresponde a la comercializadora de último recurso, es decir, la que tenemos la mayoría de los ciudadanos de a pie, la proporción de energía renovable baja del 70.9 % a un mero 13.9 %: ¿sabrá Iberdrola cómo etiquetar los electrones que recorren el tendido eléctrico para enviar unos u otros a clientes distintos?).

Pero lo grave es que descubro que el pequeño informe de Iberdrola, con sus tartas 3D y todo, responde al modelo establecido por la CNE, la Comisión Nacional de la Energía, tal y como figura en su portal de internet (sección Listado de Informes de Etiquetado de Electricidad). En los documentos PDF que uno puede descargarse allí, descubre el modelo (¿obligadorio?) de informe, que tras los párrafos

Si bien la energía eléctrica que llega a nuestros hogares es indistinguible de la que consumen nuestros vecinos u otros consumidores conectados al mismo sistema eléctrico, ahora sí es posible garantizar el origen de la producción de energía eléctrica que usted consume.
A estos efectos se proporciona el desglose de la mezcla de tecnologías de producción nacional para así comparar los porcentajes del promedio nacional con los correspondientes a la energía vendida por su Compañía Comercializadora.

pinta la siguiente calamidad:

¿Graficaca por imperativo legal?

Los economistas usan unas cosas a las que llaman variables instrumentales con las que uno apenas se tropieza fuera de contextos econométricos. El problema se plantea en el contexto de la regresión

cuando existe correlación entre X y . En tales casos, el estimador por mínimos cuadrados es

y debido a la correlación entre X y , está sesgado.

La solución que se plantea en ocasiones es el de usar variables instrumentales, es decir, variables correlacionadas con X pero no con . La siguiente simulación en R ilustra el problema:

library(MASS)
library(AER)

n <- 100

b.0 <- 1
b.1 <- 2

cov.mat <- diag(3)
cov.mat[2,1] <- cov.mat[1,2] <- sqrt(0.5)
cov.mat[3,1] <- cov.mat[1,3] <- sqrt(0.5)

foo <- function(){
	xze <- data.frame(mvrnorm(n, rep(0,3), cov.mat ))
	colnames( xze ) <- c("x", "z", "e")

	xze$y <- b.0 + b.1 * xze$x + xze$e

	c.m1 <- coefficients( lm( y ~ x, data = xze ) )
	c.m2 <- coefficients( ivreg( y ~ x | z, data = xze ) )

	data.frame( c.1.1 = c.m1[1], c.1.2 = c.m1[2],
		c.2.1 = c.m2[1], c.2.2 = c.m2[2] )
}

foo()

res <- replicate( 1000, foo(), simplify = F )
res <- do.call( rbind, res )

Lo que se hace en ella es construir 1000 conjuntos de datos con las variables y, x, z y el error e. Las tres últimas son normales y están construidas de forma que:

• La correlación entre x y el error e no es nula.
• La correlación entre x y z tampoco es nula.
• La correlación entre z y el error e es nula.

Luego, se construye y como b.0 + b.1 * x + e. Finalmente, se comparan los coeficientes obtenidos por la regresión por mínimos cuadrados tradicional con la que se obtiene usando z como variable instrumental. La comparación de los coeficientes obtenidos puede observarse gráficamente haciendo

library(ggplot2)
library(reshape)

mi.coef <- res[,c(2,4)]
colnames(mi.coef) <- c("ols", "vi")
mi.coef <- melt( mi.coef )

ggplot(mi.coef, aes(x=value, fill=variable)) + geom_density(alpha=.3)

que produce el gráfico

En esencia, lo que se ha hecho es calcular el coeficiente condicionando previamente por z, es decir, calculando

que anula el error. Por otra parte, como puede apreciarse en el gráfico, se pierde precisión, es decir, aumenta la varianza del estimador (y se ensancha la campana alrededor de su valor medio) por la pérdida de información que supone condicionar x a z.

El lector interesado hará bien en:

• Estudiar el efecto de cambiar la correlación entre x y z en la anchura de la campana.
• Leer la entrada de la Wikipedia sobre variables instrumentales para entender cómo se extiende este caso univariante al general.
• Echarle un vistazo al artículo Instrumental Variables Estimation in Political Science: A Readers’s Guide para conocer mejor el uso y el abuso de este tipo de técnicas en aplicaciones reales.

El jueves pasado, en MediaLab Prado, tuve ocasión de asistir a una presentación de los responsables de Via52,

un nuevo un semanario digital que quiere sumarse con modestia al panorama mediático. Desde hace meses venimos trabajando en este proyecto, impulsado por David Rojo (@rojovegas) y Andrés Hermosa (@andresh), y que cuenta con la colaboración de un grupo de profesionales del periodismo, el fotoperiodismo, la ilustración y la tecnología.

Lo más interesante de esta publicación para quienes siguen esta bitácora es el énfasis que hacen en el periodismo de datos. Además:

Las bases de datos que Vía52 utilice o elabore serán difundidas en un catálogo abierto a disposición de otros periodistas para que puedan utilizarlas en sus propios trabajos.

Tengo entendido, además, que próximamente van a sacar un número muy jugoso… y hasta ahí puedo leer escribir.

Óscar Perpiñán me puso el otro día al tanto del paquete osmar de R, que

proporciona la infraestructura para acceder a datos de OpenStreetMap a través de diferentes fuentes, trabajar con ellos con R de una manera unificada y aprovechando la infraestructura que proporcionan otros paquetes como, por ejemplo, sp e igraph.

Hoy voy a ilustrar el uso de este paquete adaptando un ejemplo de sus autores para encontrar la ruta óptima entre dos puntos de Zaragoza, la mercería Bell y el colegio La Salle Montemolín, ambos lugares muy vinculados a mi prehistoria. Comenzaré cargando los paquetes necesarios y los datos de OpenStreetMap correspondientes a Zaragoza:

library( igraph )
library( osmar )
 
api <- osmsource_api(url = "http://api.openstreetmap.org/api/0.6/")
box    <- corner_bbox( -0.90, 41.60, -0.85, 41.69 )
zgz <- get_osm(box, source = api)

En segundo lugar, voy a usar las funciones auxiliares del paquete para extraer las calles del objeto zgz y representarlas geométricamente:

hways_zgz <- subset(zgz, way_ids = find(zgz, way(tags(k == "highway"))))
hways <- find(hways_zgz, way(tags(k == "name")))
hways <- find_down(zgz, way(hways))
hways_zgz <- subset(zgz, ids = hways)
 
plot_ways(hways_zgz, col = "gray" )
title("Calles de Zaragoza")

El resultado es

A continuación, selecciono los puntos incial y final de mi ruta:

id <- find(zgz, node(tags(v == "Moda infantil Bell")))[1]
hway_start_node <- find_nearest_node(zgz, id, way(tags(k == "highway")))
hway_start <- subset(zgz, node(hway_start_node))
 
id <- find(zgz, node(tags(v == "La Salle Montemolín")))[1]
hway_end_node <- find_nearest_node(zgz, id, way(tags(k == "highway")))
hway_end <- subset(zgz, node(hway_end_node))

Puedo representarlos mediante

plot_nodes(hway_start, add = TRUE, col = "red", pch = 19, cex = 1)
plot_nodes(hway_end, add = TRUE, col = "blue", pch = 19, cex = 1)

para obtener

Finalmente, utilizo la infraestructura proporcionada por el paquete igraph (para el análisis de redes sociales) para calcular la ruta más corta entre ambos puntos haciendo

gr_zgz <- as_igraph(hways_zgz)
route <- get.shortest.paths(gr_zgz,
                            from = as.character(hway_start_node),
                            to = as.character(hway_end_node))[[1]]
 
route_nodes <- as.numeric(V(gr_zgz)[route]$name)
route_ids <- find_up(hways_zgz, node(route_nodes))
route_zgz <- subset(hways_zgz, ids = route_ids)
route_zgz

y puedo finalmente representarla en el mapa mediante

plot_ways(route_zgz, add = TRUE, col = "black", lwd = 2)

para obtener el resultado final

Es incluso posible (véase esto) obtener una lista de las calles y nodos que hay que seguir para ir de uno de los puntos al otro.

Interesante, ¿eh?

Notas:

• No he conseguido calcular rutas en zonas más amplias: descargar los mapas de carreteras/calles de toda España, por ejemplo, no es tan sencillo como descargar los de zonas más pequeñas (una ciudad, por ejemplo). El API por defecto de OpenStreetMap, el que he usado, sólo permite descargar zonas que ocupen menos de la cuarta parte de un grado cuadrado. Existen procedimientos para extracciones más amplias, pero no las he probado.
• Para otras ubicaciones distintas de las que he probado (y que existen como tales en el objeto zgz como, por ejemplo, el campus de la Universidad de Zaragoza) el código que he mostrado falla de manera que no he podido todavía explicar.

En España, lo de las balanzas fiscales es como lo de las manifestaciones: un número que se tiran a la cabeza y con muy mala baba tirios y troyanos. La cantinela que más se oye es la de la prensa periférica (perdón, prensa de la parte este de la periferia: existen otras periferias que callan como palabras de cuatro letras): dizque Cataluña aporta mucho más a la hacienda pública que lo que después recibe de ella por inversiones y servicios.

El interesado en este asunto puede echarle un vistazo a su página de la Wikipedia y, en particular, al documento que describe los métodos usados para calcular las balanzas fiscales.

Cabe señalar en primer lugar que tiene todo el sentido del mundo usar el plural para referirse a esta especie de ficción contable: existen seis métodos distintos con seis resultados igualmente distintos. Por ejemplo, los balances fiscales de la Comunidad de Madrid, la más perjudicada por la política fiscal española, van del -5.57 % al -9.13 % del PIB regional. Los de Melilla, en el extremo contrario, del 6.46 % al 34.27 % de su PIB. Seis valores distintos por región, en definitiva, para que cada cual elija el que más le convenga.

¿Imaginan mis lectores que en sus coches tuviesen seis velocímetros discordantes? ¿Imaginan que uno dice que circula a 55 km/h, otro que a 91 km/h y los demás, alguna otra velocidad intermedia? Pues parece que esta algarabía de cifras, gracias al nunca escaso abono del anumerismo, es clave para la gobernanza de esta peculiar nación.