0
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ 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! $\endgroup$
    – JimB
    Jan 30, 2020 at 16:51
  • $\begingroup$ 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.stackoverflow.com/questions/50078497/… $\endgroup$
    – 51sep
    Jan 31, 2020 at 15:53
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Jan 31, 2020 at 16:08
  • $\begingroup$ 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. $\endgroup$
    – 51sep
    Feb 1, 2020 at 16:33

2 Answers 2

1
$\begingroup$

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

$\endgroup$
0
$\begingroup$

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)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.