CODEWITHSUNDEEP
rglm
jcken
jcken

Reputation: 495

Recovering `init.theta` from `glm.nb()`

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

Upvotes: 0

Views: 59

Answers (1)

AkselA
AkselA

Reputation: 8846

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)

Upvotes: 0

Related Questions