Vor- und Nachteile der Protokollverknüpfung im Vergleich zur Identitätsverknüpfung für die Poisson-Regression

Ich trage eine Poisson - Regression mit dem Endziel aus zu vergleichen (und die Differenz der Einnahme) die vorhergesagten mittleren Zählungen zwischen zwei Faktorstufen in meinem , während anderes Modell Kovariaten halten (das ist alle binäre ) konstant. Ich habe mich gefragt, ob irgendjemand einen praktischen Rat geben kann, wann ein Protokolllink anstelle eines Identitätslinks zu verwenden ist. Welche Vor- und Nachteile haben diese beiden unterschiedlichen Verknüpfungsfunktionen bei der Poisson-Regression, wenn ich das Ziel habe, Unterschiede zu vergleichen?$\hat{\mu}_1-\hat{\mu}_2$

Ich habe auch das gleiche Ziel für eine logistische / binomiale Regression (Verwendung eines Logit-Links oder eines Identity-Links) im Sinn, um den Unterschied in den Anteilen zwischen zwei Faktorstufen zu vergleichen, und benötige einen ähnlichen Rat. Ich habe hier einige Posts gelesen, die sich mit diesem Thema befassen, aber keiner scheint zu erklären, warum oder wann man einen Link über den anderen wählt und was die Vor- / Nachteile sein könnten. Vielen Dank im Voraus für Ihre Hilfe!

AKTUALISIEREN:

Mir ist auch klar, dass der Hauptzweck der Verwendung bestimmter Verknüpfungsfunktionen darin besteht, den Bereich möglicher vorhergesagter Werte auf den Bereich der mittleren Antwort zu beschränken (z. B. ist der Bereich für die Logistik auf 0 bis 1 und für das Protokoll beschränkt Link, die Vorhersagen sind auf positive Zahlen beschränkt). Ich frage mich also, ob ich, wenn ich einen Identitätslink für eine logistische / binomiale Regression verwende und meine Ergebnisse innerhalb des Bereichs (0,1) liegen, wirklich eine logistische Verknüpfungsfunktion verwenden muss oder nicht Könnte ich einfach das Denken vereinfachen und einen Identitätslink verwenden?

Nachteile eines Identitätslinks bei der Poisson-Regression sind:

• Wie Sie bereits erwähnt haben, können Vorhersagen außerhalb des Bereichs erstellt werden.
• Beim Versuch, das Modell anzupassen, kann es zu seltsamen Fehlern und Warnungen kommen, da der Link zulässt, dass Lambda kleiner als 0 ist, aber die Poisson-Verteilung für solche Werte nicht definiert ist.
• Da die Poisson-Regression davon ausgeht, dass Mittelwert und Varianz gleich sind, ändern Sie beim Ändern des Links auch die Annahmen über die Varianz. Ich habe die Erfahrung gemacht, dass dieser letzte Punkt am aussagekräftigsten ist.

Letztendlich ist dies jedoch eine empirische Frage. Passen Sie beide Modelle. Führen Sie nach Belieben Überprüfungen durch. Wenn der Identitätslink einen niedrigeren AIC hat und bei allen anderen Überprüfungen mindestens genauso gut funktioniert, führen Sie ihn mit dem Identitätslink aus.

Im Fall des Logit-Modells im Vergleich zum linearen Wahrscheinlichkeitsmodell (dh, was Sie als Identitätsverknüpfung bezeichnen) ist die Situation viel einfacher. Abgesehen von einigen sehr exotischen Fällen in der Ökonometrie (die Sie bei einer Suche finden), ist das logit-Modell besser: Es macht weniger Annahmen und wird von den meisten Menschen verwendet. Die Verwendung des linearen Wahrscheinlichkeitsmodells an seiner Stelle wäre fast pervers.

Was die Interpretation der Modelle angeht , gibt es bei Verwendung von R zwei großartige Pakete, die das ganze schwere Heben erledigen: Effekte , die sehr einfach zu verwenden sind, und zelig , die schwieriger zu verwenden sind, aber großartig, wenn Sie Vorhersagen treffen möchten .

