12
$\begingroup$

I am carrying out a Poisson regression with the end goal of comparing (and taking the difference of) the predicted mean counts between two factor levels in my model: $\hat{\mu}_1-\hat{\mu}_2$, while holding other model covariates (which are all binary) constant. I was wondering if anyone could provide some practical advice on when to use a log link versus an identity link. What are the pros can cons of these two different link functions in Poisson regression, given my goal of comparing differences?

I also have the same goal in mind for a logistic/binomial regression (to use a logit link or an identity link) to compare the difference in proportions between two factor levels and need similar advice. I've read some of the posts here that touch on this issue, but none seem to explain why or when one might chose one link over the other and what the pros/cons might be. Thanks in advance for your help!

UPDATE:

I also realize that the main purpose of using certain links functions is to restrict the range range of possible predicted values to be within the range of the mean response (e.g. for logistic, the range is restricted to be between 0 and 1 and for the log link, the predictions are restricted to be positive numbers). So, I guess what I'm asking is that if I use an identity link for say a logistic/binomial regression, and my results are within the range (0,1), is there really any need to use a logistic link function or could I just make thinks simpler an use an identity link?

$\endgroup$
  • 2
    $\begingroup$ It's a good question. Given how it's stated, though, it might be useful to know that when you have just one binary factor and no other variables then it makes no difference what link you choose. $\endgroup$ – whuber Mar 23 '16 at 21:25
  • 1
    $\begingroup$ Thanks, @whuber. I've updated my question to make it clear that there are other covariates in the model. I've also added an "UPDATE" section that explains my question a bit further. $\endgroup$ – StatsStudent Mar 23 '16 at 21:54
  • 1
    $\begingroup$ For a different point of view concerning the role of link functions, see my answer to the closely related question at stats.stackexchange.com/questions/63978 . $\endgroup$ – whuber Mar 23 '16 at 22:27
  • 1
    $\begingroup$ Fascinating example @whuber! $\endgroup$ – StatsStudent Mar 23 '16 at 22:45
  • 1
    $\begingroup$ Usually I would say the choice of link function is dictated by the problem and the data at hand - see below for a concrete example... $\endgroup$ – Tom Wenseleers Apr 28 at 20:00
4
$\begingroup$

Cons of an identity link in the case of the Poisson regression are:

  • As you have mentioned, it can produce out-of-range predictions.
  • You may get weird errors and warnings when attempting to fit the model, because the link permits lambda to be less than 0, but the Poisson distribution is not defined for such values.
  • As Poisson regression assumes that the mean and variance are the same, when you change the link you are also changing assumptions about the variance. My experience has been that this last point is most telling.

But, ultimately this is an empirical question. Fit both models. Perform whatever checks you like. If the identity link has a lower AIC, and does as well or better on all your other checks, then run with the identity link.

In the case of the logit model vs the linear probability model (i.e., what you refer to as the identity link), the situation is a lot more straightforward. Except for some very exotic cases in econometrics (which you will find if you do a search), the logit model is better: it makes fewer assumptions and is what most people use. Using the linear probability model in its place would verge on being perverse.

As regards interpreting the models, if you are using R, there are two great packages that will do all the heavy lifting: effects, which is super easy to use, and zelig, which is harder to use but great if you want to make predictions.

