Prós e contras do Log Link versus Link de identidade para regressão de Poisson

Estou realizando uma regressão de Poisson com o objetivo final de comparar (e tendo a diferença de) as contagens médias previstas entre dois níveis de fator no meu $\hat{\mu}_1-\hat{\mu}_2$ , mantendo outras variáveis do modelo (que são todos binária ) constante. Eu queria saber se alguém poderia fornecer alguns conselhos práticos sobre quando usar um link de log versus um link de identidade. Quais são os prós e contras dessas duas funções diferentes de link na regressão de Poisson, dado o meu objetivo de comparar diferenças?

Também tenho o mesmo objetivo em mente para uma regressão logística / binomial (para usar um link de logit ou de identidade) para comparar a diferença de proporções entre dois níveis de fator e precisar de conselhos semelhantes. Eu li algumas das postagens aqui que abordam esse problema, mas nenhuma parece explicar por que ou quando alguém pode escolher um link em vez do outro e quais podem ser os prós / contras. Agradeço antecipadamente por sua ajuda!

ATUALIZAR:

Também percebo que o principal objetivo do uso de certas funções de links é restringir o intervalo de valores possíveis previstos para estar dentro do intervalo da resposta média (por exemplo, para logística, o intervalo é restrito para estar entre 0 e 1 e para o log link, as previsões são restritas para serem números positivos). Então, acho que o que estou perguntando é que, se eu usar um link de identidade para, por exemplo, uma regressão logística / binomial, e meus resultados estiverem dentro do intervalo (0,1), será realmente necessário usar uma função de link logístico ou eu poderia simplificar o pensamento e usar um link de identidade?

Os contras de um vínculo de identidade no caso da regressão de Poisson são:

• Como você mencionou, ele pode produzir previsões fora da faixa.
• Você pode receber erros e avisos estranhos ao tentar ajustar o modelo, porque o link permite que o lambda seja menor que 0, mas a distribuição Poisson não está definida para esses valores.
• Como a regressão de Poisson assume que a média e a variação são as mesmas, quando você altera o link, também altera as suposições sobre a variação. Minha experiência foi que esse último ponto é mais revelador.

Mas, finalmente, esta é uma questão empírica. Encaixe nos dois modelos. Faça as verificações que desejar. Se o link de identidade tiver um AIC mais baixo e for bom ou melhor em todas as suas outras verificações, execute com o link de identidade.

No caso do modelo logit x modelo de probabilidade linear (ou seja, o que você chama de link de identidade), a situação é muito mais direta. Exceto por alguns casos muito exóticos em econometria (que você encontrará se fizer uma pesquisa), o modelo de logit é melhor: faz menos suposições e é o que a maioria das pessoas usa. Usar o modelo de probabilidade linear em seu lugar quase perversa.

Com relação à interpretação dos modelos, se você estiver usando R, existem dois ótimos pacotes que farão todo o trabalho pesado: efeitos , que são super fáceis de usar, e zelig , que é mais difícil de usar, mas ótimo se você quiser fazer previsões .

No caso dos modelos de Poisson, eu também diria que o aplicativo geralmente determina se suas covariáveis ​​agiriam de maneira aditiva (o que implicaria um link de identidade) ou multiplicativamente em uma escala linear (o que implicaria um link de log). Mas os modelos de Poisson com um link de identidade também normalmente normalmente só fazem sentido e podem ser ajustados de forma estável se alguém impuser restrições de não-negatividade aos coeficientes ajustados - isso pode ser feito usando a nnpoisfunção no addregpacote R ou a nnlmfunção noNNLMpacote. Portanto, não concordo que se deva ajustar os modelos de Poisson com um link de identidade e log e ver qual deles acaba tendo a melhor AIC e inferir o melhor modelo com base em fundamentos puramente estatísticos - na verdade, na maioria dos casos, é ditado pelo estrutura subjacente do problema que se tenta resolver ou os dados disponíveis.

Por exemplo, em cromatografia (análise GC / MS), costumava-se medir o sinal sobreposto de vários picos aproximadamente em forma de Gauss e esse sinal sobreposto é medido com um multiplicador de elétrons, o que significa que o sinal medido é contado por íons e, portanto, distribuído por Poisson. Como cada um dos picos tem, por definição, uma altura positiva e age de maneira aditiva e o ruído é Poisson, um modelo não-negativo de Poisson com link de identidade seria apropriado aqui, e um modelo de Poisson com log-log seria completamente errado. Na engenharia , a perda de Kullback-Leibler é frequentemente usada como uma função de perda para esses modelos, e minimizar essa perda é equivalente a otimizar a probabilidade de um modelo de Poisson não-negativo de ligação de identidade (também existem outras medidas de divergência / perda, como divergência alfa ou beta que têm Poisson como um caso especial).

Abaixo está um exemplo numérico, incluindo uma demonstração de que um Poisson GLM de link de identidade irrestrito regular não se encaixa (devido à falta de restrições de não-negatividade) e alguns detalhes sobre como ajustar modelos de Poisson de link de identidade não-negativos usandonnpois, aqui no contexto da desconvolução de uma superposição medida de picos cromatográficos com ruído de Poisson, usando uma matriz covariada em faixas que contém cópias deslocadas da forma medida de um único pico. A não-negatividade aqui é importante por várias razões: (1) é o único modelo realista para os dados em mãos (os picos aqui não podem ter alturas negativas), (2) é a única maneira de ajustar de forma estável um modelo de Poisson com link de identidade (como caso contrário, as previsões poderiam para alguns valores covariáveis ​​serem negativas, o que não faria sentido e causaria problemas numéricos quando se tentasse avaliar a probabilidade), (3) a não-negatividade atua para regularizar o problema de regressão e ajuda muito a obter estimativas estáveis ​​(por exemplo, normalmente você não tem problemas de adaptação excessiva como na regressão irrestrita comum,restrições de não-negatividade resultam em estimativas mais esparsas, freqüentemente mais próximas da verdade básica; para o problema de desconvolução abaixo, por exemplo, o desempenho é tão bom quanto a regularização do LASSO, mas sem a necessidade de ajustar qualquer parâmetro de regularização. (A regressão penalizada L0-pseudonorm ainda apresenta um desempenho ligeiramente melhor, mas com um custo computacional maior )

# 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)