[R] accessing log likelihood of poison model

[Peter Dalgaard]

Renaud Lancelot <renaud.lancelot at cirad.fr> writes:

> >>>And that's great; but I need the log likelihood.
> >>>
> >>>Anyone know?
> >>
> >> The deviance will not suffice?  sum(dpois(nfix, fitted(freq.mod),
> >> log.p=T))
> >>
> >>should do the trick otherwise.

> > Thank you, that did the trick.  I should note that the method
> > doesn't work when I have a noninteger in Y.  I'm doing a "log
> > linear" analysis
> > of a cross table, and one of the cells of the cross table has a zero
> > value.  I put a .5 in that cell.  With the ".5 analysis" the above
> > method doesn't work.  So I had to revert to an analysis with a zero
> > cell to make it work. Which has some problems.  So I'm still left
> > looking for the log likelihood of the ".5 analysis".  I've had success
> > with Stata using this method, but not R.
> > Thanks again.
> 
> ?logLik

Right. Forgot about that. It's the same sum(dpois(....,log=T))
calculation, though (after a detour around the AIC), so that too will
generate -Inf if noninteger values are plugged in. Fair enough: The
likelihood is by definition the probability of obtaining the data, and
the Poisson distribution has probability zero of generating fractional
values. 

However, the deviance is well-defined in such cases. I'd suspect that
what Stata outputs is not the true log-likelihood, but a number
obtained by "analytic continuation" (wrong term, I know) of the
Poisson density to nonintegers, e.g.

> p <- function(x,lambda) lambda^x*exp(-lambda)/gamma(x+1)
> p(1,1)
[1] 0.3678794
> dpois(1,1)
[1] 0.3678794
> dpois(.5,1)
[1] 0
Warning message:
non-integer x = 0.500000
> p(.5,1)
[1] 0.4151075

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907