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 

He leído recientemente el artículo Social Ties and their Relevance to Churn in Mobile
Telecom Networks porque ilustra una técnica muy de moda: el análisis de redes sociales (SNA) en en ámbito de las telecomunicaciones y, en particular, la construcción de indicadores tempranos de baja (churn) de clientes de telefonía móvil. Más aún, permite rediseñar estrategias basadas en los resultados para retener clientes: al clasificarlos mejor usando técnicas de SNA, pueden diseñarse estrategias activas para aquellos que no sólo tienen una mayor predisposición a darse de baja sino, además, a arrastrar con ellos a parte de su entorno social.

El artículo, en resumen, introduce dos indicadores. El primero, p(k), es más ilustrativo que práctico: se trata de la probabilidad de que un cliente que tiene k conexiones —una conexión es alguien con quien el cliente ha hablado durante un determinado periodo— que se han dado de baja previamente se dé él mismo de baja. El gráfico siguiente muestra cómo p(k) es una función creciente de k. Sin embargo, el indicador puede no ser particularmente útil dado que, estoy seguro, el número de clientes para los que k > 1 es, casi seguro, muy pequeño.

En la segunda parte los autores construyen un modelo de propagación. Les interesa no sólo contar —y construir, de paso, probabilidades de corte frecuentista— sino explicar la dinámica y aprovecharla para construir modelos más útiles. La idea es la siguiente: un cliente que se da de baja transmite una señal a aquellos con los que se comunica. La señal puede ser del tipo esta compañía es malísima, me voy a ir a esta otra. No se sabe realmente cómo es la influencia, pero los autores la aproximan de la siguiente manera:

1. Asignan a cada cliente que se da de baja en un periodo determinado un cierto nivel de energía.
2. Un porcentaje de este nivel de energía se transmite de ellos a sus contactos en función de ciertos criterios (a mayor nivel de contacto, mayor flujo de energía). Este criterio preserva la energía: la energía total del sistema antes y después de la redistribución es la misma.
3. Los contactos que tienen un nivel de energía mayor que cero lo transmiten recursivamente a los suyos.
4. El proceso se itera hasta que alcanza un equilibrio razonable.

Al final, a muchos clientes (técnicamente, a los que pertenecen a la unión de las componentes conexas que contienen a las bajas) se les habrá asociado un nivel de energía. Y este nivel de energía es, según los autores, un indicador temprano de baja de alto valor predictivo.

¿Será?

El otro día me llegó por correo el Informe sobre el Uso del Software Libre en los Hogares Españoles 2011. Lo realiza el CENATIC, Centro Nacional de Referencia de Aplicación de las Tecnologías de Información y la Comunicación basadas en Fuentes Abiertas, por lo que uno espera, de antemano, cierto sesgo.

Una de las tablas de resultados es:

Entiendo que los porcentajes de uso se refieren al universo de la población española, extrapolados mediante un [...] muestreo por cuotas, donde se incluyen cuotas con afijación proporcional al peso real de la población objeto, obteniendo estos datos del Instituto Nacional de Estadística, en el período más actualizado.

Pero, ¿veis el porcentaje de uso de Linux? ¿Cómo puede ser que difiera tanto de estimaciones como esta o esta?

Puede que la desviación se deba al modo en que se ha realizado la encuesta:

Para la captación de la muestra, se han enviado por correo electrónico más de 10.000 invitaciones de participación en la encuesta, reforzadas con cinco recordatorios (mediante correo electrónico y teléfono) diferentes para poder maximizar la tasa de respuesta.

O quién sabe a qué. Pero en tanto no nos expliquen de manera creíble unos porcentajes de penetración tan distintos de los que aparecen publicados por doquier, mejor tomar cum grano salis los resultados del estudio.

Arranca Tolstoi en Ana Karenina con esta frase mítica: Todas las familias felices se parecen entre sí; las infelices son desgraciadas en su propia manera.

Me recuerda mucho a la contraposición entre la probabilidad que estudié en la universidad y la que regía fuera. Dentro, mis variables eran, casi indefectivamente, iid, es decir, independientes e idénticamente distribuidas. Y las variables independientes son muy parecidas entre sí. Incluso más, diría yo, que las familias felices.

Sin embargo, extramuros, no hallé independencia. Y las variables aleatorias que no son independientes, dejan de serlo, cada una, en su propia manera. Hipotecas, por ejemplo. Probabilidad de que ambas entren en mora. Una de Murcia, otra de Cádiz. Bernoulli, por supuesto. Pero, ¿igualmente distribuidas? Y, además, ¿independientes? Por supuesto que no. Pero, más aún, es imposible conocer la relación de dependencia (y por lo tanto, aproximarla: ¿aproximarla a qué?). Y es que aunque fuesen condicionalmente independientes con respecto a algunas magnitudes macroeconómicas, ni siquiera sabríamos cuáles. Ni de qué manera.

