0

I have a simple poisson model with unobserved heterogeneity eta that is normally distributed. This can be solved by integrating eta out of the log-likelihood. However when i try to do this it produces wrong coefficients and i can't see why. The code is shown below. It should return coefficients for beta and sigma close to 0.018 and 0.013 correspondingly, but it does not.

library("maxLik")
library("pracma")
library("rmutil")

n = 3000
x = rnorm(n, mean=0, sd=2)
beta1 = 0.018
sigma = 0.013
eta = rnorm(n, mean=0, sd= sigma)
mu = exp(x*beta1 + eta)
y = rpois(n,mu)


poisson <- function(y,mu) {
  a = mu^y*exp(-mu)/factorial(y)
  return(a)
}




loglik <- function(param) {
  beta = param[1]
  sigma = param[2]
  f <- function(eta) dnorm(eta, sd = sigma)*poisson(y, exp((x*beta + eta)))
  ll <- int(f, a = -100, b = 100) 
  ll[ll==0] <- 0.00000000001
  logll = log(ll)
  return((sum(logll)))
}
bet = runif(16, min=-1, max = 0)/10
ml <- maxLik(loglik, start = c(beta= 0.1, sigma = 0.1), print.level = 2)
Michael Berend
  • 51
  • 1
  • 3
  • Don't know if this has anything to do with your issue, but I'd strongly advise using your `poisson` function on the log scale. `factorial(y)` will get big enough to lose a lot of precision very quickly. Working on log scale will give you a lot more precision. (And why not just use `dpois` instead of your own function? It has a `log = TRUE` argument.) You can use the `log = TRUE` argument for `dnorm` to put that on the log scale too... – Gregor Thomas Jul 25 '19 at 15:34
  • Thank you for your comments! I didn't use the normal poisson function because i built i for specific data which contains non-integer numbers (building on other research). I do like your suggestion, because i had problems with the factorial with the original data as some observations go up to 160. – Michael Berend Jul 25 '19 at 15:45
  • If you need non-discrete factorials, the [gamma function](https://en.wikipedia.org/wiki/Gamma_function#Relation_to_other_functions) is often used. R has built-in `gamma()` and `lgamma()`, a logged version. – Gregor Thomas Jul 27 '19 at 17:25

0 Answers0