Difference between glm and optim for likelihood value on logistic regression in R

Question

I still confuse my previous question on here1 and here2. About logLik of logistic regression in the case of proportion(=yes/yes+no). I try to validate it using optim() by following program. But it was not same. (I could check the same value in the case with “weight=n”). When estimating as the proportion without “weight=n”, I can’t understand how to estimate log-likelihood . Please give me some advice.

logLik() : -1.547104

optim : 2.474444

x<-c(2,3,5,6)
yes<-c(2,1,3,4)
no<-c(3,4,2,1)
n<-yes+no
yp<-yes/n

#-----glm
modelcp<- glm(yp~x,family=binomial)
(result<-summary(modelcp))
#            Estimate Std. Error z value Pr(>|z|)
#(Intercept)  -2.0608     3.0155  -0.683    0.494
#x             0.5152     0.7038   0.732    0.464
#    Null deviance: 0.85152  on 3  degrees of freedom
#Residual deviance: 0.25523  on 2  degrees of freedom
logLik(modelcp)
#'log Lik.' -1.547104 (df=2)

#-----optim
f1<-function(para){
eta<-para[1]+para[2]*x
p<-1/(1+exp(-eta))
-sum(log(choose(1,yp))+yp*log(p)+(1-yp)*log(1-p),na.rm=TRUE)
}
(optim1<-optim(c(1,1),fn=f1,hessian=TRUE))
#$par
#[1] -2.0608361  0.5152331
#$value
#[1] 2.474444

it was same, “with weight = n”

#-----glm
modelcp<- glm(yp~x,family=binomial,weight=n)
(result<-summary(modelcp))
logLik(modelcp)
#'log Lik.' -4.548172 (df=2)

#-----optim
f1<-function(para){
eta<-para[1]+para[2]*x
p<-1/(1+exp(-eta))
-sum(log(choose(n,yes))+yes*log(p)+(n-yes)*log(1-p),na.rm=TRUE)
}
(optim1<-optim(c(1,1),fn=f1,hessian=TRUE))
#$value
#[1] 4.548172

my previous question1 :Difference between binary and count data of same data on logistic regression in R

my previous question2 :Difference between with and without “weight” option of the same data on logistic regression in R

Answer 1

You need to use the observed counts and not the fraction in the likelihood function.

It is actually a mystery to me how you can make choose(1,yp) work when yp is not an integer number.

The following code will give you the same likelihood as your manual optim function

modelcp <- glm(cbind(yes, no) ~ x , family = binomial)
print(logLik(modelcp))

Giving -4.548172

The Q&A here explains more how you can treat the data as Bernoulli or as Binomial distributed and why this gives a different value for the likelihood (it is only a difference by a constant, the functional shape is the same).

Answer 2

I checked them by SAS. If I have to say the answer to my question, only “logLik()” may be wrong but glm() is OK in R. Because SE is same all(glm, optim, SAS). But anyway, I think now "2.474444" would be answer, and I couldn't find where the value “-1.547104” comes from. I don't have confidence about this answer yet, if someone makes it, please give me some advice, thank you.

optim

(SE<-sqrt(abs(diag(solve(optim1$hessian)))))
#[1] 3.0157115 0.7037959

SAS

data dt00;
input x yes no n yp w;
cards;
2 2 3 5 0.4 0.2
3 1 4 5 0.2 0.2
5 3 2 5 0.6 0.2
6 4 1 5 0.8 0.2
;
run;

proc genmod data = dt00 descending;
 model yes/n = x / dist = binomial link=logit;
#weight w;
#weight n;
run;

#without weight
#Deviance            :1.2762
#Log Likelihood      :-12.3722
#Full Log Likelihood :-4.5482
#Intercept           :-2.0608(1.3486)
#x                   :0.5152(0.3147)

#with weight=w
#Deviance            :0.2552
#Log Likelihood      :-2.4744
#Full Log Likelihood :-0.9096
#Intercept           :-2.0608(3.0155)
#x                   : 0.5152(0.7038)

#with weight=n
#Deviance            :6.3808
#Log Likelihood      :-61.8611
#Full Log Likelihood :-22.7409
#Intercept           :-2.0608(0.6031)
#x                   :0.5152(0.1408)