glm.nb() automatically calculates an init.theta when it is not user specified.
I cannot however seem to recover the user specified init.theta. When querying the model object, I obtain init.theta as 5.028696. This coincides with the value that is used when a user does not specify an init.theta.
How can I recover the value specified by myself? In this case, it is 10^6.
set.seed(54321)
n <- 500
initial_theta <- 10^6
input_data <- data.frame(
group = sample(letters[1:4], size = n, replace = TRUE),
y = rnbinom(n, size = 5, mu = 2)
)
user_model <- MASS::glm.nb(
formula = y ~ group,
data = input_data,
init.theta = initial_theta
)
auto_model <- MASS::glm.nb(
formula = y ~ group,
data = input_data
)
user_model$call$init.theta
#> [1] 5.028696
auto_model$call$init.theta
#> [1] 5.028696
Created on 2024-07-16 with reprex v2.1.0
Looking at the source code, this does not seem to be possible, the init.theta argument is simply overwritten:
Call$init.theta <- signif(as.vector(th), 10)
Where th is the initially estimated theta, as provided by theta.ml()
init.theta is used for the initial fit, and it’s possible to trace the whole iterative process, but annoyingly, what’s returned there as Initial value for 'theta': is again th. This seems odd to me, but I’m no expert. If there are no-one else here who can tell if this is intentional, and in that case why, you might give the R-devel mailing list a shot.
user_model2 <- MASS::glm.nb(
formula = y ~ group,
data = input_data,
init.theta = initial_theta,
control = glm.control(trace=2)
)
# Initial fit:
# Deviance = 870.37355 Iterations - 1
# Deviance = 756.80908 Iterations - 2
# Deviance = 754.42813 Iterations - 3
# Deviance = 754.42665 Iterations - 4
# Deviance = 754.42665 Iterations - 5
# Initial value for 'theta': 4.929440
# Deviance = 578.94259 Iterations - 1
# Theta(1) = 4.929440, 2(Ls - Lm) = 578.943000
(You can set trace=3 for even more verbosity, but it doesn’t help us here)