3
$\begingroup$

I have a Poisson model that is performing well. Now we need to put it into Java code and release it to the world. What is the equation that I plug the Poisson coefficients into?

Similar to this question: Find the equation from generalized linear model output

Here is some sample code from R that shows what I am doing:

d = data.frame(y=(10:100)^-1*100, x=10:100)
m = glm(y~x, data=d, family=poisson(link="log"))
plot(d$x, predict.glm(m, type="response"))
points(d$x, d$y, col="blue")
summary(m)

output:

Call:
glm(formula = y ~ x, family = poisson(link = "log"), data = d)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.34781  -0.25933  -0.05075   0.21286   1.20265  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.155676   0.126778  17.004   <2e-16 ***
x           -0.025912   0.002819  -9.191   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 104.5820  on 90  degrees of freedom
Residual deviance:   8.5623  on 89  degrees of freedom
AIC: Inf

Number of Fisher Scoring iterations: 4

enter image description here

> cor(d$y, predict.glm(m, type="response"))
[1] 0.9461024

That is pretty good! :)

$\endgroup$
6
$\begingroup$

The equation is

$$\log(\mu_i) = \beta_0 + \beta_1 x_i$$

where $\mu_i$ is the conditional expectation of $y_i$, $E(y | x)$, $\beta_0$ is the coefficient marked Intercept and $\beta_1$ the coefficient marked x. The $\log$ bit is the link function you specified. Hence to get actual predictions on the scale of your response data $y$, you need to apply the inverse of the link function (anti-log) to the both sides of the equation:

$$\mu_i = \exp(\beta_0 + \beta_1 x_i)$$

$\mu_i$ is then the predicted mean count given the value of x.

Print out the coefficients if you want them:

coef(m)

They may be a little more precise than the ones in the summary() output and so your Java code will be closer to the predictions given by predict().

The model doesn't look fantastic, despite the high linear correlation. There is bias throughout the range of x. Without knowing anything about the data, did you consider a model with $x$ and $x^2$?

$\endgroup$
  • 1
    $\begingroup$ Oh that tiny model is a bogus example. :) $\endgroup$ – Chris Jun 18 '15 at 16:57
  • $\begingroup$ @Chris Ooops; now I see the made up data. $\endgroup$ – Gavin Simpson Jun 18 '15 at 16:59
  • $\begingroup$ So this really is a log-level regression. $\endgroup$ – Chris Jun 19 '15 at 0:37
  • 1
    $\begingroup$ Hmm, I'm not sure. This is not just a model with a log-transformed $y$. The errors are Poisson. $\endgroup$ – Gavin Simpson Jun 19 '15 at 2:34
  • 3
    $\begingroup$ @Chris, you are not log transforming your response data & running a regression. The idea behind the generalized linear model is that you are transforming the parameter that governs the response distribution. For the Poisson distribution, that parameter is often called $\lambda$, but as that is the conditional mean, $\mu_i$ is equally correct. It may help you to read my answer here: Difference between logit and probit models. $\endgroup$ – gung Jun 19 '15 at 3:34
0
$\begingroup$

This is great, and I have used it to verify part of a plot in which I am modeling the number of trees counted at certain elevations. However, I am using a zero-inflated Poisson model through zeroinfl() in pscl and the predicted values of zeroinfl() begin to deviate from the equation provided by Gavin after about 1000m in elevation. I assume this is because the probability of zero counts increases after about an elevation of 1000m, and so the zeroinfl() model begins to account for this increased zero count probability. Thats great and it does a good job of reflecting the actual observations, but I also need to know the equation for this line. What would be the zero-inflated Poisson version of this model? Something involving probabilities? Below is a plot in which the dots are the actual observation, the red is the predicted values based on zeroinfl() and the blue is the predicted valued based on the equation provided by Gavin and the coefficients of the zeroinfl() model. I am using a logarithmic representation of the y-axis to clarify the deviation. Included is a general representation of the code.

LogMega[,1] <- yjPower(LogMega[,1],0) #This is a log transform of the average density
m1 <- zeroinfl(Avg_dens~Z, data = LogMega) #the zero-inflated model
m2 <- 2.718^(coef(m1)[1]+(coef(m1)[2]*Z), data=LogMega) #the Poisson equation model
newdata1 <- expand.grid(LogMega[,n+3]) #create a table of the zeroinfl() predicted values
colnames(newdata1)<-"pred"
newdata1$resp <- predict(m1,newdata1)
p <- ggplot(LogMega,aes_string(x=names(LogMega)[n+3],y=LogMega[1]))
theme(
  geom_line(data=newdata1, aes(x=pred,y=resp), size=1, color="red")+
  geom_line(data=LogMega, aes(x=LogMega[,n+3], y=m2), size=1, color="blue"))

enter image description here

$\endgroup$
  • $\begingroup$ I seem to have found it through this link. It is a combination of both the Poisson and the zero inflated probability function. A bit long and complicated when its all drawn out but it seems to work. Thanks for the head start. $\endgroup$ – Benjie Jul 21 '15 at 14:51

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.