Method to calculate the error around multinomial predictions using R

Question

First, sorry for reposting but I added an example code to explain my request. Hope this is clearer.

After fitting a multinomial model to my data with the "multinom" function (package nnet), I want to show the effect of selected variables controlling for other variable values. I know that the "effects" package do mainly what I want, but I want to be able to calculate the prediction error (confidence interval) by myself. Does someone could tell me the methodology and if possible the R code? I think we should use the delta method, but I'm not sure how to apply it in this case.

Here is a small example code (based on data available in the effects package)

library(nnet)
library(effects)
mod <- multinom(vote ~ age + gender, data=BEPS)
summary(mod)

# Call:
# multinom(formula = vote ~ age + gender, data = BEPS)

# Coefficients:
#                  (Intercept)         age gendermale
# Labour             1.2241862 -0.01562320  0.1682676
# Liberal Democrat   0.4979706 -0.01551381  0.1240998

# Std. Errors:
#                  (Intercept)         age gendermale
# Labour             0.2277826 0.003830006  0.1204621
# Liberal Democrat   0.2694373 0.004578836  0.1436882

# Residual Deviance: 3186.266 
# AIC: 3198.266 

plot(allEffects(mod))

The only thing I need is to be able to calculate the values of errors shown in this graph !

Thank you in advance,