如何使用 optim() 为广义线性模型生成系数估计？

Question

我在R写了如下模型：

model1_b <- glm(bupacts ~ sex + couples + women_alone, data = risky, family = poisson(link="log"))

我想使用optim 查找此模型的系数估计值。也就是说，我想达到上述模型产生的系数与使用optim 产生的系数相匹配的程度。下面是似然函数的代码：

log.like <- function(par) {
    xb <- par[1] + par[2] * risky$sex + par[3] * risky$couples + par[4] * risky$women_alone
   lambda <- exp(xb)
    ll <- sum(log(exp(-lambda) * (lambda ^ risky$bupacts) / factorial(risky$bupacts)))
    return(-ll)}
result <- optim(c(0, 0, 0, 0), log.like, hessian = TRUE, method = "BFGS")

我不断收到错误消息：initial value in 'vmmin' is not finite。有人可以帮我吗？我不确定发生了什么。

Answer 1

使用如下所示的负对数似然函数。

library(foreign)

# input
u <- file.path("http://www.stat.columbia.edu", 
  "~gelman/arm/examples/risky.behavior/risky_behaviors.dta")
risky <-  read.dta(u, convert.factors = TRUE)

# glm
model1_b <- glm(bupacts ~ sex + couples + women_alone, data = risky, 
  family = poisson(link = "log"))
coef(model1_b)
##   (Intercept)        sexman       couples   women_alone 
##  3.1743188102  0.0474967008  0.1448180037 -0.0007948884 

## glm using optim
neg.LL <- function(p) with(risky, {
   eta <- p[1] + p[2] * (sex == "man") + p[3] * couples +  p[4] * women_alone
  -sum(dpois(bupacts, lambda = exp(eta), log = TRUE))
})
fm <- lm(bupacts ~ sex + couples + women_alone, risky)
optim(coef(fm), neg.LL, method = "BFGS")

给予：

$par
  (Intercept)        sexman       couples   women_alone 
 3.1744020275  0.0474937915  0.1447499657 -0.0009630441 

$value
[1] 7083.872

$counts
function gradient 
     178       37 

$convergence
[1] 0

$message
NULL