Estadística

Una vez escribí al respecto. Y cuanto más lo repienso y lo reeleo, menos clara tengo mi interpretación. De hecho, estoy planteándome retractar esa entrada.

Y reconozco que llevo tiempo buscando en ratos libres algún artículo serio (no extraído del recetario de algún script kiddie de Kaggle) que justifique el uso del procedimiento. Es decir, que lo eleve de técnica a categoría. Sin éxito.

He hecho probaturas y experimentos mentales en casos extremos (p.e., cuando todos los niveles de la variable categórica son distintos, cuando son iguales, etc.) con los decepcionantes resultados que cabe esperar. Lo cual contradice las presuntas virtudes casi taumatúrgicas del procedimiento.

David Hand, en Classifier Technology and the Illusion of Progress, resume el asunto así:

A great many tools have been developed for supervised classification, ranging from early methods such as linear discriminant analysis through to modern developments such as neural networks and support vector machines. A large number of comparative studies have been conducted in attempts to establish the relative superiority of these methods. This paper argues that these comparisons often fail to take into account important aspects of real problems, so that the apparent superiority of more sophisticated methods may be something of an illusion. In particular, simple methods typically yield performance almost as good as more sophisticated methods, to the extent that the difference in performance may be swamped by other sources of uncertainty that generally are not considered in the classical supervised classification paradigm.

Primero, una simulación:

n <- 100
delta <- 0.2
n.iter <- 10000

p_valores <- function(n, delta){
  tmp <- replicate(n.iter, {
    x <- rnorm(n)
    y <- rnorm(n, mean = delta)
    t.test(x, y)$p.value
  })

  res <- tmp[tmp < 0.05]

  hist(res, freq = FALSE, xlab = "p value", ylab = "", col = "gray", main = "histograma de p-valores publicables")

  res
}

null_effect_p_values <- p_valores(n, 0)
some_effect_p_values <- p_valores(n, delta)

Lo que simula son n.iter experimentos en los que se comparan n valores N(0,1) con otros n valores N(delta, 1) y se extrae el correspondiente p-valor. Luego se grafican los publicables (<0.05).

Alguno de mis lectores, supongo, estará metido en ese mundo de ir escribiendo cosas y cosechado méritos, impactos y anecosas para salir del precariado y pillar moscosos. Que dejen de leer. Es una orden.

A aquellos que tengan tiempo y talento los invito a escribir el artículo titulado Temperaturas umbrales de disparo de la mortalidad atribuible al frío y al calor en España en el periodo 2007-2017.

Se trata, esencialmente, de aggiornar metodológica, gráfica y sintácticamente esta cosa viejuna y manifiestamente mejorable en todas las dimensiones concebibles.

El objeto único de la estadística es informar decisiones.

V.g, si conceder un préstamo, proceder a la quimio, construir una línea de AVE entre Calatayud y Soria o permitir aparcar mañana en el centro de Madrid a los de Móstoles.

Pero quienes toman decisiones y quienes analizan datos suelen ser personas distintas. Típicamente, ni se conocen. Lo cual es tanto pésimo como tema para otra entrada distinta de esta.

Lo fundamental es que estas personas se comunican a través de, metafóricamente, APIs. Unas de las más usadas son los p-valores. Que son tan pésismos como tema para otra entrada distinta de esta.

He usado el vídeo

en un curso de estadística básica para ilustrar a través de experimentos se construyen histogramas y estos convergen a y, en última instancia, justifican el uso de distribuciones de probabilidad.

Es decir,

experimentos -> histogramas -> funciones de distribución.

Y de ahí, el resto.

Envalentonado por el comentario de Iñaki Úcar a mi entrada del otro día, que me remitía a este artículo, decidí rizar el rizo y crear intervalos de confianza no ya discontinuos sino con otra propiedad topológica imposible: homeomorfos con un toro.

Y aquí está:

El modelo, el código y demás,

library(rstan)
library(ggplot2)

n <- 100

a1 <- 1
a2 <- 1
sigma <- 0.4

datos <- data.frame(x1 = rnorm(n, 2, 0.1),
                    x2 = rnorm(n, 2, 0.1))

datos$y <- a1^datos$x1 + a2^datos$x2 + rnorm(n, 0, sigma)

codigo <- "
data {
  int<lower=1> N;
  real y[N];
  real x1[N];
  real x2[N];
}

parameters {
  real<lower=-3, upper="3"> a1;
  real<lower=-3, upper="3"> a2;
  real<lower=0, upper="3"> sigma;
}

model {
  for (n in 1:N)
    y[n] ~ normal(fabs(a1)^x1[n] +
      fabs(a2)^x2[n], sigma);
}"

fit <- stan(model_code = codigo,
    data = list(N = length(datos$y), y = datos$y,
                x1 = datos$x1, x2 = datos$x2),
    iter=40000, warmup=2000,
    chains=1, thin=10)

res <- as.data.frame(fit)

ggplot(res, aes(x = a1, y = a2)) + geom_point(alpha = 0.1)

De nuevo, no son intervalos propiamente dichos, lo convengo. Pero son configuraciones más fieles al espíritu de lo que un intervalo de confianza es y representa que su(s) letra(s) I N T E R V A L O.

Me sorprende haber leído tantos de los mejores 100 libros del siglo XX según Le Monde. Sobre todo porque no leo ficción casi en lo que va de siglo y porque, carajo, los libros estupendos que he leído de tapa, como el Análisis Real de Folland o la Introducción a la Teoría de la Probabilidad de Feller parece que no cualifican para esa listeja de textos sin una mala integral preparada por gentecilla de letras.

Es el de b:

(A ver cuál es el primero de mis excolegas que protesta que pinto la unión de dos intervalos de confianza y no un intervalo propiamente dicho).

Ahora un poco más en serio: esta entrada se me ocurrió mientras pensaba en las distintas opciones existentes para crear intervalos de confianza, desde las canónicas (simétricos, de longitud mínima) a cualquier otra elección de algo que contenga la debida cantidad de probabilidad.

El gráfico

ha estado dando vueltas por el ciberespacio. Lo vi en Twitter de mano de alguien que lo usaba para justificar que la distribución de la renta no es tan desigual en España al fin y al cabo. Está comentado desde el punto de vista de la interpretación y tufneado en términos de la forma

aquí.

Pero lo que no he visto comentar es que las variaciones reflejan más cómo es el tamaño de las provincias (o regiones, estados, o las divisiones administrativas que se haya considerado) en cada uno de los países que si la renta está mejor o peor repartida.