Stata:双向固定效应模型中是否要控制公司年龄?

发布时间:2022-06-19 阅读 110

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc, ihelp, rdbalance, gitee, installpkg

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

PDF下载 - 推文合集

刘欣妍 (香港中文大学),liuxinyan@link.cuhk.edu.hk
史  柯 (中央财经大学),shike2231128@gmail.com

编者按:本文主要摘译自下文,特此致谢!
Source:Controlling for Log Age -Link-


目录


1. 引言

在实证分析中,我们通常需要控制某些变量来缓解遗漏变量问题,例如公司年龄 (age) 和公司规模 (size)。不过,当我们同时控制了公司固定效应和时间固定效应后,公司年龄与固定效应会存在共线性问题。为缓解上述问题,学者往往会对公司年龄取对数。那么这一方法是否可以发挥作用呢?

通过下文的分析,我们认为在实证研究中,需要将公司年龄的对数作为控制变量纳入到模型中。

2. R 中的实现

2.1 模拟数据

首先我们来模拟一个面板数据

# Install Packages
install.packages("magrittr") # package installations are only needed the first time you use it
install.packages("dplyr")    # alternative installation of the %>%
install.packages("fastDummies") 
library(magrittr) # needs to be run every time you start R and want to use %>%
library(dplyr) 
library(tidyr)
library(fastDummies)

# simulate data
# set seed 
set.seed(74792766)
# Generate data - 20,000 firms are placed in each group. Groups 3 and 4 are 
# treated in 1998, Groups 1 and 2 are untreated
make_data <- function(...) {
  # Fixed Effects
  # get a list of 4,000 units that are randomly treated sometime in the panel 
  treated_units <- sample(1:20000, 4000)
  # unit fixed effects
  unit <- tibble(
    unit = 1:20000, 
    unit_fe = rnorm(20000, 0, 1),
    treat = if_else(unit %in% treated_units, 1, 0)) %>% 
    # make first and last year per unit, and treat year if treated
    rowwise() %>% 
    mutate(first_year = sample(seq(1981, 2010), 1),
           # pull last year as a randomly selected date bw first and 2010
           last_year = if_else(first_year < 2010, sample(seq(first_year, 2010), 1), 
                               as.integer(2010)),
           # pull treat year as randomly selected year bw first and last if treated
           treat_year = if_else(treat == 1,
                                if_else(first_year != last_year,
                                        sample(first_year:last_year, 1), as.integer(first_year)),
                                as.integer(0))) %>% 
    ungroup()
  # year fixed effects 
  year <- tibble(
    year = 1981:2010,
    year_fe = rnorm(30, 0, 1))
  # make panel
  crossing(unit, year) %>% 
    arrange(unit, year) %>% 
    # keep only if year between first and last year 
    rowwise() %>% 
    filter(year %>% between(first_year, last_year)) %>% 
    ungroup() %>%
    # make error term, treat term and log age term
    mutate(error = rnorm(nrow(.), 0, 1),
           posttreat = if_else(treat == 1 & year >= treat_year, 1, 0),
           rel_year = if_else(treat == 1, year - treat_year, as.integer(NA)),
           tau = if_else(posttreat == 1, .2, 0),
           firm_age = year - first_year,
           log_age = log(firm_age + 1)) %>% 
    # make cumulative treatment effects
    group_by(unit) %>% 
    mutate(cumtau = cumsum(tau)) %>% 
    ungroup() %>% 
    # finally, make dummy variables
    bind_cols(dummy_cols(tibble(lag = .$rel_year), select_columns = "lag", 
                         ignore_na = TRUE, remove_selected_columns = TRUE)) %>%
    # replace equal to 0 for all lead lag columns
    mutate_at(vars(starts_with("lag_")), ~replace_na(., 0))
}

# make data
data <- make_data()

得到的模拟数据的情况如下:

2.2 估计模型

在得到模拟数据后,我们假设公司年龄变量的对数 (log_age,即 year-first_year 的对数) 与结果变量和处理变量都不相关。进而,我们可以通过估计一个简单的双向固定效应模型 (Two-way fixed effect),来画出真实的处理效应与估计的处理效应。具体模型如下所示:

#Install packages
install.packages("tidyverse")
install.packages("ggplot2")
install.packages("lfe")
library(ggplot2)
library(tidyverse) 
library(lfe)
#Draw plot
theme_set(theme_clean() + theme(plot.background = element_blank()))

# first make a plot where log age has no impact on the outcome variable
# get var names in a vector - need to drop the most negative lag (lag_1) and 
min <- min(data$rel_year, na.rm = TRUE) + 1
max <- max(data$rel_year, na.rm = TRUE)

# make string for rel time indicators 
indics <- paste0(c(paste0("`lag_", seq(min, -2), "`"),
                   paste0("lag_", seq(0, max))),
                 collapse = " + ")

