Esta entrada responde y complementa Malditas proporciones pequeñas I y II_ _trayendo a colación un artículo que ya mencioné en su día y que cuelgo de nuevo: On the Near Impossibility of Measuring the Returns to Advertising. ¡Atención al teorema de la imposibilidad de la Super Bowl!

Y el resumen breve: cada vez estamos abocados a medir efectos más y más pequeños. La fruta que cuelga a la altura de la mano ya está en la fragoneta del rumano. Solo nos queda la morralla y cada vez va a costar más separar grano y paja.

Tienes dos variables aleatorias positivamente correlacionadas, $latex X$ y $latex Y$ y una muestra de $latex n$ parejas de ellas $latex (x_i, y_i)$.

La esperanza de $latex X$, $latex E(X)$, es conocida y la de $latex Y$ no. Obviamente, la puedes estimar haciendo

$$ E(Y) \sim \frac{1}{n} \sum_i y_i.$$

Sin embargo, la varianza del estimador

$$ E(Y) \sim E(X) \frac{\sum y_i}{\sum x_i}$$

es menor.

Tengo una explicación de la intuición de por qué eso es cierto en lugar de no serlo. Pero como no sé si es suficientemente buena, dejo que alguien proponga la suya en los comentarios.

Era como

y se ha convertido en

¡Qué horror!

Coda: En otra página de la Wikipedia en la que he caído después por azar he leído la siguiente frase (que por algún motivo encuentro relevante insertar aquí):

Los ríos arrastran sedimentos que consiguen colmatar y rellenar de lodo los lagos. Además, la proliferación de ciertas plantas, como el lirio acuático, los obstruye por completo.

Si ejecutas

import numpy as np
from sklearn.linear_model import LinearRegression

n = 1000
X = np.random.rand(n, 2)

Y = np.dot(X, np.array([1, 2])) + 1 + np.random.randn(n) / 2
reg = LinearRegression().fit(X, Y)

reg.intercept_
reg.coef_

se obtiene más o menos lo esperado. Pero si añades una columna linealmente dependiente,

X = np.column_stack((X, 1 * X[:,1]))

ocurren cosas de la más calamitosa especie:

Y = np.dot(X, np.array([1, 2, 1])) + 1 + np.random.randn(n) / 2
reg = LinearRegression().fit(X, Y)
reg.coef_
# array([ 9.89633991e-01, -1.63740303e+14,  1.63740303e+14])

Comentarios:

• Diríase que la implementación del modelo lineal en scikit-learn no es la que se estudia por doquier (la prima, la inversa, etc.); sospecho que convierte el error cuadrático en una función que depende de los coeficientes y se la pasan a un optimizador (más o menos) genérico.
• Supongo que la implementación actual pasa todos las pruebas unitarias.
• Y sospecho, además, que las pruebas unitarias no las ha planteado un estadístico.
• Así que tal vez los de scikit-learn no saben que tienen problemas de colinealidad; y si alguien se lo ha comentado, igual no han comprendido el issue.
• Para la predicción, el problema arriba apuntado no es tal. Aun con coeficientes desaforados y siempre que no haya problemas de precisión numérica, tanto da hacer las cosas como todo el mundo o implementando ocurrencias como la anterior.
• Pero para todo lo demás (p.e., impacto de variables, etc.), la implementación es de traca y no vale para maldita de Dios la cosa.
• Aunque a estas alturas de siglo, ¿quién en su sano juicio usa el modelo lineal básico?
• Además, en la práctica, no hay problemas de colinealidad (ni aproximada o, mucho menos, exacta). Hummm…
• ¿O sí? Mi colega Cañadas ha escrito una entrada en su blog sobre la codificación de variables categóricas donde denuncia cómo en Python las funciones más habituales crean por defecto columnas linealmente dependientes por defecto (por no omitir el primer nivel). O sea, en Python, si no andas con cuidado, acabas con la suela llena de kk de perro.

Voy a abundar sobre la entrada de hace unos días, ¿Informática o matemáticas?, una pregunta muy mal planteada, mostrando simplemente un ejemplo del tipo de cosas que se espera de los matemáticos y/o estadísticos cuando trabajan en ciencia de datos y para las cuales los informáticos no están particularmente mejor entrenados (de serie) que otras especies faunísticas.

