The difference of parameters between “glm” and “optim” in R

user_12389311 0
Viewed 342

I’d like to know the difference of parameters(intercept, slopes) between “glm” and “optim” in R. I think those predictions would be fine, but I can’t understand why those parameters are different. If it’s a misinterpretation, please give me some advice.

glm; -23.36, 46.72, 46.72

optim; -73.99506, 330.09424, 122.50453

#data
x1<-c(0,0,1,1)
x2<-c(0,1,0,1)
y<-c(0,1,1,1)

#glm
model<-glm(y~x1+x2,family=binomial(link=logit))
summary(model)
#             Estimate Std. Error z value Pr(>|z|)
#(Intercept)    -23.36   71664.47       0        1
#x1              46.72  101348.81       0        1
#x2              46.72  101348.81       0        1

round(fitted(model))
#0 1 1 1


#optim
f<-function(par){
eta<-par[1]+par[2]*x1+par[3]*x2
p<-1/(1+exp(-eta))
-sum(log(choose(1,y))+y*log(p)+(1-y)*log(1-p),na.rm=TRUE)
}
(optim<-optim(c(1,1,1),f))
$par
#-73.99506   330.09424  122.50453

round(1/(1+exp(-(optim$par[1]+optim$par[2]*x1+optim$par[3]*x2))))
#0 1 1 1
glm r optimization
1 Answers

I'm going to answer this in a more general context. You have set up a logistic regression with complete separation (you can read about this elsewhere); there is a linear combination of parameters that perfectly separates all-zero from all-one outcomes, which means that the maximum likelihood estimates are actually infinite in this case. This has several consequences:

  • you will have seen a warning glm.fit: fitted probabilities numerically 0 or 1 occurred (this will usually, but not always, happen in this case)
  • the standard deviations, which are based on the local curvature and depend on the assumption that the log-likelihood surface is quadratic, are ridiculous
  • since the log-likelihood surface becomes flatter and flatter as you go toward extreme values of the parameter (slope approaching zero in the infinite limit), different optimizers will get essentially arbitrarily different answers, depending on the point where the surface happens to be flat enough for the algorithm to conclude that it has gotten sufficiently close to an optimum (= zero gradient). Not only different optimization methods (IRLS as in glm vs. the different method options in optim), but even the same optimization method on different operating systems, or even possibly built with different compilers, will give different answers.
  • the brglm2 package has methods for diagnosing complete separation (?brglm2::detect_separation), as well as methods for computing bias-reduced or penalized fits that will avoid the problem (by changing the objective function in a sensible way).
user_1902770
Related