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))
}
我的目的是使用 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))
}