$\endgroup$
  • 1
    $\begingroup$ You mention that linear probability models are "exotic" but from my interactions with economists (I'm a statistician myself) it seems there are two camps, one of which argues that linear probability is better because it involves less assumptions and directly models the expectation, which is what one normally cares about. $\endgroup$ – zipzapboing Feb 26 at 19:16
  • 1
    $\begingroup$ I did caveat my answer by referring to exotic cases in economics. Having said that, the trouble with the linear probability model is that if you estimate it via OLS its assumptions are commonly violated. The assumption that the model is linear in parameters is not plausible in many cases (i.e., when estimated using OLS you can get probabilities outside of 0 and 1). And, the residuals cannot be remotely close to normal, so you then need to use a sandwich estimator or something. $\endgroup$ – Tim Feb 28 at 6:55
  • $\begingroup$ In the case of Poisson models I would also say that the application often dictates whether your covariates would act additively (which would then imply an identity link) or multiplicatively on a linear scale (which would then imply a log link). But Poisson models with an identity link also normally only make sense and can only stably be fit if one imposes nonnegativity constraints on the fitted coefficients - this can be done using the nnpois function in the R addreg package or using the nnlm function in the NNLM package. $\endgroup$ – Tom Wenseleers Apr 26 at 20:52
0
$\begingroup$

In the case of Poisson models I would also say that the application often dictates whether your covariates would act additively (which would then imply an identity link) or multiplicatively on a linear scale (which would then imply a log link). But Poisson models with an identity link also normally only make sense and can only stably be fit if one imposes nonnegativity constraints on the fitted coefficients - this can be done using the nnpois function in the R addreg package or using the nnlm function in the NNLM package. So I don't agree that one should fit Poisson models with both an identity and log link and see which one ends up having the best AIC and infer the best model based on purely statistical grounds - rather, in most cases it is dictated by the underlying structure of the problem that one tries to solve or the data at hand.

For example, in chromatography (GC/MS analysis) one would often measure the superimposed signal of several approx Gaussian shaped peaks and this superimposed signal is measured with an electron multiplier, which means that the measured signal are ion counts and therefore Poisson distributed. Since each of the peaks have by definition a positive height and act additively and the noise is Poisson, a nonnegative Poisson model with identity link would be appropriate here, and a log link Poisson model would be plain wrong. In engineering Kullback-Leibler loss is often used as a loss function for such models, and minimizing this loss is equivalent to optimizing the likelihood of a nonnegative identity-link Poisson model (there is also other divergence/loss measures like alpha or beta divergence that have Poisson as a special case).

Below is a numerical example, including a demonstration that a regular unconstrained identity link Poisson GLM does not fit (because of the lack of nonnegativity constraints) and some details on how to fit nonnegative identity-link Poisson models using nnpois, here in the context of deconvoluting a measured superposition of chromatographic peaks with Poisson noise on them using a banded covariate matrix that contains shifted copies of the measured shape of a single peak. Nonnegativity here is important for several reasons: (1) it is the only realistic model for the data at hand (peaks here cannot have negative heights), (2) it is the only way to stably fit a Poisson model with identity link (as otherwise predictions could for some covariate values go negative, which wouldn't make sense and would give numerical problems when one would try to evaluate the likelihood), (3) nonnegativity acts to regularize the regression problem and greatly helps to obtain stable estimates (e.g. you typically don't get the overfitting problems as with ordinary unconstrained regression, where overfitting typically results in coefficients with very large absolute magnitude and opposite signs - you see this effect in the OLS estimate below) and (4) nonnegativity constraints result in sparser estimates that are frequently closer to the ground truth; for the deconvolution problem below e.g. performance is about as good as LASSO regularization, but without requiring one to tune any regularization parameter. (L0-pseudonorm penalized regression still performs slightly better but at a greater computational cost)

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

$\endgroup$
  • 1
    $\begingroup$ I don't understand the need or even the validity of nonnegativity constraints, which apparently are intended to force all fitted coefficients to be non-negative. After all, if you were simply to reverse the coding of the response (convert $Y$ to $1-Y$), then all coefficients would be negated, demonstrating that the sign of any one of them is without meaning. Am I perhaps misunderstanding what you mean by such constraints? $\endgroup$ – whuber Apr 26 at 22:11
  • 1
    $\begingroup$ @whuber I now added a concrete example to more clearly make my point! Any thoughts on my use of weighted nonnegative least squares to approximate an actual nonnegative identity-link Poisson model would also be welcome! $\endgroup$ – Tom Wenseleers Apr 28 at 0:20
  • $\begingroup$ Btw - the weighted nnls I use to approximate a nonnegative identity-link Poisson GLM in fact corresponds to using a single iteration of an iteratively reweighted nonnegative least squares algo to fit a nonnegative Poisson GLM (R itself uses 1/(y+0.1) in Poisson GLM fits as initialization) $\endgroup$ – Tom Wenseleers May 14 at 11:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.