Discrepancies in the standard errors calculated in different ways I noticed discrepancies in the standard errors (SE) calculated in different ways and I would like to know the reason for it. This question is a follow-up of this one.
Here is an example based on a dataset of Poisson-distributed count data (but Normally distributed data also show discrepancies - see my footnote):
set.seed(100)
dd <- data.frame(
                group = c(rep("group1", 100),rep("group2", 100)),
                values = c(rpois(n=100, lambda=2), rpois(n=100, lambda=7))
                )

I computed the following summary statistics:
library(doBy)
logmean <- function(x){
  log(mean(x))
  }
ss <- summaryBy(values ~ group, data=dd, 
            FUN=c(length, logmean, mean, var, sd)
                )
ss$SE.normal <- sqrt(ss$values.var/(ss$values.length-1))
ss$sd.Poiss <- sqrt(ss$values.mean)
ss$se.Poiss <- sqrt(ss$values.mean/ss$values.length)
ss <- cbind(ss[,"group"], round(ss[,-1],3))
names(ss) <- c("group", "N", "log.mean", "mean", "variance", "sd", "SE","sd.Poiss","se.Poiss")

Thus obtaining object ss:
   group   N log.mean mean variance    sd    SE sd.Poiss se.Poiss
1 group1 100    0.723 2.06    1.653 1.286 0.129    1.435    0.144
2 group2 100    1.937 6.94    8.340 2.888 0.290    2.634    0.263

My goal was to compare the SE I estimated above with those estimated by the model:
glm1 <- glm(values ~ group, data=dd, family="poisson"(link = "log"))

I opted for comparing the Confidence Intervals (CI) associated to SE calculated in different ways. I computed CI using: 


*

*ss$se.Poiss,

*SE provided by summary(glm1),

*SE provided by summary(glm1), back-transformed using
Taylor expansion,

*SE estimates provided by predict(...,se.fit=T).


The code for the latter is:
preddata <- data.frame(group=unique(dd$group))
preds <- predict(glm1, newdata=preddata, type="response", se.fit=T)

I computed the four sets of CI is as follows:
# CI based on the predict() output
fit <- preds$fit
lwr <- fit - 1.96*preds$se.fit
upr <- fit + 1.96*preds$se.fit

# CI based on the SE calculated using the formula for Poisson distribution:
fit.ss <- ss$mean
lwr.ss <- fit.ss - 1.96*ss$se.Poiss
upr.ss <- fit.ss + 1.96*ss$se.Poiss

# CI obtained from summary(glm1) by estimating the CI on the link scale and back-transforming it:
fit.glm.link <- c(summary(glm1)$coefficients[1,1],
                  summary(glm1)$coefficients[1,1]+
                  summary(glm1)$coefficients[2,1])
fit.glm <- exp(fit.glm.link)
lwr.glm <- exp(fit.glm.link - 1.96*summary(glm1)$coefficients[,2])
upr.glm <- exp(fit.glm.link + 1.96*summary(glm1)$coefficients[,2])

# CI calculated using SE from summary(glm1), back-transformed using Taylor expansion formula:
fit.taylor <- exp(fit.glm.link)
# SE on the log scale are transformed by multiplying them by the associated exp(log(parameter)):
se.taylor <- exp(fit.glm.link)*summary(glm1)$coefficients[,2]
lwr.taylor <- fit.taylor - 1.96*se.taylor
upr.taylor <- fit.taylor + 1.96*se.taylor

Here's is a plot of the outcome:

[Estimates are off-set along the x axis to avoid overlap].
Confidence intervals based on back-transformations of estimates provided by summary(glm1), either by computing the CI on the link scale and back-transforming them or by using back-transformed SE, are wider than the CI obtained with the other two methods. Why is that?
Furthermore, here is a comparison of the SE estimates:
        se.Taylor se.Poiss se.predict()
Group 1     0.046    0.046        0.046
Group 2     0.176    0.084        0.084

The SE estimated using Taylor's expansion for Group 2 is different from those estimated otherwise. Why is that?

Note: I produced an example using data following a Normal distribution and comparing SE from summary(lm1), SE=sd/sqrt(N), and SE from predict() (for brevity I am not showing the code). They also show discrepancies:
        SE from summary(lm1) SE=sd/sqrt(N) SE from predict()
Group 1                0.275         0.306             0.275
Group 2                0.388         0.239             0.275


 A: I do not have a full answer about all the ways to obtain standard errors, but  the standard error of a group mean is not the same concept as the standard error of the difference between that group and a given reference. Here, the standard error of the coefficient 'group2' in your model is not the standard error of the mean of group2, because the estimates of the intercept and of the difference 'group2-group1' are correlated (negatively). You can see it with vcov(glm1).
The predict functions take the correlation into account, and you should do it too if you want to match their result in your "back-transformed" and "Taylor" CI. 
In your case, the standard errors for the two groups on the latent scale are sqrt(vcov(glm1)[1,1]) and sqrt(vcov(glm1)[2,2]+vcov(glm1)[1,2]). The data-scale confidence interval for group 2 can be calculated as exp(sum(coef(glm1)) + 1.96*sqrt(vcov(glm1)[2,2]+vcov(glm1)[1,2])); exp(sum(coef(glm1)) - 1.96*sqrt(vcov(glm1)[2,2]+vcov(glm1)[1,2])) and should match closely what you obtained from predict.
Cheers
