如何使用optim()为广义线性模型生成系数估计？
如何使用optim()为广义线性模型生成系数估计？
提问于 2021-04-01 22:57:11
我用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。有人能帮帮我吗？我不知道是怎么回事。
回答 1
Stack Overflow用户
回答已采纳
发布于 2021-04-02 00:05:47
使用如下所示的负对数似然函数。
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页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/66906859
复制相关文章
相似问题
腾讯云开发者
