You have discovered an intimate, but generic, property of GLMs fit by maximum likelihood. The result drops out once one considers the simplest case of all: Fitting a single parameter to a single observation!
One sentence answer: If all we care about is fitting separate means to disjoint subsets of our sample, then GLMs will always yield $\hat\mu_j = \bar y_j$ for each subset $j$, so the actual error structure and parametrization of the density both become irrelevant to the (point) estimation!
A bit more: Fitting orthogonal categorical factors by maximum likelihood is equivalent to fitting separate means to disjoint subsets of our sample, so this explains why Poisson and negative binomial GLMs yield the same parameter estimates. Indeed, the same happens whether we use Poisson, negbin, Gaussian, inverse Gaussian or Gamma regression (see below). In the Poisson and negbin case, the default link function is the $\log$ link, but that is a red herring; while this yields the same raw parameter estimates, we'll see below that this property really has nothing to do with the link function at all.
When we are interested in a parametrization with more structure, or that depends on continuous predictors, then the assumed error structure becomes relevant due to the mean-variance relationship of the distribution as it relates to the parameters and the nonlinear function used for modeling the conditional means.
GLMs and exponential dispersion families: Crash course
An exponential dispersion family in natural form is one such that the log density is of the form
$$
\log f(y;\,\theta,\nu) = \frac{\theta y - b(\theta)}{\nu} + a(y,\nu) \>.
$$
Here $\theta$ is the natural parameter and $\nu$ is the dispersion parameter. If $\nu$ were known, this would just be a standard one-parameter exponential family. All the GLMs considered below assume an error model from this family.
Consider a sample of a single observation from this family. If we fit $\theta$ by maximum likelihood, we get that $y = b'(\hat\theta)$, irrespective of the value of $\nu$. This readily extends to the case of an iid sample since the log likelihoods add, yielding $\bar y = b'(\hat\theta)$.
But, we also know, due to the nice regularity of the log density as a function of $\theta$, that
$$
\frac{\partial}{\partial \theta} \mathbb E \log f(Y;\theta,\nu) = \mathbb E \frac{\partial}{\partial \theta} \log f(Y;\theta,\nu) = 0 \>.
$$
So, in fact $b'(\theta) = \mathbb E Y = \mu$.
Since maximum likelihood estimates are invariant under transformations, this means that
$
\bar y = \hat\mu
$
for this family of densities.
Now, in a GLM, we model $\mu_i$ as $\mu_i = g^{-1}(\mathbf x_i^T \beta)$ where $g$ is the link function. But if $\mathbf x_i$ is a vector of all zeros except for a single 1 in position $j$, then $\mu_i = g(\beta_j)$. The likelihood of the GLM then factorizes according to the $\beta_j$'s and we proceed as above. This is precisely the case of orthogonal factors.
What's so different about continuous predictors?
When the predictors are continuous or they are categorical, but cannot be reduced to an orthogonal form, then the likelihood no longer factors into individual terms with a separate mean depending on a separate parameter. At this point, the error structure and link function do come into play.
If one cranks through the (tedious) algebra, the likelihood equations become
$$
\sum_{i=1}^n \frac{(y_i - \mu_i)x_{ij}}{\sigma_i^2}\frac{\partial \mu_i}{\partial \lambda_i} = 0\>,
$$
for all $j = 1,\ldots,p$ where $\lambda_i = \mathbf x_i^T \beta$. Here, the $\beta$ and $\nu$ parameters enter implicitly through the link relationship $\mu_i = g(\lambda_i) = g(\mathbf x_i^T \beta)$ and variance $\sigma_i^2$.
In this way, the link function and assumed error model become relevant to the estimation.
Example: The error model (almost) doesn't matter
In the example below, we generate negative binomial random data depending on three categorical factors. Each observation comes from a single category and the same dispersion parameter ($k = 6$) is used.
We then fit to these data using five different GLMs, each with a $\log$ link: (a) negative binomial, (b) Poisson, (c) Gaussian, (d) Inverse Gaussian and (e) Gamma GLMs. All of these are examples of exponential dispersion families.
From the table, we can see that the parameter estimates are identical, even though some of these GLMs are for discrete data and others are for continuous,and some are for nonnegative data while others are not.
negbin poisson gaussian invgauss gamma
XX1 4.234107 4.234107 4.234107 4.234107 4.234107
XX2 4.790820 4.790820 4.790820 4.790820 4.790820
XX3 4.841033 4.841033 4.841033 4.841033 4.841033
The caveat in the heading comes from the fact that the fitting procedure will fail if the observations don't fall within the domain of the particular density. For example, if we had $0$ counts randomly generated in the data above, then the Gamma GLM would fail to converge since Gamma GLMs require strictly positive data.
Example: The link function (almost) doesn't matter
Using the same data, we repeat the procedure fitting the data with a Poisson GLM with three different link functions: (a) $\log$ link, (b) identity link and (c) square-root link. The table below shows the coefficient estimates after converting back to the log parameterization. (So, the second column showns $\log(\hat \beta)$ and the third shows $\log(\hat \beta^2)$ using the raw $\hat\beta$ from each of the fits). Again, the estimates are identical.
> coefs.po
log id sqrt
XX1 4.234107 4.234107 4.234107
XX2 4.790820 4.790820 4.790820
XX3 4.841033 4.841033 4.841033
The caveat in the heading simply refers to the fact that the raw estimates will vary with the link function, but the implied mean-parameter estimates will not.
R code
# Warning! This code is a bit simplified for compactness.
library(MASS)
n <- 5
m <- 3
set.seed(17)
b <- exp(5+rnorm(m))
k <- 6
# Random negbin data; orthogonal factors
y <- rnbinom(m*n, size=k, mu=rep(b,each=n))
X <- factor(paste("X",rep(1:m,each=n),sep=""))
# Fit a bunch of GLMs with a log link
con <- glm.control(maxit=100)
mnb <- glm(y~X+0, family=negative.binomial(theta=2))
mpo <- glm(y~X+0, family="poisson")
mga <- glm(y~X+0, family=gaussian(link=log), start=rep(1,m), control=con)
miv <- glm(y~X+0, family=inverse.gaussian(link=log), start=rep(2,m), control=con)
mgm <- glm(y~X+0, family=Gamma(link=log), start=rep(1,m), control=con)
coefs <- cbind(negbin=mnb$coef, poisson=mpo$coef, gaussian=mga$coef
invgauss=miv$coef, gamma=mgm$coef)
# Fit a bunch of Poisson GLMs with different links.
mpo.log <- glm(y~X+0, family=poisson(link="log"))
mpo.id <- glm(y~X+0, family=poisson(link="identity"))
mpo.sqrt <- glm(y~X+0, family=poisson(link="sqrt"))
coefs.po <- cbind(log=mpo$coef, id=log(mpo.id$coef), sqrt=log(mpo.sqrt$coef^2))
Best Answer
It depends on what you mean by "the R output." It's not conceptually different for generalized linear models or ordinary least squares models. Maybe it's simplest to start thinking about this in the context of ordinary least squares regression. Consider the example from the
lm()help page:If you then type
lm.D9at the prompt, you get the coefficient estimatesIf that's all you have to work with you can't get the standard errors of the coefficients even in this simple case. The standard errors of the coefficients are the square roots of the diagonal elements of the coefficient variance-covariance matrix. As this page (among many others) shows, that matrix is $\hat \sigma^2 (X'X)^{-1}$, where $\hat \sigma^2 $ is the estimated error variance and $X$ is the design matrix that submits the data to the model. Even for an intercept-only model the standard error of the mean isn't $\sqrt {\hat \sigma^2}$ by itself, as one might infer from the question as stated. It's $\sqrt {\hat \sigma^2/n}$, where $n$ is the number of observations.
R models typically allow for a
vcov()function to report that matrix:For many model types (like this one) you can get more information simply by asking for a
summary()that directly reports the standard errors along with the coefficient estimates (not shown).The principle is the same for generalized linear models like Poisson that are fit by maximum likelihood estimation (MLE). As this answer says:
That is, the variance-covariance matrix of the coefficient estimates is the inverse of the matrix of second derivatives of the log-likelihood with respect to the covariate values at the final solution. For ordinary linear regression under the classic assumption of independent, normally distributed errors having constant variance, that ends up identical to the formula above. For other generalized linear models, however, there is not typically a closed-form solution and the entire variance-covariance matrix must be estimated from the combination of the model and the data.
A Poisson model assumes that the variance around any estimate equals the mean. So the formula for the likelihood used to fit the model, where it includes a parameter $\lambda$ for the mean, inherently also includes the variance $\lambda$. For the Poisson example of the
glm()help page:Again in this case,
summary(glm.D93)will produce those values along with additional information about the model.In this Poisson case, it's mainly how the maximum-likelihood estimate $\hat \lambda$ of the mean also estimates the variance. For an intercept-only model with $n$ observations, this page shows that the standard error of the estimate of the mean is $\sqrt{\hat \lambda/n}$. This answer shows a similar intercept-only result for a negative binomial (in a different parameterization from yours). In both cases you have the inverse relationship of the coefficient standard error to $\sqrt{n}$.
For a single binary predictor and independent observations, as in your model, you essentially have 2 such intercept-only models in the Poisson case, or for the negative binomial when you pre-specify a known $\theta$. For more complicated regressions, in which you need to estimate the negative-binomial $\theta$ or have more predictors whose coefficients need to be estimated together, the standard errors also have to do with how precisely the model fits the data under the modeling assumptions.