Computing the standard error when dividing coefficients of different regressions in R

Question

Consider the following two regressions from the same dataset mtcars.

#load the data
data(mtcars)


# Run the regression
model1<-lm(mpg~cyl+gear+drat, data = mtcars)
model2<-lm(wt~cyl+gear+drat, data = mtcars)

summary(model1)
summary(model2)

#Calculate ratio of coefficients
g<-model1$coefficients[2] / model2$coefficients[2]

#calculate clustered standard errors
vcov<-cluster.vcov(model1, mtcars$vs)
coeftest(model1, vcov)
vcov<-cluster.vcov(model2, mtcars$vs)
coeftest(model2, vcov)

We observe that the ratio of the cyl variable in the two regressions is equal to -8.16. Now I would like to calculate the standard error that corresponds with this ratio (the clustering of my standard errors here does not make much sense, but it just provides us with a variance covariance matrix for both models, which we may need). Stata has a command called "nlcom" that can do this, but I cannot find a similar command in R. Does anyone of you know whether it exists? If not, then how should I do this? I appreciate any help.

Answer 1

As @MDEWITT suggests, you could use the delta method in the article, though I think you would probably need to estimate the model differently - you would probably need a single multivariate model rather than two independent regression models because you need the covariance of the two coefficients, which doesn't exist unless you estimate a joint model.

Joint model with deltamethod for computing SE:

library(msm)
data(mtcars)
# Run the regression
model<-lm(cbind(mpg, wt)~cyl+gear+drat, data = mtcars)
b <- c(coef(model))
v <- vcov(model)
## calcualte se
est <- b[2]/b[6]
se <- deltamethod(g = ~ x2/x6, b, vcov(model))
> est 
# [1] -8.160363
> se
# [1] 1.770336

There are other methods, too.

A non-parametric bootstap could be used:

## write a function to calculate the statistic of interest
boot.fun <- function(data, inds){
  m <- lm(cbind(mpg, wt) ~ cyl + gear + drat, data=data[inds, ])
  # return the appropraite ratio
  coef(m)[2,1]/coef(m)[2,2]
}

library(boot)
## bootstrap the function
out <- boot(mtcars, statistic=boot.fun, R=5000)
## calculate confidence intervals
boot.ci(out, type=c("perc", "bca"))
# BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
# Based on 5000 bootstrap replicates
# 
# CALL : 
#   boot.ci(boot.out = out, type = c("perc", "bca"))
# 
# Intervals : 
#   Level     Percentile            BCa          
# 95%   (-12.993,  -5.152 )   (-12.330,  -4.838 )  
# Calculations and Intervals on Original Scale

You could also use a parametric boostrap:

## draw coefficients from implied sampling distributions
B <- MASS::mvrnorm(5000, c(coef(model)), vcov(model))
## calculate the ratio for each draw
rat <- B[,2]/B[,6]
## calculate the confidence itnerval
round(c(mean(rat), quantile(rat, c(.025,.975))), 3)
#           2.5%   97.5% 
# -8.559 -14.418  -5.554 

The three methods generate these three confidence intervals:

q1 <- est + qt(c(.025,.975), df=model$df.residual)*se
q2 <- boot.ci(out, type="perc")$percent[1,4:5]
q3 <- quantile(rat, c(.025,.975))
cis <- rbind(q1, q2, q3)
colnames(cis) <- c("lwr", "upr")
rownames(cis) <- c("Delta Method", "BS (non-parametric)", "BS (parametric)")
cis
#                           lwr       upr
# Delta Method        -11.78673 -4.533994
# BS (non-parametric) -12.99338 -5.151545
# BS (parametric)     -14.41812 -5.553773

Both the parametric bootstrap distributions (these are for the middle 99% of the data) are skewed left, so a normal-theory confidence interval based on the estimate and standard error may well be inappropriate.

par(mfrow=c(1,2))
hist(out$t[which(out$t > quantile(out$t, .005) & out$t < quantile(out$t, .995))], main="Non-parametric Bootstrap", xlab="ratio")
hist(rat[which(rat < quantile(rat, .995) & rat > quantile(rat, .005))], main="Parametric Boostrap", xlab="ratio")