0

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

51sep
  • 227
  • 1
  • 5
  • Why would you not use "weight=n"? The warning message you get when you run your `glm` model should cause some concern: Warning message: In eval(family$initialize) : non-integer #successes in a binomial glm! – JimB Jan 30 '20 at 16:51
  • Thank you for your advice. Yes, you are right. I didn’t notice the Warning message, sorry. And I found it would be the similar question.https://stackoverflow.com/questions/50078497/weighted-logistic-regression-in-r – 51sep Jan 31 '20 at 15:53
  • While the estimates of the parameters are the same (because all of the sample sizes are equal), you do get the wrong standard errors with you leave off `weight=n`. I am a little surprised that the "warning" isn't a show-stopping error message. – JimB Jan 31 '20 at 16:08
  • Hmm, yes, this question was riddle, so, I did cross validation using SAS, following answer. If you are still interested in the matter, see it. – 51sep Feb 01 '20 at 16:33

2 Answers2

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

Sextus Empiricus
  • 43,080
  • 1
  • 72
  • 161
0

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)
51sep
  • 227
  • 1
  • 5