Delta Method Standard Errors for Odds Ratios

Alternatively, we can use the SE for the odds ratio to determine a normal (and symmetric) approximation for the 95% CI. But what is the SE for the odds ratio?

An approach known as the delta method is used frequently to come up with standard errors for nonlinear transformations of model parameters.

It is based on computing the variance for a Taylor series linearization of the function.

A Taylor Series rewrites a function at a given location \(a\) as a (possibly infinite) sum of the function’s derivatives.

\[ f(x) = f(x) + f'(a)(x - a) + \frac{f''(a)}{2!}(x - a)^2 + \frac{f'''(a)}{3!}(x - a)^3 + \ldots \]

A Taylor series approximation chops off all but the first one or two derivatives. A linear approximation to \(f(x)\) at \(a\) is thus

\[ f(x) = f(x) + f'(a)(x - a) \]

Take a linear transformation of a random variable \(x\).

\[ f(x) = a + bx \]

The variance of \(f(x)\) is known to be

\[ \mbox{Var}(a + bx) = b^2\mbox{Var}(x) \] or

\[ \mbox{Var}(a + bx) = b\mbox{Var}(x)b \] or

\[ \mbox{Var}(a + bx) = f'(x)\mbox{Var}(x)f'(x) \]

Generalizing to any (univariate) differentiable function of a random varaible \(x\) with mean \(\mu\), we can approximate a function of \(x\) at \(\mu\) with

\[ f(x) \approx f(\mu) + f'(x)(x - \mu) \]

meaning that the variance of the function is

\[ \mbox{Var}\left(f(x)\right) = f'(\mu)\mbox{Var}(x)f'(\mu) \]

This generalizes to functions of multiple variables. Simply replace the derivate with the gradient vector and the variance with the variance-covariance matrix.

\[ \mbox{Var}\left(f(\mathbf x)\right) = \nabla(\boldsymbol\mu )^{T}\mbox{Cov}(\mathbf x)\nabla(\boldsymbol\mu) \]

Going back to logistic regression, our random variables are the sample estimates \(\widehat{\beta}_k\), and our function is \(f(\beta_k) = e^{\beta_k}\). The maximum likelihood estimates are the values for the vector \(\boldsymbol \mu\). The covariance matrix is the covariance matrix of the estimates. Both are easily recovered in R from a glm object.

coef(trump_model)
vcov(trump_model)

The function in which we are interested is

\[ f(\beta_k) = e^{\beta_k} \] The delta-method variance is

\[ f(\beta_k) = f'(\beta_k)\mbox{Var}(\widehat{\beta_k})f'(\beta_k) \]

This turns out to be a pretty simple problem, given that

\[ \frac{d e^x}{dx} = e^x \]

What would the variance be for the odds ratio on the noise term?

• The estimate \(\beta_{\mbox{noise}}\) was 0.0804698.
• The variance was 0.6725694

\[ \mbox{Var}\left(\mbox{exp}\left(\beta_{\mbox{noise}}\right)\right) = e^{\beta}\mbox{Var}\left(\hat{\beta}\right)e^{\beta} = 0.79 \]

A normal-approximated 95% confidence interval is found as

\[ 95\% \mbox{ CI} = \beta_k \pm 1.96 \times \mbox{SE}_{\beta_k} \]

where \(\mbox{SE}\) is the square root of the variance.

A function to return odds ratios and confidence intervals based on normal approximations.

OR_SE <- function(glm_object, variable){
  
  d <- exp(coef(glm_object)[variable])
  v <- vcov(glm_object)[variable, variable]
  
  output_df <- as_tibble(list(Estimate = variable,
                              Beta     = coef(glm_object)[variable],
                              OR       = d,
                              `OR SE`  = sqrt(d*v*d))) %>%
    mutate(`Lower CI` = OR - 1.96*`OR SE`,
           `Upper CI` = OR + 1.96*`OR SE`)

  output_df
  
}

Map over all estimates and reduce to tibble.

or_tibble <- map(names(coef(trump_model)), function(i) OR_SE(trump_model, i)) %>%
  reduce(bind_rows)
kable(or_tibble %>% var_mapping(Estimate), align = c("l", rep("c", 5)))

cis <- exp(confint(trump_model)) %>%
  as.data.frame %>%
  rownames_to_column("Variable") %>%
  var_mapping(Variable)
## Waiting for profiling to be done...
graph_data_2 <- tidy(trump_model) %>%
  dplyr::rename(Variable = term, Beta = estimate, SE = std.error, z = statistic, p = p.value) %>%
  mutate_if(is.numeric, funs(round(.,3))) %>%
  var_mapping(Variable) %>%
  mutate(OR = exp(Beta)) %>%
  left_join(cis) %>%
  mutate(Source = "exp(confint(.))") %>%
  dplyr::select(Source, Variable, OR, `Lower CI` = `2.5 %`, `Upper CI` = `97.5 %`) %>%
  bind_rows(or_tibble %>% 
              mutate(Source = "Delta Method") %>%
              dplyr::select(Source, Variable = Estimate, OR, `Lower CI`, `Upper CI`) %>%
              var_mapping(Variable)) %>%
  mutate(Source = factor(Source, levels = c("exp(confint(.))", "Delta Method"))) %>%
  filter(Variable != "Constant")
## Joining, by = "Variable"
ggplot(data = graph_data_2, aes(x = Variable, y = OR, ymin = `Lower CI`, ymax = `Upper CI`)) +
  geom_pointrange() +
  geom_hline(yintercept = 1, lty = 2) +
  coord_flip() +
  xlab("Variable") + ylab("Odds Ratio with 95% CI") + 
  facet_grid(~ Source)

Which is the correct method?

• Delta method provides a standard error for the odds ratio, which can be used to create a normal-approximated (i.e. symmetric) confidence interval.
• But delta method confidence intervals can also extend into negative territory.

What does Stata do?

• Stata reports standard errors for odds ratios determined by the delta method.
• But its 95% confidence intervals around the odds ratios are based on \(\mbox{exp}(\beta \pm 1.96*\mbox{SE}_{\beta})\).

That is, the standard error is the delta method, but the confidence intervals are equal to Rs exp(confint(model_object))!

Let’s export our data to Stata and take a look.

haven::write_dta(anes, "anes-with-noise.dta")
. logistic vote i.gender i.educ age noise

Logistic regression                             Number of obs     =      2,368
                                                LR chi2(7)        =     147.10
                                                Prob > chi2       =     0.0000
Log likelihood = -1565.3592                     Pseudo R2         =     0.0449

----------------------------------------------------------------------------------------
                  vote | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-----------------------+----------------------------------------------------------------
                gender |
               Female  |   .6881121     .05868    -4.38   0.000      .582199    .8132929
                       |
                  educ |
         Completed HS  |    1.92396   .4436975     2.84   0.005     1.224319    3.023412
    College < 4 Years  |   2.006222   .4372836     3.19   0.001     1.308723     3.07546
College 4 Year Degree  |   1.508235   .3344887     1.85   0.064     .9765486    2.329401
      Advanced Degree  |   .6541158   .1500084    -1.85   0.064     .4173001    1.025323
                       |
                   age |   1.015176   .0025378     6.03   0.000     1.010214    1.020162
                 noise |   1.083796   .8888249     0.10   0.922     .2172073    5.407803
                 _cons |    .349481   .0890344    -4.13   0.000     .2121143    .5758072
----------------------------------------------------------------------------------------

Comparing the standard errors from Stata to our delta method SEs, we find the two match.

But the delta method CIs do not match. In fact, Stata’s confidence intervals are close to R’s exp(confint()) results.

Understanding Stata’s output for logit models with results reported as odds ratios:

• The CIs are exponentiated \(\widehat{\beta_k} \pm 1.96 \times \mbox{SE}_{\widehat{\beta_k}}\).
• The p-values reported for the odds ratio come from \(z = \frac{\widehat{\beta_k}}{\mbox{SE}_{\widehat{\beta_k}}}\). To see this, compare the p-values with and without the or option to the .logit command.

If the sampling distribution is asymmetric, then the confidence interval should be asymmetric. We can use the nonparametric bootstrap to get a sense of the shape of the sampling distribution. The R code applied to the ANES data is:

or_boot_func <- function(form, data, indices) {
  
  data <- data[indices,]
  ors  <- exp(coef(glm(form, data = data, family = "binomial")))
  
  ors
}  
  
or_boot_values <- boot(data = anes, statistic = or_boot_func,
                       R = 1000,
                       form = vote ~ gender + educ + age + noise) 

boot_ci <- apply(or_boot_values$t, 2, function(i) {
  
  as_tibble(list(Estimate = mean(i),
                 Lower    = quantile(i, .025),
                 Upper    = quantile(i, .975)))
}) %>%
  reduce(bind_rows)
boot_ci %>%
  mutate(Variable = c("Constant", graph_data$Variable)) %>%
  filter(Variable != "Constant") %>%
  ggplot(aes(x = Variable, y = `Estimate`, ymin = `Lower`, ymax = `Upper`)) +
  geom_pointrange() +
  geom_hline(yintercept = 1, lty = 2) +
  coord_flip() +
  xlab("Variable") + ylab("Bootrapped Mean Odds Ratio and 2.5-97.5 Pecentiles") 

Delta Method Standard Errors for Predicted Probabilities

It makes most sense to stick with exp(confint(glm_object)) for confidence intervals around ORs, since this is a trivial task. But what about other quantities that are functions of model parameters, such as predicted probabilities?

Predicted probabilities are obtained from the results of a logit model by:

\[ \begin{array} \mbox{Pr}(y=1|x) &= F(x\beta) \\ & = \frac{\mbox{exp}(x\beta)}{1 + \mbox{exp}(x\beta)} \end{array} \]

What is the prediction for a 55-year-old male who finished high school but did not go to college (with average noise)? First, get \(x\beta\).

estimates <- coef(trump_model)

xb        <- estimates[1] + estimates[2]*0 + estimates[3]*1 + estimates[4]*0 + estimates[5]*0 + 
             estimates[6]*0 + estimates[7]*55 + estimates[8]*0
(xb <- as.numeric(round(xb, 3)))
## [1] 0.431

\[ \begin{array} \mbox{Pr}(y=1|x) &= \frac{\mbox{exp}(x\beta)}{1 + \mbox{exp}(x\beta)} \\ & = \frac{\mbox{exp}(0.431)}{1 + \mbox{exp}(0.431)} \\ & = 0.606 \end{array} \]

Equivalently.

round(plogis(xb),3)
## [1] 0.606

The predicted probabilty is a function of all parameter estimates, so we need to use the matrix version of the Taylor series approximation to get our SE.

\[ \mbox{Var}\left(f(\mathbf x)\right) = \nabla(\boldsymbol\mu )^{T}\mbox{Cov}(\mathbf x)\nabla(\boldsymbol\mu) \]

First, determine the gradient vector for the \(\beta_k\).

We’ve used \(F()\) to represent the standard logistic cdf. Let \(f()\) be the standard logistic pdf.

By the chain rule, and by the fact that the first derivative of a cdf is its pdf,

\[ \begin{array} \frac{\partial F(x'\beta)}{\partial \beta} &= \frac{\partial F((x' \beta))}{\partial x' \beta} \\ &= f(x'\beta)x \end{array} \]

Define the \(x\) vector for our subject of interest along with the beta vector from the model results.

x    <- c(1,0,1,0,0,0,55,0)
beta <- coef(trump_model)

cbind(beta,x)
##                                  beta  x
## (Intercept)               -1.05130623  1
## genderFemale              -0.37380350  0
## educCompleted HS           0.65438538  1
## educCollege < 4 Years      0.69625321  0
## educCollege 4 Year Degree  0.41094016  0
## educAdvanced Degree       -0.42447089  0
## age                        0.01506207 55
## noise                      0.08046982  0

And \(x'\beta\) is

xb <- t(x)%*%beta

The variance is thus:

\[ \begin{array} \mbox{Var}(p(\mbox{Vote} = \mbox{Trump})) &= \nabla(\boldsymbol\mu )^{T}\mbox{Cov}(\mathbf x)\nabla(\boldsymbol\mu) \\ &= f(x'\beta)x'\mbox{Cov}(\widehat{\beta})xf(x'\beta) \end{array} \]

Calculating “by hand” in R:

dlogis(xb)%*%t(x)%*%vcov(trump_model)%*%x%*%dlogis(xb) %>%
  sqrt
##            [,1]
## [1,] 0.02747443

The deltamethod function in the msm package will also do this for you. The arguments are

1. The function expressed as a formula
2. A vector of means, here \(\mu = \beta\)
3. The covariance matrix

deltamethod(
  ~ exp(x1 + x2*0 + x3*1 + x4*0 + x5*0 + x6*0 + x7*55 + x8*0)/
          (1 + exp(x1 + x2*0 + x3*1 + x4*0 + x5*0 + x6*0 + x7*55 + x8*0)), 
  coef(trump_model), vcov(trump_model))
## [1] 0.02747443

This is the standard error for our predicted probability.

One catch: because the derivatives are found symbolically, you are limited to what stats::deriv knows. For example, the following would be simpler to code.

deltamethod(~ plogis(x1 + x2*0 + x3*1 + x4*0 + x5*0 + x6*0 + x7*55 + x8*0), 
            coef(trump_model), vcov(trump_model))

But R throws an error, Function plogis is not in the derivatives table.

By the way, this is where Stata truly shines and why I still turn to it despite protestations that “R can do everything.” Here are all of the predicted probabilities, with standard errors, for all combinations of education and gender, for 55 year-olds.

. margins educ, over(gender) at(age = 55 noise = 0)


---------------------------------------------------------------------------------------
                      |            Delta-method
                      |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
----------------------+----------------------------------------------------------------
          gender#educ |
  Male#HS Incomplete  |   .4445065   .0514349     8.64   0.000     .3436959     .545317
    Male#HS Complete  |   .6062301   .0274744    22.07   0.000     .5523812     .660079
   Male#Some College  |   .6161789   .0205271    30.02   0.000     .5759466    .6564112
          Male#BA/BS  |   .5468739   .0231694    23.60   0.000     .5014627    .5922851
         Male#Higher  |    .343584   .0251349    13.67   0.000     .2943204    .3928475
Female#HS Incomplete  |      .3551   .0482039     7.37   0.000     .2606221    .4495779
  Female#HS Complete  |   .5144184   .0278629    18.46   0.000     .4598081    .5690286
 Female#Some College  |   .5248688   .0201584    26.04   0.000      .485359    .5643786
        Female#BA/BS  |   .4536941   .0231244    19.62   0.000     .4083711    .4990171
       Female#Higher  |   .2648002   .0215839    12.27   0.000     .2224966    .3071038
---------------------------------------------------------------------------------------

Type one more word, and this becomes a pretty graph.

. marginsplot

Again, the delta method gives us an approximation that may not be accurate.

Predicted probabilities are bounded by zero and one, yet a delta method CI may extend below or above these limits.

Alternatives:

1. Boostrapping
2. Simulations by repeated draws from \(\boldsymbol \beta \sim \mathcal{N}\left(\widehat{\boldsymbol \beta}, \mbox{Cov}(\widehat{\boldsymbol \beta})\right)\).

Both are computationally intensive and may be less attractive for biger data sets than the delta method.