Search code examples
rglmestimationmle

R: Different result from glm and mle2 package in R


So I want to find the estimate parameter using GLM and compare it with mle2 package. Here's my code for GLM

d <- read.delim("http://dnett.github.io/S510/Disease.txt")

d$disease=factor(d$disease)
d$ses=factor(d$ses)
d$sector=factor(d$sector)
str(d)
glm2 <- glm(disease~ses+sector, family=binomial(link=logit), data=d)
summary(glm2)

And my code for mle2()

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

library(bbmle)
nlldbin=function(A,B,C,D){
  eta<-A+B*(x3==2)+C*(x3==3)+D*(x2==2)
  p<-1/(1+exp(-eta))
  joint.pdf= (p^y)*((1-p)^(1-y))
  -sum(joint.pdf, log=TRUE ,na.rm=TRUE)
}
st <- list(A=0.0001,B=0.0001,C=0.0001,D=0.0001)
est_mle2<-mle2(start=st,nlldbin,hessian=TRUE)
summary(est_mle2)

But the result is quiet different. Please help me to fix this, thank you!

> summary(est_mle2)
Maximum likelihood estimation

Call:
mle2(minuslogl = nlldbin, start = st, hessian.opts = TRUE)

Coefficients:
     Estimate  Std. Error z value  Pr(z)
A    -20.4999   5775.1484 -0.0035 0.9972
B     -5.2499 120578.9515  0.0000 1.0000
C     -7.9999 722637.2670  0.0000 1.0000
D     -2.2499  39746.6639 -0.0001 1.0000

> summary(glm2)

Call:
glm(formula = disease ~ ses + sector, family = binomial(link = logit), 
    data = d) 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.52001    0.33514  -4.535 5.75e-06 ***
ses2        -0.08525    0.41744  -0.204 0.838177    
ses3         0.16086    0.39261   0.410 0.682019    
sector2      1.28098    0.34140   3.752 0.000175 ***

Solution

  • I'm not sure your definition of eta is correct. I would use the model matrix.

    X <- model.matrix(~ ses + sector, data = d)
    
    nlldbin <- function(A,B,C,D){
      eta <- X %*% c(A, B, C, D)
      p <- 1/(1+exp(-eta))
      logpdf <- y*log(p) + (1-y)*log(1-p)
      -sum(logpdf)
    }