Mean Differences in Poisson Regression I performed a Poisson regression.  My dependent variable is number of days and my independent variables are all binary indicators (e.g. facility = facility1 facility2, facility3 and severity = extreme, major, moderate, minor and treatment level=high, medium, low).  So my model looks like this:
$\log(\mu)=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_3+\beta_4X_4+\beta_5X_5+\beta_6X_6+\beta_7X_7+\beta_8X_8+\beta_9X_9+\beta_{10}X_{10}$
where $X_1-X_3$ are the binary facility indicators, $X_4-X_7$ are the binary severity indicators and $X_8-X_{10}$ are the binary treatment indicators.  
I am interested in obtaining the mean predicted differences for a given facility and severity between each of the treatment levels.  It is straight forward for me to obtain a point estimate.  For example, to compare high to low treatments in facility 1 and extreme severity, I simply calculate:
$\hat{Z}=\hat{\mu}_{1,4,8}-\hat{\mu}_{1,4,10}=\exp(\hat{\beta}_0+\hat{\beta}_1+\hat{\beta}_4+\hat{\beta}_8)-\exp(\hat{\beta}_0+\hat{\beta}_1+\hat{\beta}_4+\hat{\beta}_{10})$
Now, I need to calculate the variance of this estimate (in order to obtain confidence intervals).  I understand that I can get this using the Taylor Series/Delta Method approach by linearizing $Z.$  I have been able to obtain these variance estimates using SAS' nlmixed for a single difference (high-low) (see here for a SAS note that explains something very similar), but I was wondering if these is any way to conveniently obtain these variance estimates using without having to specify every single estimate statement.  I'm looking for a way to compare all pairwise treatments (similar to using an lsmeans statement in SAS).  I'm well versed in SAS, R, and dabble in STATA, and can use any of these if there are suggestions.  I just want to have to avoid specifying every single parameter to obtain all combinations of my estimates (which in my actual model contains over 60 parameters -- this was just a simplified example).
 A: The key is to apply a couple statistical principles. You can estimate the mean difference pretty easily, just use the "predict" command. But getting an approximate confidence interval is difficult. You could bootstrap the data to get an efficient robust 95% exact interval (very attractive indeed), or you can use a normal approximation to the difference. 
The LSMeans command does the former, however if you were fitting a mixed effects model (as would be done with nlmixed), there are many more considerations (too many to discuss). It's beyond the scope of the answer. Therefore, I'll assume you fit a plain old Poisson regression model using myglm.fit <- glm(..., family=poisson) as you've said.
The normal approximation is easily computed using the delta method. The covariance matrix for the parameter estimates $\Sigma$ is obtained by computing Sigma <- vcov(myglm.fit).
Note $\hat{\mu}_1 = \exp \left( X_1 \hat{\beta} \right)$ and similarly for $\hat{\mu}_2$. Then $f(u, v) = \exp(u) - \exp(v)$ and $\nabla f(u, v) = [f(u, v), -f(u, v)]$. Using your notation $Z = f(X_1 \hat{\beta}, X_2 \hat{\beta})$. We can computer the variance of those components, it's just $X_1 \Sigma X_1^T$ and $X_2 \Sigma X_2^T$, and the covariance is $X_1 \Sigma X_2^T$. The $2\times 2$ covariance matrix is $\Xi = [X_1, X_2] \Sigma [X_1, X_2]^T$.
The Delta method, then, is simply applied: 
$$\mbox{Var} (Z) = \nabla f \Xi \nabla f^T$$
and that is an approximate, asymptotically correct 95% confidence interval.
Example:
set.seed(1)

## from ?glm
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
d.AD <- data.frame(treatment, outcome, counts)
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())

## create the proper levels x_1, x_2 to create linear combos
new.data <- data.frame(outcome = factor(c(1, 3), levels=1:3), 
                       treatment=factor(c(1,3), levels=1:3))


## how to create a lsmeans statement using delta method
lincom <- model.matrix(formula(glm.D93)[-2], data=new.data)
mu12 <- exp(lincom %*% beta)
nabla <- matrix(c(mu12[1], -mu12[2]), c(1,2))
Z <- diff(mu12)
varZ <- nabla %*% lincom %*% Sigma %*% t(lincom) %*% t(nabla)

## CI:
Z + c(-1.96 * sqrt(varZ), 1.96 * sqrt(varZ))

## ex bootstrap
bs <- replicate(1000, {
  d.AD[, 'counts'] <- rmultinom(n = 1, size=sum(counts), prob = counts/sum(counts))
  glm.D93 <- glm(counts ~ outcome + treatment, data=d.AD, family = poisson())
  diff(predict(glm.D93, newdata=new.data, type='response'))
})

quantile(bs, c(0.25, 0.975))

Giving disparate results for CIs:
> Z + c(-1.96 * sqrt(varZ), 1.96 * sqrt(varZ))
[1] -15.281464   4.614797

> quantile(bs, c(0.25, 0.975))
    25%   97.5% 
-8.6350  3.8875 