Es este.

¿Cosas sobre las que podría hacer comentarios? Por ejemplo:

• Tampoco sé si el matemático o estadístico promedio podría desenvolverse con mediana soltura con ese tipo de modelos. Y sí, cuando la sal se vuelve sosa, no es de extrañar que la tiren fuera y que la pise la gente.
• Ese tipo de modelos no se usan y no porque no sean [más] adecuados [que otros]. No se usan, principalmente, por motivos que mi colega José Luis Cañadas expone en párrafos de su blog que suelen contener la palabra ingenazi.

He construido

que, obviamente no es la gran maravilla, basándome en Rectangular Statistical Cartograms in R: The recmap Package y usando

library(rgdal)
library(pxR)
library(recmap)

provs <- readOGR(dsn = "provincias/",
    layer = "Provincias")

pobl <- as.data.frame(read.px("2852.px",
    encoding = "latin1"), use.codes = T)
pobl2 <- as.data.frame(read.px("2852.px",
    encoding = "latin1"))

pobl$nombre <- pobl2$Provincias

pobl <- pobl[, c("Provincias", "nombre", "value")]
colnames(pobl) <- c("COD_PROV", "nombre", "poblacion")
pobl <- pobl[pobl$COD_PROV != "null",]

pobl <- pobl[!pobl$COD_PROV %in%
    c("51", "52", "38", "07", "35"),]


dat <- merge(provs, pobl,
    by = "COD_PROV", all.x = FALSE)
dat@data$NOM_PROV <- NULL
dat$z <- dat$poblacion

tmp <- as.recmap(dat)

tmp$name <- dat@data$nombre
tmp$ccaa <- dat@data$COD_CCAA

res <- recmapGA(tmp, popSize = 300,
    maxiter = 30, run = 10)

cartogram <- res$Cartogram

ccaa <- tmp[, c("name", "ccaa")]
ccaa$ccaa <- as.numeric(factor(ccaa$ccaa))
cartogram <- merge(cartogram, ccaa)

plot.recmap(cartogram, col.text = "black",
    main = "cartograma -- población\n  españa peninsular",
    col = cartogram$ccaa)

Como los datos los he bajado de por ahí y no recuerdo dónde, dejo como referencia el objeto arriba llamado tmp aquí.

Continúo con esto.

(El pico de primeros de agosto corresponde a la ola de calor).

Resumen:

• Acusada variación intrasemanal (con un intevalo de variación del 5%, que es mucho).
• Los intervalos de confianza de los nuevos modelos contienen mucho mejor las observaciones reales; muchos picos observados dejan de parecer anómalos para quedar dentro de las franjas de variación esperadas.

Dos noticias para la entrada de hoy:

• En Economipedia han tenido a bien entrevistarme.
• El jueves y viernes de la semana que viene impartiré sendas conferencias en Alicante dentro del ciclo Matemáticos en la Sociedad. Más información, aquí.

Es sabido que nacen menos niños en domingo, efecto, parece de la planificación de partos. Tengo cierto indicio, que voy a ver cuándo (y si) puedo corroborar, de que también hay menos fallecimientos. Mírese

que es una gráfica de mortalidad diaria en España y en la que, contumazmente, los días 2 de junio (domingo), 9 de junio (domingo) y 16 de junio (domingo), las observaciones quieren salirse de las bandas que para que no trotaran libérrimamente por el plano cartesiano construyó el bueno de Monsieur Poisson.

En vigilancia epidemiológica contamos eventos (p.e., muertes o casos de determinadas enfermedades). Lo que pasa es que el caso ocurrido en el día 0 puede notificarse con un retraso de 1, 2, 3… o incluso más días. En algunas aplicaciones, incluso semanas.

¿Cómo estimar el número de casos ocurridos el día 0 el día, p.e., 5?

Se puede aplicar el análisis de la supervivencia donde el evento muerte se reinterpreta como notificación. El el día 0 todos los sujetos están vivos y, poco a poco, van cayendo. Como en los consabidos modelos/gráficos de Kaplan-Meier,