Граничные эффекты / графики взаимодействия для объекта регрессии lfe felm - PullRequest
Новая
       57

Граничные эффекты / графики взаимодействия для объекта регрессии lfe felm

0 голосов
/ 23 марта 2020

Мне нужно создать график взаимодействия / предельных эффектов для модели с фиксированными эффектами, включая кластерные стандартные ошибки, сгенерированные с помощью команды 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")

Пример вывода выглядит следующим образом: enter image description here

Рад любым советам о том, как это сделать. лучшая, «хорошо закодированная», быстрая функция, чтобы впоследствии она была более полезной для других. Тем не менее, я в основном пытаюсь подтвердить, является ли он «правильным» в первую очередь.

Кроме того, я хотел бы проверить у сообщества вопрос о некоторых оставшихся вопросах, в частности:

  • Могу ли я использовать стандартные ошибки, сгенерированные мной, для средних предельных эффектов для случаев «высокого» и «низкого» лечения, или мне нужно генерировать разные стандартные ошибки для этих случаев? Если да, то как?

  • Вместо использования аналитически выведенных стандартных ошибок, я мог бы также рассчитать загрузочные стандартные ошибки, создав множество оценок коэффициентов на основе повторных выборок данных. Как бы я генерировал стандартные погрешности при загрузке для верхнего / нижнего регистра?

  • Есть ли что-то в моделях с фиксированными эффектами или моделях с фиксированными эффектами с кластерными стандартными ошибками, которые создают графики предельных эффектов или что-то еще, что я неужели в коде принципиально недопустимо?

PS .: Вышеприведенные функции и вопросы являются своего рода расширением Как изобразить предельный эффект взаимодействия после функции felm ()

...