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](https://rpubs.com/fhernanb/deviance_glm)) 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: 
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))

```
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`.

```r
get_null_deviance <- function(fit) {

    # Get data from fit object and create new `.idx` column for
    # saturated model
    data <- transform(fit$model, .idx = factor(1:nrow(fit$model)))

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

    # Formula for saturated model and fit model
    # In this case, the saturated model is `y ~ .idx`
    fm_saturated <- reformulate(".idx", response = resp)
    fit_saturated <- 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 <- reformulate("1", response = resp)
    fit_null <- 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

```