用置信区间绘制回归系数

Question

假设我有 2 个数据框，一个用于 2015 年，一个用于 2016 年。我想为每个数据框运行回归，并绘制每个回归的系数之一及其各自的置信区间。例如：

set.seed(1020022316)
library(dplyr)
library(stargazer)

df16 <- data.frame(
  x1 = rnorm(1000, 0, 2),
  t = sample(c(0, 1), 1000, T),
  e = rnorm(1000, 0, 10)
) %>% mutate(y = 0.5 * x1 + 2 * t + e) %>%
  select(-e)

df15 <- data.frame(
  x1 = rnorm(1000, 0, 2),
  t = sample(c(0, 1), 1000, T),
  e = rnorm(1000, 0, 10)
) %>% mutate(y = 0.75 * x1 + 2.5 * t + e) %>%
  select(-e)

lm16 <- lm(y ~ x1 + t, data = df16)

lm15 <- lm(y ~ x1 + t, data = df15)

stargazer(lm15, lm16, type="text", style = "aer", ci = TRUE, ci.level = 0.95)

我想用各自的 .95 CI 绘制 t=1.558, x=2015 和 t=2.797, x=2016。这样做的最佳方法是什么？

我可以“手动”完成，但我希望有更好的方法。

library(ggplot2)
df.plot <-
  data.frame(
    y = c(lm15$coefficients[['t']], lm16$coefficients[['t']]),
    x = c(2015, 2016),
    lb = c(
      confint(lm15, 't', level = 0.95)[1],
      confint(lm16, 't', level = 0.95)[1]
    ),
    ub = c(
      confint(lm15, 't', level = 0.95)[2],
      confint(lm16, 't', level = 0.95)[2]
    )
  )
df.plot %>% ggplot(aes(x, y)) + geom_point() +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.1) + 
  geom_hline(aes(yintercept=0), linetype="dashed")

---

---

最佳：图形质量（看起来不错）、代码优雅、易于扩展（超过 2 个回归）

Answer 1

评论有点太长了，所以我把它作为部分答案发布。

从您的帖子中不清楚您的主要问题是将数据转换成正确的形状，还是绘图本身。但只是为了跟进其中一个 cmets，让我向您展示如何使用 dplyr 和 broom 运行多个模型，这使得绘图变得容易。考虑mtcars-dataset：

 library(dplyr)
 library(broom)
 models <- mtcars %>% group_by(cyl) %>% 
           do(data.frame(tidy(lm(mpg ~ disp, data = .),conf.int=T )))

 head(models) # I have abbreviated the following output a bit

    cyl        term estimate std.error statistic   p.value conf.low conf.high
  (dbl)       (chr)    (dbl)     (dbl)     (dbl)     (dbl)    (dbl)     (dbl)
     4 (Intercept)  40.8720    3.5896     11.39 0.0000012   32.752  48.99221
     4        disp  -0.1351    0.0332     -4.07 0.0027828   -0.210  -0.06010
     6 (Intercept)  19.0820    2.9140      6.55 0.0012440   11.591  26.57264
     6        disp   0.0036    0.0156      0.23 0.8259297   -0.036   0.04360

您会发现，这会在一个漂亮的数据框中为您提供所有系数和置信区间，这使得使用ggplot 进行绘图更容易。例如，如果您的数据集具有相同的内容，您可以为它们添加一个年份标识符（例如df1$year <- 2000; df2$year <- 2001 等），然后将它们绑定在一起（例如使用bind_rows，您可以使用bind_rows 的.id选项）。那么你可以在上面的例子中使用年份标识符而不是cyl。

那么绘图就很简单了。要再次使用mtcars 数据，让我们仅绘制disp 的系数（尽管您也可以使用faceting、grouping 等）：

 ggplot(filter(models, term=="disp"), aes(x=cyl, y=estimate)) + 
          geom_point() + geom_errorbar(aes(ymin=conf.low, ymax=conf.high))

使用您的数据：

 df <- bind_rows(df16, df15, .id = "years")

 models <- df %>% group_by(years) %>% 
           do(data.frame(tidy(lm(y ~ x1+t, data = .),conf.int=T ))) %>%
           filter(term == "t") %>% 
           ggplot(aes(x=years, y=estimate)) + geom_point() + 
           geom_errorbar(aes(ymin=conf.low, ymax=conf.high)) 

请注意，您只需将越来越多的数据绑定到主数据框即可轻松添加越来越多的模型。如果要绘制多个系数，也可以轻松使用faceting、grouping 或 position-dodgeing 来调整相应绘图的外观。

Answer 2

这是我现在的解决方案：

gen_df_plot <- function(reg, coef_name){
  df <- data.frame(y = reg$coefficients[[coef_name]],
                   lb = confint(reg, coef_name, level = 0.95)[1],
                   ub = confint(reg, coef_name, level = 0.95)[2])
  return(df)
}

df.plot <- lapply(list(lm15,lm16), gen_df_plot, coef_name = 't')

df.plot <- data.table::rbindlist(df.plot)

df.plot$x <- as.factor(c(2015, 2016))

df.plot %>% ggplot(aes(x, y)) + geom_point(size=4) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.1, linetype="dotted") + 
  geom_hline(aes(yintercept=0), linetype="dashed") + theme_bw()

我不喜欢它，但它有效。

Answer 3

这可能是通用代码。我对“x”的定义方式进行了更改，这样您就不必担心因子的字母顺序重新排序。

#
# Paul Gronke and Paul Manson
# Early Voting Information Center at Reed College
#
# August 27, 2019
#
#
# Code to plot a single coefficient from multiple models, provided
# as an easier alternative to "coefplot" and "dotwhisker". Some users
# may find those packages more capable
#
# Code adapted from https://stackoverflow.com/questions/35582052/plot-regression-coefficient-with-confidence-intervals


# gen_df_plot function will create a tidy data frame for your plot
#   Currently set up to display 95% confidence intervals

gen_df_plot <- function(reg, coef_name){
  df <- data.frame(y = reg$coefficients[[coef_name]],
                   lb = confint(reg, coef_name, level = 0.95)[1],
                   ub = confint(reg, coef_name, level = 0.95)[2])
  return(df)
}

# Populate the data frame with a list of your model results.

df.plot <- lapply(list(model1,      # List your models here
                       model2), 
                  gen_df_plot, 
                  coef_name = 'x1') # Coefficient name

  # Convert the list to a tidy data frame

df.plot <- data.table::rbindlist(df.plot)

# Provide the coefficient or regression labels below, in the
# order that you want them to appear. The "levels=unique(.)" parameter
# overrides R's desire to order the factor alphabetically

df.plot$x <- c("Group 1", 
               "Group 2") %>%
  factor(., levels = unique(.),
         ordered = TRUE)

# Create your plot

df.plot %>% ggplot(aes(x, y)) + 
  geom_point(size=4) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.1, linetype="dotted") + 
  geom_hline(aes(yintercept=0), linetype="dashed") + 
  theme_bw() +
  ggtitle("Comparing Coefficients") +
  ylab("Coefficient Value")```