Glm

Hace años tomé el curso de NLP de M. Collings en Coursera (¡muy recomendable!), uno de cuyos capítulos trataba de los llamados modelos loglineales. En esto, Collings sigue una nomenclatura un tanto personal porque la mayor parte de la gente se refiere con ese nombre a algo que no es exactamente lo mismo (y dentro del mundo de las tablas de contingencia).

El otro día, sin embargo, me pensé que los modelos loglineales à la Collings me serían muy útiles para un problema de clasificación en el que estamos trabajando. Y repasándolos… me di cuenta de que eran versiones de algo ya conocido: GLMs multinomiales con regularización. Sí, como estos.

Observo en The deadly effects of losing health insurance cómo el efecto de interés, 15% sobre una población se convierte en efectos del 16%, 23% y 30% en sus tres subpoblaciones (útimas columnas de la tabla que ocupa la página 25). Es raro que el efecto combinado no esté cerca de la media ponderada (por población) de cada uno de sus subcomponentes.

Podría ser Simpson, pero hay motivos para pensar que hayan cambiado las proporciones de las poblaciones subyacentes (demasiado). Habría un efecto Simpson, por ejemplo, si se hubiese incrementado sustancialmente la proporción del grupo con el efecto (no confundir con la variación del efecto) globalmente más pequeño antes y después del tratamiento. Pero dudo que sea el caso.

Liberados del estrecho ámbito de nuestra original mentira sugerente gracias a la relación que descubrimos entre residuos y gradientes cuando las pérdidas son cuadráticas podemos adentrarnos en ámbitos más extensos.

Lo que discutimos del gradiente tiene una interpretación fácilmente inteligible en el caso de pérdidas cuadráticas. Pero ni la pérdida de interpretabilidad nos impide extender el razonamiento de la entrada anterior a funciones de pérdida distintas de la cuadrática siempre que podamos calcular un gradiente.

Estos días he estado haciendo de campaña promoviendo el uso de nuevas técnicas de análisis de datos en ámbitos como, p.e., el riesgo de crédito, uno de esos campos sujetos al parecer de un regulador (el Banco de España, en este caso).

La gente con la que he debatido al respecto tiende a aplicar esa forma cuasiperfecta de censura que es la autocensura previa. La autocensura previa ni siquiera requiere la acción explícita del censor: es el potencial censurado el que la aplica de mejor o peor gana automáticamente… por si las moscas.

Preguntaba el otro día Emilio Torres esto en R-help-es. Resumo la pregunta. Se trata de una simulación de unos datos y su ajuste mediante una regresión logística para ver si los coeficientes obtenidos son o no los esperados (teóricamente y por construcción).

El código de Emilio (cuyos resultados no podemos reproducir porque no nos ha contado qué similla usa) es

logisticsimulation <- function(n){
  dat <- data.frame(x1=sample(0:1, n,replace=TRUE),
                    x2=sample(0:1, n,replace=TRUE))
  odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 *dat$x1* dat$x2 )
  pr <- odds/(1+odds)
  res <- replicate(100, {
    dat$y <- rbinom(n,1,pr)
    coef(glm(y ~ x1*x2, data = dat, family = binomial()))
  })
  t(res)
}

res <- logisticsimulation(100)
apply(res,2,median)
## (Intercept)          x1          x2       x1:x2
## -1.0986123 -18.4674562  20.4823593  -0.0512933

Efectivamente, los coeficientes están lejos de los esperados, i.e., -1, -4, 7 y 1.

Creo que todos sabéis la historia de las admisiones de la Universidad de Berkeley y la paradoja de Simpson. Con palabras, muchas palabras, está contado, por ejemplo, aquí. Y si buscáis ubc admissions simpson en Google la encontraréis también en modo --verbose en muchos más sitios.

En R puede resumirse en

library(reshape2)
library(plyr)

data(UCBAdmissions)

raw <- as.data.frame(UCBAdmissions)

dat <- dcast(raw, Gender + Dept ~ <a href="http://inside-r.org/packages/cran/AdMit">Admit)

mod.0 <- glm(cbind(Admitted, Rejected) ~ Gender, data = dat, family = binomial)
mod.1 <- glm(cbind(Admitted, Rejected) ~ Gender + Dept, data = dat, family = binomial)

Echad un vistazo a los coeficientes de Gender en ambos modelos y veréis.

Tuve que saltarme una diapositiva en el DataBeers de Madrid del pasado jueves.

(A propósito, aquí están las 1+20 diapositivas.)

La decimonona, de la que trata la entrada, viene a hablar de lo siguiente. Tenemos una base de datos con sujetos (ids) que hacen cosas en determinados momentos. No es inhabitual calcular la frecuencia de esos sujetos así:

select id, count(*) as freq
from mytabla
where fecha between current_date - 7 and current_date
group by id
;

Esa variable se utiliza frecuentemente ya sea como descriptor de los sujetos o como alimento de otros modelos.

Lo que escribí hace un par de días sobre procesos puntuales, ahora me doy cuenta, podía haberse resuelto con nuestro viejo amigo glm.

