Pros y contras de Log Link Versus Identity Link para la regresión de Poisson

Estoy llevando a cabo una regresión de Poisson con el objetivo final de comparar (y tomar la diferencia de) los recuentos medios pronosticados entre dos niveles de factores en mi modelo: , mientras otras covariables modelo (que son todas binarias) constantes. Me preguntaba si alguien podría proporcionar algún consejo práctico sobre cuándo usar un enlace de registro frente a un enlace de identidad. ¿Cuáles son las ventajas y desventajas de estas dos funciones de enlace diferentes en la regresión de Poisson, dado mi objetivo de comparar diferencias?$\hat{\mu}_1-\hat{\mu}_2$

También tengo el mismo objetivo en mente para una regresión logística / binomial (usar un enlace logit o un enlace de identidad) para comparar la diferencia de proporciones entre dos niveles de factores y necesito un consejo similar. He leído algunas de las publicaciones aquí que abordan este tema, pero ninguna parece explicar por qué o cuándo uno podría elegir un enlace sobre el otro y cuáles podrían ser los pros / contras. ¡Gracias de antemano por tu ayuda!

ACTUALIZAR:

También me doy cuenta de que el propósito principal de usar ciertas funciones de enlaces es restringir el rango de rango de posibles valores pronosticados para que esté dentro del rango de la respuesta media (por ejemplo, para logística, el rango está restringido entre 0 y 1 y para el registro enlace, las predicciones están restringidas a números positivos). Entonces, supongo que lo que pregunto es que si uso un enlace de identidad para decir una regresión logística / binomial, y mis resultados están dentro del rango (0,1), ¿realmente hay alguna necesidad de usar una función de enlace logístico o ¿podría simplificar las ideas y utilizar un enlace de identidad?

Los contras de un enlace de identidad en el caso de la regresión de Poisson son:

• Como ha mencionado, puede producir predicciones fuera de rango.
• Puede obtener errores y advertencias extraños al intentar ajustar el modelo, porque el enlace permite que lambda sea menor que 0, pero la distribución de Poisson no está definida para dichos valores.
• Como la regresión de Poisson supone que la media y la varianza son las mismas, cuando cambia el enlace también está cambiando los supuestos sobre la varianza. Mi experiencia ha sido que este último punto es muy revelador.

Pero, en última instancia, esta es una pregunta empírica. Montar ambos modelos. Realice los controles que desee. Si el enlace de identidad tiene un AIC más bajo y funciona igual o mejor en todas sus otras comprobaciones, ejecute con el enlace de identidad.

En el caso del modelo logit frente al modelo de probabilidad lineal (es decir, a lo que se refiere como el enlace de identidad), la situación es mucho más sencilla. A excepción de algunos casos muy exóticos en econometría (que encontrará si realiza una búsqueda), el modelo logit es mejor: hace menos suposiciones y es lo que usa la mayoría de la gente. Usar el modelo de probabilidad lineal en su lugar estaría a punto de ser perverso.

Con respecto a la interpretación de los modelos, si está utilizando R, hay dos paquetes geniales que harán todo el trabajo pesado: efectos , que son súper fáciles de usar, y zelig , que es más difícil de usar pero genial si desea hacer predicciones .

En el caso de los modelos de Poisson, también diría que la aplicación a menudo dicta si sus covariables actuarían de forma aditiva (lo que implicaría un enlace de identidad) o multiplicativamente en una escala lineal (lo que implicaría un enlace logarítmico). Pero los modelos de Poisson con un enlace de identidad también normalmente tienen sentido y solo pueden ajustarse de manera estable si se imponen restricciones de no negatividad a los coeficientes ajustados; esto se puede hacer usando la nnpoisfunción en el addregpaquete R o la nnlmfunción en elNNLMpaquete. Por lo tanto, no estoy de acuerdo en que uno deba ajustar los modelos de Poisson con una identidad y un enlace de registro y ver cuál termina teniendo el mejor AIC e inferir el mejor modelo basado en motivos puramente estadísticos, más bien, en la mayoría de los casos, lo dicta estructura subyacente del problema que se intenta resolver o los datos disponibles.

