2
$\begingroup$

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 ] $$

and the gradient being

$$ \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

enter image description here

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.

enter image description here

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

def grad(beta, xs, ys):
    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))
$\endgroup$
  • 1
    $\begingroup$ Perhaps just $x$ and $exp(β^Tx)$ would do. $\endgroup$ – StatsStudent May 30 at 18:46
  • 4
    $\begingroup$ 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. $\endgroup$ – Reinstate Monica - G. Simpson May 30 at 18:47
  • 1
    $\begingroup$ You aren't using the link function when you generate the data, but the model assumes the link function in fitting the data. $\endgroup$ – gung - Reinstate Monica May 30 at 19:06
  • 2
    $\begingroup$ 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. $\endgroup$ – gung - Reinstate Monica May 30 at 19:18
  • 3
    $\begingroup$ 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. $\endgroup$ – Reinstate Monica - G. Simpson May 30 at 20:18
1
$\begingroup$

As some of the comments have already pointed out, you didn't fit the same model you used to generate the data.

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)$.

You if you want data from a log-link Poisson you need to change the way you generate the data (exponentiate in the step where you compute means).

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.

$\endgroup$

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.