Ejecuto el código del otro día y obtengo (para un caso nuevo)

          mu       alfa verosimilitud delta
    1  0.4493158 0.50000000      340.6141     1
    2  0.2675349 0.40457418      307.3939     2
    3  0.1894562 0.28917407      293.4696     3
    4  0.1495654 0.22237707      287.0784     4
    5  0.1243791 0.18079703      281.3900     5
    6  0.1142837 0.14913172      284.9227     6
    7  0.1217504 0.12150745      288.5448     7
    8  0.1214365 0.10424818      289.3282     8
    9  0.1204605 0.09148817      290.9081     9
    10 0.1315896 0.07857330      295.3935    10</code>

que significa que el parámetro óptimo es delta = 5, mu = 0.124 y alfa = 0.18.

Ahora hago

    cuantos.previos <- function(i, muestra, delta){
      indices <- Filter(function(x) x < i & x > i - delta, 1:n)
      cuantos <- sum(muestra[indices])
    }

    fit.glm <- function(delta){
      prev <- sapply(1:length(muestra),
                     cuantos.previos, muestra, delta)
      dat  <- data.frame(muestra = muestra, prev = prev)

      res.glm <- glm(muestra ~ prev, data = dat,
                     family = poisson(link = "identity"))
      c(delta, res.glm$coefficients, summary(res.glm)$aic)
    }

    res.glm <- sapply(1:10, fit.glm)
    res.glm <- as.data.frame(t(res.glm))
    colnames(res.glm) <- c("delta", "mu", "alfa", "aic")

y obtengo

Supongamos que queremos construir un modelo para predecir quién ganará un determinado partido de baloncesto basándonos en datos diversos. Y en un histórico, por supuesto.

Podemos utilizar una regresión logística así:

set.seed(1234)

my.coefs <- -2:2
n <- 200
train.n <- floor(2*n/3)

test.error.glm <- function(){
  X <- matrix(rnorm(n*5), n, 5)
  Y <- (0.2 + X %*% my.coefs + rnorm(n)) > 0

  train <- sample(1:n, train.n)

  X <- as.data.frame(X)
  X$Y <- Y

  mod.glm <- glm(Y ~ ., data = X[train,],
    family = binomial)

  glm.pred <- predict(mod.glm, X[-train,],
    type = "response")

  error <- length(glm.pred) -
    sum(diag(table(glm.pred > 0.5, Y[-train,])))
}

errores.glm <- replicate(1000, test.error.glm())

El código anterior hace lo siguiente:

No conocía el paquete gbm. Pero como ahora ando rodeado de data scientists que no son estadísticos…

Bueno, la cuestión es que había que ajustar un modelo para el que yo habría hecho algo parecido a

dat <- read.csv("http://www.ats.ucla.edu/stat/data/poisson_sim.csv")
summary(m.glm <- glm(num_awards ~ prog + math, family = "poisson", data = dat))
# Call:
#   glm(formula = num_awards ~ prog + math, family = "poisson", data = dat)
#
# Deviance Residuals:
#   Min       1Q   Median       3Q      Max
# -2.1840  -0.9003  -0.5891   0.3948   2.9539
#
# Coefficients:
#   Estimate Std. Error z value Pr(>|z|)
# (Intercept) -5.578057   0.676823  -8.242   <2e-16 ***
#   prog         0.123273   0.163261   0.755     0.45
# math         0.086121   0.009586   8.984   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# (Dispersion parameter for poisson family taken to be 1)
#
# Null deviance: 287.67  on 199  degrees of freedom
# Residual deviance: 203.45  on 197  degrees of freedom
# AIC: 385.51
#
# Number of Fisher Scoring iterations: 6

como en esta página.

Fueron problemas planteados por Frank Harrell, recopilados aquí y ahora traducidos por mí para mi bitácora.

Problemas de la regresión paso a paso:

• La R-cuadrado obtenida está muy sesgada hacia arriba.
• Los test F y chi-cuadrado que aparecen al lado de las variables no siguen dichas distribuciones.
• Los intervalos de confianza son demasiado (e incorrectamente) estrechos.
• Los p-valores obtenidos no tienen el significado esperado y el de corregirlos adecuadamente es un problema muy difícil.
• Proporciona coeficientes sesgados y excesivamente grandes.
• Tiene problemas serios en caso de colinealidad en las variables.
• Está basado en métodos que fueron pensados para probar hipótesis preestablecidas.
• Incrementar el número de muestras no corrige los problemas anteriores.
• Nos permite no tener que pensar sobre el problema.
• Consume mucho papel.

Algunas conclusiones:

He tropezado con una extensión curiosa y que no conocía del modelo logístico que lo emparenta un tanto con los modelos de supervivencia. Es un problema que aparece en los modelos de los actuarios, por ejemplo, y en la supervivencia de nidos (sí, nidos de bichos alados), parece.

Es el siguiente: supongamos que unos sujetos están expuestos a un cierto suceso cuya probabilidad, $latex p_i$, depende del sujeto a través del esquema habitual de la regresión logística (es decir, depende de algunas variables como el sexo, etc., a través de una fórmula lineal cuyos coeficientes interesa estimar).