[R] How to obtain individual log-likelihood value from glm?

[peter dalgaard]

Briefly, you shouldn't. One way of seeing it is if you switch the model to y~1, you still get logLik==0.

The root cause is the rounding in binomial()$aic:

> binomial()$aic
function (y, n, mu, wt, dev) 
{
    m <- if (any(n > 1)) 
        n
    else wt
    -2 * sum(ifelse(m > 0, (wt/m), 0) * dbinom(round(m * y), 
        round(m), mu, log = TRUE))
}

which, if wt is small enough ends up calculating dbinom(0, 0, p, log=TRUE) which is zero. 

(Not rounding gives you NaN, because you're trying to fit a model with a non-integer number of observations.)

-pd

> On 29 Aug 2020, at 03:28 , John Smith <jswhct using gmail.com> wrote:
> 
> If the weights < 1, then we have different values! See an example below. How  should I interpret logLik value then?
> 
> set.seed(135)
>  y <- c(rep(0, 50), rep(1, 50))
>  x <- rnorm(100)
>  data <- data.frame(cbind(x, y))
>  weights <- c(rep(1, 50), rep(2, 50))
>  fit <- glm(y~x, data, family=binomial(), weights/10)
>  res.dev <- residuals(fit, type="deviance")
>  res2 <- -0.5*res.dev^2
>  cat("loglikelihood value", logLik(fit), sum(res2), "\n")
> 
> On Tue, Aug 25, 2020 at 11:40 AM peter dalgaard <pdalgd using gmail.com> wrote:
> If you don't worry too much about an additive constant, then half the negative squared deviance residuals should do. (Not quite sure how weights factor in. Looks like they are accounted for.)
> 
> -pd
> 
> > On 25 Aug 2020, at 17:33 , John Smith <jswhct using gmail.com> wrote:
> > 
> > Dear R-help,
> > 
> > The function logLik can be used to obtain the maximum log-likelihood value
> > from a glm object. This is an aggregated value, a summation of individual
> > log-likelihood values. How do I obtain individual values? In the following
> > example, I would expect 9 numbers since the response has length 9. I could
> > write a function to compute the values, but there are lots of
> > family members in glm, and I am trying not to reinvent wheels. Thanks!
> > 
> > counts <- c(18,17,15,20,10,20,25,13,12)
> >     outcome <- gl(3,1,9)
> >     treatment <- gl(3,3)
> >     data.frame(treatment, outcome, counts) # showing data
> >     glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
> >     (ll <- logLik(glm.D93))
> > 
> >       [[alternative HTML version deleted]]
> > 
> > ______________________________________________
> > R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> 
> -- 
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Office: A 4.23
> Email: pd.mes using cbs.dk  Priv: PDalgd using gmail.com
> 
> 
> 
> 
> 
> 
> 
> 
> 

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes using cbs.dk  Priv: PDalgd using gmail.com