Ya, ya sé que algunos estudian una cosa que se llama AR(1) donde las observaciones no son independientes. Pero tienen un mecanismo de dependencia simple, explícito e igualmente alejado de prácticamente toda realidad que no sea la simulada por uno mismo en condiciones de laboratorio. Suponen, además, un avance minúsculo desde el puerto seguro de la independencia hacia el océano inmenso e insondable de la falta de ella: solo es uno de los infinitos, incontables, mecanismos que introducen dependencia, cada cual en su propia manera, entre variables. ¡Y eso que ni siquiera es posible saber de manera indubitable que una serie real es AR(1)!

Esto lo sé. De mercados eficientes e ineficientes, particularmente de los segundos, sé infinitamente menos. Pero cuando leo a Jesús Fernández Villaverde escribir:

[...] con repeticiones de banalidades de Krugman [...] como la que los economistas estamos enamorados de modelos donde los mercados son perfectos. Esto de nuevo son noticias para mi [sic] ya que yo tenía la impresión de que casi todos mis papers son modelos con imperfecciones de mercado. Bueno debe de ser que mis papers solo existen en un proyecto existencial distinto. Pero quien lo deberá llevar peor son mis estudiantes de doctorado, a los que acabo de meter 15 horas de clase bien intensas (y 3 más que les quedan) de modelos con imperfecciones y que se deberán estar preguntado los pobres que en qué dimensión existencial se han ido esas horas.

Desde mi ignorancia —pero apoyándome en la analogía—, sospecho que la relación entre los modelos con imperfecciones de mercado de JFV y los que los suponen eficientes deben de parecerse muchísimo a los existentes entre los modelos de dependencia de las series autorregresivas y el ideal universo iid. Pero apostaría a que luego, los mercados ineficientes —igual que las familias infelices y las variables aleatorias no independientes— lo son, también, cada uno de ellos, a su propia manera.

Y no sé si JFV ha zanjado definitivamente viejos problemas epistemológicos (sobre el problemático encaje entre los constructos teorico-matemáticos y la siempre rebelde realidad) o si le falta una dosis de sana humildad socrática. Pero si tuviese que apostar…

Hoy, por motivos evidentes, e igual que en septiembre de 2010, voy a hablar de huelgas. De la misma fuente que entonces he descargado este fichero. Y he ejecutado

library( pxR )
library( reshape )
library( ggplot2 )

dat <- read.px( "pcaxis-623612450.px" )
dat <- as.data.frame( dat )

dat.mes <- cast( dat, Periodo ~ series )
colnames(dat.mes) <- c( "mes", "n.huelgas", "n.trabajadores", "n.jornadas" )

p <- ggplot( data = dat.mes ) + geom_line( aes( x = mes, y = n.huelgas, group = rep(1, nrow(dat)) ) )
p
ggsave( "huelgas_por_mes.png" )

dat.anno <- dat

tmp <- strsplit( as.character(dat.anno$Periodo), "M" )
dat.anno$Periodo <- sapply( tmp, function(x) x[1] )

dat.anno <- cast( dat.anno, Periodo ~ series, fun.aggregate = sum )
colnames(dat.anno) <- c( "anno", "n.huelgas", "n.trabajadores", "n.jornadas" )

p <- ggplot( data = dat.anno, aes( x = anno, y = n.huelgas, group = rep(1, nrow(dat.anno)) ) ) + geom_line()
p <- p + geom_point( aes(size = n.jornadas ) )
p <- p + scale_x_discrete( "año" ) + scale_y_continuous( "número de huelgas" )
p
ggsave( "huelgas_por_anno.png" )

p <- ggplot( data = dat.anno, aes( x = anno, y = n.trabajadores/n.huelgas, group = rep(1, nrow(dat.anno)) ) ) + geom_line()
p <- p + scale_x_discrete( "año" ) + scale_y_continuous( "número de trabajadores por huelga" )
p
ggsave( "trabajadores_huelga_por_anno.png" )

p <- ggplot( data = dat.anno, aes( x = anno, y = n.jornadas /n.huelgas, group = rep(1, nrow(dat.anno)) ) ) + geom_line()
p <- p + scale_x_discrete( "año" ) + scale_y_continuous( "número de jornadas por huelga" )
p
ggsave( "jornadas_huelga_anno.png" )

para obtener, por un lado, el número de huelgas por mes desde enero de 1995 a noviembre de 2011:

