R:使用 MaxLik() 包估计 mle 时出错

R: error on estimating mle with MaxLik() package

我的目的是使用 Newton Raphson 算法找到最大似然估计并将解与 glm() 进行比较。所以我尝试在 R 中使用 maxLik() 。结果出错,我以前没有使用过这个包,请修复这个错误,谢谢!!

d <- read.delim("http://dnett.github.io/S510/Disease.txt")
    
d$disease=as.factor(d$disease)
d$ses=as.factor(d$ses)
d$sector=as.factor(d$sector)
str(d)

y<-as.numeric(as.character(d$disease))
x1<-as.numeric(as.character(d$age))
x2<-as.numeric(as.character(d$sector))
x3<-as.numeric(as.character(d$ses))

oreduced <- glm(disease ~ age + sector, family=binomial(link = logit), data=d)
summary(oreduced)
    
    
nlldbin=function(param){
  eta<-param[1]+param[2]*x1+param[3]*(x2==2)
  p<-1/(1+exp(-eta))
  -sum(y*log(p)+(1-y)*log(1-p),na.rm=TRUE)
}
est_nr<-maxLik(nlldbin,start=c(0.01,0.01,0.01))
summary(est_nr)

结果是

Iteration 1 
Parameter:
[1]   9841290 377928533   4325584
Gradient (first 30 components):
     [,1] [,2] [,3]
[1,]  NaN  NaN  NaN
Error in maxNRCompute(fn = function (theta, fnOrig, gradOrig = NULL, hessOrig = NULL,  : 
  NA in gradient

我们正在尝试最大化对数似然,但由于对总和应用负数,您的函数正在最小化。因此只需删除否定,给出:

nlldbin <- function(param){
  eta <- param[1] + param[2]*x1 + param[3]*(x2==2)
  p <- 1/(1+exp(-eta))  
  sum(y*log(p) + (1-y) * log(1-p), na.rm=TRUE)
}

这与其他优化器不同,例如 optim,后者通常默认最小化。这就是为什么你会像你一样否定总和。

ps 你可以使用内置函数编写你的函数,这可能会更稳定一点(并且输入更少):

nlldbin2 <- function(param){
  eta <- cbind(1, x1, x2 == 2) %*% param
  p <- plogis(eta)
  sum(dbinom(y, 1, p, log=TRUE))
}