# make the formula we want to run
form <- as.formula(paste0("y ~ ", indics, "| unit + year | 0 | unit"))

# get true taus to merge in
true_taus <- tibble(
  time = seq(-10, 10),
  true_tau = c(rep(0, length(-10:-1)), .2*(0:10 + 1))
)

# estimate the model and make the plot
data %>% 
  # generate the outcome variable
  mutate(y = unit_fe + year_fe + cumtau + error) %>% 
  # run the model and tidy up
  do(broom::tidy(felm(form, data = .), conf.int = TRUE)) %>%
  filter(term != "firm_age") %>% 
  # make the relative time variable and keep just what we need
  mutate(time = c(min:(-2), 0:max)) %>% 
  select(time, estimate, conf.low, conf.high) %>% 
  # bring back in missing indicator for t = -1
  bind_rows(tibble(
    time = -1, estimate = 0, conf.low = 0, conf.high = 0
  )) %>% 
  # keep just -10 to 10
  filter(time %>% between(-10, 10)) %>% 
  left_join(true_taus) %>% 
  # split the error bands by pre-post
  mutate(band_groups = case_when(
    time < -1 ~ "Pre",
    time >= 0 ~ "Post",
    time == -1 ~ ""
  )) %>%
  # plot it
  ggplot(aes(x = time, y = estimate)) + 
  geom_line(aes(x = time, y = true_tau, color = "True Effect"), size = 1.5, linetype = "dashed") + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, group = band_groups),
              color = "lightgrey", alpha = 1/4) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high, color = "Estimated Effect"), 
                  show.legend = FALSE) + 
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = -0.5, linetype = "dashed") + 
  scale_x_continuous(breaks = seq(-10, 10, by = 2)) + 
  labs(x = "Relative Time", y = "Estimate") +
  scale_color_brewer(palette = 'Set1') + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        axis.title = element_text(size = 18),
        axis.text = element_text(size = 16))

2.3 加入公司年龄的对数

当公司年龄对数与结果变量和处理变量均不相关时,加入公司年龄的对数并不会影响我们的估计结果。具体如下图所示:

theme_set(theme_clean() + theme(plot.background = element_blank()))

# make string for rel time indicators 
indics <- paste0(c(paste0("`lag_", seq(min, -2), "`"),
                   paste0("lag_", seq(0, max))),
                 collapse = " + ")

# make the formula we want to run
form <- as.formula(paste0("y ~ ", indics, " + log_age | unit + year | 0 | unit"))

# estimate the model and make the plot
data %>% 
  # generate the outcome variable
  mutate(y = unit_fe + year_fe + cumtau + error) %>% 
  # run the model and tidy up
  do(broom::tidy(felm(form, data = .), conf.int = TRUE)) %>%
  filter(term != "log_age") %>% 
  # make the relative time variable and keep just what we need
  mutate(time = c(min:(-2), 0:max)) %>% 
  select(time, estimate, conf.low, conf.high) %>% 
  # bring back in missing indicator for t = -1
  bind_rows(tibble(
    time = -1, estimate = 0, conf.low = 0, conf.high = 0
  )) %>% 
  # keep just -10 to 10
  filter(time %>% between(-10, 10)) %>% 
  left_join(true_taus) %>% 
  # split the error bands by pre-post
  mutate(band_groups = case_when(
    time < -1 ~ "Pre",
    time >= 0 ~ "Post",
    time == -1 ~ ""
  )) %>%
  # plot it
  ggplot(aes(x = time, y = estimate)) + 
  geom_line(aes(x = time, y = true_tau, color = "True Effect"), size = 1.5, linetype = "dashed") + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, group = band_groups),
              color = "lightgrey", alpha = 1/4) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high, color = "Estimated Effect"), 
                  show.legend = FALSE) + 
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = -0.5, linetype = "dashed") + 
  scale_x_continuous(breaks = seq(-10, 10, by = 2)) + 
  labs(x = "Relative Time", y = "Estimate") +
  scale_color_brewer(palette = 'Set1') + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        axis.title = element_text(size = 18),
        axis.text = element_text(size = 16))

当公司年龄对数与结果变量相关,而与处理变量不相关时,加入公司年龄的对数仍然不会影响我们的估计结果。需要注意的是,这里我们其实并不需要控制 log_age,因为它不是一个混杂变量 (只影响结果变量)。

theme_set(theme_clean() + theme(plot.background = element_blank()))

# next, assume that there is an effect on the outcome variable but no effect on the treatment outcome
# make string for rel time indicators 
indics <- paste0(c(paste0("`lag_", seq(min, -2), "`"),
                   paste0("lag_", seq(0, max))),
                 collapse = " + ")

form <- as.formula(paste0("y ~ ", indics, "  + log_age| unit + year | 0 | unit"))

