Why null deviance is different from my manual calculations? 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?
 A: 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))

A: 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 <- 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

