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.

  • $\begingroup$ 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? $\endgroup$ – kjetil b halvorsen Nov 22 '17 at 10:25
  • 1
    $\begingroup$ Your offset should be log(N) rather than N $\endgroup$ – Glen_b Nov 23 '17 at 6:10
  • 2
    $\begingroup$ 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 $\endgroup$ – Glen_b Nov 23 '17 at 6:24
  • $\begingroup$ 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. $\endgroup$ – Brian Nov 24 '17 at 21:59

Your Answer

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

Browse other questions tagged or ask your own question.