(Reposting after changing some things and the title: standard error for glm: R and theory arent agreeing)
I was trying to work out by hand the standard errors for an exponential glm (gamma distn with dispersion=1) with a log link ($\log\mu=\beta_0 + \beta_1 x$).
I'm getting some contradictions. by hand I do this:
exponential has prob function $f = (1/\mu) \exp(-y/\mu)$, so log-likelihood is $\ell = -\sum \log\mu + (y/\mu)$. That means $d\ell/d\mu = \sum (y-\mu)/\mu^2$.
- $\partial\ell/\partial\beta_0 = \sum (y-\mu)/\mu$ and
- $\partial\ell/\partial\beta_0 = \sum x(y-\mu)/\mu$.
so using the chain rule the information matrix bits are:
- $\partial^2\ell/\partial\beta_0^2 = \sum y/\mu$
- $\partial^2\ell/\partial\beta_1^2 = \sum x^2y/\mu$
- $\partial^2\ell/\partial\beta_0\beta_1 = \sum xy/\mu$
I tried to see if that worked in a numerical example in R.
data(trees) # comes with R y <- trees$Volume x <- trees$Girth ### in r: mm <- glm(y ~ x, family=Gamma(link=log)) InfMat.R <- solve(summary(mm, dispersion=1)$cov.unscaled) ### by hand lb0b0 <- sum( y / fitted(mm)) lb0b1 <- sum( y * x / fitted(mm)) lb1b1 <- sum( y * x^2 / fitted(mm)) InfMat.me <- array( dim=c(2,2)) InfMat.me[1,1] <- lb0b0 InfMat.me[1, 2] <- InfMat.me[2, 1] <- lb0b1 InfMat.me[2, 2] <- lb1b1 InfMat.R # what r gives me InfMat.me # when I do it ny hand
the two information matrices are so close to being the same:
> InfMat.R (Intercept) x (Intercept) 31.0 410.70 x 410.7 5736.55 > InfMat.me [,1] [,2] [1,] 31.0 410.700 [2,] 410.7 5718.696
So it is only element (2,2) that differs. and the rest are right. It seems odd that I can get all the bits right except one.
Can anyone explain what I'm doing wrong?