How to calculate standard errors for GLMs fitted values "by-hand", without using predicted() in R?

Question

I am running a negative binomial model without creating a glm object in R. I found an answer at StackExchange on how to get standard errors "by-hand" here, but it shows only how to get values on a linear scale, not in the scale of the response, just as by default for type argument in predict(). To go from a linear scale to the scale of the fitted values for a negative binomial is just exp(fit), but how to get standard errors?

Answer 1

To get the standard errors, you could either approximate them with the delta method, or just use simulation! You can probably use the simulate() generic in R.

For help on how to use the delta method, have a look at the answer to Calculating the Variance using Delta Method

Following, an example on how to get the standard errors from the model object by simulation. I give an Poisson example, where we know the standard error from the mean, so as to be able to check on the simulations. Code:

set.seed(7*11*13)  
n  <-  100
N  <-  1000 # number of replications
y  <-  rpois(n, lambda=exp(5+x))
mod  <-  glm(y~x, family=poisson())
sims  <-  simulate(mod, nsim=N)  # this is sims from the estimated predictive distribution
preds  <-  matrix(0,n,N)
for (i in seq(along=sims)) {
    preds[,i]  <-  sims[[i]]
}
means  <-  apply(preds,1,mean)
vars   <-  apply(preds,1,var)
plot(means,vars)

the plot showing the expected linear relation.