Why null deviance is different from my manual calculations?

Question

Let's consider this very simple example with Poisson regression:

y <- c(1, 2, 3, 2, 1, 4, 5, 6, 7, 1, 2, 3)
x <- rnorm(length(y))
z <- runif(length(y))
mod <- glm(y ~ x + z, family = poisson)

As we can observe, the null deviance is equal to:

mod$null.deviance
14.10487

However, my manual calculations (based for example on RPubs) shows something different:

# Estimate null model
mod_null <- glm(y~1, family=poisson)
2 * (logLik(mod) - logLik(mod_null))
'log Lik.' 5.812325 (df=3)

Can you please explain to me what am I doing wrong?

Answer 1

There is already a good answer but I thought we could still clarify some points.

Also the Note at the end shows the input we will use. It is the same as in the question except we added a seed to make the data and models reproducible.

Let $L$ and $D$ be the likelihood and the deviance of our model, respectively, and use a subscript of $0$ or $s$ to describe the corresponding quantities for the null or saturated models. Then the code in the question calculates the null deviance, i.e. the deviance of the null model, as $2 (L - L_0)$ but the correct formula is $D_0 = -2 (L_0 - L_s)$

In terms of R code we have the following where, in R, gl(n,1) means a factor with n levels assuming n is the length of y. In the link given in the question it is written as as.factor(1:length(y)) .

n <- length(y)
identical(as.factor(1:length(y)), gl(n, 1))
## [1] TRUE

Now calculating the saturated model and using the correct formula from above we have:

mod <- glm(y ~ x + z, family = poisson)   # model - same as in question
mod_null <- glm(y~1, family = poisson)    # null model - same as in question
n <- length(y)
mod_sat <- glm(y ~ gl(n, 1), family = poisson) # saturated model
mod$null.deviance
[1] 14.10487
D0 <- c(- 2 * (logLik(mod_null) - logLik(mod_sat))); D0  # null deviance
[1] 14.10487

Note

set.seed(123)
y <- c(1, 2, 3, 2, 1, 4, 5, 6, 7, 1, 2, 3)
x <- rnorm(length(y))
z <- runif(length(y))

Answer 2

Expanding from my comment above, you are not comparing the right quantities. The null deviance is given by 2(LL(Saturated Model) - LL(Null Model)).

Let's calculate the null deviance manually and compare with mod$null.deviance.

get_null_deviance <- function(fit) {
# Get data from fit object and create new `.idx` column for
# saturated model
data &lt;- transform(fit<span class="math-container">$model, .idx = factor(1:nrow(fit$</span>model)))

# Get response variable from fit object
# In this case, this returns `&quot;y&quot;`
resp &lt;- all.vars(attr(fit$model, &quot;terms&quot;))[1]

# Formula for saturated model and fit model
# In this case, the saturated model is `y ~ .idx`
fm_saturated &lt;- reformulate(&quot;.idx&quot;, response = resp)
fit_saturated &lt;- glm(fm_saturated, data = data, family = poisson)

# Formula for null model and fit model
# In this case, the null model is `y ~ 1`
fm_null &lt;- reformulate(&quot;1&quot;, response = resp)
fit_null &lt;- glm(fm_null, data = data, family = poisson)

# Return null deviance
2 * (logLik(fit_saturated) - logLik(fit_null))


}
get_null_deviance(mod)
#'log Lik.' 14.10487 (df=12)
mod$null.deviance
#[1] 14.10487