Formato de entrada para la respuesta en glm binomial en R

Question

En R, hay tres métodos para formatear los datos de entrada para una regresión logística usando la función glm:

1. Los datos pueden estar en un formato "binario" para cada observación (por ejemplo, y = 0 o 1 para cada observación);
2. Los datos pueden estar en formato "Wilkinson-Rogers" (por ejemplo, y = cbind(success, failure)) con cada fila representando un tratamiento; o
3. Los datos pueden estar en un formato ponderado para cada observación (por ejemplo, y = 0.3, pesos = 10).

Los tres enfoques producen los mismos coeficientes estimados, pero difieren en los grados de libertad y los valores resultantes de desviación y puntuaciones AIC. Los dos últimos métodos tienen menos observaciones (y por lo tanto grados de libertad) porque utilizan cada tratamiento para el número de observaciones, mientras que el primero utiliza cada observación para el número de observaciones.

Mi pregunta: ¿Existen ventajas numéricas o estadísticas al utilizar un formato de entrada sobre otro? La única ventaja que veo es no tener que reformatear los datos en R para usarlos con el modelo.

He revisado la documentación de glm, he buscado en la web y en este sitio y encontré un post tangencialmente relacionado, pero no encontré orientación sobre este tema.

Aquí hay un ejemplo simulado que demuestra este comportamiento:

# Escribir una función para ayudar a simular datos
drc4 <- function(x, b =1.0, c = 0, d = 1, e = 0){
    (d - c)/ (1 + exp(-b * (log(x)  - log(e))))
}
# simular la forma larga del conjunto de datos
nReps = 20
dfLong <- data.frame(dosis = rep(seq(0, 10, by = 2), each = nReps))
dfLong$mortalidad <-rbinom(n = dim(dfLong)[1], size = 1,
                               prob = drc4(dfLong$dosis, b = 2, e = 5))

# agregar para crear la forma corta del conjunto de datos
dfShort <- aggregate(dfLong$mortalidad, by = list(dfLong$dosis), 
                     FUN = sum)
colnames(dfShort) <- c("dosis", "mortalidad")
dfShort$supervivencia <- nReps - dfShort$mortalidad 
dfShort$nReps <- nReps
dfShort$mortalidadP <- dfShort$mortalidad / dfShort$nReps

fitShort <- glm( cbind(mortalidad, supervivencia) ~ dosis, 
                 data = dfShort, 
                 family = "binomial")
summary(fitShort)

fitShortP <- glm( mortalidadP ~ dosis, data = dfShort, 
                  weights = nReps,     
                  family = "binomial")
summary(fitShortP)

fitLong <- glm( mortalidad ~ dosis, data = dfLong, 
                family = "binomial")
summary(fitLong)

Answer 1

No hay una razón estadística para preferir uno sobre el otro, además de la claridad conceptual. Aunque los valores de desviación reportados son diferentes, estas diferencias se deben completamente al modelo saturado. Por lo tanto, cualquier comparación utilizando la desviación relativa entre modelos no se ve afectada, ya que la verosimilitud logarítmica del modelo saturado se cancela.

Creo que es útil pasar por el cálculo explícito de la desviación.

La desviación de un modelo es 2*(LL(Modelo Saturado) - LL(Modelo)). Supongamos que tienes $i$ celdas diferentes, donde $n_i$ es el número de observaciones en la celda $i$, $p_i$ es la predicción del modelo para todas las observaciones en la celda $i$, y $y_{ij}$ es el valor observado (0 o 1) para la $j$-ésima observación en la celda $i$.

Forma Larga

La verosimilitud logarítmica del modelo (propuesto o nulo) es $$\sum_i\sum_j\left(\log(p_i)y_{ij} + \log(1 - p_i)(1 - y_{ij})\right)$$

y la verosimilitud logarítmica del modelo saturado es $$\sum_i\sum_j \left(\log(y_{ij})y_{ij} + \log(1 - y_{ij})(1-y_{ij})\right).$$ Esto es igual a 0, porque $y_{ij}$ es 0 o 1. Nota que $\log(0)$ es indefinido, pero por conveniencia lee $0\log(0)$ como una abreviatura de $\lim_{x \to 0^+}x\log(x)$, que es 0.

Forma Corta (ponderada)

Nota que una distribución binomial en realidad no puede tomar valores no enteros, pero de todos modos podemos calcular una "verosimilitud logarítmica" utilizando la fracción de éxitos observados en cada celda como la respuesta, y ponderando cada término de la suma en el cálculo de la verosimilitud logarítmica por el número de observaciones en esa celda.

$$\sum_in_i \left(\log(p_i)\sum_jy_{ij}/n_i + \log(1 - p_i)(1 - \sum_j(y_{ij}/n_i)\right)$$

Esto es exactamente igual a la desviación del modelo que calculamos anteriormente, lo cual puedes ver extrayendo la suma sobre $j$ en la ecuación de forma larga lo más posible.

Mientras tanto, la desviación saturada es diferente. Dado que ya no tenemos respuestas de 0-1, incluso con un parámetro por observación no podemos llegar exactamente a 0. En cambio, la verosimilitud logarítmica del modelo saturado es

$$\sum_i n_i\left(\log(\sum_jy_{ij}/n_i)\sum_jy_{ij}/n_i + \log(1 - \sum_jy_{ij}/n_i)(1-\sum_jy_{ij}/n_i)\right).$$

En tu ejemplo, puedes verificar que dos veces esta cantidad es la diferencia entre los valores de desviación nula y residual reportados para ambos modelos.

ni = dfShort$nReps
yavg = dfShort$mortalityP
sum.terms <-ni*(log(yavg)*yavg + log(1 - yavg)*(1 - yavg))
# Necesitamos manejar NaN cuando yavg es exactamente 0
sum.terms[1] <- log(1 - yavg[1])*(1 - yavg[1])

2*sum(sum.terms)
fitShortP$deviance - fitLong$deviance