Bei Poisson-Modellen würde ich auch sagen, dass die Anwendung häufig vorschreibt, ob Ihre Kovariaten additiv (was dann eine Identitätsverknüpfung impliziert) oder multiplikativ auf einer linearen Skala (was dann eine Protokollverknüpfung impliziert) wirken. Aber Poisson - Modelle mit einer Identitätsverknüpfung sind normalerweise nur sinnvoll und können nur dann stabil angepasst werden, wenn man den angepassten Koeffizienten Nicht - Negativitätsbeschränkungen auferlegt - dies kann über die nnpoisFunktion im R - addregPaket oder über die nnlmFunktion im R - Paket erfolgenNNLMPaket. Daher stimme ich nicht zu, dass man Poisson-Modelle sowohl mit einer Identität als auch mit einem Log-Link ausstatten sollte und sieht, welches das beste AIC hat und welches aus rein statistischen Gründen das beste Modell ableitet Grundstruktur des Problems, das man zu lösen versucht, oder der vorliegenden Daten.

Beispielsweise würde man in der Chromatographie (GC / MS-Analyse) häufig das überlagerte Signal mehrerer ungefähr Gauß-förmiger Peaks messen und dieses überlagerte Signal wird mit einem Elektronenvervielfacher gemessen, was bedeutet, dass die gemessenen Signale Ionenzahlen und daher Poisson-verteilt sind. Da jeder der Peaks per Definition eine positive Höhe hat und additiv wirkt und das Rauschen Poisson ist, wäre hier ein nichtnegatives Poisson-Modell mit Identitätsverknüpfung angebracht, und ein logarithmisches Poisson-Modell wäre einfach falsch. In der Technik wird Kullback-Leibler-Verlust häufig als Verlustfunktion für solche Modelle verwendet, und die Minimierung dieses Verlusts entspricht der Optimierung der Wahrscheinlichkeit eines nichtnegativen Poisson-Modells mit Identitätsverknüpfung (es gibt auch andere Divergenz- / Verlustmaße wie Alpha- oder Betadivergenz) die Poisson als Sonderfall haben).

Nachfolgend finden Sie ein numerisches Beispiel, einschließlich einer Demonstration, dass ein reguläres, nicht eingeschränktes Identitätslink-Poisson-GLM nicht passt (aufgrund des Fehlens von Nicht-Negativitätsbeschränkungen), sowie einige Details zur Anpassung nicht negativer Identitätslink-Poisson-Modellennpoishier im Zusammenhang mit der Entfaltung einer gemessenen Überlagerung von chromatographischen Peaks mit Poisson-Rauschen unter Verwendung einer bandierten Kovariatenmatrix, die verschobene Kopien der gemessenen Form eines einzelnen Peaks enthält. Nicht-Negativität ist hier aus mehreren Gründen wichtig: (1) Es ist das einzige realistische Modell für die vorliegenden Daten (Spitzen können hier keine negativen Höhen haben). (2) Es ist die einzige Möglichkeit, ein Poisson-Modell mit Identitätsverknüpfung (as) stabil anzupassen andernfalls könnten Vorhersagen für einige kovariate Werte negativ werden, was keinen Sinn macht und numerische Probleme ergeben würde, wenn man versuchen würde, die Wahrscheinlichkeit zu bewerten In der Regel treten bei Ihnen keine Überanpassungsprobleme auf, wie dies bei einer normalen, nicht eingeschränkten Regression der Fall ist.Nicht-Negativitätsbeschränkungen führen zu spärlicheren Schätzungen, die häufig näher an der Grundwahrheit liegen; B. ist die Leistung ungefähr so ​​gut wie die LASSO-Regularisierung, ohne dass ein Regularisierungsparameter eingestellt werden muss. (Die mit L0-Pseudonorm bestrafte Regression schneidet zwar immer noch etwas besser ab, ist jedoch mit einem höheren Rechenaufwand verbunden. )

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