R FAQ:
How can I estimate the standard error of a transformed parameter in R using the delta method?

Often we do not wish to directly report parameters fit by a model but rather some transformation of these parameters.  The transformation can generate the point estimates of our desired values, but the standard errors of these point estimates are not so easily calculated.  They can, however, be well approximated using the delta method. This is based on the first two terms of the Taylor expansion of the transformation function.  Let f be the transformation function and u be the mean of a random variable X.  The first two terms of the Taylor expansion are then an approximation for f(X),

f(X) = f(u) + (X - u)f'(u)

where f' is the first derivative of f. We can then take the variance of this approximation to find the variance of f(X) to find the standard error of a transformed parameter.

For example, we may fit a multinomial logit model and want to report odds ratios.  However, R fits and displays the linear coefficients. We can use the deltamethod in the msm package to find the standard errors of the odds ratios. Below we load the needed libraries and read in the dataset.

library(mlogit)
library(msm)

mydata <- read.csv("http://www.ats.ucla.edu/stat/r/dae/mlogit.csv")
attach(mydata)

To use the mlogit command in the mlogit package, we need to prepare our dataset using mlogit.data. While our starting dataset had one observation for each subject with their choice of brands 1, 2, or 3 indicated in the brand variable, we reshape the data so that for each subject, we have an observation for each of the possible brands (alt) and the brand variable is now true or false, indicating if the given brand was selected.

mydata$brand<-as.factor(mydata$brand)
mldata<-mlogit.data(mydata, varying=NULL, choice="brand", shape="wide")

mydata[1:10,]

   brand female age
1      1      0  24
2      1      0  26
3      1      0  26
4      1      1  27
5      1      1  27
6      3      1  27
7      1      0  27
8      1      0  27
9      1      1  27
10     1      0  27


mldata[1:10,]

    chid alt brand female age
1.1    1   1  TRUE      0  24
1.2    1   2 FALSE      0  24
1.3    1   3 FALSE      0  24
2.1    2   1  TRUE      0  26
2.2    2   2 FALSE      0  26
2.3    2   3 FALSE      0  26
3.1    3   1  TRUE      0  26
3.2    3   2 FALSE      0  26
3.3    3   3 FALSE      0  26
4.1    4   1  TRUE      1  27

We can now run our model predicting brand with female and age as predictors and brand 1 as our reference category.

mlogit.model<- mlogit(brand~1|female+age, data = mldata, reflevel="1")
summary(mlogit.model)

Call:
mlogit(formula = brand ~ 1 | female + age, data = mldata, reflevel = "1")

Frequencies of alternatives:
      1       2       3 
0.28163 0.41769 0.30068 

Newton-Raphson maximisation 
gradient close to zero. May be a solution 
5 iterations, 0h:0m:0s 
g'(-H)^-1g = 7.65E-19 

Coefficients :
              Estimate Std. Error  t-value  Pr(>|t|)    
alt2        -11.774655   1.774612  -6.6351 3.244e-11 ***
alt3        -22.721396   2.058028 -11.0404 < 2.2e-16 ***
alt2:female   0.523814   0.194247   2.6966  0.007004 ** 
alt3:female   0.465941   0.226090   2.0609  0.039315 *  
alt2:age      0.368206   0.055003   6.6943 2.167e-11 ***
alt3:age      0.685908   0.062627  10.9524 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Log-Likelihood: -702.97
McFadden R^2:  0.70981 
Likelihood ratio test : chisq = 185.85 (p.value=< 2.22e-16)

Next, we can pull out of the model the pieces of information needed to calculate the standard errors: the coefficients (the means of the parameters' distributions) and the variance-covariance matrix of the coefficients. We can calculate the odds ratios by exponentiating the coefficients. Thus our transformation function is f(X) = exp(X). 

estmean<-coef(mlogit.model)
var<-vcov(mlogit.model)
exp(estmean)

        alt2         alt3  alt2:female  alt3:female     alt2:age     alt3:age 
7.697192e-06 1.355886e-10 1.688456e+00 1.593514e+00 1.445140e+00 1.985574e+00

To use the deltamethod command, we need to provide the transformation function, the mean(s), and the variance-covariance matrix.  We are not interested in the odds ratios from the model intercepts, so we first use the delta method to find the standard error of the transformed third term in our estmean vector by indicating x3 in the first argument.  The deltamethod is capable of handling vectors of means and variance-covariance matrices.  We show a second implementation where we provide a single mean and variance. 

deltamethod (~ exp(x3), estmean, var) 
0.3279769

deltamethod (~ exp(x1), estmean[3], var[3,3])
0.3279769

We can similarly find the standard errors for our other RRRs.

deltamethod (~ exp(x4), estmean, var) 
0.3602768

deltamethod (~ exp(x5), estmean, var) 
0.07948726

deltamethod (~ exp(x6), estmean, var)
0.1243496

While the transformation used above is quite trivial, deltamethod works well for more elaborate transformations.  Also, a transformation can include multiple values as seen in the example below.

deltamethod (~ exp(x5)/(exp(x4) + x3), estmean, var)
0.1631436