Fallar en optimizar el modelo binomial negativo utilizando optim.

Estoy intentando optimizar manualmente un modelo de regresión binomial negativa utilizando el paquete optim en R para predecir una variable de conteo y utilizando una matriz de factores X utilizando el siguiente código:

“`{r}

<h1>generando datos ficticios</h1>

n <- 1000

X <- matrix(NA, ncol = 5, nrow = n)
X[,1] <- 1
X[,2] <- sample(size = n, x = c(0,1), replace = TRUE)
X[,3] <- sample(size = n, x = c(0,1), replace = TRUE)
X[,4] <- sample(size = n, x = c(0,1), replace = TRUE)
X[,5] <- sample(size = n, x = c(0,1), replace = TRUE)

beta0 <- 3
beta1 <- -2
beta2 <- -2
beta3 <- -4
beta4 <- -0.9
k <- 0.9

<h1>extrayendo de la distribución binomial negativa</h1>

mu <- exp(beta0 + beta1 * X[,2] + beta2 * X[,3] + beta3 * X[,4] + beta4 * X[,5])
theta <- mu + mu^2 / k

<h1>variable dependiente</h1>

y <- rnegbin(n, mu = mu, theta = theta)

<h1>función a optimizar</h1>

negbin_ll <- function(y, X, theta) {
beta <- theta[1:ncol(X)]
alpha <- theta[ncol(X) + 1]

logll <- y * log(alpha) + y %<em>% t(X) – (y + (1 / alpha)) * log(1 + alpha * exp(beta %</em>% t(X))) + lgamma(y + (1 / alpha)) – lgamma(y + 1) – lgamma(1 / alpha)

logll <- sum(logll)

return(logll)
}

stval <- rep(0, ncol(X) + 1)

res <- optim(
stval,
negbin_ll,
y = y,
X = X,
control = list(fnscale = -1),
hessian = TRUE,
method = "BFGS"
)
“`

El código debería producir estimaciones de punto del proceso de optimización, pero en su lugar falla al ejecutar la función optim con el mensaje de error “initial value in ‘vmmin’ is not finite”.

Ya intenté cambiar log(gamma(…)) a lgamma(…) en la función de verosimilitud y probé muchas otras formas, pero no logro obtener estimaciones.

Cambiar los valores de inicio de optim tampoco ayuda.

¿Tienes alguna idea si hay alguna particularidad en la función de verosimilitud que cause que los valores se traten de manera extraña?

Se agradecería la ayuda.