# estimate the model and make the plot
data %>% 
  # generate the outcome variable
  mutate(y = unit_fe + year_fe + cumtau + -.85*log_age + error) %>% 
  # run the model and tidy up
  do(broom::tidy(felm(form, data = .), conf.int = TRUE)) %>% 
  # make the relative time variable and keep just what we need
  filter(term != "log_age") %>% 
  mutate(time = c(min:(-2), 0:max)) %>% 
  select(time, estimate, conf.low, conf.high) %>% 
  # bring back in missing indicator for t = -1
  # bring back in missing indicator for t = -1
  bind_rows(tibble(
    time = -1, estimate = 0, conf.low = 0, conf.high = 0
  )) %>% 
  # keep just -10 to 10
  filter(time %>% between(-10, 10)) %>% 
  left_join(true_taus) %>% 
  # split the error bands by pre-post
  mutate(band_groups = case_when(
    time < -1 ~ "Pre",
    time >= 0 ~ "Post",
    time == -1 ~ ""
  )) %>%
  # plot it
  ggplot(aes(x = time, y = estimate)) + 
  geom_line(aes(x = time, y = true_tau, color = "True Effect"), size = 1.5, linetype = "dashed") + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, group = band_groups),
              color = "lightgrey", alpha = 1/4) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high, color = "Estimated Effect"), 
                  show.legend = FALSE) + 
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = -0.5, linetype = "dashed") + 
  scale_x_continuous(breaks = seq(-10, 10, by = 2)) + 
  labs(x = "Relative Time", y = "Estimate") +
  scale_color_brewer(palette = 'Set1') + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        axis.title = element_text(size = 18),
        axis.text = element_text(size = 16)) 

当公司年龄对数与结果变量和处理变量都相关时,公司年龄对数就是一个混杂因素,此时是否需要将其纳入回归呢?我们通过以下示例可以看出,不加入公司年龄对数将得到一个有偏的估计。

theme_set(theme_clean() + theme(plot.background = element_blank()))

# Now assume that the log age variable is correlated in some way with the treatment assignment decision
# in particular assume that younger firms are more likely to get targeted 

# Generate data - 20,000 firms are placed in each group. Groups 3 and 4 are 
# treated in 1998, Groups 1 and 2 are untreated

make_data <- function(...) {
  # Fixed Effects ------------------------------------------------
  # unit fixed effects
  unit <- tibble(
    unit = 1:20000, 
    unit_fe = rnorm(20000, 0, 1)) %>%
    # make first and last year per unit, and treat year if treated
    rowwise() %>% 
    mutate(first_year = sample(seq(1981, 2010), 1),
           # pull last year as a randomly selected date bw first and 2010
           last_year = if_else(first_year < 2010, sample(seq(first_year, 2010), 1), 
                               as.integer(2010))) %>% 
    ungroup()
  
  # year fixed effects 
  year <- tibble(
    year = 1981:2010,
    year_fe = rnorm(30, 0, 1))
  
  # make panel
  data <- crossing(unit, year) %>% 
    arrange(unit, year) %>% 
    # keep only if year between first and last year 
    rowwise() %>% 
    filter(year %>% between(first_year, last_year)) %>% 
    ungroup() %>% 
    mutate(firm_age = year - first_year)
  
    # make an average age data frame 
  avg_age <- data %>% 
    group_by(unit) %>% 
    summarize(avg_age = mean(firm_age)) %>% 
    mutate(weight = 1 / (avg_age + 1))
  
  # sample 4,000 firms to get treatment, weighted by average age
  treated_units <- sample_n(avg_age, 4000, replace = FALSE, weight = weight) %>%
    mutate(treat = 1) %>% 
    select(unit, treat)
  
  # merge treat back into the data 
  treat_data <- data %>% 
    select(unit, first_year, last_year) %>% 
    distinct() %>% 
    left_join(treated_units) %>% 
    replace_na(list(treat = 0)) %>% 
    rowwise() %>% 
    mutate(treat_year = if_else(treat == 1,
                                if_else(first_year != last_year,
                                        sample(first_year:last_year, 1), as.integer(first_year)),
                                as.integer(0)))
  
  # merge back into main data
  data <- data %>% 
    left_join(treat_data) %>% 
    # make error term, treat term and log age term
    mutate(error = rnorm(nrow(.), 0, 1),
           posttreat = if_else(treat == 1 & year >= treat_year, 1, 0),
           rel_year = if_else(treat == 1, year - treat_year, as.integer(NA)),
           tau = if_else(posttreat == 1, .2, 0),
           firm_age = year - first_year,
           log_age = log(firm_age + 1)) %>% 
    # make cumulative treatment effects
    group_by(unit) %>% 
    mutate(cumtau = cumsum(tau)) %>% 
    ungroup() %>% 
    # finally, make dummy variables
    bind_cols(dummy_cols(tibble(lag = .$rel_year), select_columns = "lag", 
                         ignore_na = TRUE, remove_selected_columns = TRUE)) %>%
    # replace equal to 0 for all lead lag columns
    mutate_at(vars(starts_with("lag_")), ~replace_na(., 0))
}

