Non-linear Poisson Regression

I am trying to fit a count regression (rate model) of the form
$y_{r,k}$ ~ $Poisson(N_k t_r^{-\alpha}exp(X_k\beta) + \delta)$.
This looks complicated, so let me explain.

(1) $y_{r,k}$ are the counts that we observe for group k at time r.
(2) $N_r$ is the total "population" in the group r. (all groups have obeservations from time 1 to T, these do not sum to $N_r$)
(3) $t_r$ is the observation time for $y_{r,k}$
(4) $X_k$ is a matrix of covariates and $\beta$ is the coefficients for these covariates
(5) The coefficients that I want to find are $\alpha, \beta, \delta$ (where $\beta$ is a vector)

In the absence of $\delta$ the model is the usual count regression. It can be fit in R as below:
glm(y ~ offset(N) + log(t) + as.factor(X1) + as.factor(X2) + as.factor(X3), family = "poisson")

The reason why $\delta$ is added to the model because when t is very big $N_k t_r^{-\alpha}X_k\beta$ is very close to zero, in which case the pearson residuals will become very big if there is even one count. ($\delta$ is expected to be around 0.02)

However, fitting a model with $\delta$ is not easy, since this is no longer a linear model. More specifically, $log(E(y))$ is not longer linear.

My question is, how can i fit this model?

So far, I have some preliminary ideas. (1) use gnm (i am not sure how to write the model formula for this) (2) fit the model without $\delta$ and then iterate through possible $\delta$ values that improve the deviance. This can be done by simply adding $\delta$ to the predict(model,"type = response") and calculating the deviance for a poisson model. (3) try to find the MLE's simultaneously by building a user defined function and use R's optim function (i suspect this is what gnm is doing)

I have already done (2) but am not very satisfied with this approach. How can I do (1) or (3) or some other method? I understand that my post may not be very clear and am willing to provide more details if needed.

• From your explanation it seems the problems are caused by bad fit when t islarge. Thos m8ght be becaus3 of nonlinearities. Maybe investigate that rather than adhoc modify the model? – kjetil b halvorsen Nov 22 '17 at 10:25
• Your offset should be log(N) rather than N – Glen_b Nov 23 '17 at 6:10
• reading the gnm overview vignette, I think you can write the formula for gnm with an identity link as follows: X1f <- as.factor(X1); X2f <- as.factor(X2); X3f <- as.factor(X3) (to make the formula compact enough to read) -- then the identity link version of the formula should be -- Mult(Exp(offset(log(N))+log(t)), X1f+X2f+X3f) +1 – Glen_b Nov 23 '17 at 6:24
• Tks glen you are right about the offest part. and tks a lot for the gnm, I was getting confused on how to write this formula. – Brian Nov 24 '17 at 21:59