1

I'm trying to follow this answer here and the linked paper. But I can't seem to get the same results as the confint gives. Sometimes I get something very simmilar, but sometimes my CI is more conservative or less conservative than confint.

I want to know also if the general idea of the implementation is correct - say your model is $a+ b*x$, and you wish to have a CI for $b$ - are you supposed to find, for every value of $b$ the maximum log-likelihood of $a$, and compute these log-likelihoods, and find the $b$'s whose log-likelihood is bigger than $ll(b^{MLE}) - \chi^2_{1,0.975}/2$ ?

Here's the code:

# Generate the data
rpt = 4
x = rep(seq(1, 10, 0.5), times=rpt)
mu = 2.3*x # + rnorm(length(x), sd=0.5)
y = rpois(length(mu), mu)
plot(x, y)

fit = glm(y ~ x, family=poisson(link="identity"))
fit$coefficients

# profile log likelihood function
loglik = function(b) {
  func = function(a) {
    mu.hat = a + b*x
    ll = 0
    for (i in 1:length(x)) {
      ll = ll + dpois(y[i], mu.hat[i], log=T)
    }
    return(ll)
  }
  fit1 = optimize(func, c(0, 5), maximum=T)
  return(fit1$objective)
}

logLik(fit) # built-in function that computes the log-likelihood
loglik(fit$coefficients[2]) # profile ll - should be equal for MLE

beta = seq(1.9, 2.7, 0.0001)
ll = rep(0, length(beta))
for (i in 1:length(beta)){
  ll[i] = loglik(beta[i])
}
plot(beta, ll)

wch = which(ll > logLik(fit)-qchisq(0.975, 1)/2)
(ind.lo = wch[1])
(beta.lo = beta[ind.lo])
(ind.hi = wch[length(wch)])
(beta.hi = beta[ind.hi])
confint(fit) # compare to confint

(beta.MLE = fit$coefficients[2])
pchisq(2*(loglik(beta.MLE)-loglik(beta.lo)), 1)
pchisq(2*(loglik(beta.MLE)-loglik(beta.hi)), 1)

UPDATE:

Using this glm seems to give exact results:

loglik = function(b) {
  fit1 = glm(y ~ 1, offset=b*x, family=poisson(link="identity"))
  return(logLik(fit1))
}
Maverick Meerkat
  • 2,147
  • 14
  • 27
  • why was this closed? I'm wondering if this is even the correct way to achieve the profile likelihood CI? – Maverick Meerkat Aug 12 '20 at 16:59
  • 1
    Hi @DavidRefaeli, it seems like you are trying to come up with another function to calculate the loglikelihood which involves fitting again. So it's a lot of code without a clear indication of what you want to achieve. Hence I believed it was closed – StupidWolf Aug 13 '20 at 10:00
  • 1
    My suggestion is, if you have already obtained a fit through glm, just extract the loglik from it. Don't need to reinvent the wheel – StupidWolf Aug 13 '20 at 10:01
  • @StupidWolf I was trying to calculate the profile likelihood. And understand how it was calculated. Specifically how is `confint` calculates CI for parameters in GLM fit. In the end it seems my logic was correct, but my implementation of the profile log likelihood function was probably faulty - if I use the GLM fit and extract the likelihood from it, the CI seems to be correct. – Maverick Meerkat Aug 13 '20 at 14:43
  • you have a glm object, so it calls stats:::confint.glm and if you look at this code, it calls MASS:::confint.glm ... i might not have time to check the difference between your code and that implemented.. hopefully this gets you started – StupidWolf Aug 14 '20 at 17:12

0 Answers0