Avantages et inconvénients du lien de log contre le lien d'identité pour la régression de Poisson

J'effectue une régression de Poisson dans le but final de comparer (et de prendre la différence de) les comptes moyens prévus entre deux niveaux de facteurs dans mon modèle: , tout en maintenant d'autres covariables du modèle (qui sont toutes binaires) constantes. Je me demandais si quelqu'un pouvait fournir des conseils pratiques sur le moment d'utiliser un lien de journal par rapport à un lien d'identité. Quels sont les avantages de ces deux fonctions de liaison différentes dans la régression de Poisson, étant donné mon objectif de comparer les différences?$\hat{\mu}_1-\hat{\mu}_2$

J'ai également le même objectif en tête pour une régression logistique / binomiale (utiliser un lien logit ou un lien d'identité) pour comparer la différence de proportions entre deux niveaux de facteurs et j'ai besoin de conseils similaires. J'ai lu certains des articles ici qui traitent de ce problème, mais aucun ne semble expliquer pourquoi ou quand l'un pourrait choisir un lien plutôt que l'autre et quels pourraient être les avantages / inconvénients. Merci d'avance pour votre aide!

MISE À JOUR:

Je me rends également compte que le but principal de l'utilisation de certaines fonctions de liens est de restreindre la plage de plage des valeurs prédites possibles à la plage de la réponse moyenne (par exemple, pour la logistique, la plage est limitée entre 0 et 1 et pour le journal lien, les prévisions sont limitées à des nombres positifs). Donc, je suppose que ce que je demande, c'est que si j'utilise un lien d'identité, par exemple une régression logistique / binomiale, et que mes résultats se situent dans la plage (0,1), est-il vraiment nécessaire d'utiliser une fonction de lien logistique ou pourrais-je simplement simplifier la réflexion et utiliser un lien d'identité?

Les inconvénients d'un lien d'identité dans le cas de la régression de Poisson sont:

• Comme vous l'avez mentionné, il peut produire des prédictions hors de portée.
• Vous pouvez obtenir des erreurs et des avertissements étranges lorsque vous tentez d'ajuster le modèle, car le lien permet à lambda d'être inférieur à 0, mais la distribution de Poisson n'est pas définie pour ces valeurs.
• Comme la régression de Poisson suppose que la moyenne et la variance sont les mêmes, lorsque vous modifiez le lien, vous modifiez également les hypothèses sur la variance. D'après mon expérience, ce dernier point est le plus révélateur.

Mais, finalement, c'est une question empirique. Convient aux deux modèles. Effectuez les vérifications que vous souhaitez. Si le lien d'identité a un AIC inférieur et fait aussi bien ou mieux sur tous vos autres contrôles, exécutez-le avec le lien d'identité.

Dans le cas du modèle logit par rapport au modèle de probabilité linéaire (c'est-à-dire ce que vous appelez le lien d'identité), la situation est beaucoup plus simple. À l'exception de quelques cas très exotiques en économétrie (que vous trouverez si vous effectuez une recherche), le modèle logit est meilleur: il fait moins d'hypothèses et est ce que la plupart des gens utilisent. Utiliser le modèle de probabilité linéaire à sa place reviendrait à être pervers.

En ce qui concerne l'interprétation des modèles, si vous utilisez R, il y a deux super packages qui feront tout le gros du travail: les effets , qui est super facile à utiliser, et zelig , qui est plus difficile à utiliser mais génial si vous voulez faire des prédictions .

Dans le cas des modèles de Poisson, je dirais également que l'application dicte souvent si vos covariables agiraient de manière additive (ce qui impliquerait alors un lien d'identité) ou multiplicative sur une échelle linéaire (ce qui impliquerait alors un lien logarithmique). Mais les modèles de Poisson avec un lien d'identité n'ont normalement de sens et ne peuvent être ajustés de manière stable que si l'on impose des contraintes de non négativité aux coefficients ajustés - cela peut être fait en utilisant la nnpoisfonction dans le addregpackage R ou en utilisant la nnlmfonction dans leNNLMpaquet. Je ne suis donc pas d'accord pour dire que l'on devrait adapter les modèles de Poisson à la fois à une identité et à un lien de log et voir lequel finit par avoir le meilleur AIC et inférer le meilleur modèle basé sur des motifs purement statistiques - plutôt, dans la plupart des cas, il est dicté par le structure sous-jacente du problème que l'on essaie de résoudre ou des données disponibles.

Par exemple, en chromatographie (analyse GC / MS), on mesure souvent le signal superposé de plusieurs pics de forme gaussienne approximative et ce signal superposé est mesuré avec un multiplicateur d'électrons, ce qui signifie que le signal mesuré est le nombre d'ions et donc la distribution de Poisson. Étant donné que chacun des pics a par définition une hauteur positive et agit de manière additive et que le bruit est Poisson, un modèle de Poisson non négatif avec lien d'identité serait approprié ici, et un modèle de Poisson à lien log serait tout à fait faux. En ingénierie, la perte de Kullback-Leibler est souvent utilisée comme fonction de perte pour de tels modèles, et minimiser cette perte équivaut à optimiser la probabilité d'un modèle de Poisson à lien d'identité non négatif (il existe également d'autres mesures de divergence / perte comme la divergence alpha ou bêta qui ont Poisson comme cas particulier).

Vous trouverez ci-dessous un exemple numérique, comprenant une démonstration qu'un lien d'identité non contraint régulier Poisson GLM ne correspond pas (en raison du manque de contraintes de non-négativité) et quelques détails sur la façon d'adapter les modèles de Poisson à lien d'identité non négatifs en utilisantnnpois, ici dans le contexte de la déconvolution d'une superposition mesurée de pics chromatographiques avec du bruit de Poisson sur eux en utilisant une matrice de covariables en bandes qui contient des copies décalées de la forme mesurée d'un seul pic. La non négativité ici est importante pour plusieurs raisons: (1) c'est le seul modèle réaliste pour les données disponibles (les pics ici ne peuvent pas avoir des hauteurs négatives), (2) c'est le seul moyen d'ajuster de manière stable un modèle de Poisson avec un lien d'identité (comme sinon, les prédictions pourraient devenir négatives pour certaines valeurs de covariables, ce qui n'aurait pas de sens et poserait des problèmes numériques lorsque l'on tenterait d'évaluer la probabilité), (3) la non négativité agit pour régulariser le problème de régression et aide grandement à obtenir des estimations stables (par exemple vous n'obtenez généralement pas les problèmes de surajustement comme avec la régression ordinaire sans contrainte,les contraintes de non négativité entraînent des estimations plus clairsemées qui sont souvent plus proches de la vérité du terrain; pour le problème de déconvolution ci-dessous, par exemple, les performances sont à peu près aussi bonnes que la régularisation LASSO, mais sans qu'il soit nécessaire de régler un paramètre de régularisation. (La régression pénalisée L0-pseudonorm fonctionne toujours légèrement mieux mais à un coût de calcul plus élevé )

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