Плюсы и минусы лог-ссылки в сравнении с идентификационной ссылкой для регрессии Пуассона

Я выполняю регрессию Пуассона с конечной целью сравнения (и взятия разности) предсказанного среднего значения между двумя уровнями фактора в моей модели: , удерживая другие модельные ковариаты (которые являются двоичными) постоянными. Мне было интересно, если кто-нибудь может дать несколько практических советов о том, когда использовать ссылку на журнал, а не ссылку на личность. Каковы плюсы и минусы этих двух различных функций связи в регрессии Пуассона, учитывая мою цель сравнения различий?$\hat{\mu}_1-\hat{\mu}_2$

Я также имею в виду ту же цель для логистической / биномиальной регрессии (использовать ссылку логита или идентификационную ссылку) для сравнения различий в пропорциях между двумя уровнями фактора и нуждаюсь в аналогичном совете. Я читал некоторые из постов здесь, которые касаются этой проблемы, но ни один, кажется, не объясняет, почему или когда можно было бы выбрать одну ссылку над другой и каковы могут быть плюсы / минусы. Заранее спасибо за помощь!

ОБНОВИТЬ:

Я также понимаю, что основная цель использования определенных функций связей состоит в том, чтобы ограничить диапазон возможных прогнозируемых значений диапазоном среднего отклика (например, для логистики диапазон ограничен значениями от 0 до 1 и для журнала). ссылка, прогнозы ограничены, чтобы быть положительными числами). Итак, я предполагаю, что я спрашиваю, что если я использую идентификационную ссылку, скажем, для логистической / биномиальной регрессии, и мои результаты находятся в диапазоне (0,1), действительно ли есть необходимость использовать функцию логистической ссылки или Могу ли я просто сделать мысли проще использовать идентификационную ссылку?

Минусы идентификационной ссылки в случае регрессии Пуассона:

• Как вы упомянули, он может давать прогнозы вне диапазона.
• Вы можете получить странные ошибки и предупреждения при попытке подгонки модели, потому что ссылка позволяет лямбда быть меньше 0, но распределение Пуассона не определено для таких значений.
• Поскольку пуассоновская регрессия предполагает, что среднее значение и дисперсия одинаковы, при изменении ссылки вы также меняете предположения относительно дисперсии. Мой опыт показывает, что этот последний момент наиболее показателен.

Но, в конечном счете, это эмпирический вопрос. Подходят обе модели. Выполните все проверки, которые вам нравятся. Если идентификационная ссылка имеет более низкий AIC, а также хорошо или лучше подходит для всех остальных ваших проверок, запустите идентификационную ссылку.

В случае логит-модели против линейной вероятностной модели (то есть того, что вы называете идентификационной связью), ситуация намного проще. За исключением некоторых очень экзотических случаев в эконометрике (которые вы найдете, если выполните поиск), модель логита лучше: она делает меньше предположений и используется большинством людей. Использование линейной вероятностной модели на ее месте могло бы оказаться извращенным.

Что касается интерпретации моделей, если вы используете R, есть два замечательных пакета, которые сделают всю тяжелую работу: эффекты , которые очень просты в использовании, и zelig , которые сложнее в использовании, но хороши, если вы хотите делать прогнозы ,

В случае моделей Пуассона я бы также сказал, что приложение часто диктует, будут ли ваши ковариаты действовать аддитивно (что затем будет означать идентификационную ссылку) или мультипликативно в линейном масштабе (что затем будет подразумевать лог-ссылку). Но модели Пуассона с тождественной связью также, как правило, имеют смысл и могут стабильно подходить только в том случае, если на навязанные коэффициенты накладываются ограничения неотрицательности - это можно сделать с помощью nnpoisфункции в addregпакете R или с помощью nnlmфункции вNNLMпакет. Поэтому я не согласен с тем, что модели Пуассона следует согласовывать как с идентификационной, так и с журнальной связью, и посмотреть, какая из них в итоге будет иметь лучший AIC, и вывести лучшую модель на основе чисто статистических соображений - скорее, в большинстве случаев это диктуется основная структура проблемы, которую каждый пытается решить, или данные под рукой.

Например, в хроматографии (анализ ГХ / МС) часто измеряют наложенный сигнал из нескольких пиков приблизительно гауссовой формы, и этот наложенный сигнал измеряется с помощью умножителя электронов, что означает, что измеренный сигнал представляет собой счетчик ионов и, следовательно, распределение Пуассона. Поскольку каждый из пиков по определению имеет положительную высоту и действует аддитивно, а шумом является Пуассон, здесь будет уместна неотрицательная модель Пуассона с тождественной связью, а логарифмическая связь - модель Пуассона будет совершенно неверной. В разработке потеря Кулбека-Лейблера часто используется в качестве функции потерь для таких моделей, и минимизация этой потери эквивалентна оптимизации вероятности неотрицательной модели Пуассона с тождественным звеном (есть также другие меры расхождения / потери, такие как расхождение альфа или бета что есть пуассон как частный случай).

Ниже приведен числовой пример, включающий демонстрацию того, что обычная неограниченная тождественная ссылка Пуассона GLM не подходит (из-за отсутствия ограничений неотрицательности), и некоторые подробности о том, как подобрать неотрицательные модели Пуассона с тождественной связью, используяnnpoisздесь, в контексте деконволюции измеренной суперпозиции хроматографических пиков с пуассоновским шумом на них, используя полосчатую ковариатную матрицу, которая содержит сдвинутые копии измеренной формы одного пика. Неотрицательность здесь важна по нескольким причинам: (1) это единственная реалистичная модель для имеющихся данных (пики здесь не могут иметь отрицательную высоту), (2) это единственный способ стабильно согласовать модель Пуассона с тождественной связью (как в противном случае предсказания для некоторых ковариатных значений могут стать отрицательными, что не имеет смысла и приведет к численным проблемам, когда кто-то попытается оценить вероятность), (3) неотрицательность действует для регуляризации проблемы регрессии и значительно помогает получить стабильные оценки (например, Вы, как правило, не получаете проблем с переоснащением, как с обычной неограниченной регрессией,ограничения неотрицательности приводят к более редким оценкам, которые часто ближе к основной истине; для проблемы деконволюции, приведенной ниже, например, производительность примерно такая же, как и для регуляризации LASSO, но без необходимости настройки какого-либо параметра регуляризации. ( Штрафная регрессия L0-псевдонорм все еще работает немного лучше, но с большими вычислительными затратами )

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