The difference of Standard Error between glm(y~x, family=poisson(link=identity)) and optim() in R

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

• 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. – mlofton Dec 14 '19 at 19:23
• 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. – Sycorax says Reinstate Monica Dec 14 '19 at 19:27
• What are the error tolerances your software is using for the solutions? – whuber Dec 14 '19 at 20:56

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.