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.

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.

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.

No sé si esto que voy a contar me obliga a tragarme mis propias palabras. Porque siempre he pensado que era poco menso que imposible. Pero hace unos pocos días escribí sobre el asunto y hoy traigo otro similar a colación.

La variable más importante a la hora de construir un modelo es, precisamente, la que se quiere predecir. Casi todos los textos asumen que se conoce sin ningún género de dudas en, al menos, una determinada muestra que, además, corresponde más o menos a la población subyacente: si el paciente sobrevive o no; si la hipoteca entra en mora o no; si el cliente responde a la oferta o no, etc.

Pero hay muchos problemas famosos y relativamente urgentes en los que la situación es distinta. En la entrada a la que hago referencia más arriba, sólo se conocía el valor que predecir para un conjunto de casos, los positivos. Pero era desconocido para la gran masa.

El problema aparece también en el riesgo de crédito: el banco sólo tiene información sobre la situación crediticia de aquellos clientes que no fueron rechazados previamente. Pero es necesario crear un mecanismo de medición del riesgo para todos los clientes. Y todos es una población distinta de aceptados. ¡Y qué peligrosas son las extrapolaciones!

En otro contexto en el que aparece es en el de la determinación de lo que llaman share of wallet, el porcentaje de, por ejemplo, la cesta de la compra que realiza un consumidor en una determinada cadena de supermercados (desconociéndose las compras que realiza en los de la competencia) o la cantidad de transacciones financieras que realiza en una determinada entidad (cuando existe la posibilidad de que tenga también cuentas activas en otras).

Ahora me encuentro con la corrección de Heckman, que le valió al susodicho el premio Nobel de economía. Traduzco de la Wikipedia:

Supóngase que un investigador quiere estimar cuáles son los determinantes de las ofertas de salarios pero sólo tiene acceso a los salarios de aquellos que trabajan. Dado que quienes trabajan forman una muestra no aleatoria de la población, estudiar estos determinantes sobre este subconjunto introduciría un sesgo.

Información adicional sobre esta corrección puede encontrarse en este enlace.

Nota: El modelo puede implementarse en R usando el paquete sampleSelection

He tropezado con una extensión curiosa y que no conocía del modelo logístico que lo emparenta un tanto con los modelos de supervivencia. Es un problema que aparece en los modelos de los actuarios, por ejemplo, y en la supervivencia de nidos (sí, nidos de bichos alados), parece.

Es el siguiente: supongamos que unos sujetos están expuestos a un cierto suceso cuya probabilidad, , depende del sujeto a través del esquema habitual de la regresión logística (es decir, depende de algunas variables como el sexo, etc., a través de una fórmula lineal cuyos coeficientes interesa estimar).

Supongamos también, y esta es la novedad, que no todos los sujetos están expuestos al factor de riesgo el mismo tiempo. Pero si el sujeto i está expuesto al riesgo durante periodos, la probabilidad de que no le ocurra cierta calamidad es .

Planteemos el problema en R:

# número de sujetos
n <- 10000
 
# coeficientes "verdaderos"
beta.0 <- 1
beta.1 <- -2
 
# un factor que afecta a la probabilidad
x <- factor( sample( c("a", "b"), n, replace = T ) )
 
# número de días de "exposición"
days <- sample( c( 1,2,3,4 ), n, replace = T )
 
# probabilidad de que *no* ocurra el suceso
p <- plogis( beta.0 + beta.1 * (x == "b") )^days
 
# variable objetivo (1 significa que *no* ocurre el suceso)
y <- rbinom( n, 1, p )
 
# ¡gráfico!
mosaicplot( table( days, x, y ) )
 
 
# modelo logístico "normal"
dat <- data.frame( x = x, y = y, days = days )
m.0 <- glm( y ~ x, family = binomial(), data = dat )
m.0$coefficients

Adaptando código ajeno a nuestro contexto, podemos escribir:

logexp <- function(days)
{
    linkfun <- function(mu) qlogis(mu^(1/days))
    linkinv <- function(eta) plogis(eta)^(days)
    mu.eta <- function(eta) days * plogis(eta)^(days-1) * .Call("logit_mu_eta", eta, PACKAGE = "stats")
    valideta <- function(eta) TRUE
    link <- paste("logexp(", days, ")", sep="")
    structure(list(linkfun = linkfun, linkinv = linkinv,
                   mu.eta = mu.eta, valideta = valideta, name = link),
              class = "link-glm")
}

Esta definición de la función de enlace es bastante peculiar: depende del sujeto no sólo a través de mu_i, como es habitual y nos enseñan los libros, sino también a través de la variable tiempo, que depende del sujeto.

Et voilá:

m.1 <- glm( y ~ x, family=binomial(logexp(days=dat$days)), data=dat )
m.1$coefficients

Prácticamente (y, en gran medida, gracias a que el número de observaciones es grande), recobramos el valor original de los coeficientes (cosa que como habrán comprobado los más diligentes de mis lectores, no ocurre con el primer modelo).

