Estadística

Un recordatorio: MOMOCalor está "up and running"

Por desgracia, MoMo ya no exige presentación. Pero con los termómetros acariciando los 40º no está mal recordar la existencia de MoMoCalor, su hermanito, que trata atribuir mortalidad a los excesos de temperaturas.

¿Por qué es particularmente importante MoMoCalor hoy? Recuérdese que MoMo estima, simplemente, desviaciones de mortalidad con respecto a la que sería la normal en una fecha determinada. Cuando hay una epidemia o una ola de calor, la mortalidad crece y MoMo lo detecta. Pero cuando hay una epidemia y una ola de calor simultáneas, MoMo es incapaz de atribuir muertos las causas anómalas subyacentes. Pero MoMoCalor sí.

Por supuesto que tengo más variables que observaciones... ¿y?

He intentado replicar los resultados de la entrada de ayer con GAM (vía mgcv) así (véase el enlace anterior para la definición de los datos):

library(mgcv)
modelo_gam <- gam(
    y ~ x + s(id, bs = "re"),
    data = datos,
    method = "REML",
    family = "poisson")

Y nada:

Error in gam(y ~ x + s(id, bs = "re"), data = datos, method = "REML", : Model has more coefficients than data

Sí, ya sé que tengo más variables que observaciones. Pero, ¿no es para eso que estoy usando efectos aleatorios?

Aún más sobre la presunta sobredispersión en modelos de Poisson

[Esta entrada continúa el ciclo al que he dedicado esta y esta otra entradas durante los últimos días.]

Las dos entradas anteriores de la serie se resumen en que:

  • el modelo de Poisson no recoge todas las fuentes de error que pueden existir en los datos y que
  • las soluciones al uso (como, p.e., usar modelos quasi-Poisson) son puros remiendos.

Si el error en el modelo de Poisson entra (también) en el término lineal, podemos modelar ese error explícitamente. Podría haber implementado la solución INLA o Stan del problema, pero me conformaré con la lme4. Primero, generaré los datos (igual que en las entradas anteriores) y añadiré una variable categórica que identifique cada registro:

Análisis de arquetipos

De eso trata un artículo de los noventa de Breiman. Es decir, de encontrar dentro de conjuntos de datos conjuntos finitos de sujetos puros que permiten representar cualquier otro como una mezcla (o combinación convexa) de ellos.

Ideas a vuelapluma:

  • Cuando leo sobre el asunto, la palabra que no deja de aparecérseme es outlier. Curiosamente, la busco en el texto y se resiste a aparecer. Pero me aterra la posibilidad de estar caracterizando a los sujetos normales (¿aún se puede usar la expresión?) como combinación convexa de raritos.
  • La técnica podía competir muy favorablemente con el clústering tanto conceptualmente (resuelve el problema de la heterogeneidad de los clústers) como operativamente (se podrían extraer para algún fin los sujetos que participasen en una proporción determinada de un cierto arquetipo).
  • En el fondo, se solapa con otras técnicas bien establecidas y que hacen cosas parecidas como LDA (con D de Dirichlet) o NMF (factorización no negativa de matrices).

Más sobre la presunta sobredispersión en el modelo de Poisson

[Esta entrada abunda sobre la de ayer y sin la cual no se entiende.]

Generemos unos datos, las x:

n <- 1000
sigma <- .5
x <- rep(-2:2, each = n)
x_real <- -1 + .5 * x + rnorm(length(x), 0, sigma)

En el bloque anterior hemos creado una/la variable observada, x, el término lineal que operará en el modelo de Poisson, -1 + .5 * x, y el real, -1 + .5 * x + rnorm(length(x), 0, sigma), que agrega al anterior el impacto de otras variables no tenidas en cuenta a través de un error normal al uso.

No, tus datos no "tienen sobredispersión": es que el gato de Nelder se ha merendado la epsilon

El modelo de Poisson viene a decir que si y es una variable con valores 0, 1,… y x1,…, xn son variables explicativas tiene cierto sentido en algunos casos plantear un modelo de la forma

$$ y | x_i \sim \text{Pois}(\exp(a_0 + \sum_i a_i x_i) ),$$

Es decir , para cada combinación de las xi, el modelo proporciona el parámetro de una distribución de Poisson de la que y es una realización. Hay una incertidumbre (o un error irreductible) que reside en que de y solo conocemos la distribución.

Sobre el efecto medio

Traduzco de aquí:

En estadística y econometría se habla a menudo del efecto medio de un tratamiento. A menudo, he sido [Gelman] escéptico con respecto al efecto medio por la sencilla razón de que, si se trata de un efecto medio, se está reconociendo la posibilidad de variación; y si hay una variación importante (tanto como para hablar del efecto medio y no solo del efecto) es que nos preocupa tanto que deberíamos estudiarla directamente en lugar de reducirla a su promedio.

Sobre la curva ROC como medida de bondad de clasificadores

Esta entrada se entiende mal sin esta otra donde se daba noticia de un clasificador que era mucho mejor o peor (de acuerdo con ciertas métricas) según la tasa de prevalencia de la clase relevante a pesar de que tanto su sensibilidad como su especificidad no eran particularmente malas. Efectivamente, con lo del coronavirus hemos reaprendido a darle la vuelta a las probabilidades condicionales y aplicar el teorema de Bayes para ver qué cabía esperar de un clasificador cuyas bondades se predican en términos de la sensibilidad y la especificidad.

¿Qué queda de la "estadística robusta" clásica?

Estos días estoy muy atento a todo lo que tiene que ver con estadística robusta. El motivo es doble:

  • Estoy involucrado en un proyecto donde quieren ajustar ciertos modelos usando funciones de pérdida robustas (Huber, Tukey, etc.).
  • Hay una $latex 1 > p > 0$ de que me toque meter mano a MOMO y sus derivados para que lo del coronavirus no joda los contrafactuales de 2021 y sucesivos (¿bastará con eliminar unos cuantos meses de 2020?).

Así las cosas, ha aterrizado en mi tableta The Changing History of Robustness, donde, el autor, Stigler:

Un artículo muy raro, raro, raro

Hoy voy a comentar un artículo muy raro que me ha llegado recientemente y que se titula nada menos que Bayesian Estimation with Informative Priors is Indistinguishable from Data Falsification.

Argumenta el artículo alrededor de lo siguiente (que creo que ya sabemos todos: son ejercicios matemáticos básicos de un curso introductorio de probabilidad):

  • Que la inferencia bayesiana con prioris planas (degeneradas, de ser necesario) es equivalente a la inferencia frecuentista.
  • Que para tres ejemplos concretos (binomial, Poisson y normal), de usarse prioris a través de las distribuciones conjugadas, el resultado de la inferencia bayesiana es equivalente a haber añadido datos (problemas de redondeo aparte) a los originales.

Luego añade unos experimentos numéricos para dejar constancia de que no se ha equivocado en las cuentas y mostrar que, efectivamente, sustituyendo las letras por números y operando se obtienen los resultados que anuncian las matemáticas con su árido simbolismo.