1
$\begingroup$

I try to understand how is calculated residual deviance after a glm with binomial distribution and logit link: I am not able to reproduce the value that is reported by R (I do not blame R; I am sure that the output is right!). Thanks if someone can help me:

A <- c(10, 12, 15, 0, 1, 2, 3)
B <- c(2, 5, 3, 14, 15, 20, 30)
S <- c(25, 26, 27, 28, 29, 30, 31)

g <- glm(cbind(A=A, B=B) ~ S, 
         family = binomial(link = "logit"))
g$deviance
g$null.deviance

The values are:

> g$deviance
[1] 22.16312
> g$null.deviance
[1] 77.99713

Residual deviance being -2*(LnL - LnLSat), I estimate LnL and LnLSat:

invlogit <- function (n) 1/(1 + exp(-n))
LnL <- sum(pbinom(q=A, size = A+B, 
           prob = invlogit(predict(g)), log.p = TRUE))
LnLSat <- sum(pbinom(q=A, size = A+B, 
           prob = A/(A+B), log.p = TRUE))
-2*(LnL - LnLSat)

Error here; it was obviously dbinom() not pbinom() that should be used !

The correct code is :

invlogit <- function (n) 1/(1 + exp(-n))
LnL <- sum(dbinom(x=A, size = A+B, 
           prob = invlogit(predict(g)), log = TRUE))
LnLSat <- sum(dbinom(x=A, size = A+B, 
           prob = A/(A+B), log = TRUE))
-2*(LnL - LnLSat)

All is ok now

$\endgroup$
2
  • $\begingroup$ It's usually helpful to look at the source code for the whichever function you're using. $\endgroup$ Aug 27, 2019 at 13:46
  • $\begingroup$ I found !!! it was just stupid error !I used pbinom() whereas it should be dbinom()! And it works perfectly now... sorry. I was stuck during 3 hours about this. $\endgroup$
    – MarcG
    Aug 27, 2019 at 14:18

1 Answer 1

0
$\begingroup$

I made a mistake: pbinom() is for the cumulative distribution function whereas I should use dbinom() for the density/mass function.

And all is ok now !

I do not delete the question just in case someone did the same mistake as me !

$\endgroup$
1
  • $\begingroup$ Don't beat yourself up about this. Why not develop a full answer that explains (conceptually) how the deviances are calculated, & illustrates that w/ R code that's made as transparent as possible (not everyone will use R, or be able to read it well)? $\endgroup$ Aug 27, 2019 at 15:13

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.