未能使用 optim 优化负二项式模型

Question

我正在尝试使用 R 中的 optim 包手动优化负二项式回归模型，尝试使用以下代码使用因子矩阵 X 预测计数变量 y：


# generating some fake data

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 

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

# dependent variable

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



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

  logll <- y * log(alpha) + y *( beta %*% t(X) ) - (y + (1 / alpha ) ) * log( 1 + alpha * exp(beta %*% 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"
  )

代码应该从优化过程中产生点估计，但在使用 error in optim(stval, negbin_ll, y = y, X = X, control = list(fnscale = -1), : initial value in 'vmmin' is not finite.

执行优化函数时却失败了

我已经尝试将似然函数中的 log(gamma(...)) 更改为 lgamma(...) 并尝试了许多其他方法，但我无法获得估计值。

更改 optim 的起始值也无济于事。

您是否知道似然函数是否有任何特殊性导致以任何奇怪的方式处理值？

不胜感激。

Answer 1

optim 尝试了几个点以达到最小值，在您的情况下，它在日志中的参数中遇到了一些非正值。一种方法是通过 return 一个负数（在你的情况下）大数，如 -lenght(series)*10^6 来丢弃问题函数中 return 任何非正数的值。重新制作了对数似然函数，就像这样它有点管用：

negbin_ll <- function(y, X, theta){
  
  beta <- theta[1:ncol(X)]
  alpha <- theta[ncol(X) + 1]
  
  if(any(alpha<=0)) return(-length(y)*10^6)
  if(any(1 + alpha * exp(beta %*% t(X))<=0)) return(-length(y)*10^6)
  
  logll <- y * log(alpha) + y *( beta %*% t(X) ) - (y + (1 / alpha ) ) * log( 1 + alpha * exp(beta %*% t(X))) + lgamma(y + (1 / alpha)) - lgamma ( y + 1)  - lgamma ( 1 / alpha)
  
  
  logll <- sum( logll  )
  
  return(logll)
  
}

未能使用 optim 优化负二项式模型

Failing to optimise negative binomial model using optim

optimization

r

glm