In R, How can I calculate the elasticity of Y with respect to X, when Y is binary?

Question

I have a dataset that I am doing in R and I need to calculate elasticities in it. To simplify my model, I have Y = XB + u, and I need to find the elasticity of Y with respect to X.

My investigation led me to the margins package in R, which seems to emulate the margins command in Stata. However, from what I understand of the Stata description of margins, dy dx and ey ex measure two different things - the former measures the marginal effect, whereas the latter measures the elasticity.

What I want to do is calculate the elasticity of a linear probability model, so I cannot just directly take logs on both sides and look at the coefficient, since the LHS is a binary variable. In Stata, again, my understanding is that this is achieved with the Margins package.

So I guess, at the end of the day, my question is: how do I calculate the elasticity of Y with respect to X in a linear probability model and/or a probit, using R?

dimitriy · Answer 1 · 2018-11-01T15:15:42.047

The short answer is that I don't think there is a way to do what you want with the margins package in R.

However, below I will:

Demonstrate how Stata is calculating the AME elasticity with logit and OLS/LPM
Replicate the AME point estimates in R for logit and LPM by hand
Use the bootstrap to get the standard error of the logit AME

I am using the rsource command from SSC to run R from within Stata, but you can just cut the R code out and paste it directly into R after exporting the CSV, though installing rsource is a great investment that future you will appreciate.

Here's the output (code and some interpretation below):

. /* Stata Code */
. sysuse auto, clear
(1978 Automobile Data)

. export delimited mpg foreign using "~/Desktop/cars.csv", replace
file /Users/dimitriy/Desktop/cars.csv saved

. 
. /* Logit and Elasticity */
. qui logit foreign mpg, nolog

. margins, eyex(mpg) // elasticity

Average marginal effects                        Number of obs     =         74
Model VCE    : OIM

Expression   : Pr(foreign), predict()
ey/ex w.r.t. : mpg

------------------------------------------------------------------------------
             |            Delta-method
             |      ey/ex   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         mpg |   2.225875    .713342     3.12   0.002     .8277506       3.624
------------------------------------------------------------------------------

. 
. predict double phat // same elasticity by hand
(option pr assumed; Pr(foreign))

. gen double ME  = _b[mpg]*phat*(1-phat)*mpg/phat

. gen double ME2 = _b[mpg]*(1-phat)*mpg // same as above since phat cancels

. sum ME ME2

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
          ME |         74    2.225875    .3252027   .6704181   2.472696
         ME2 |         74    2.225875    .3252027   .6704181   2.472696

. drop ME ME2 phat

. 
. /* LPM and Elasticity */
. qui reg foreign mpg, robust
. margins, eyex(mpg)      // elasticity

Average marginal effects                        Number of obs     =         74
Model VCE    : Robust

Expression   : Linear prediction, predict()
ey/ex w.r.t. : mpg

------------------------------------------------------------------------------
             |            Delta-method
             |      ey/ex   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         mpg |   4.284932    21.2706     0.20   0.841    -38.11723     46.6871
------------------------------------------------------------------------------

. 
. predict double phat // same elasticity by hand
(option xb assumed; fitted values)

. gen double ME = _b[mpg]*mpg/phat

. sum ME

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
          ME |         74    4.284932    9.259253   1.403936    58.9375

. drop ME phat

. 
. 
. /* R equivalent for AMEs */
. rsource, terminator(END_OF_R)
Unknown #command
Unknown #command
Assumed R program path: "/usr/local/bin/R"

Beginning of R output
> 
> library(foreign)
> cars<-read.csv("~/Desktop/cars.csv")
> logit<-glm(foreign~mpg,data=cars, family="binomial")
> (AME_logit<-mean( coef(logit)["mpg"] * logit$fit * (1-logit$fit) * cars$mpg/logit$fit ))
[1] 2.225875
> 
> cars$foreign <- as.numeric(cars$foreign)-1
> lpm<-glm(foreign~mpg,data=cars, family="gaussian")
> (AME_lpm<-mean( coef(lpm)["mpg"] *cars$mpg/lpm$fit))
[1] 4.284932
> 
> 
> set.seed(10312018)
> B <- 10000
> AME <- c()
> for (b in 1:B){
+   samp_b = sample.int(nrow(cars), replace=TRUE)
+   logit<-glm(foreign~mpg,data=cars[samp_b,], family="binomial")
+   AME_logit<-mean(coef(logit)["mpg"] * logit$fit * (1-logit$fit) * cars[samp_b,]$mpg/logit$fit )
+   AME <- rbind(AME, AME_logit)
+ }
> 
> AME_logit                 # logit AME from above
[1] 2.037947
> sd(AME)                   # BS SE of logit AME (compare to Stata margins)
[1] 0.8071795
> quantile(AME, c(.05,.95)) # BS 95% CI for logit AME (compare to Stata margins)
      5%      95% 
