Glm

Según la teoría de la relatividad, las velocidades (lineales) se suman así:

v1 <- 100000
v2 <- 100000
velocidad_luz <- 300000

suma_relativista <- function(x,y){
  (x + y) / (1 + x * y / velocidad_luz^2)
}

suma_relativista(v1, v2)
# 180000

Lo que es todavía menos conocido es que esa operación es equivalente a la suma ordinaria de velocidades a través de una transformación de ida y vuelta vía la arcotangente hiperbólica (véase esto). En concreto:

f1 <- function(x) {
  atanh(x / velocidad_luz)
}

f2 <- function(x) {
  velocidad_luz * tanh(x)
}

f2(f1(v1) + f1(v2))
# 180000

Ahora imaginemos un universo donde la velocidad máxima no es la de la luz, sino que solo están permitidas las velocidades entre 0 y 1:

Una de las consencuencias del coronavirus es que vamos a tener que replantearnos lo que significa ajustar series temporales. Es decir, comenzar a ajustar series temporales y no repetir la consabida teoría que subyace a los modelos ARIMA simplemente porque es guay.

También tendremos que replantearnos qué hacer con los outliers que la pandemia va dejando tras de sí. Y tratar de hacerlo más elegantemente que cierta gente, por supuesto. En particular, habrá que ver cuál y cómo es el efecto de los outliers en determinados modelos. En particular, en esos en los que yo más trabajo últimamente, que son los de Poisson.

[Quod si sal evanuerit in quo sallietur ad nihilum valet ultra nisi ut mittatur foras et conculcetur ab hominibus.]

Vuelvo con mi monotema de los últimos días: cómo hacer GLMs de Poisson robustos. Encuentro la tesis Robust Inference for Generalized Linear Models: Binary and Poisson Regression y pienso: ajá, será cuestión de copipegar.

Nada más lejos de la realidad. El método propuesto en la tesis está basado en asignaciones de pesos a las observaciones usando kernels con centros y anchuras basadas respectivamente en

[Este es un extracto, una píldora atómica, de mi charla del otro día sobre el modelo de Poisson y al sobredispersión.]

Aunque me guste expresar el modelo lineal de la forma

$$ y_i \sim N(a_0 + \sum_j a_j x_{ij}, \sigma_i)$$

hoy, para lo que sigue, es más conveniente la representación tradicional

$$ y_i = a_0 + \sum_j a_j x_{ij} + \epsilon_i$$

donde si no sabes lo que es cada cosa, más vale que no sigas leyendo.

Este es un anuncio de una charla que daré este viernes (2020-09-18) dentro del congreso virtual EncuentRos en la fase R. Ni que decir tiene que los detalles logísticos pueden consultarse en el enlace anterior.

Hablaré de cuestiones relativas al modelo de Possion (gran parte de las cuales pueden trasladarse también al logístico) de las que se habla poco y sobre las que la teoría que uno tropieza por ahí no es del todo clara pero que se manifiestan claramente en datos como los de la monitorización de la mortalidad, que será discutida también de pasada.

Le he dado muchas vueltas en estos últimos tiempos al asunto de la sobredispersión, particularmente en dos tipos de modelos: Poisson y logístico. Así que, aunque solo sea por proximidad semántica, se me quedan pegados ejemplos y casos de ese fenómeno mucho menos frecuente que es el de la infradispersión.

Un ejemplo ilustrativo del fenómeno que se me ocurrió el otro día era

pero hace nada, ese señor lleno de paz y amor que es Putin, nos ha regalado otro:

[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:

[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.

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.

Hay cosas tan obvias que ni se plantea la alternativa. Pero luego va R. Gomila y escribe Logistic or Linear? Estimating Causal Effects of Treatments on Binary Outcomes Using Regression Analysis que se resume en lo siguiente: cuando te interese la explicación y no la predicción, aunque tu y sea binaria, usa regresión lineal y pasa de la logística.

Nota: La sección 4.2 de An Introduction to Statistical Learning de se titula precisamente Why Not Linear Regression?

Sin embargo, basta con mirar la foto

leer la entrada de hace unos días, que se refiere a algo muy parecido (y que, en particular, describe los datos usados en el modelo que representa) y, en el peor de los casos, esto, para hacerse idea de su utilidad y relevancia.

Pensé que ya había escrito sobre el asunto porque tropecé con él en un proyecto hace un tiempo. Pero mi menoria se había confundido con otra entrada, Sobre la peculiarisima implementacion del modelo lineal en (pseudo-)Scikit-learn, donde se discute, precisamente, un problema similar si se lo mira de cierta manera o diametralmente opuesto si se ve con otra perspectiva.

Allí el problema era que Scikit-learn gestionaba muy sui generis el insidioso problema de la colinealidad. Precisamente, porque utiliza un optimizador ad hoc y no estándar para ajustar el modelo lineal.

• Alguien quería un glm forzando determinados coeficientes >0.
• Una solución 100% bayesiana no era una opción.

Hay varias opciones por ahí. Pero me ha sorprendido que la opción esté disponible en glmnet::glmnet:

Filosóficamente, es un tanto sorprendente: de alguna manera, glmnet es glm con prioris alrededor del cero. Los límites superiores e inferiores permiten introducir información a priori adicional no necesariamente compatible con la anterior.

Desde el punto de vista de la implementación, tiene sentido que estas opciones estén disponibles. glmnet usa coordinate descent como algoritmo de minimización e introducir restricciones en ese tipo de algoritmos es una trivialidad.