Delta Method in R

Suppose you want the standard error of a nonlinear function of the parameters of an estimated model. Let’s work with the model \[y_{t} = \alpha + \beta x_{t} + \gamma z_{t} + \varepsilon_{t}\] We want the standard error of the product \(\widehat{\beta} \widehat{\gamma}\). R package msm provides a function to do that calculation using the delta method.

Generate some data:

set.seed(100)
y <- rnorm(100)
x <- rnorm(100)
z <- rnorm(100)

Estimate the model:

fit <- lm(y ~ x+z)

Pull out the relevant coefficients. In this case, it is the coefficients on x and z. There are two equivalent ways to do this.

b <- coefficients(fit)
theta <- b[c("x", "z")]
theta <- b[2:3]

The first approach, where we have referred to coefficients by name, is recommended in most cases. The second approach, where we refer to coefficients by numerical index, is only useful if you need to pull a large block of coefficients. It is too easy to make mistakes and too hard to spot mistakes when you refer to coefficients by numerical index. Even worse, when you make a change to your model, you have to remember to come back and change this line of code to pull the correct coefficients. You might forget to do that after adding more control variables to your model to satisfy a referee 18 months down the road.

Pull out the relevant block from the covariance matrix. Once again, we prefer to do this by name rather than by numerical index, but it can be done either way:

v <- vcov(fit)[c("x", "z"), c("x", "z")]
v <- vcov(fit)[2:3, 2:3]

Load the msm package and calculate the standard error:

std.err <- deltamethod(~(x1*x2), theta, v)

The first argument to deltamethod is a formula to calculate the nonlinear transformation of the parameters. What are x1 and x2? x1 refers to the first element of theta, x2 to the second element, and so on. You can think of the first argument as shorthand for a function, with arguments x1, x2, and so on.

If you want to calculate the standard errors of multiple nonlinear functions of the parameters (\(\beta \gamma\), \(\beta + \gamma\), \(\beta - \gamma\)), you can send a list of transformations as the first argument:

std.err <- deltamethod(list(~(x1*x2), ~(x1+x2), ~(x1-x2)), theta, v)

Sums and cumulative sums of coefficients

If you want to the standard error of a sum of coefficients, there is a function in the tstools package to do that calculation. As above, you can index coefficients by name (preferred) or by numerical index:

library(tstools)
std.err <- se.sum(fit, c("x", "z"))
std.err <- se.sum(fit, 2:3)

There is also a function for the standard error of the cumulative sum of the parameters:

std.errors <- se.cumsum(fit, c("x", "z"))