Por ejemplo, en la cromatografía (análisis GC / MS) a menudo se mediría la señal superpuesta de varios picos con forma de Gauss aproximadamente y esta señal superpuesta se mide con un multiplicador de electrones, lo que significa que la señal medida es un recuento de iones y, por lo tanto, distribución de Poisson. Dado que cada uno de los picos tiene, por definición, una altura positiva y actúa de manera aditiva y el ruido es Poisson, un modelo de Poisson no negativo con enlace de identidad sería apropiado aquí, y un modelo de Poisson de enlace de registro sería completamente incorrecto. En ingeniería , la pérdida de Kullback-Leibler a menudo se usa como una función de pérdida para tales modelos, y minimizar esta pérdida es equivalente a optimizar la probabilidad de un modelo de Poisson de enlace de identidad no negativo (también hay otras medidas de divergencia / pérdida como la divergencia alfa o beta que tienen a Poisson como un caso especial).

A continuación se muestra un ejemplo numérico, que incluye una demostración de que un enlace de identidad regular sin restricciones Poisson GLM no encaja (debido a la falta de restricciones de no negatividad) y algunos detalles sobre cómo ajustar modelos de Poisson de enlace de identidad no negativos utilizandonnpois, aquí en el contexto de desconvolucionar una superposición medida de picos cromatográficos con ruido de Poisson sobre ellos usando una matriz covariable con bandas que contiene copias desplazadas de la forma medida de un solo pico. La no negatividad aquí es importante por varias razones: (1) es el único modelo realista para los datos disponibles (los picos aquí no pueden tener alturas negativas), (2) es la única forma de ajustar de manera estable un modelo de Poisson con enlace de identidad (como de lo contrario, las predicciones podrían ser negativas para algunos valores covariables, lo que no tendría sentido y daría problemas numéricos cuando uno intentara evaluar la probabilidad), (3) los actos de no negatividad para regularizar el problema de regresión y ayudan en gran medida a obtener estimaciones estables (p. ej. por lo general, no obtienes los problemas de sobreajuste como con la regresión sin restricciones ordinaria,las restricciones de no negatividad dan como resultado estimaciones más dispersas que con frecuencia están más cerca de la verdad básica; para el siguiente problema de desconvolución, por ejemplo, el rendimiento es casi tan bueno como la regularización de LASSO, pero sin requerir que uno ajuste ningún parámetro de regularización. (La regresión penalizada con pseudonorma L0 aún funciona un poco mejor pero a un costo computacional mayor )

# we first simulate some data
require(Matrix)
n = 200
x = 1:n
npeaks = 20
set.seed(123)
u = sample(x, npeaks, replace=FALSE) # unkown peak locations
peakhrange = c(10,1E3) # peak height range
h = 10^runif(npeaks, min=log10(min(peakhrange)), max=log10(max(peakhrange))) # unknown peak heights
a = rep(0, n) # locations of spikes of simulated spike train, which are assumed to be unknown here, and which needs to be estimated from the measured total signal
a[u] = h
gauspeak = function(x, u, w, h=1) h*exp(((x-u)^2)/(-2*(w^2))) # peak shape function
bM = do.call(cbind, lapply(1:n, function (u) gauspeak(x, u=u, w=5, h=1) )) # banded matrix with peak shape measured beforehand
y_nonoise = as.vector(bM %*% a) # noiseless simulated signal = linear convolution of spike train with peak shape function
y = rpois(n, y_nonoise) # simulated signal with random poisson noise on it - this is the actual signal as it is recorded
par(mfrow=c(1,1))
plot(y, type="l", ylab="Signal", xlab="x", main="Simulated spike train (red) to be estimated given known blur kernel & with Poisson noise")
lines(a, type="h", col="red")

# let's now deconvolute the measured signal y with the banded covariate matrix containing shifted copied of the known blur kernel/peak shape bM

# first observe that regular OLS regression without nonnegativity constraints would return very bad nonsensical estimates
weights <- 1/(y+1) # let's use 1/variance = 1/(y+eps) observation weights to take into heteroscedasticity caused by Poisson noise
a_ols <- lm.fit(x=bM*sqrt(weights), y=y*sqrt(weights))$coefficients # weighted OLS
plot(x, y, type="l", main="Ground truth (red), unconstrained OLS estimate (blue)", ylab="Peak shape", xlab="x", ylim=c(-max(y),max(y)))
lines(x,-y)
lines(a, type="h", col="red", lwd=2)
lines(-a_ols, type="h", col="blue", lwd=2)

# now we use weighted nonnegative least squares with 1/variance obs weights as an approximation of nonnegative Poisson regression
# this gives very good estimates & is very fast
library(nnls)
library(microbenchmark)
microbenchmark(a_wnnls <- nnls(A=bM*sqrt(weights),b=y*sqrt(weights))$x) # 7 ms
plot(x, y, type="l", main="Ground truth (red), weighted nnls estimate (blue)", ylab="Signal (black) & peaks (red & blue)", xlab="Time", ylim=c(-max(y),max(y)))
lines(x,-y)
lines(a, type="h", col="red", lwd=2)
lines(-a_wnnls, type="h", col="blue", lwd=2)
# note that this weighted least square estimate in almost identical to  the nonnegative Poisson estimate below and that it fits way faster!!!

# an unconstrained identity-link Poisson GLM will not fit:
glmfit = glm.fit(x=as.matrix(bM), y=y, family=poisson(link=identity), intercept=FALSE)
# returns Error: no valid set of coefficients has been found: please supply starting values

# so let's try a nonnegativity constrained identity-link Poisson GLM, fit using bbmle (using port algo, ie Quasi Newton BFGS):
library(bbmle)
XM=as.matrix(bM)
colnames(XM)=paste0("v",as.character(1:n))
yv=as.vector(y)
LL_poisidlink <- function(beta, X=XM, y=yv){ # neg log-likelihood function
  -sum(stats::dpois(y, lambda = X %*% beta, log = TRUE)) # PS regular log-link Poisson would have exp(X %*% beta)
}
parnames(LL_poisidlink) <- colnames(XM)
system.time(fit <- mle2(
  minuslogl = LL_poisidlink ,
  start = setNames(a_wnnls+1E-10, colnames(XM)), # we initialise with weighted nnls estimates, with approx 1/variance obs weights
  lower = rep(0,n),
  vecpar = TRUE,
  optimizer = "nlminb"
)) # very slow though - takes 145s 
summary(fit)
a_nnpoisbbmle = coef(fit)
plot(x, y, type="l", main="Ground truth (red), nonnegative Poisson bbmle ML estimate (blue)", ylab="Signal (black) & peaks (red & blue)", xlab="Time", ylim=c(-max(y),max(y)))
lines(x,-y)
lines(a, type="h", col="red", lwd=2)
lines(-a_nnpoisbbmle, type="h", col="blue", lwd=2)

# much faster is to fit nonnegative Poisson regression using nnpois using an accelerated EM algorithm:
library(addreg)
microbenchmark(a_nnpois <- nnpois(y=y,
                                  x=as.matrix(bM),
                                  standard=rep(1,n),
                                  offset=0,
                                  start=a_wnnls+1.1E-4, # we start from weighted nnls estimates 
                                  control = addreg.control(bound.tol = 1e-04, epsilon = 1e-5),
                                  accelerate="squarem")$coefficients) # 100 ms
plot(x, y, type="l", main="Ground truth (red), nonnegative Poisson nnpois estimate (blue)", ylab="Signal (black) & peaks (red & blue)", xlab="Time", ylim=c(-max(y),max(y)))
lines(x,-y)
lines(a, type="h", col="red", lwd=2)
lines(-a_nnpois, type="h", col="blue", lwd=2)

# or to fit nonnegative Poisson regression using nnlm with Kullback-Leibler loss using a coordinate descent algorithm:
library(NNLM)
system.time(a_nnpoisnnlm <- nnlm(x=as.matrix(rbind(bM)),
                                 y=as.matrix(y, ncol=1),
                                 loss="mkl", method="scd",
                                 init=as.matrix(a_wnnls, ncol=1),
                                 check.x=FALSE, rel.tol=1E-4)$coefficients) # 3s
plot(x, y, type="l", main="Ground truth (red), nonnegative Poisson nnlm estimate (blue)", ylab="Signal (black) & peaks (red & blue)", xlab="Time", ylim=c(-max(y),max(y)))
lines(x,-y)
lines(a, type="h", col="red", lwd=2)
lines(-a_nnpoisnnlm, type="h", col="blue", lwd=2)