如何在使用felm()函数后绘制交互作用的边际效应

Question

如何在使用felm()函数后绘制交互作用的边际效应

3

我基于一个包含大量单位固定效应的“巨型”面板数据运行了回归。因此，我使用了“lfe”包中的“felm()”函数。此外，在回归中有两个连续变量的交互项。但是，当绘制x对y的边际效应如何随x2变化时，使用“felm()”生成的对象似乎与大多数绘图函数（例如“ggplot”、“interplot()”和“meplot”）不兼容。但是我必须使用“felm()”，因为我需要控制大量的单位固定效应（就像Stata中的“reghdfe”一样）。那么，我该如何在R中解决这些问题？请随意让我知道一些方法。谢谢！

以下是interplot()无法与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)

然后出现错误：

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

我也尝试了meplot()，但仍无法正常工作：

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

我遇到了一个错误：

Error in meplot(model = reg1, var1 = "x", var2 = "x2", int = "x:x2", vcov = vcov(reg1),  : 
  (list) object cannot be coerced to type 'double'

- Ronald

1

看一下 help(lfe, pac=lfe)，我看到了一个实例。你能否使用它来创建一个可以重现你的问题的 data 对象？目前你的问题无法重现，因为只有你能看到 data。顺便说一句，在 R 中，名称 data 是基本函数的名称，因此将其用作对象名称是不被赞同的。 - IRTFM

谢谢，42。我已经编辑了上面的代码，并使用了你提到的工作示例。但效果并不太好。我认为问题仍然存在，因为felm()函数与许多绘图函数不兼容。 - Ronald

1

呃，你这样写代码：data.frame(cbind(y, x, x2, id, firm, year))。这会把所有的因子转换成数字类型。我不惊讶你遇到了麻烦。data.frame(cbind(...)) 这种策略经常是 R 初学者问题的根源。扔掉你从中复制的书吧。 - IRTFM

即使我不使用 data.frame，只是使用像 y、x、x2 等创建的对象，也会发生同样的事情。 - Ronald

非常感谢 @42。我已经使用 ggplot2 实现了自己想要的，但是手动绘制它。 - Ronald

显示剩余3条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ronald · Accepted Answer

我使用了 ggplot2、coef() 和 vcov() 来实现我想要的结果，通过手动绘制边际效应图。

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)