# make data
data2 <- make_data()

# get var names in a vector - need to drop the most negative lag (lag_1) and 
min <- min(data2$rel_year, na.rm = TRUE) + 1
max <- max(data2$rel_year, na.rm = TRUE)

# make string for rel time indicators 
indics <- paste0(c(paste0("`lag_", seq(min, -2), "`"),
                   paste0("lag_", seq(0, max))),
                 collapse = " + ")

# make the formula we want to run
form1 <- as.formula(paste0("y ~ ", indics, "| unit + year | 0 | unit"))
form2 <- as.formula(paste0("y ~ ", indics, " + log_age| unit + year | 0 | unit"))

# estimate the model and make the plot
data2 %>% 
  # generate the outcome variable
  mutate(y = unit_fe + year_fe + cumtau + -.85*log_age + error) %>% 
  # run the model and tidy up
  do(bind_rows(
    broom::tidy(felm(form1, data = .), conf.int = TRUE) %>% mutate(mod = "No Control Log Age"),
    broom::tidy(felm(form2, data = .), conf.int = TRUE) %>% mutate(mod = "Control Log Age"))) %>% 
  # make the relative time variable and keep just what we need
  filter(term != "log_age") %>% 
  mutate(time = rep(c(min:(-2), 0:max), 2)) %>% 
  select(time, estimate, conf.low, conf.high, mod) %>% 
  # bring back in missing indicator for t = -1
  bind_rows(tibble(
    time = rep(-1, 2), estimate = rep(0, 2), conf.low = rep(0, 2),
    conf.high = rep(0, 2), mod = c("No Control Log Age", "Control Log Age")
  )) %>% 
  # keep just -10 to 10
  filter(time %>% between(-10, 10)) %>% 
  left_join(true_taus) %>% 
  # split the error bands by pre-post
  mutate(band_groups = case_when(
    time < -1 ~ "Pre",
    time >= 0 ~ "Post",
    time == -1 ~ ""
  )) %>%
  # plot it
  ggplot(aes(x = time, y = estimate)) + 
  geom_line(aes(x = time, y = true_tau, color = "True Effect"), size = 1.5, linetype = "dashed") + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, group = band_groups),
              color = "lightgrey", alpha = 1/4) + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high, color = "Estimated Effect"), 
                  show.legend = FALSE) + 
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = -0.5, linetype = "dashed") + 
  scale_x_continuous(breaks = seq(-10, 10, by = 2)) + 
  labs(x = "Relative Time", y = "Estimate") +
  scale_color_brewer(palette = 'Set1') + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        axis.title = element_text(size = 18),
        axis.text = element_text(size = 16),
        panel.spacing = unit(2, "lines")) + 
  facet_wrap(~mod)

3. Stata 中的实现

下面本文将通过 Stata 复现 R 中的操作,由于篇幅所限,本文仅复现第一个图像。为了保证结果的统一,本文使用 R 中生成的数据进行相关回归分析及图像绘制。具体代码如下:

* 回归分析
reghdfe y `x', absorb(unit year)

* 保存回归结果
matrix a2=r(table)
matrix bandse1=a2[1..2,1..9]
matlist bandse1
matrix bandse2=a2[1..2,24..34]
matlist bandse2
matrix  bandse0=[0,0]'    // 补充-1期
mat bandse = (bandse0,bandse1,bandse2)'
mat list bandse
clear
svmat bandse             // 将矩阵转化为变量

* 修改变量名
rename bandse1 estimated
rename bandse2 se

* 生成时间变量
gen reltime= (-1)*_n in 1/10
replace reltime=_n-11 in 11/21

* 生成 coef.min 及 coef.max (95%的置信区间) 
gen estmin=estimated-se*1.96
gen estmax=estimated+se*1.96

* 生成 real effect
gen realeffect= (reltime+1)*0.2 in 11/21
replace realeffect= 0 in 1/10

* plot it
sort reltime
graph twoway (rarea estmax estmin reltime, color(gs12))         ///
    (scatter estimated reltime, msize(small) mcolor(maroon))    ///
    (line realeffect reltime,lpattern(dash) lcolor(navy) lwidth(medium)), ///
    xtitle("Relative Time")  ytitle("Estimate") xticks(-10(2)10)          ///
    tline(-0.5, noextend lcolor(clack) lpattern(dash)) ///
    yline(0, noextend lcolor(clack))                   ///
    legend(label(1 "95% confidence interval")          ///
    label(2 "Estimated Effect") label(3 "True Effect"))  

4. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh 控制变量, m
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh