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

12

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?μ^1μ^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é?

StatsStudent
la source
2
C'est une bonne question. Compte tenu de la façon dont il est énoncé, cependant, il pourrait être utile de savoir que lorsque vous avez un seul facteur binaire et aucune autre variable, le lien que vous choisissez ne fait aucune différence.
whuber
1
Merci, @whuber. J'ai mis à jour ma question pour indiquer clairement qu'il existe d'autres covariables dans le modèle. J'ai également ajouté une section "MISE À JOUR" qui explique un peu plus ma question.
StatsStudent
1
Pour un point de vue différent concernant le rôle des fonctions de lien, voir ma réponse à la question étroitement liée à stats.stackexchange.com/questions/63978 .
whuber
1
Exemple fascinant @whuber!
StatsStudent
1
Habituellement, je dirais que le choix de la fonction de liaison est dicté par le problème et les données disponibles - voir ci-dessous pour un exemple concret ...
Tom Wenseleers

Réponses:

4

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 .

Tim
la source
1
Vous mentionnez que les modèles de probabilité linéaires sont "exotiques" mais d'après mes interactions avec les économistes (je suis moi-même statisticien), il semble qu'il y ait deux camps, dont l'un soutient que la probabilité linéaire est meilleure car elle implique moins d'hypothèses et modélise directement les attentes , c'est ce qui compte normalement.
zipzapboing le
1
J'ai mis ma réponse en garde en faisant référence à des cas exotiques en économie. Cela dit, le problème avec le modèle de probabilité linéaire est que si vous l'estimez via OLS, ses hypothèses sont généralement violées. L'hypothèse selon laquelle le modèle est linéaire dans les paramètres n'est pas plausible dans de nombreux cas (c.-à-d., Lorsqu'il est estimé en utilisant OLS, vous pouvez obtenir des probabilités en dehors de 0 et 1). Et, les résidus ne peuvent pas être à distance proches de la normale, vous devez donc utiliser un estimateur sandwich ou quelque chose.
Tim
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 fonction nnpois dans le package R addreg ou en utilisant la fonction nnlm dans le package NNLM .
Tom Wenseleers
0

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

entrez la description de l'image ici

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

entrez la description de l'image ici

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

entrez la description de l'image ici

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

entrez la description de l'image ici

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

entrez la description de l'image ici

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

entrez la description de l'image ici

Tom Wenseleers
la source
1
Y1Y
1
@whuber J'ai maintenant ajouté un exemple concret pour mieux faire comprendre mon point! Toute réflexion sur mon utilisation des moindres carrés non négatifs pondérés pour approximer un modèle de Poisson à lien d'identité non négatif réel serait également la bienvenue!
Tom Wenseleers
Btw - les nnls pondérés que j'utilise pour approximer un GLM de Poisson à lien d'identité non négatif correspondent en fait à l'utilisation d'une seule itération d'un algo de moindres carrés non négatifs repondéré de manière itérative pour s'adapter à un GLM de Poisson non négatif (R lui-même utilise 1 / (y + 0,1) dans Poisson GLM correspond à l'initialisation)
Tom Wenseleers