The delta method for dummies

Reporting variances (or standard deviations) for all estimates is normally regarded as a basic requirement for scientific publications. In my work as consultant I am very frequently asked such questions: “I have the measure Q with variance q and the measure W with variance w. What is the variance of exp(Q + W)?”

I recently wrote a post relating to linear transformations (see here), but I feel that being able to menage nonlinear transformations is much more useful. Frequentist statistics offer an approximate solution to this problem, that is normally known as the 'delta' method.

The basic idea is to create a linear approximation to our nonlinear function and use the former to approximate the variance, by using the same methods as given in my above mentioned post. Indeed, every nonlinear function G (close to linear, differentiable) of X can be locally approximated by the tangent line on a point \(\mu\) (in the close vicinity of X), with slope equal to the first derivative G'(mu) (Taylor series approximation to first order):

\[G(X) \simeq G(\mu ) + G'(\mu )(X - \mu )\]

Having transformed a nonlinear function into a linear function, we can approximate the variance by using:

\[Var[G(X)] \simeq [G'(\mu)]^2 Var(X) \]

In words: if I have a nonlinear (close to linear) function G(X), where X is a random variable with small variance and normal errors, the variance of G(X) is approximated by knowing the first derivative G' and the variance of X

Example 1: exp(X)

To calculate the derivatives, we can use an R function:

D(expression(G(x), “x”)

In this case:

D(expression(exp(b0)), "b0")

## exp(b0)

If we have:

X <- 3.5
VarX <- 0.15
exp(X)

## [1] 33.12

VarX * exp(X)^2

## [1] 164.5

In R there is a shortcut function to calculate the delta method, from the car package. The syntax is:

library(car)
deltaMethod(object = c(X = X), g = "exp(X)", vcov = matrix(VarX, 1, 1))

##        Estimate    SE
## exp(X)    33.12 12.83

By using this function we obtain the standard error, i.e. the standard deviation of the sampling distribution. By squaring the above SE (12.82556) we see that the results are obviously the same.

Example 2: 1/X

The first derivative is:

D(expression(1/x), "x")

## -(1/x^2)

Delta standard errors may be obtained as:

X <- 3.5
VarX <- 0.15
1/X

## [1] 0.2857

VarX * ((-1/X)^2)^2

## [1] 0.0009996

With R:

library(car)
deltaMethod(object = c(X = X), g = "1/X", vcov = matrix(VarX, 1, 1))

##     Estimate      SE
## 1/X   0.2857 0.03162

Example 3: log(X)

The first derivative is:

D(expression(log(x)), "x")

## 1/x

Delta standard errors may be obtained as:

X <- 3.5
VarX <- 0.15
log(X)

## [1] 1.253

VarX * (1/X)^2

## [1] 0.01224

Example 4: Inverse logit function

First derivative is:

D(expression(exp(X)/(1 + exp(X))), "X")

## exp(X)/(1 + exp(X)) - exp(X) * exp(X)/(1 + exp(X))^2

Delta standard errors may be obtained as:

X <- 3.5
VarX <- 1.5
exp(X)/(1 + exp(X))

## [1] 0.9707

VarX * (exp(X)/(1 + exp(X)) - exp(X) * exp(X)/(1 + exp(X))^2)^2

## [1] 0.001214

---

The extension to the case of two random variables is straightforward: we need the partial derivatives in each variable, plus all covariances

\[ var[g(x, y,..)] = D S D^T \]

where D is the matrix of partial derivatives with respect to all the variables, S is the variance-covariance matrix of all variables.

Example 5: X Y

First derivatives are:

D(expression(X * Y), "X")

## Y

D(expression(X * Y), "Y")

## X

Delta standard errors may be obtained as:

X <- 3.5
Y <- 4.2
VarX <- 1.5
VarY <- 1.2
CoVar <- 0.8
X * Y

## [1] 14.7

(Y^2) * VarX + (X^2) * VarY + 2 * X * Y * CoVar

## [1] 64.68

With R:

library(car)
sigma <- matrix(c(VarX, CoVar, CoVar, VarY), 2, 2, byrow = T)
deltaMethod(object = c(X = X, Y = Y), g = "X*Y", vcov = sigma)

##       Estimate    SE
## X * Y     14.7 8.042

Example 6: X/Y

First partial derivatives are:

D(expression(X/Y), "X")

## 1/Y

D(expression(X/Y), "Y")

## -(X/Y^2)

Delta standard errors may be calculated as:

X <- 3.5
Y <- 4.5
VarX <- 0.9
VarY <- 0.8
CoVar <- 0.5
X/Y

## [1] 0.7778

(1/Y)^2 * VarX + (-(X/(Y^2)))^2 * VarY + 2 * (1/Y) * (-X/(Y^2)) * CoVar

## [1] 0.02993

With R:

library(car)
sigma <- matrix(c(VarX, CoVar, CoVar, VarY), 2, 2, byrow = T)
deltaMethod(object = c(X = X, Y = Y), g = "X/Y", vcov = sigma)

##     Estimate    SE
## X/Y   0.7778 0.173

Example 7: exp(X + Y)

Partial derivatives are:

D(expression(exp(X + Y)), "X")

## exp(X + Y)

D(expression(exp(X + Y)), "Y")

## exp(X + Y)

Delta standard errors may be obtained as:

X <- 3.5
Y <- 0.5
VarX <- 0.15
VarY <- 0.08
CoVar <- 0.05
exp(X + Y)

## [1] 54.6

(exp(X + Y)^2) * VarX + (exp(X + Y)^2) * VarY + 2 * exp(X + Y) * exp(X + Y) * 
    CoVar

## [1] 983.7

With R:

library(car)
sigma <- matrix(c(VarX, CoVar, CoVar, VarY), 2, 2, byrow = T)
deltaMethod(object = c(X = X, Y = Y), g = "exp(X+Y)", vcov = sigma)

##            Estimate    SE
## exp(X + Y)     54.6 31.36

Two final notes relating to nonlinear transformations are very important:

1. The estimated variance is always approximate.
2. If the original variables are gaussian, the transformed variable is normally not gaussian.