r - lfe felm 回归对象的边际效应/交互图

Question

我需要为固定效应模型创建一个交互/边际效应图，包括使用 lfe“felm”命令生成的聚集标准误差。

我已经创建了一个实现这一点的函数。但是，在我开始使用它之前，我想仔细检查这个函数是否被正确指定。请在下面找到函数和可重现的示例。

library(lfe)

### defining function
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.05, se = NULL) {
        library(ggplot2)
        library(ggthemes)
        library(gridExtra)

        ### defining function to get average marginal effects
        getmfx <- function(betas, data, treatment, moderator) {
                betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
        }

        ### defining function to get marginal effects for specific levels of the treatment variable
        getmfx_high_low <- function(betas, data, treatment, moderator, treatment_val) {
                betas[treatment] * treatment_val + betas[paste0(treatment, ":", moderator)] * data[, moderator] * treatment_val
        }

        ### Defining function to analytically derive standard error for marginal effects
        getvarmfx <- function(my_vcov, data, treatment, moderator) {
                my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
        }

        ### constraining data to relevant variables
        data <- data[, c(treatment, moderator)]

        ### getting marginal effects
        data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)

        ### getting  marginal effects for high and low cases of treatment variable
        data[, "marginal_effects_treatment_low"] <- getmfx_high_low(coef(regression_model), data, treatment, moderator, quantile(data[,treatment], 0.05))
        data[, "marginal_effects_treatment_high"] <- getmfx_high_low(coef(regression_model), data, treatment, moderator, quantile(data[,treatment], 0.95))

        ### getting robust SEs
        if (is.null(se)) {
                data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
        } else if (se == "clustered") {
                data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
        } else if (se == "robust") {
                data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
        }

        ### Getting CI bounds
        data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
        data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)

        ### plotting marginal effects plot
        p_1 <- ggplot(data, aes_string(x = moderator)) +
                geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
                geom_line(aes(y = marginal_effects)) +
                theme_fivethirtyeight() +
                theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
                geom_rug() + 
                xlab(moderator_translation) +
                ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
                ggtitle("Average marginal effects")

        p_2 <- ggplot(data, aes_string(x = moderator)) +
                geom_line(aes(y = marginal_effects_treatment_high, color = paste0("High ",treatment_translation))) +
                geom_line(aes(y = marginal_effects_treatment_low, color = paste0("Low ",treatment_translation))) +
                theme_fivethirtyeight() + 
                theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10), axis.title.y = element_blank(), legend.justification = c(0.95, 0.95), legend.position = c(1, 1), legend.direction = "vertical") +
                geom_rug() + 
                xlab(moderator_translation) +
                ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
                ggtitle("Marginal effects at high / low levels of treatment") +
                scale_color_manual(name = NULL, values = c(rgb(229, 93, 89, maxColorValue = 255), rgb(75, 180, 184, maxColorValue = 255)), labels=c(paste0("High ",treatment_translation), paste0("Low ",treatment_translation)))

        ### exporting plots as combined grob
        return(grid.arrange(p_1, p_2, ncol = 2))
}

### example:
# example model (just for demonstration, fixed effects and cluster variables make little sense here)
model <- felm(mpg ~ cyl + am + cyl:am | carb | 0 | cyl, data = mtcars)

# creating marginal effects plot
felm_marginal_effects(regression_model = model, data = mtcars, treatment = "cyl", moderator = "am", treatment_translation = "Number of cylinders", moderator_translation = "Transmission", dependent_variable_translation = "Miles per (US) gallon")

示例输出如下所示：

很高兴有任何关于如何使它成为更好、“编码良好”、快速的功能的建议，以便以后对其他人更有用。但是，我主要是想首先确认它是否“正确”。

此外，我想与社区再次核对一些剩余的问题，特别是：

• 我可以使用我为“高”和“低”处理案例的平均边际效应生成的标准误差，还是我需要为这些案例生成不同的标准误差？如果有怎么办？

• 除了使用分析得出的标准误差，我还可以通过基于重复的数据子样本创建许多系数估计来计算自举标准误差。我将如何为高/低情况生成自举标准错误？

• 固定效应模型或具有聚集标准误差的固定效应模型是否存在使边际效应图或我在代码中所做的任何其他事情从根本上不可接受的东西？

PS.：上述函数和问题是如何在 felm() 函数之后绘制交互的边际效应的扩展

标签： rlfe