刘欣妍 (香港中文大学),liuxinyan@link.cuhk.edu.hk
史 柯 (中央财经大学),shike2231128@gmail.com
编者按:本文主要摘译自下文,特此致谢!
Source:Controlling for Log Age -Link-
温馨提示: 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码,直达原文:
目录
-
1. 引言
-
2. R 中的实现
-
2.1 模拟数据
-
2.2 估计模型
-
2.3 加入公司年龄的对数
-
3. Stata 中的实现
-
4. 相关推文
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
-
专题:Stata命令 -
Stata:控制变量组合的筛选-tuples -
Stata新命令-pdslasso:众多控制变量和工具变量如何挑选? -
专题:回归分析 -
调节效应是否需要考虑对控制变量交乘? -
控制变量!控制变量! -
不用太关心控制变量,真的! -
加入控制变量后结果悲催了! -
专题:IV-GMM -
Lasso一下:再多的控制变量和工具变量我也不怕-T217 -
专题:断点回归RDD -
RDD:断点回归可以加入控制变量吗? -
Stata:RDD-中可以加入控制变量 -
专题:其它 -
锚定情境法(一):有效控制变量自评偏差
#artContent h1{font-size:16px;font-weight: 400;}#artContent p img{float:none !important;}#artContent table{width:100% !important;}