# The standard error of the ratio of two regression coefficients

One question which I find occurs repeatedly around the department is on finding confidence intervals for various transformations of coefficients from a regression. One of the most common, (and problematic), transformation asked about is the ratio of two regression coefficients. This is, in fact, a very old problem, and, as with many things in regression, exact form solutions are impossible without resort to either the bootstrap or assumptions of normality.

If one is willing to assume normality, then I believe the following logic should apply.
Consider the standard regression set up:

$y_i=\beta_0+x_i\beta_1+...+\beta_k+\epsilon_i$

Written more compactly as $y=X\beta+epsilon$, if $\epsilon_i \sim N(0,sigma^2)$, then we have following the standard regression derivations, that

$\hat\beta \sim N(0,\sigma^2_{\epsilon}(X'X)^{-1})$

Let us suppose our ratio of interest is $\gamma=\frac{\beta_1}{\beta_2}$, and let us suppose we are interested in knowing the the distribution of the plug-in estimator $\hat\gamma=\frac{\hat\beta_1}{\hat\beta_2}$
All we are asking here, as it turns out, is “what is the ratio of two standard normal random variables with known mean, variance, and correlation?” This question was solved in 1969 by David Hinkley, who would later go on to be an important contributor in the development of the bootstrap, in his appropriately named paper

Unsurprisingly, perhaps, given the complexity of the distribution, I have had little fortune finding actual implementations of the algorithm, barring some old C code from a retired medical researcher who solved the same problem in a different way at around the same time. Modern statisticians no longer need mathematics in the same way our forbearers did, and most practitioners to today would opt for computing the confidence intervals via monte carlo simulation. Here is a simple R function which I find performs the estimation quite well

It takes an object returned by lm() as an output as well as the name of the numerator and denominator variables of interest as its primary arguments.



RatioCI <- function(lmoutput,numerator,denominator,replications=1000, confidence.level=.95){
#A simple monte-carlo algorithm for the CI of two regression coefficients.
#Args:
#lmoutput: an object return from running an OLS regression with LM
# numerator: the name of the numerator variable (as string, put "(Intercept)" for intercept)
#denominator: the name of the denominator variable (as string, put "(Intercept)" for intercept)
#replications:  The number of monte-carlo replications to be used in estimating the confidence interval
#confidence.level: the size of the confidence interval to be computed (the computed interval should capture the true ratio confidence.level precent of the time)
#Returns:
#CIout: an estimate of our confidence interval (numeric vector, two elements)
#
require("MASS")
variance1 <- vcov(lmoutput)[numerator,numerator]
variance2 <- vcov(lmoutput)[denominator,denominator]
covar <- vcov(lmoutput)[numerator,denominator]
mean.plugin <- lmoutput[["coefficients"]][numerator] /  lmoutput[["coefficients"]][denominator]
Sigma.plugin <- matrix(c(variance1,covar,covar,variance2),nrow=2)
distribution <-mvrnorm(replications, c(lmoutput[["coefficients"]][numerator], lmoutput[["coefficients"]][denominator] ), Sigma.plugin)
ratio <- sort(distribution[,1] / distribution[,2])
confidence.band <- c( floor( replications * (1 - confidence.level) / 2) + 1, ceiling( replications * (1 + confidence.level) / 2) )
CIout <- ratio[ confidence.band ]
return(CIout)

}


As a quick consistency check, monte-carlo simulations run on a group of correlated independent variables show a 94.5% coverage probability for a 95% confidence interval with 1000 replications.