# Poisson regression predictions

I've started looking at GLMs and I've worked out point estimates for Poisson regression using the canonical exponential link function. So the likelihood being

$$\ell(y_i \vert x_i, \beta) = \sum_{i=1}^N \big [ y_i\beta^Tx_i - \exp(\beta^Tx_i) - \log \Gamma (x_i + 1) \big ]$$

$$\frac{\partial\ell}{\partial \beta} = \sum_{i=1}^N \big [ y_i - \exp(\beta^Tx_i) \big] x_i$$

I generated some data using the following process

xs = np.random.uniform(low, high, size=100)
means = intercept + slope * xs
ys = np.array([stats.poisson(mu).rvs() for mu in means])


and plotting that on a graph gives the following

Very straightforward so far. So the further left we move on the $$x$$ axis, the larger the variance.

However, when do maximum likelihood estimation using the previously derived gradients, I get the following picture.

The blue dots are the true generated data and the green dots represent the predictions of the model i.e. $$\exp(\beta^T x_i)$$. It definitely seems to be fitting the data properly, but what threw me off what the exponential curve of the predictions. Given the canonical link function being the exponential, I suppose it makes sense, but it's not really what I expected.

I expected the predictions to follow the straight line, be it a bit shifted up or down to accommodate the sampling noise. Was I wrong to assume this? Is there a way to achieve the model following a straight line? I am aware that in this particular case, Gaussian linear regression would probably do just fine, but that's not my question. Have I generated the data in such a way that the model couldn't fit it?

Code used:

slope = 2
intercept = 10
low, high = 0, 500
boundaries = np.array([low, high])

xs = np.random.uniform(low, high, size=100)
means = intercept + slope * xs
ys = np.array([stats.poisson(mu).rvs() for mu in means])

xs_ = np.vstack((np.ones_like(xs), xs))

def neg_log_likelihood(beta, xs, ys):
return -np.sum(ys * (beta @ xs) - np.exp(beta @ xs))

return -np.sum((ys - np.exp(beta @ xs)) * xs, axis=1)

beta = np.random.normal(0, 0.01, 2)

beta_prime = optimize.fmin_bfgs(neg_log_likelihood, beta, fprime=grad, args=(xs_, ys))

• Perhaps just $x$ and $exp(β^Tx)$ would do. May 30 '19 at 18:46
• The canonical link function in a Poisson GLM is the log link. If the relationships is linear, then it won't make sense to use a non-linear link function, where the $\exp$ in your model post comes from the inverse of the $\log$ link function. In which case, you might look at the identity link function. May 30 '19 at 18:47
• You aren't using the link function when you generate the data, but the model assumes the link function in fitting the data. May 30 '19 at 19:06
• Note that a "linear model" (more correctly a Poisson GLM with an identity link) would necessarily go out of bounds at some point (you can't have counts <0). Do you really need to fit this manually? Plenty of software already exist to do this easily. There are functions available in R and Python, eg. May 30 '19 at 19:18
• That set-up doesn't sound like a Poisson GLM; there isn't a separate noise term in that model. It sounds more like a heteroscedastic linear model problem, which could be estimated using Generalised least squares, a Gaussian location scale model, quantile regression, or sandwich estimators with a linear regression. May 30 '19 at 20:18

The glm you used to generate data has $$E(Y|x) = \beta_0+\beta_1 x$$ but the one you fitted has $$E(Y|x) = \exp(\beta_0+\beta_1 x)$$.
If you want to fit an identity-link Poisson (i.e. you want the straight-line relationship between $$E(Y|x)$$ and $$x$$), then you must specify an identity-link Poisson model to the glm-fitting function.