Convergence comparison
For each of 200 replications we draw a fresh design matrix and
response, then fit with both fastglm_nb() and
MASS::glm.nb(). We record whether each method converges
without error or warning, and when both converge, we compare the
maximized log-likelihoods.
results <- data.frame(
seed = integer(n_reps),
mass_conv = logical(n_reps),
fg_conv = logical(n_reps),
mass_loglik = numeric(n_reps),
fg_loglik = numeric(n_reps),
mass_theta = numeric(n_reps),
fg_theta = numeric(n_reps)
)
for (i in seq_len(n_reps)) {
set.seed(i)
X <- cbind(1, matrix(rnorm(n * (p - 1)), n, p - 1))
mu <- exp(X %*% beta_true)
y <- MASS::rnegbin(n, mu = mu, theta = theta_true)
## MASS::glm.nb
mass_ok <- FALSE
mass_ll <- NA_real_
mass_th <- NA_real_
tryCatch(
withCallingHandlers(
{
m <- MASS::glm.nb(y ~ X[, -1])
if (m$converged) {
mass_ok <- TRUE
mass_ll <- as.numeric(logLik(m))
mass_th <- m$theta
}
},
warning = function(w) invokeRestart("muffleWarning")
),
error = function(e) NULL
)
## fastglm_nb
fg_ok <- FALSE
fg_ll <- NA_real_
fg_th <- NA_real_
tryCatch({
f <- fastglm_nb(X, y, link = "log")
if (f$converged && !any(is.nan(coef(f)))) {
fg_ok <- TRUE
fg_ll <- as.numeric(logLik(f))
fg_th <- f$theta
}
}, error = function(e) NULL)
results$seed[i] <- i
results$mass_conv[i] <- mass_ok
results$fg_conv[i] <- fg_ok
results$mass_loglik[i] <- mass_ll
results$fg_loglik[i] <- fg_ll
results$mass_theta[i] <- mass_th
results$fg_theta[i] <- fg_th
}Summary
tab <- with(results, data.frame(
Method = c("MASS::glm.nb", "fastglm_nb"),
Converged = c(sum(mass_conv), sum(fg_conv)),
Failed = c(sum(!mass_conv), sum(!fg_conv)),
Rate = sprintf("%.1f%%", 100 * c(mean(mass_conv), mean(fg_conv)))
))
knitr::kable(tab, align = "lrrl",
caption = sprintf("Convergence over %d replications", n_reps))| Method | Converged | Failed | Rate |
|---|---|---|---|
| MASS::glm.nb | 165 | 35 | 82.5% |
| fastglm_nb | 200 | 0 | 100.0% |
fastglm_nb() converges on 200 of 200 replications, and
MASS::glm.nb() fails on 35.
Log-likelihood comparison
When both methods converge they should find the same MLE, so their
log-likelihoods should agree. In a small number of replications
MASS::glm.nb() reports convergence with an extreme theta
estimate with \(\hat\theta > 10^6\)
and very poor log-likelihood. We exclude these (identified by
fastglm achieving a strictly higher log-likelihood by more
than 1 unit) from the following agreement plot. These exclusions mean
that the plot describes only the remaining subset, and the excluded fits
must still be counted when evaluating numerical reliability:
both <- results$mass_conv & results$fg_conv
ll_diff <- results$fg_loglik[both] - results$mass_loglik[both]
n_degenerate <- sum(ll_diff > 1)
agree <- both
agree[both] <- ll_diff <= 1
ll_good <- results$fg_loglik[agree] - results$mass_loglik[agree]Both methods converge on 165 of 200 replications.
MASS::glm.nb() reaches a degenerate fit in 3 of these,
which leaves 162 agreeing fits. Among them, the largest absolute
log-likelihood difference is 1.5e-06, and the mean is 1.9e-08.
library(ggplot2)
agreement <- results[agree, ]
ggplot(agreement, aes(mass_loglik, fg_loglik)) +
geom_abline(intercept = 0, slope = 1, colour = "grey40", linetype = 2) +
geom_point(colour = "#0B5CAD", alpha = 0.5, size = 1.2) +
labs(x = "MASS log-likelihood", y = "fastglm log-likelihood",
title = "Log-likelihood agreement") +
theme_bw(base_size = 11) + theme(panel.grid.minor = element_blank())
Theta estimates
ggplot(agreement, aes(mass_theta, fg_theta)) +
geom_abline(intercept = 0, slope = 1, colour = "grey40", linetype = 2) +
geom_point(colour = "#0B5CAD", alpha = 0.5, size = 1.2) +
labs(x = expression(hat(theta) ~ "(MASS)"), y = expression(hat(theta) ~ "(fastglm)"),
title = expression("Dispersion" ~ hat(theta) ~ "agreement")) +
theme_bw(base_size = 11) + theme(panel.grid.minor = element_blank())
The largest absolute theta difference among these fits is 4.3e-07.