3

I ran a regression based on a "giant" panel data with a bunch of unit fixed effects. So I employed function "felm()" from package "lfe". In addition, I have an interaction term of two continuous variables in the regression. But when plotting how the marginal effects of x on y vary with x2, it seems that the objects produced by "felm()" are often incompatible to most plotting functions like "ggplot", "interplot()" and "meplot". But I have to use "felm()" because I need to control for a large amount of unit fixed effects (like people do by "reghdfe" in Stata). So, how could I address this issues in R? Feel free to let me know some ways out. Thanks!

Here is an example about how interplot() does not work with felm():

# An example data:
library(lfe)
library(ggplot2)
library(interplot)
oldopts <- options(lfe.threads=1)
x <- rnorm(1000)
x2 <- rnorm(length(x))
id <- factor(sample(10,length(x),replace=TRUE))
firm <- factor(sample(3,length(x),replace=TRUE,prob=c(2,1.5,1)))
year <- factor(sample(10,length(x),replace=TRUE,prob=c(2,1.5,rep(1,8))))
id.eff <- rnorm(nlevels(id))
firm.eff <- rnorm(nlevels(firm))
year.eff <- rnorm(nlevels(year))
y <- x + 0.25*x2 + id.eff[id] + firm.eff[firm] +
  year.eff[year] + rnorm(length(x))
mydata <- data.frame(cbind(y, x, x2, id, firm, year))

# Regression using felm():
reg1 <- felm(y ~ x + x2 + x:x2|id+firm+year|0|id, data=mydata)
summary(reg1)

# Using interplot() to plot marginal effects
interplot(m=reg1, var1="x", var2="x2", ci=0.9)

Then errors appear:

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘sim’ for signature ‘"felm"’

Also I tried meplot() but it still does not work:

# Using meplot() to plot marginal effects
library(evir)
meplot(model=reg1, var1="x", var2="x2", int="x:x2", vcov=vcov(reg1), data=mydata)

I got an error:

Error in meplot(model = reg1, var1 = "x", var2 = "x2", int = "x:x2", vcov = vcov(reg1),  : 
  (list) object cannot be coerced to type 'double'
jay.sf
  • 60,139
  • 8
  • 53
  • 110
Ronald
  • 61
  • 6
  • 1
    Looking at `help(lfe, pac=lfe)` I see a worked example. Are you able to use that to make a `data` object that reproduces your issues? At the moment your issues are not reproducible sinc eonly you can see `data`. BTW in R the name `data` refers to a base function so it's use as an object name is frowned upon. – IRTFM Dec 07 '19 at 17:41
  • Thanks, 42. I have edited my codes above and used the worked example you mentioned. But it does not work as well. I think the problem is still there as felm() is not compatible to many plotting functions. – Ronald Dec 07 '19 at 20:40
  • 1
    Ewww. You did: `data.frame(cbind(y, x, x2, id, firm, year))`. That coerced all your factors to numeric. I'm not surprised you had trouble. The `data.frame(cbind(...))` strategy is often the cause of R-noob problems. Throw away the book you copied it from. – IRTFM Dec 07 '19 at 21:49
  • Even though I do not use data.frame and just use created objects like y, x, x2, and etc, the same thing happens. – Ronald Dec 07 '19 at 22:58
  • I would have not used `cbind`. The data.frame function will accept a set of vectors and give the columns appropriate names. Using cbind is what is causing the most damage because matrices in R are not capable of holding varying classes. However, using just `data.frame` does not solve the problem with `interplot` or for that matter with `sjPlot` which I also tried) since they have no methods for models of class "felm". I think you may just need to use coef() and do your own plotting. – IRTFM Dec 08 '19 at 00:14
  • @42 Yes, they are not compatible. Just plotting the marginal effects by hands using `coef()` is a good idea, even though it might not look that fancy like `ggplot`, it also demonstrates exactly what I want. Thanks so much for your time and help! – Ronald Dec 08 '19 at 00:35
  • I'm sure you could coerce ggplot using a dataframe input and it's annotation capacities, but it might take a bit of fiddling. You might want to do a search on: `"[r] [ggplot] interaction plot"` – IRTFM Dec 08 '19 at 01:46
  • Thank you ssssoo much @42. I have realized what I wanted using `ggplot2` but plotting it by hands. – Ronald Dec 09 '19 at 03:55

1 Answers1

3

I have used ggplot2, coef() and vcov() to realize what I want, plotting the marginal effects by hands.

library(ggplot2)
beta.hat <- coef(reg1)
vcov1 <- vcov(reg1)
z0 <- seq(min(x2), max(x2), length.out = 1000)
dy.dx <- beta.hat["x"] + beta.hat["x:x2"]*z0
se.dy.dx <- sqrt(vcov1["x", "x"] + (z0^2)*vcov1["x1:x2", "x1:x2"] + 2*z0*vcov1["x", "x1:x2"])
upr <- dy.dx + 1.96*se.dy.dx
lwr <- dy.dx - 1.96*se.dy.dx
ggplot(data=NULL, aes(x=z0, y=dy.dx)) +
  labs(x="x2", y="Marginal Effects",
       title=paste("Marginal Effects of x on y vary with x2"), 
       cex=4) +
  geom_line(aes(z0, dy.dx),size = 1) +
  geom_line(aes(z0, lwr), size = 1, linetype = 2, color="blue") +
  geom_line(aes(z0, upr), size = 1, linetype = 2, color="blue") +
  geom_hline(yintercept=0, size = 1, linetype=3) +
  geom_ribbon(aes(ymin=lwr, ymax=upr), alpha=0.3)
Ronald
  • 61
  • 6
  • Hi @Ronald - I faced the same problems as you did. I tried to copy your code for plotting marginal effects. But I see some errors. Below is my code - `model.1_AL.I <- felm(AUDITLAG ~ Ov_deg_pct*AGE_CEO_D +CONTROLS|factor(sic2)+factor(year)|0|TICKER, data = CEO_DEMO_FINAL)` Your code works ok except this line - `se.dy.dx <- sqrt(vcov1["Ov_deg_pct", "Ov_deg_pct"] + (z0^2)*vcov1["AGE_CEO_D", "AGE_CEO_D"] + 2*z0*vcov1["Ov_deg_pct", "AGE_CEO_D"])` The error is - `Error in vcov1["AGE_CEO_D", "AGE_CEO_D"] ` Do you have any idea how can I solve it? – Sharif May 23 '20 at 21:04
  • Hi @Sharif - I'm not sure where it goes wrong, since I don't have your data. I usually use ":" to add interaction terms. For example, I would use `model.1_AL.I <- felm(AUDITLAG ~ Ov_deg_pct + AGE_CEO_D + Ov_deg_pct:AGE_CEO_D + CONTROLS | factor(sic2)+factor(year) | 0 | TICKER, data = CEO_DEMO_FINAL)`. Have you tried this? Sorry I've been a long time absent in StackOverflow..... – Ronald Nov 28 '21 at 19:50
  • Hi again @Sharif - I think I know what's going wrong. See my new codes above. – Ronald Nov 28 '21 at 20:40