Imaginemos que queremos predecir y, que toma valores 0 y 1 a partir de indicios (o variables) x mediante una función (un clasificador) f. Podemos visualizar el error de clasificación usando la matriz

Efectivamente, el error es (perdónenme la notación) B+C. Como estadísticos estamos habituados a tratar de minimizar C mientras mantenemos B lo suficientemente pequeño. Un test es tanto más potente cuanto menor es C.

Pero podemos reescribir B+C como

B+C = (A+B) - (A+C) + 2C

Identifiquemos estas partes:

• A+C es fijo y corresponde a la probabilidad (o frecuencia) de que y=1.
• A+B corresponde a los casos en que f(x)=1.
• C es el error que se produce cuando y=1 y f(x)=0, es decir, la incapacidad de detectar el valor 1.

Esta representación es inhabitual. Pero permite controlar el error tratando de minimizar A+B por una parte y 2C por otra; es decir, por una parte, tratando de asignar el menor número posible de valores 1. Y por la otra, tratando de reducir el número total de 0 erróneos sobre la colección de los y=1.

¿Por qué es útil esta representación? Pues porque permite afrontar un problema habitual en muchos ámbitos: tratar de predecir un fenómeno cuando éste no es observable en toda la población, pero existe un número suficiente de marcas y=1 conocido. Un ejemplo: buscamos un clasificador de artículos científicos que distinga los de matemáticas de los de otras disciplinas. Tenemos un corpus amplio de artículos sin clasificar de todo tipo y un conjunto de artículos de matemáticas (nuestros y=1).

Otro ejemplo: tenemos una población inmensa de “clientes” entre los que hay “defraudadores” pero no sabemos cuáles son. Pero tenemos una colección de “defraudadores” previamente identificados.

Esta aproximación al problema viene no sin ciertos caveats presenta serios peligros: ¿qué pasaría con los artículos de estadística? ¿Cómo podríamos, a ciegas, saber si son de matemáticas o, como sucede realmente, de una disciplina semántica y sintácticamente próxima?

Sobre estos asuntos podrá averiguar más quien lea Partially Supervised Classification of Text Documents.

Hablé el otro día con Emilio Torres y comentamos de pasada la situación política en Asturias, donde vive, después de las últimas elecciones. El escaño obtenido por UPyD otorgaba a tal partido un poder en exceso del tamaño de su representación porque era clave para formar el futuro gobierno del principado. Pero, ¿cuánto poder realmente supone ese escaño en esas condiciones? ¿Puede cuantificarse?

Porque se habla mucho en periodo electoral de la ley D’Hondt pero, una vez asignados los escaños, cambia el juego.

Existe un método, el índice de poder de Banzhaf. El índice de Banzhaf para un determinado partido político mide su poder en términos del porcentaje de las posibles alianzas mínimas ganadoras en las que participa dentro de su universo total. Una alizanza es ganadora cuando reúne más de la mitad de los votos. Y es mínima cuando todos sus integrantes son necesarios para que sea ganadora; excluye, por ejemplo, la alianza trivial formada por todos los partidos.

Veamos cómo calcular este índice con R y lo utilizaremos para cuantificar el valor de ese escaño:

escannos <- c(17,12,10,5,1)
names(escannos) <- c( "psoe", "fac", "pp", "iu", "upyd")
 
banzhaf <- function(x){
  x <- -sort(-x)
  x <- x/sum(x)
 
  foo <- function(a,b,p){
    if(p>1/2)
      return(list(a))
 
    if (length(b)==0)
      return(NULL)
 
    b.prima <- b[-1]
    delta <- b[1]
    p.delta <- x[delta]
 
    return(c( foo(c(a,delta), b.prima, p+p.delta), foo(a,b.prima,p)) )
  }
 
  res <- foo( NULL, names(x), 0)
 
  sort( table(unlist(res)) / length(res) )
 
}
 
banzhaf( escannos )

El resultado es:

  iu upyd  fac   pp psoe
 0.4  0.4  0.6  0.6  0.6

Es decir, el escaño de UPyD le concede el mismo poder (según Banzhaf) que los cinco de IU. Y la diferencia de siete escaños entre el PP y el PSOE no le concede a este último cuota de poder adicional alguna.

Existen limitaciones obvias a este indicador que resultarán evidentes a quien piense mínimamente en él o lo aplique a casos conocidos y concretos. Así que no abundaré en ellos.

Pero sí que dejaré, por referencia, otra aplicación de este índice al resultado de las últimas elecciones generales:

generales <- c(186,110,16,11,7,5,5,3,2,2,1,1,1)
names(generales) <- c("pp", "psoe", "ciu", "iu", "amaiur", "upyd",
   "pnv", "esquerra", "bng", "cc", "compromis", "fac", "gbai")
 
banzhaf( generales )
# pp 
#  1