R: error on estimating mle with MaxLik() package

Question

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

Does `nlldbin` have NA? It seems `x3` is unused? Have you walked through the `maxLik()` documentation example and does that work? https://www.rdocumentation.org/packages/maxLik/versions/1.4-8/topics/maxLik — M.Viking, Apr 29 '21 at 00:57
the example does work, but when I apply the package to the model, it wont work — Jasmine Helen, Apr 29 '21 at 01:21
I can't see why this doesn't work. Your function works in `optim(c(0,0,0), nlldbin)` . I tried changing the loss to `-sum(dbinom(y, 1, p, log=TRUE))` in case there was a stability issue but same error. — user20650, Apr 29 '21 at 01:44
I have tried using optim package and it works. but in this case I want to compare the estimation with maxLik package where I can use newton-raphson algorithm to find the estimation. Do you have any idea to fix this? — Jasmine Helen, Apr 29 '21 at 02:03
yes, my second comment (which may have been unclear ); use `sum(y*log(p)+(1-y)*log(1-p))` or `sum(dbinom(y, 1, p, log=TRUE))` i.e. no negative — user20650, Apr 29 '21 at 02:17

score 2 · Accepted Answer · answered Apr 29 '21 at 10:12

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

R: error on estimating mle with MaxLik() package

1 Answers1