8
$\begingroup$

I'm executing the following program glm(y~x, family=poisson(link=identity)). I can't understand the difference of SE from glm(y~x, family=poisson(link=identity)) and optim(). Please give me some advice.

optim()

x<-c(1,2,3,4)
y<-c(2,3,5,4)
f<-function(arg){
  a=arg[1]
  b=arg[2]
  mu=a+b*x
  -sum(y*log(mu)-mu-log(factorial(y)))
}
(reso<-optim(c(1,1),f,hessian=TRUE))
#$par
#[1] 1.2785717 0.8886162
#$hessian
#[1,] 1.239716 2.717462
#[2,] 2.717462 7.343171
sqrt(diag(solve(reso$hessian)))#SE from optim()
#[1] 2.0669196 0.8492641

glm()

resg<-glm(y~x,family=poisson(link=identity))
summary(resg)
#            Estimate Std. Error z value Pr(>|z|)
#(Intercept)   1.2784     1.9766   0.647    0.518
#x             0.8887     0.8141   1.092    0.275
diag(sqrt(vcov(resg)))#SE from glm(poisson, identity)
#(Intercept)           x 
#   1.976575    0.814139
$\endgroup$
  • $\begingroup$ Hi: it depends on the algorithm that glm uses because that determines how the covariance is calculated. In the first case, you're maximizing a likelihood so what the gradient is ( close to zero but maybe not zero ) when the algorithm stops is crucial, since the hessian is the derivative of the gradient. So, even though you get ~ the same coefficient estimates, ( again, I'm not sure what glm does. it may run IRWLS or maximize a likelhood ? ), this doesn't mean that the covariance estimates are calculated in the same way. I would look at John Fox's CAR which covers glm R code nicely. $\endgroup$ – mlofton Dec 14 '19 at 19:23
  • $\begingroup$ The first clue that there’s a difference isn’t in the standard error estimates but in the estimated coefficients. This tells us that there are different solutions found by each method. I also suspect that the numerical approximation of the Hessian differs from a direct computation of it. $\endgroup$ – Sycorax says Reinstate Monica Dec 14 '19 at 19:27
  • 1
    $\begingroup$ What are the error tolerances your software is using for the solutions? $\endgroup$ – whuber Dec 14 '19 at 20:56
6
$\begingroup$

In statistical likelihood theory, minus the second derivative of the log-likelihood function is called the observed information. We might write this as $$ I = -\ddot \ell(y; \theta) $$ where the dots indicate differentiation with respect to $\theta$. The expected value of the observed information $$ {\cal I} = E(I) $$ is called Fisher information or expected information.

Observed and expected information are asymptotically equivalent (by the law of large numbers) under the same regularity conditions that guarantee that the maximum likelihood estimators are consistent. This implies that the observed and expected information will usually be close, in relative as well as absolute terms, when the standard errors are small.

In glm theory, Fisher information is preferred over observed information because it (1) has a much simpler analytic form, (2) is guaranteed to be positive definite (which observed information is not) and (3) is the same as the Cramer-Rao Lower Bound for the variance of unbiased estimators.

If the glm model has a canonical link, then the distinction is less important because observed and expected information are in that case identical when computed at the maximum likelihood estimator of $\theta$.

If you use optim to maximumize the log-likelihood of a glm model, then the maximum likelihood estimates returned by glm and optim will be the same apart from rounding errors if both algorithms are run to convergence. The standard errors from glm however will generally differ from those from optim because

  1. glm returns Fisher information whereas optim computes observed information and

  2. glm uses an exact analytic formula for Fisher information whereas optim approximates the Hessian numerically from second differences of the log-likelihood.

Had you had used a log-link instead of identity link for your example, then item 1 would no longer cause any difference because the log-link is canonical and observed and expected information would then become identical at convergence.

Note on terminology

Some authors use "observed Fisher information" as a synonym for "observed information" and "expected Fisher information" as a synonym for "Fisher information". I think this terminology probably originates from Efron & Hinkley (1978):

B. Efron and D.V. Hinkley (1978). Assessing the accuracy of the maximum likelihood estimator: Observed versus expected Fisher information. Biometrika 65(3), 457–483.

$\endgroup$
  • $\begingroup$ Very nice answer. Thanks. $\endgroup$ – mlofton Dec 15 '19 at 1:48
  • 1
    $\begingroup$ @mlofton Thanks, but an upvote would be more in the spirit of the site, see stats.stackexchange.com/help/someone-answers $\endgroup$ – Gordon Smyth Dec 16 '19 at 11:58
  • $\begingroup$ Thank you very much for the detail to my question. I'll check it and have to study more.. $\endgroup$ – 51sep Dec 17 '19 at 17:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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