Plusy i minusy połączenia logarytmicznego a identyfikacyjnego dla regresji Poissona

Ja przeprowadzenia regresji Poissona z celem końcowym porównywania (i przy różnicy) prognozowane średnie zliczenie pomiędzy dwoma poziomami czynnika w moim , trzymając inny model współzmiennych (które są binarne ) stała. Zastanawiałem się, czy ktokolwiek mógłby udzielić praktycznych porad, kiedy użyć linku dziennika zamiast linku tożsamości. Jakie są zalety tych dwóch różnych funkcji łącza w regresji Poissona, biorąc pod uwagę mój cel polegający na porównaniu różnic?$\hat{\mu}_1-\hat{\mu}_2$

Mam również na uwadze ten sam cel regresji logistycznej / dwumianowej (użycie linku logit lub linku tożsamości) w celu porównania różnicy proporcji między dwoma poziomami czynników i potrzebuję podobnych porad. Przeczytałem niektóre posty tutaj, które dotyczą tego problemu, ale żaden nie wydaje się wyjaśniać, dlaczego lub kiedy można wybrać jeden link nad drugim i jakie mogą być zalety / wady. Z góry dziękuje za twoją pomoc!

AKTUALIZACJA:

Zdaję sobie również sprawę, że głównym celem korzystania z niektórych funkcji linków jest ograniczenie zakresu możliwych przewidywanych wartości, aby mieściły się w zakresie średniej odpowiedzi (np. W przypadku logistyki zakres jest ograniczony od 0 do 1, a dla logu link, prognozy są ograniczone do liczb dodatnich). Myślę więc, że pytam o to, że jeśli użyję linku tożsamości do powiedzenia regresji logistycznej / dwumianowej, a moje wyniki mieszczą się w zakresie (0,1), czy naprawdę istnieje potrzeba użycia funkcji linku logistycznego lub czy mogę po prostu uprościć myślenie i użyć linku tożsamości?

Wady łącza tożsamości w przypadku regresji Poissona to:

• Jak już wspomniałeś, może generować prognozy spoza zakresu.
• Podczas próby dopasowania modelu mogą pojawić się dziwne błędy i ostrzeżenia, ponieważ łącze pozwala, aby lambda była mniejsza niż 0, ale rozkład Poissona nie jest zdefiniowany dla takich wartości.
• Ponieważ regresja Poissona zakłada, że ​​średnia i wariancja są takie same, po zmianie łącza zmieniasz również założenia dotyczące wariancji. Z mojego doświadczenia wynika, że ​​ten ostatni punkt jest najbardziej wymowny.

Ale ostatecznie jest to pytanie empiryczne. Dopasuj oba modele. Wykonuj dowolne kontrole. Jeśli link tożsamości ma niższy AIC i działa równie dobrze lub lepiej na wszystkich innych czekach, uruchom go z linkiem tożsamości.

W przypadku modelu logit vs liniowy model prawdopodobieństwa (tj. Tak zwany link tożsamości) sytuacja jest o wiele prostsza. Poza kilkoma bardzo egzotycznymi przypadkami w ekonometrii (które znajdziesz, jeśli przeprowadzasz wyszukiwanie), model logit jest lepszy: robi mniej założeń i jest tym, czego używa większość ludzi. Zastosowanie liniowego modelu prawdopodobieństwa na jego miejscu byłby bliski przewrotności.

Jeśli chodzi o interpretację modeli, jeśli używasz R, istnieją dwa świetne pakiety , które wykonają wszystkie ciężkie podnoszenie: efekty , które są bardzo łatwe w użyciu i zelig , który jest trudniejszy w użyciu, ale świetny, jeśli chcesz przewidywać .

W przypadku modeli Poissona powiedziałbym również, że aplikacja często dyktuje, czy twoje zmienne towarzyszące działałyby addytywnie (co oznaczałoby wówczas link tożsamości), czy multiplikacyjnie w skali liniowej (co oznaczałoby link logu). Ale modele Poissona z łączem tożsamości również zwykle mają sens i można je stabilnie dopasować tylko wtedy, gdy nałoży się ograniczenia nieujemności na dopasowane współczynniki - można tego dokonać za pomocą nnpoisfunkcji w addregpakiecie R lub nnlmfunkcji wNNLMpakiet. Nie zgadzam się więc, że należy dopasować modele Poissona zarówno z identyfikatorem, jak i logiem i zobaczyć, który z nich ma najlepszy AIC i wywnioskować najlepszy model na podstawie czysto statystycznych podstaw - raczej w większości przypadków jest to podyktowane podstawowa struktura problemu, który próbuje się rozwiązać, lub dostępne dane.

Na przykład w chromatografii (analiza GC / MS) często mierzy się nałożony sygnał kilku pików w przybliżeniu Gaussa i ten nałożony sygnał jest mierzony za pomocą multiplikatora elektronów, co oznacza, że ​​mierzony sygnał to liczba jonów, a zatem rozkład Poissona. Ponieważ każdy z pików ma z definicji wysokość dodatnią i działa addytywnie, a hałasem jest Poisson, odpowiedni byłby nieujemny model Poissona z łączem tożsamości, a logarytmiczny model Poissona byłby po prostu błędny. W inżynierii strata Kullbacka-Leiblera jest często stosowana jako funkcja straty dla takich modeli, a minimalizacja tej straty jest równoważna z optymalizacją prawdopodobieństwa nieujemnego modelu Poissona powiązanego z tożsamością (istnieją również inne miary dywergencji / straty, takie jak dywergencja alfa lub beta które mają specjalny przypadek Poissona).

Poniżej znajduje się numeryczny przykład, w tym demonstracja, że ​​zwykły nieograniczony link tożsamości Poisson GLM nie pasuje (z powodu braku ograniczeń nieujemności) oraz kilka szczegółów na temat dopasowania nieujemnych modeli tożsamości Poison za pomocą łączannpois, tutaj w kontekście dekonwolacji zmierzonej superpozycji pików chromatograficznych z szumem Poissona na nich za pomocą pasmowej macierzy kowariancyjnej, która zawiera przesunięte kopie zmierzonego kształtu pojedynczego piku. Nieegatywność jest tutaj ważna z kilku powodów: (1) jest to jedyny realistyczny model dostępnych danych (szczyty tutaj nie mogą mieć wysokości ujemnych), (2) jest to jedyny sposób na stabilne dopasowanie modelu Poissona z łączem tożsamości (ponieważ w przeciwnym razie prognozy dla niektórych wartości zmiennych towarzyszących mogą być ujemne, co nie miałoby sensu i dawałoby problemy numeryczne, gdy ktoś próbowałby ocenić prawdopodobieństwo), (3) nieujemność działa na rzecz uregulowania problemu regresji i znacznie pomaga uzyskać stabilne szacunki (np. zazwyczaj nie występują problemy z nadmiernym dopasowaniem, jak w przypadku zwykłej regresji nieograniczonej,ograniczenia nieujemności skutkują mniejszymi szacunkami, które często są bliższe podstawowej prawdzie; dla poniższego problemu dekonwolucji np. wydajność jest prawie tak dobra jak regularyzacja LASSO, ale bez konieczności dostrajania jakiegokolwiek parametru regularyzacji. ( Regresja karna pseudonorma L0 nadal działa nieco lepiej, ale przy wyższych kosztach obliczeniowych )

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