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