Glm

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:

Rob Tibshirani ha seleccionado recientemente una lista de nueve artículos de estadística publicados desde 1970. Son estos:

• Regression models and life tables (with discussion) (Cox 1972). Según Tibshirani, David Cox merece el Nobel de Medicina por él.
• Generalized linear models (Nelder and Wedderburn 1972). Es la base del paquete glm de R.
• Maximum Likelihood from Incomplete Data via the {EM} Algorithm (with discussion) (Dempster, Laird, and Rubin 1977).
• Bootstrap methods: another look at the jackknife (Efron 1979).
• Classification and regression trees (Breiman, Friedman, Olshen and Stone 1984).
• How biased is the error rate of a prediction rule? (Efron 1986).
• Sampling based approaches to calculating marginal densities (Gelfand and Smith 1990).
• Controlling the false discovery rate: a practical and powerful approach to multiple testing (Benjamini and Hochberg 1995).
• A decision-theoretic generalization of online learning and an application to boosting (Freund and Schapire 1995).

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, $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).