1.211233 3.744386 
> 
> 
End of R output

We are using the standard elasticity formula

$$ \epsilon = \frac{\partial E[y \vert x]}{\partial x} \cdot \frac{x}{E[y \vert x]}$$

With the logit, the first term for a continuous covariate $x_k$ is

$$\Lambda(X'\beta)\cdot \left[1-\Lambda(X'\beta)\right]\cdot\beta_k,$$ where $$\Lambda(z)=\frac{\exp{z}}{1+\exp{z}}=\hat p.$$

LPM/OLS is much simpler because of the linearity, so you just get back the coefficient when you calculate $\frac{\partial E[y \vert x]}{\partial x}=\beta_k$.

The R bootstrap gives us a standard error of 0.8071795 for the logit AME, compared to Stata's 0.713342, and a 95% CI of [1.211233,3.744386] compared to Stata's [.8277506, 3.624]. Stata is using the delta method to calculate the SE/CI, which is the asymptotic approximation to the finite sample distribution that the bootstrap is producing. Here with only 74 observations the difference is still noticeable, but hopefully you have more data to work with.

The interpretations of the logit point estimate is that for a 1% increase in MPG, you get a 4.3% bump in probability that the car is foreign. The LPM estimate is about half that. With continuous regressors, you will often get fairly different estimates between the two types of model. Note that these are not percentage point changes in probability that margins, dydx would calculate, which would be much smaller here.

Stata With Embedded R Code:

/* Stata Code */
sysuse auto, clear
export delimited mpg foreign using "~/Desktop/cars.csv", replace

/* Logit and Elasticity */
qui logit foreign mpg, nolog
margins, eyex(mpg) // elasticity

predict double phat // same elasticity by hand
gen double ME  = _b[mpg]*phat*(1-phat)*mpg/phat
gen double ME2 = _b[mpg]*(1-phat)*mpg // same as above since phat cancels
sum ME ME2
drop ME ME2 phat

/* LPM and Elasticity */
qui reg foreign mpg, robust
margins, eyex(mpg)  // elasticity

predict double phat // same elasticity by hand
gen double ME = _b[mpg]*mpg/phat
sum ME
drop ME phat


/* R equivalent for AMEs */
rsource, terminator(END_OF_R)

library(foreign)
cars<-read.csv("~/Desktop/cars.csv")
logit<-glm(foreign~mpg,data=cars, family="binomial")
(AME_logit<-mean( coef(logit)["mpg"] * logit$fit * (1-logit$fit) * cars$mpg/logit$fit ))

cars$foreign <- as.numeric(cars$foreign)-1
lpm<-glm(foreign~mpg,data=cars, family="gaussian")
(AME_lpm<-mean( coef(lpm)["mpg"] *cars$mpg/lpm$fit))


# Bootstrap the AME to get the SE and CI for AME
# This is finite sample version of what margins, eyex() is doing using the delta-method in Stata  
set.seed(10312018)
B <- 10000
AME <- c()
for (b in 1:B){
  samp_b = sample.int(nrow(cars), replace=TRUE) 
  logit<-glm(foreign~mpg,data=cars[samp_b,], family="binomial")
  AME_logit<-mean(coef(logit)["mpg"] * logit$fit * (1-logit$fit) * cars[samp_b,]$mpg/logit$fit )
  AME <- rbind(AME, AME_logit)
}  

AME_logit                 # logit AME from above
sd(AME)                   # BS SE of logit AME (compare to Stata margins)
quantile(AME, c(.05,.95)) # BS 95% CI for logit AME (compare to Stata margins)

END_OF_R

In R, How can I calculate the elasticity of Y with respect to X, when Y is binary?

1 Answers1