Estadística

En defensa de Simón: variaciones diarias de la mortalidad

Qué cafres tenéis que ser para que tenga que salir yo —precisamente yo, que tantas cosas no buenas tengo para decir del buen hombre— en defensa de Simón. Tiene delito que de todo lo que se le pueda echar en cara os hayáis fijado en una intervención en la que os trataba de desasnar para que no le anduviéseis buscando tres pies a la varianza.

Es un tema que vengo tratando de antiguo en estas páginas y de ello dan fe:

Consensus clustering

No hay nada tan corrosivo para la fe en el clústering que probar una y otra vez k-medias (por ejemplo) sobre los mismos datos y ver cómo los resultados cambian drásticamente de ejecución en ejecución.

Pero eso viene a ser, esencialmente, lo que hay detrás del consensus clústering (CC), una técnica que puede ser usada, entre otros fines, para determinar el número óptimo de grupos.

La idea fundamental de la cosa es que observaciones que merezcan ser agrupadas juntas lo serán muy frecuentemente aunque cambien ligeramente las condiciones iniciales (por ejemplo, se tome una submuestra de los datos o cambien las condiciones iniciales de k-medias, por ejemplo). Si uno altera esas condiciones iniciales repetidas veces puede contar la proporción de las veces que las observaciones i y j fueron emparejadas juntas y crear la correspondiente matriz (simétrica, para más señas) $latex C(i,j)$.

lme4 + simulate

Esta entrada es casi una referencia para mí. Cada vez tiro más de lme4 en mis modelos y en uno en concreto que tengo entre manos toca simular escenarios. Para lo cual, simulate.merMod.

Véamoslo en funcionamiento. Primero, datos (ANOVA-style) y el modelo que piden a gritos:

library(plyr)
library(lme4)

a <- c(0,0,0, -1, -1, 1, 1, -2, 2)
factors <- letters[1:length(a)]

datos <- ldply(1:100, function(i){
    data.frame(x = factors, y = a + rnorm(length(a)))
})
modelo <- lmer(y ~ (1 | x), data = datos)

El resumen del modelo está niquelado:

summary(modelo)

# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ (1 | x)
# Data: datos
#
# REML criterion at convergence: 2560.3
#
# Scaled residuals:
#     Min      1Q  Median      3Q     Max
# -3.6798 -0.6442 -0.0288  0.6446  3.3582
#
# Random effects:
#     Groups   Name        Variance Std.Dev.
# x        (Intercept) 1.5197   1.2328
# Residual             0.9582   0.9789
# Number of obs: 900, groups:  x, 9
#
# Fixed effects:
#     Estimate Std. Error t value
# (Intercept) -0.009334   0.412212  -0.023

En particular,

Intervalos de confianza, intervalos de predicción

Contexto:

modelo <- lm(dist ~ speed, data = cars)

Intervalos de confianza:

head(predict(modelo, interval = "confidence"))
#        fit        lwr       upr
#1 -1.849460 -12.329543  8.630624
#2 -1.849460 -12.329543  8.630624
#3  9.947766   1.678977 18.216556
#4  9.947766   1.678977 18.216556
#5 13.880175   6.307527 21.452823
#6 17.812584  10.905120 24.720047

Intervalos de predicción:

head(predict(modelo, interval = "prediction"))
#        fit       lwr      upr
#1 -1.849460 -34.49984 30.80092
#2 -1.849460 -34.49984 30.80092
#3  9.947766 -22.06142 41.95696
#4  9.947766 -22.06142 41.95696
#5 13.880175 -17.95629 45.71664
#6 17.812584 -13.87225 49.49741

Creo que la diferencia (y el significado) es claro. Para todos los demás, esto.

Sobre los peligros del "Tukey biweight"

Sigo con ajustes robustos. Y cosas que como matemático, me ponen muy nervioso.

Una de las maneras de hacer ajustes robustos es la de sustituir la función cuadrática por la biweight. Es decir, utilizar la función que aparece la derecha en

en lugar de la de la izquierda. O, dicho de otra manera, en lugar de tratar de minimizar

$$ \sum_i \rho(y_i - f_\alpha(x_i))$$

usando $latex \rho(x) = x^2$, que es la función que se representa a la izquierda y a la que estamos acostumbrados, usar la de la derecha. Que es la función biweight de Tukey.

La probabilidad de que el parámetro esté en el intervalo de confianza es .95

Si dices lo anterior, corres el riesgo de que un estadístico gruñón frunza mucho el ceño.

Hace muchos, muchos años, las gentes ávidas de saber más acudieron al tabernáculo donde se congregaban los sapientísimos estadísticos frecuentistas implorándoles una herramienta con que estimar el error de sus estimaciones puntuales. Estos cavilaron luengamente y décadas después entregaron a los representantes de los hombres, reunidos en el ágora, unas tablas de piedra que tenían grabadas a cincel la teoría de los intervalos de confianza. Pero, les advirtieron, los intervalos de confianza no son lo que vosotros queréis sino otra cosa y a quien ose interpretarlos torcidamente le pasará lo que a aquella señora que comió la manzana inadecuada: será expulsado del paraíso de la teoría como Dios manda.

WoE,... pero ¿y las interacciones?

Esto del WoE he tenido que aplicarlo (de manera no estándar, además) en alguna ocasión. Pero forzado por las circunstancias (que, concretamente, eran el misteriosísimo y no siempre conforme a lo que cabría esperar que hace ranger de las variables categóricas). Digamos que a veces toca, pero no es tampoco algo de lo que enorgullecerse.

Pero cuando escucho o leo a los apologetas del WoE, siempre me pregunto mucho por lo que tendrán que decir sobre la pérdida de información en términos abstractos y, en otros más concretos, qué ocurre con las interacciones.

Comparación y selección de modelos bayesianos

En el mundo bayesiano existen, cuando menos, dos escuelas:

  • La flowerpower, que sostiene que los modelos bayesianos son subjetivos y, por lo tanto, inasequibles a la confrontación con la realidad objetiva.
  • La de los que tienen un jefe que les paga un salario, al que le da igual si los modelos son bayesianos o no pero a quien le interesa por encima de todo saber si representan razonablemente el proceso subyacente.

Los segundos cuentan con referencias como Comparison of Bayesian predictive methods for model selection. Es un artículo, en cierto modo, desasosegadoramente antibayesiano: miradlo y encontraréis en él cosas que se parecen demasiado a la validación cruzada, al RMSE, etc.

GoF para modelos bayesianos

Existe una muy perezosa escuela de pensamiento que sostiene que dado que las probabilidades son subjetivas, cualquier modelo y, en particular, los bayesianos, como expresión de la subjetividad de sus autores, no necesita ser contrastado con la realidad. Porque, de hecho, la realidad no existe y es una construcción que cada cual hace a su manera, deberían añadir.

Existe, por supuesto, una escuela realista tan mayoritaria que ni siquiera es consciente de que lo es. Basta leer la primera página de Statistical Modeling: The Two Cultures para hacerse una idea muy clara de a lo que me refiero.