Skip to content

How to obtain individual log-likelihood value from glm?

8 messages · Bert Gunter, Peter Dalgaard, John Smith +1 more

John Smith
# 
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))
Bert Gunter
# 
If you look at

stats:::logLik.glm  #3 ":" because it's unexported, as is true of most
methods

it should be obvious.

Cheers,
Bert

Bert Gunter

"The trouble with having an open mind is that people keep coming along and
sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
On Tue, Aug 25, 2020 at 8:34 AM John Smith <jswhct at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Peter Dalgaard
    

    

      
      #
    

  


  

      
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

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    
2 days later

    

      
      

        


  

        

      John Smith
    

    

      
      #
    

  


  

      
Thanks Peter for a very promising tip.

        
On Tue, Aug 25, 2020 at 11:40 AM peter dalgaard <pdalgd at gmail.com> wrote:

        

            

      
      
      
    

        

      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      John Smith
    

    

      
      #
    

  


  

      
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 at gmail.com> wrote:

        

            

      
      
      
    

        

      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Mark Leeds
    

    

      
      #
    

  


  

      
Hii: It's been a long time but John Fox's "Companion to Appied Regression"
book has the expressions
for the likelihood of the binomial glm. ( and probably the others also ).
Just running logLik is not so useful
because it could be leaving out multiplicative factors. If you can get your
hands on any edition of John's book,
I remember it being quite helpful in terms of providing all of the gory
details and  for understanding GLM's in
general.

        
On Fri, Aug 28, 2020 at 9:28 PM John Smith <jswhct at gmail.com> wrote:

        

            

      
      
      
    

        

      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Peter Dalgaard
    

    

      
      #
    

  


  

          

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      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:

            

      
      
      
    

        
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