El número de huelgas por año, con el número de jornadas afectadas por las huelgas indicado por el diámetro de punto:

El número de trabajadores que siguieron cada huelga en cada año:

Finalmente, el número de jornadas afectadas por las huelgas en cada año.

Es relativamente fácil ubicar la huelga general del 2002. La que hubo en septiembre de 2010… no lo tengo tan claro.

Finalmente, ¿podría decirse que aunque parece aumentar el número de huelgas, estas afectan cada vez a menos trabajadores y restan mejos jornadas de trabajo? ¿Se ve esa tendencia? ¿Es una tendencia o una mera ilusión sensorial?

De ser tendencia, (y sólo de ser tendencia) ¿podría aventurarse alguna causa? ¿Huelgas más locales? ¿O menos interés y seguimiento por parte de los trabajadores?

Finalmente, ¿nos satisface la información sobre huelgas que publica el INE? ¿Existen fuentes más completas? ¿Alguien las tiene ubicadas? ¿Son accesibles? ¿Permitiría la nueva ley de transparencia acceder a ellas?

Las leyes son un carajal. Últimamente he tenido que enfrentarme a algunas y me doy cuenta de que es un error que los abogados (y los legisladores) no hayan seguido nunca un buen curso de geometría euclídea.

Pongo un ejemplo. La Constitución Española (artículo 22) dice que las asociaciones deberán inscribirse en un registro a los solos efectos de publicidad. Pues bien, parece ser que una asociación, y nos ocurrió con la Comunidad R Hispano no es legal (para poder abrir una cuenta corriente, por ejemplo) de no inscribirse en el registro. ¿Pero no era a los solos efectos de publicidad? ¿Son lo mismo publicidad y legalidad?

En fin.

Hoy me toca volver a hablar de leyes y de cómo afectan a la labor estadística. En particular, a los estadísicos que nos movemos (en términos amplios) por libre los que usamos software libre y a los que nos encantan los datos libres. A los que nos gusta poder enriquecernos de las aportaciones de otros libremente y ponemos también libremente al servicio de la comunidad las propias.

Desde el punto de vista del software, disponemos de R. Podemos acceder a él y podemos reutilizarlo, copiarlo, cederlo, regalarlo. También cedemos, distribuimos y permitimos que otros reutilicen nuestras propias aportaciones. Acceso y reutilización son dos caras de la misma moneda.

¿Y respecto a los datos? Y en particular, ¿de los datos públicos? En España, está regulada por ley su redistribución (otra cosa es que los organismos titulares de dichos datos, como el INE, no estén al corriente de ella). Y ¿el acceso (que se supone que debería ser, en pura lógica, previo)? Pues hubo un borrador, del que hablé en otra ocasión, y que los avatares electorales convirtieron en papel mojado.

Pero ahora el gobierno ha sacado definitivamente a la luz (desde las cinco de la tarde de ayer) un nuevo borrador de la Ley de Transparencia con la particularidad adicional de que ha puesto a disposición de la ciudadanía un portal para que esta pueda realizar alegaciones y propuestas de mejora.

Por eso, desde estas páginas, quiero invitar a aquellos interesados en los datos públicos y en garantizar el acceso a los datos públicos tanto para él como para las generaciones posteriores, que lean el borrador, evalúen en qué medida salvaguardan su derecho de acceso a esos ficheros csv tan jugosos estadísticamente y que, de observar cualquier indicio de mejora, hablen ahora o callen para siempre (y no den el coñazo, con perdón, después).

Nota: el día 26 acudí al acto de inauguración del portal tuderechoasaber.es que facilita a los ciudadanos la solicitud de información al amparo de esta ley en ciernes. Invito a los lectores de esta bitácora a visitarlo y aprovecharlo y, a la vez, quiero felicitar (y a la vez agradecer) a David Cabo (@dcabo) y el resto del equipo que ha puesto en marcha dicho proyecto por su capacidad de iniciativa, trabajo y vocación de servicio público.

En 2013 celebraremos el Año Internacional de la Estadística. Diversas organizaciones estadísticas, entre las que no veo al INE, se han unido para dar a conocer “la importancia de la estadística en la comunidad científica, el mundo de los negocios, la administración pública, los medios de comunicación, las empresas, los estudiantes y el público en general”.

Pueden consultarse las actividades programadas y la lista de organizaciones patrocinadoras, a la que es posible que es sume la Comunidad R Hispano.

Mucha gente se hace un lío con los errores de tipo I y II. Para ellos, esta regla nemotécnica:

• Tipo I, un false: I falsely think hypothesis is true
• Tipo II, dos falses: I falsely think hypothesis is false