R: error on estimating mle with MaxLik() package

user_15544553 0
Viewed 179

My purpose is to find maximum likelihood estimator using Newton Raphson algorithm and compare the solution with glm(). so I tried to use maxLik() in R. And it turns out error, I have not use this package before, please fix this error, thank you!!

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)

The result is

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
mle glm r
1 Answers

We are trying to maximise the log-likelihood, but due to the negative applied to the sum, your function is minimising. Therefore just remove the negative, giving:

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)
}

This is different from other optimisers, like optim, which often minimise by default. Which is why you would negate the sum, as you did.

ps you could write your function using in built functions, which may be a bit more stable (and less typing):

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