I am fitting a series of generalized linear models in R using glm. With a simple regression model (family = gaussian), how can I compute the log likelihood of the model on a new dataset? A similar question has been asked here, but in that question a logistic regression model was fitted where the outcome is already a probability which doesn't directly transfer to other model families.

Here's a minimal reproducible example:

set.seed(1)
n <- 200
p <- 120
X <- matrix(rnorm(n * p), ncol = p)
colnames(X) <- paste0("X", 1:p)
beta <- sample((-5):5, p, replace = TRUE)
dat <- data.frame(X,
                  y = 5 + X %*% beta)
model <- glm(y ~ ., family = gaussian(), data = dat)

logLik only returns the log likelihood of the training data. I tried to compute the log likelihood on the training data with

custom_loglik <- function(mod, newdata) {
  sum(dnorm(
    newdata$y,
    predict(mod, newdata, type = "response"),
    sqrt(summary(mod)$dispersion),
    log = TRUE
  ))
}

However,

> logLik(model)
'log Lik.' 5886.64 (df=122)

and

> custom_loglik(model, dat)
[1] 5854.253

After some experiments, the difference seems to get larger with larger p, or more specifically with a larger p/n ratio. I assume that me using the dispersion isn't correct, but I feel like I'm missing something obvious, since I imagine out of sample likelihood to be a rather common thing to need.