16 固定效应
16.1 方法原理
16.1.1 控制未观测混杂变量
本书大量篇幅讨论了通过控制变量来解答研究问题的方法。若已构建因果图，我们可据此确定需要控制的变量，并参照\(第13章和第14章\)的方法进行控制。
然而我们始终面临两大难题：其一，若无法构建正确的因果图来确定需要控制的变量怎么办？其二，若无力测量所有需控制的变量又当如何？这两种情况下，仅靠控制变量法均无法识别真实因果效应。
不过固定效应模型真的能实现这一目标吗？该方法的本质在于：只要变量在特定层级类别内保持恒定，无论其是否可观测，都能通过控制该类别变量来间接控制所有组内不随时间变化的混杂因素。其实现原理非常简单——只需在模型中加入分组虚拟变量或进行组内离差变换，即可自动控制所有组内恒定的特征12。
“变量在更高层级类别中保持恒定”具体是指该变量在给定的分组维度内不随时间或个体变化。以评估农村地区通电对生产率的影响为例：
地理因素构成了一个典型的混杂变量——位于山区的乡镇（更高层级类别为”县”）往往具有以下特征：
通电难度更大（因地形复杂导致基建成本升高）
生产率基础值不同（受限于土壤贫瘠、交通不便等地理约束）
这种层级恒定性的计量特征是：
在县级层面：山地/平原的地理属性是固定不变的
在乡镇层面：通电状态与生产率可能随时间变化
形成混杂的机制：县级地理特征同时影响通电进度（处理变量）和农业生产效率（结果变量）
但地理特征本身具有高度时间稳定性：对同一城镇进行多年观测时，其地理属性基本保持不变。那么，如果在许多不同的年份观察同一个城镇，我们控制城镇的影响会怎么样呢3？如果我们这样做，我们将消除城镇所解释的任何变异。而且，由于地理位置的任何变异都由城镇解释，因此在控制城镇的影响之后，我们就控制了地理位置的影响——地理位置将不再有任何变异！
在这种情况下，我们或许本可以对地理因素进行测量和控制，这没什么大不了的。但对于那些我们无法控制的因素又该怎么办呢？
举一个非常典型的例子，在许多社会科学的研究情境中，我们希望对“个人背景”进行控制，但这是一个高维的概念，包含了很多我们根本无法衡量的方面。它属于不可观测因素。但如果加入个体固定效应呢？哇塞！这个问题就得到控制了4。
固定效应可以被看作是，把某一类别内保持不变的、存在于后门路径上的一长串变量整合起来，将它们归为这一个类别，然后对该类别进行控制。
举个例子，假设我们想了解德国总理访问某国对该国与德国贸易水平的影响，如\(图16.2\)所示。
在这幅图中，与德国开展贸易有诸多影响因素，其中很多因素难以追踪，甚至难以衡量。然而，我们也可以注意到，这些变量中有几个——该国的地理位置、民众的文化以及该国与德国的历史渊源——在一个国家内部是恒定不变的，或者至少（比如民众的文化）在我们所掌握的任何数据涵盖的时间段内不太可能发生太大变化5 。因此，我们可以像\(图16.2\)那样重新绘制这幅图。
现在，有了我们简化后的图表，我们无需对\(图16.2\)中的每个变量进行控制就能确定影响。我们只需对国家进行控制即可6。
请记住，我们并没有将所有因素都归结到国家层面。任何可能随时间规律变化的因素，比如一个国家当前的政治状况，都无法通过固定效应来解决。所以，如果我们想识别出这种效应，除了对国家设置固定效应之外，我们仍然需要控制“当前政治状况”这一因素。
16.1.2 组内变异
从本质上讲，固定效应是对个体进行控制，这里的“个体”在你的研究情境中可以指”个人”“公司”“学校”或者”国家”等等。
个体内变异\(（Within \space variation）\)：同一个体在不同时间段（通常）内某变量的变化。
个体间变异\(（Between \space variation）\)：不同个体之间某变量的差异，通常在同一时间段内比较，或基于跨时间平均值的对比。
这意味着它消除了个体 之间 的任何差异。我们之所以知道这一点，是因为这就是”控制个体因素”的含义。我们正在去除数据中由个体因素所解释的所有差异（即个体之间的所有差异）。
剩下的是个体内部的差异。什么是个体内部的差异呢？例如，假设你是一个相当高的人，我们在跟踪你随时间变化的身高情况。你比其他人高这一事实属于个体间差异。将你与其他人比较，总体而言你更高。然而，我们也能观察到，当你从儿童成长为成年人时，你的身高随时间发生了变化。将现在的你与更矮、更年幼的自己进行比较，这就是个体内部差异。个体内部发生了某些变化，而这正是我们所关注的。
固定效应的核心思想在于：通过消除所有个体间的变异，我们实际上控制了所有不随时间变化的变量7。例如，无论我一生中去往何处，我的出生城市始终不变。同理适用于所有人8。因此，“出生城市”的所有变异都属于个体间变异。如果我们消除所有个体间变异，就等于消除了出生城市的所有变异——从而实现了对该变量的控制。
需要特别注意的是：固定效应模型确实完全消除了个体间变异。因此，如果我们真正感兴趣的是个体间变异（比如想研究”出生城市”的影响效应），那么固定效应模型就无能为力了9。
library(tidyverse) ## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(vtable) ## Loading required package: kableExtra
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(cowplot) ##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamplibrary(Cairo)
library(extrafont) ## Registering fonts with Rlibrary(gapminder)
library(lfe) ## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(Cairo)
library(ggpubr) ## Registered S3 method overwritten by 'broom':
## method from
## nobs.felm lfe
##
## Attaching package: 'ggpubr'
##
## The following object is masked from 'package:cowplot':
##
## get_legendlibrary(modelsummary)
library(splines)
library(stargazer)##
## Please cite as:
##
## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer让我们通过两个技术示例来说明这一点：一个用表格展示，另一个用图表呈现。假设我们想研究运动对每年感冒次数的影响，研究对象是你和我。我们推测存在某些不随时间变化的变量（如”遗传因素”），可以通过固定效应模型来控制。现在，先来看看我们收集的运动数据：
tb <- tibble(Individual = c('You','You','Me','Me'), Year = c('2019','2020','2019','2020'), Exercise = c(5,7,4,3))
tb## # A tibble: 4 × 3
## Individual Year Exercise
## <chr> <chr> <dbl>
## 1 You 2019 5
## 2 You 2020 7
## 3 Me 2019 4
## 4 Me 2020 3\(表16.1\)展示了我们的原始数据（运动量以每周小时数衡量）。要分析个体内变异和个体间变异，我们需要先计算每个个体的运动量均值，如\(表16.2\)所示。
\(表16.2\)展示了我们如何计算每个个体的平均运动时长，并减去该均值。个体间均值的差异即为个体间变异。因此，运动的个体间变异来源于比较你每周\(6\)小时的平均运动量与我每周\(3.5\)小时的平均运动量之间的差异。
# Table example of between/within
tb %>% group_by(Individual) %>%
mutate(MeanExercise = mean(Exercise)) %>%
mutate(WithinExercise = Exercise - MeanExercise) ## # A tibble: 4 × 5
## # Groups: Individual [2]
## Individual Year Exercise MeanExercise WithinExercise
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 You 2019 5 6 -1
## 2 You 2020 7 6 1
## 3 Me 2019 4 3.5 0.5
## 4 Me 2020 3 3.5 -0.5若从原始数据中减去这些个体均值，剩下的就是各时间段内个体相对于自身平均值的波动情况——这就是个体内变异，反映的是个体内部的变化特征。以我的运动量为例：\(2019年\)我比平均运动量每周多锻炼半小时，而\(2020年\)则比平均值每周少锻炼半小时，这两个偏离均值的差异就构成了我的个体内变异。
我们也可以通过图形化方式直观呈现这一点。\(图16.3\)左上角的子图展示了当我们增加多年数据采集后，数据可能呈现的分布形态。
\[Figure \space 16.3: \space Between \space and \space Within \space Variation \space of \space Exercise \space and \space Colds \space for \space You \space and \space Me\]
# Graphically
set.seed(2000)
tb <- tibble(Individual = c(rep('You',8),rep('Me',8))) %>%
mutate(ExercisePerWeek = 4*runif(16) + 3*(Individual == "You")) %>%
mutate(ColdsPerYear = 5+ 5*runif(16) + 5*(Individual == 'You') - ExercisePerWeek) %>%
group_by(Individual) %>%
mutate(ExMean = mean(ExercisePerWeek), ColdsMean = mean(ColdsPerYear))
p1 <- ggplot(tb, aes(x = ExercisePerWeek,y=ColdsPerYear, shape = Individual, color = Individual)) +
geom_point() +
scale_color_manual(values = c('Me' = 'black', 'You' = '#777777')) +
theme_pubr() +
theme(legend.position.inside = c(.85,.8),
legend.background = element_rect(color = 'black'),
text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family= "sans")) +
labs(x = "Exercise Hours Per Week", y = "Colds Per Year", title = "Raw Data") +
scale_y_continuous(limits = c(3,11)) + scale_x_continuous(limits = c(0,8))
p2 <- ggplot(tb, aes(x = ExercisePerWeek,y=ColdsPerYear, shape = Individual, color = Individual)) +
geom_point() +
geom_point(aes(x = ExMean, y = ColdsMean, color = Individual),shape = 3, size = 40) +
scale_color_manual(values = c('Me' = 'black', 'You' = '#777777')) +
theme_pubr() +
guides(shape = FALSE, color = 'none') +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family= "sans")) +
labs(x = "Exercise Hours Per Week", y = "Colds Per Year", title = "Raw Data with Individual Means")+
scale_y_continuous(limits = c(3,11)) +
scale_x_continuous(limits = c(0,8)) ## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.p3 <- ggplot(tb, aes(x = ExMean,y=ColdsMean, shape = Individual, color = Individual)) +
geom_point() + geom_point(aes(x = ExMean, y = ColdsMean, color = Individual),shape = 3, size = 40) +
scale_color_manual(values = c('Me' = 'black', 'You' = '#777777')) +
theme_pubr() +
guides(shape = FALSE, color = 'none') +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family= "sans")) +
labs(x = "Exercise Hours Per Week", y = "Colds Per Year", title = "Just Between Variation")+
scale_y_continuous(limits = c(3,11)) +
scale_x_continuous(limits = c(0,8))
p4 <- ggplot(tb %>% dplyr::filter(Individual == 'Me'),
aes(x = ExercisePerWeek - ExMean, y=ColdsPerYear - ColdsMean, shape = Individual, color = Individual)) +
geom_point() +
geom_hline(aes(yintercept = 0)) +
geom_vline(aes(xintercept = 0)) +
scale_color_manual(values = c('Me' = 'black', 'You' = '#777777')) +
theme_pubr() +
guides(shape = FALSE, color = 'none') +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family= "sans")) +
labs(x = "Exercise Hours Per Week", y = "Colds Per Year", title = "Just Within Variation Just for Me")
plot_grid(p1,p2,p3,p4, nrow = 2) 
ggsave('FixedEffects/between_within_graphs.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)\(图16.3\)右上角的子图在原始数据基础上叠加了个体均值基准线——相当于为每个研究对象建立了独立坐标系。具体操作是：先计算你运动时长和感冒次数的均值，将该均值点作为你个人坐标系的\(原点（0，0）\)；同理以我的均值点建立我的坐标系。当然，在原始共享坐标系中，这些点并非真正的\(（0，0）\)，而是我们在\(X\)轴（运动量）和\(Y\)轴（感冒次数）上的个体均值位置。
\(图16.3\)左下角的子图则聚焦于这些坐标系本身——我们将每个坐标系的中心点标出，并仅保留这些个体均值数据。此时若通过连线方式比较不同坐标系中心点（例如将你的均值点与我的均值点相连），我们观察到的就是个体间变异。 在\(图16.3\)右下角的子图中，我们暂时忽略你的数据，仅聚焦于我的坐标系10。此时将我的均值点设为新原点\(（0，0）\)，所有数据点都将基于我的平均值重新定位。以左下角的数据点为例：在右上子图中，该点坐标约为每周运动\(0.3\)小时/每年感冒\(5.5\)次；而在当前坐标系下，该点表示每周运动量比我的平均值少\(1.3\)小时，每年感冒次数比我的平均值少\(0.8\)次。
要实现固定效应模型，我们需要将你的坐标系整体平移，与我的坐标系完全重合！下面具体说明这个操作过程
16.1.3 固定效应的实际应用
固定效应模型究竟如何处理数据？其本质原理与\(第13章\)所述的控制分类变量方法相同——核心在于提取个体内变异。具体而言：首先计算组别（个体）间的均值差异，然后将这些差异从数据中剔除。这就是固定效应模型的根本原理。
我们首先通过\(图16.4\)展示原始数据，需要特别注意的是数据中包含多个个体。本研究主要关注健康饮食提醒对实际饮食行为的影响效应。\(4\)名受试者各自下载了一款应用程序，该程序会随机发送健康饮食提醒。虽然每位受试者可自主设置提醒频率，但无法控制具体接收提醒的日期。我们记录了每位受试者接收提醒的频率，并对其实际饮食健康程度进行了评分。
因此，我们的数据呈现出两个关键维度的变化：首先是每位个体接收健康提醒的密集程度随时间的变化，其次是他们实际饮食健康水平的变化。值得注意的是，可能存在个体层面的混杂因素（即统计中的”后门路径”）——这些因素既影响了个体对提醒频率的选择，也与其基础健康饮食水平相关。正因如此，我们需要运用固定效应模型来有效阻断这些潜在的混杂路径。
\(图16.4\)显示，健康提醒与实际饮食健康程度之间存在微弱正相关性\(（r=0.111）\)——接收更多提醒似乎意味着饮食更健康！虽然效应量很小，但这可能暗示提醒消息确实有效11。不过，由于我们已意识到存在个体层面的混杂因素（后门路径），这一结果尚不能作为因果结论。因此，我们启动固定效应估计：首先分别计算每个个体的平均提醒接收频率和平均饮食健康评分。
set.seed(2000)
tb <- tibble(Individual = factor(c(rep('You',8),rep('Me',8), rep('Shamma',8),rep('Liqing',8)),
levels = c('Me','You','Liqing','Shamma')),
IndNo = sort(rep(1:4,8))) %>%
mutate(IntensityOfReminders = runif(32)*5 + IndNo) %>%
# the relationships of HealthyEatingScore and IntensityOfReminders
mutate(HealthyEatingScore = runif(32)*10 + IntensityOfReminders - 2*IndNo) %>%
mutate(HealthyEatingScore = case_when(HealthyEatingScore < 0 ~ 0, TRUE ~ HealthyEatingScore)) %>%
group_by(IndNo) %>%
mutate(MeanIntensity = mean(IntensityOfReminders), MeanScore = mean(HealthyEatingScore)) %>%
# demean
mutate(Intensity.R = IntensityOfReminders - MeanIntensity, Score.R = HealthyEatingScore - MeanScore)
cor(tb$IntensityOfReminders,tb$HealthyEatingScore)## [1] 0.1117119tb## # A tibble: 32 × 8
## # Groups: IndNo [4]
## Individual IndNo IntensityOfReminders HealthyEatingScore MeanIntensity
## <fct> <int> <dbl> <dbl> <dbl>
## 1 You 1 1.98 1.20 3.87
## 2 You 1 4.58 11.7 3.87
## 3 You 1 2.81 6.57 3.87
## 4 You 1 2.96 7.26 3.87
## 5 You 1 5.07 4.28 3.87
## 6 You 1 3.14 10.9 3.87
## 7 You 1 5.80 7.98 3.87
## 8 You 1 4.63 5.92 3.87
## 9 Me 2 6.10 8.86 4.01
## 10 Me 2 5.06 1.76 4.01
## # ℹ 22 more rows
## # ℹ 3 more variables: MeanScore <dbl>, Intensity.R <dbl>, Score.R <dbl>
# Walk through FE
ggplot(tb, aes(x = IntensityOfReminders, y = HealthyEatingScore, shape = Individual)) +
geom_point() +
theme_pubr() +
labs(x = "Intensity of Reminders", y = "Healthy Eating Score")+
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family= "sans")) 
ggsave('FixedEffects/fe_action_raw.pdf', width = 6, height = 5,device=cairo_pdf) \(图16.5\)展示了各受试者的均值分布情况。从个体均值来看，健康提醒与饮食健康程度呈现明显的负相关关系——若仅观察图中的大十字标记（代表个体均值），可见其呈现下降趋势。这种组间关系为负向：选择接收更频繁提醒消息的个体，其平均饮食健康水平反而更低。组间相关系数达到\(-0.44\)。
cor(tb$MeanIntensity,tb$MeanScore)## [1] -0.4395609
ggplot(tb, aes(x = IntensityOfReminders, y = HealthyEatingScore, shape = Individual)) +
geom_point() +
geom_point(aes(x = MeanIntensity, y = MeanScore), shape = 3, size = 40) +
theme_pubr() +
labs(x = "Intensity of Reminders", y = "Healthy Eating Score") +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family= "sans")) 
ggsave('FixedEffects/fe_action_means.pdf', width = 6, height = 5, device = cairo_pdf) 我们接下来的任务是消除所有个体间变异。这意味着，在统计操作上，我们将把四个”+“标记（代表个体均值）完全重叠。首先，让我们通过\(图16.6\)单独观察四位受试者的个体内变异情况
\[Figure \space 16.6: \space Healthy \space Eating \space for \space Four \space Individuals\]
shapes <- c(16,17,15,3)
p <- 1:4 %>%
map(function(x)
ggplot(tb %>% dplyr::filter(Individual == levels(tb$Individual)[x]),aes(x = Intensity.R, y = Score.R)) +
geom_point(shape = shapes[x]) +
geom_vline(aes(xintercept = 0)) +
geom_hline(aes(yintercept = 0)) +
theme_pubr() +
labs(x = 'Within Variation in Intensity',
y = 'Within Variation in Score',
title = levels(tb$Individual)[x]) +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family= "sans"))
)
plot_grid(plotlist = p, ncol = 2) 
ggsave('FixedEffects/fe_action_individual.pdf',width = 7, height = 6, units = 'in', device = cairo_pdf)请注意，在\(图16.6\)的四个子图中，所有数据点的\(X\)轴（提醒频率）和\(Y\)轴（饮食健康评分）都以 \(0\) 值为中心进行了标准化处理。因此，当我们将这四个子图叠加时，所有数据点将集中在图形的同一区域分布。
cor(tb$Intensity.R,tb$Score.R)## [1] 0.3632897如\(图16.7\)所示，我们将所有数据叠加呈现——通过将所有”+“标记（个体均值点）精确对齐，实现坐标系统一。
ggplot(tb, aes(x = Intensity.R, y = Score.R, shape = Individual)) +
geom_point() + geom_vline(aes(xintercept = 0)) +
geom_hline(aes(yintercept = 0)) + theme_pubr() +
labs(x = "Within Variation in Intensity", y = "Within Variation in Score") +
theme(text = element_text(size = 10, family="sans"),
axis.title.x = element_text(size = 10, family="sans"),
axis.title.y = element_text(size = 10, family= "sans")) 
ggsave('FixedEffects/fe_action_after.pdf', width = 6, height = 5, device = cairo_pdf) 由于我亲自设定了这个数据生成过程，因此可以确认：健康提醒与饮食改善之间确实存在正向因果关系。
16.2 固定效应模型的实现方法
16.2.1 回归估计量
固定效应回归通常始于如下回归方程：\(Y_{it} = \beta_i + \beta_1 X_{it} + \varepsilon_{it} \quad\quad \text{(16.1)}\)
该方程与\(第13章\)的标准回归方程形式相似，但存在两个关键区别：
双下标变量：\(X_{it}\)中的下标”\(it\)“表明，该数据同时具备个体间差异 \(i\) 和时序变化 \(t\)
个体截距项：\(\beta_i\)表示每个个体独有的基线效应12
第二个区别在于：截距项 \(\beta\) 的下标从 \(0\) 变为 \(i\)，记作 \(\beta_i\)。将回归理解为直线拟合时，这意味着：
统一斜率约束：所有个体必须共享相同斜率（\(\beta_1\) 无 \(i\) 下标）
可变截距特性：每个个体可拥有独特的截距值
由此自然引出三个关键问题：
\(问题1\)：允许截距项变化为何能获得仅使用个体内变异的固定效应估计？
这个问题可以从两个角度直观理解。首先是操作层面：允许截距变化相当于为每个研究对象单独添加控制变量。以国家研究为例，如果样本包含印度，那么模型中就会存在一个专属印度的截距项 \(\beta_{India}\)——这个截距仅适用于印度数据，不参与其他国家的计算。
这非常类似于在回归模型中添加一个”是否为印度”的二值控制变量，而 \(\beta_{India}\) 就是这个控制变量的系数。我们的数据中每个国家都对应这样一个控制变量——本质上，我们是在为每个国家进行单独控制。正如计量经济学基本原理所示：当我们控制国家变量时，所得的就是国家内部的变异。
理解这一机制的第二种直观方式是通过图形化分析。\(图16.8\) 采用了 \(Gapminder\) 研究所提供的 \(人均GDP\) 与 \(预期寿命\) 数据——我们特意选取了印度和巴西这两个国家，展示它们自\(1952年\)以来的数据轨迹。
Gapminder Institute. Gapminder. https://www.gapminder.org/,2020. Accessed: 2020-03-09.
library(ggplot2)
library(ggpubr)
# 优化后的绘图代码
ggplot(gapminder, aes(x = gdpPercap, y = lifeExp, group = country, color = country)) +
# 1. 基础几何对象
geom_point(alpha = 0.6) + # 添加透明度提高可读性
# 2. 文本标签优化（按需显示）
ggrepel::geom_text_repel(aes(label = ifelse(year %% 1 == 0, year, "")),size = 3,show.legend = FALSE,max.overlaps = 20) +
# 3. 比例尺调整
scale_x_log10(breaks = c(500, 1000, 5000, 10000, 50000),labels = scales::dollar) +
# 4. 颜色定制
scale_color_manual(values = c('Brazil' = '#2ca02c', 'India' = '#d62728'), guide = "none") +
# 5. 注释优化
geom_label(
data = data.frame(x = c(1250, 7000),y = c(65, 57),label = c("India", "Brazil"),color = c('#d62728', '#2ca02c')),
aes(x = x, y = y, label = label),
inherit.aes = FALSE,
size = 3.5
) +
# 6. 趋势线
geom_smooth(method = 'lm',formula = y ~ x,color = 'red', linetype = 'dashed',se = FALSE,aes(group = 1)) +
# 7. 主题与标签
labs(x = 'GDP per capita (PPP-adjusted, log scale)',y = 'Life expectancy (years)',caption = "Source: Gapminder dataset") +
theme_pubr() +
theme(
text = element_text(size = 10, family = "sans"),
axis.title = element_text(size = 10),
panel.grid.major = element_line(color = "grey90", linewidth = 0.2)
)## Warning: ggrepel: 1626 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
ggsave('FixedEffects/fe_intercept_shared.pdf', width = 6, height = 5, device = cairo_pdf)## Warning: ggrepel: 1641 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps# Varying Intercept
data(gapminder)
gapminder <- gapminder %>% dplyr::filter(country %in% c('India','Brazil'), year > 1940)
# Additional data that spans the box
singleline <- coef(lm(lifeExp ~ log(gdpPercap), data=gapminder))
doubleline <- coef(lm(lifeExp ~ -1 + factor(country) + log(gdpPercap), data=gapminder))
singleline## (Intercept) log(gdpPercap)
## -0.8176992 7.5888733doubleline## factor(country)Brazil factor(country)India log(gdpPercap)
## -69.31981 -52.10183 15.35575从\(图16.8\)中可以明确得出三点关键发现：首先，两个国家的各项指标均呈现出随时间推移而显著改善的趋势。然而，图中也显示出明显的国家间差异——这些正是我们需要通过模型加以控制的系统性差异。值得注意的是，当前采用单一截距项的回归线明显低估了巴西和印度各自增长轨迹的真实斜率，这种设定未能捕捉两国在发展路径上的异质性特征。
通过分解截距项，我们可以改用两条斜率相同但截距不同的平行回归线，通过上下平移使其分别拟合两国数据。这种方法实现了双重优势：在纵轴维度上，我们无需关注各国数据点的整体高低分布（通过截距调整即可精准捕捉各国\(Y\)变量的组内变异）；在横轴维度上，由于斜率固定，无论数据点沿\(X\)轴如何分布（即各国经济水平差异），平移后的平行线均能有效拟合（从而提取\(X\)变量的组内变异）。这正是固定效应模型提取组内变异的数学本质——通过截距的灵活性剥离组间差异，仅保留核心解释变量的纯净时序影响。
library(ggplot2)
library(ggpubr)
library(scales)##
## Attaching package: 'scales'## The following object is masked from 'package:purrr':
##
## discard## The following object is masked from 'package:readr':
##
## col_factorlibrary(ggrepel)
# 优化后的可视化代码
ggplot(gapminder %>% dplyr::filter(country %in% c("Brazil", "India")), aes(x = gdpPercap, y = lifeExp, color = country)) +
# 1. 数据点层（增加透明度避免重叠）
geom_point(size=1, alpha = 0.5) +
# 2. 智能标签
ggrepel::geom_text_repel(aes(label = year),size = 3.0,show.legend = FALSE) +
# 3. 回归线（替代stat_function的更精确方法）
geom_smooth(method = "lm",formula = y ~ log(x),se = FALSE,linewidth = 0.5) +
# 4. 坐标轴优化
scale_x_log10(breaks = c(500, 1000, 2500, 5000, 10000, 25000),
labels = dollar_format(accuracy = 1),
name = "GDP per Capita (log scale, PPP-adjusted)") +
scale_y_continuous(name = "Life Expectancy (years)",breaks = seq(40, 80, by= 5)) +
# 5. 颜色和图例设置
scale_color_manual(values = c('Brazil' = '#1F77B4', 'India' = '#FF7F0E'), name = NULL) +
# 6. 注释优化
geom_label(
data = data.frame(x = c(1250, 7000),y = c(65, 57),label = c("India", "Brazil"),color = c('#FF7F0E', '#1F77B4')),
aes(x = x, y = y, label = label),
inherit.aes = FALSE,
size = 3,
label.padding = unit(0.4, "lines")
) +
# 7. 主题和标签
labs(title = "GDP and Life Expectancy Relationship (1960-2015)",
subtitle = "Country-specific trajectories with common slope",
caption = "Data source: Gapminder Project") +
theme_pubr(base_size = 10) +
theme(
text = element_text(family = "sans"),
legend.position = "top",
plot.title = element_text(face = "bold", size = 10),
plot.subtitle = element_text(size = 10),
plot.caption = element_text(size = 10, color = "gray50"),
panel.grid.major.y = element_line(color = "grey90", size = 0.2)
) ## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggsave('FixedEffects/fe_intercept_separate.pdf', width = 6, height = 5, device = cairo_pdf)我们在\(图16.9\)中精确实现了这一设定——保持斜率一致但采用差异化截距。如图所示，调整后的回归线显著贴近实际数据点，相较于单一回归线，这些平行线能更准确地捕捉两国数据中真实存在的关联性。
\(问题2\)：如何估计包含个体截距项的回归模型？
实际上，实现这一目标的方法有多种，其中部分将在本章后续章节详细探讨13。在这些方法中，有两种标准方法特别值得关注——它们不仅极其简单，而且仅需使用我们已经掌握的回归技术即可实现。
第一种方法是通过手动提取个体内变异来实现的14。具体步骤如下：
对每个个体计算因变量的时间序列均值\(\bar{Y}_i\)
构造去中心化的因变量：\(Y_{it}-\bar{Y}_i\)
对每个自变量重复相同操作：\(X_{it} - \bar{X}_i\)
最后估计以下回归模型：\(Y_{it}-\bar{Y}_i=\beta_0 + \beta_1 (X_{it}-\bar{X}_i) + \epsilon_{it} \quad\quad (16.2)\)
（注：原方程中的\(\beta_0\)项在去均值处理后理论上应为零，实际操作中可省略）
这种方法在计量经济学中被称为”固定效应吸收法15“。我们将再次运用 \(Gapminder\) 数据集来演示该方法的实际操作。
library(tidyverse)
library(modelsummary)
gm <- causaldata::gapminder
gm <- gm %>%
# Put GDP per capita in log format since it's very skewed
mutate(log_GDPperCap = log(gdpPercap)) %>%
# Perform each calculation by group
group_by(country) %>%
# Get within variation by subtracting out the mean
mutate(lifeExp_within = lifeExp - mean(lifeExp),log_GDPperCap_within = log_GDPperCap - mean(log_GDPperCap)) %>%
# We no longer need the grouping
ungroup()
# Analyze the within variation
m1 <- lm(lifeExp_within ~ log_GDPperCap_within, data = gm)
msummary(m1, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | -0.000 |
| (0.123) | |
| log_GDPperCap_within | 9.769*** |
| (0.284) | |
| Num.Obs. | 1704 |
| R2 | 0.410 |
| R2 Adj. | 0.410 |
| AIC | 10367.6 |
| BIC | 10383.9 |
| Log.Lik. | -5180.806 |
| F | 1183.326 |
| RMSE | 5.06 |
import numpy as np
import statsmodels.formula.api as sm
from causaldata import gapminder
gm = gapminder.load_pandas().data
# Put GDP per capita in log format since it's very skewed
gm['logGDPpercap'] = gm['gdpPercap'].apply('log')
# Use groupby to perform calculations by group
# Then use transform to subtract each variable's within -group mean to get within variation
gm[['logGDPpercap_within','lifeExp_within']] = (gm.groupby('country')[['logGDPpercap','lifeExp']].transform(lambda x: x - np.mean(x)))
# Analyze the within variation
m1 = sm.ols(formula = 'lifeExp_within ~ logGDPpercap_within',data = gm).fit()
m1.summary()| Dep. Variable: | lifeExp_within | R-squared: | 0.410 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.410 |
| Method: | Least Squares | F-statistic: | 1183. |
| Date: | 周六, 23 8月 2025 | Prob (F-statistic): | 2.51e-197 |
| Time: | 23:50:23 | Log-Likelihood: | -5180.8 |
| No. Observations: | 1704 | AIC: | 1.037e+04 |
| Df Residuals: | 1702 | BIC: | 1.038e+04 |
| Df Model: | 1 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 7.216e-16 | 0.123 | 5.88e-15 | 1.000 | -0.241 | 0.241 |
| logGDPpercap_within | 9.7690 | 0.284 | 34.400 | 0.000 | 9.212 | 10.326 |
| Omnibus: | 23.105 | Durbin-Watson: | 0.626 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 33.240 |
| Skew: | -0.143 | Prob(JB): | 6.05e-08 |
| Kurtosis: | 3.622 | Cond. No. | 2.32 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
第二种简单的方法是为每个国家（个体）添加一个二元控制变量。更准确地说，是除了一个国家之外的每个国家——如果我们的模型中还有一个截距/常数项，那么包含所有国家将会使模型无法估计，所以我们要剔除一个国家，或者剔除截距项16。那个被剔除的国家仍然在分析中；它只是没有得到自己的系数，而是被卡在了那个糟糕的常数项里。谢天谢地，大多数软件都会自动处理这种情况。
需要注意的是，这种做法通常并非理想选择。当模型中仅包含少量分类变量时，引入一组二元指标是可行的，但若需为每个个体设置控制变量，则可能涉及数百甚至数千个参数项。这将导致计算效率急剧下降——说真的，难道您还打算逐一解释所有这些系数吗？采用这种方法的常见理由之一，是通过联合显著性检验来验证是否存在截距项的异质性证据（参见\(第13章\)关于”联合\(F\)检验”的论述）。尽管如此，这绝非估计固定效应的常规方法。不过出于演示目的，我们将运行这个版本来观察结果。
library(tidyverse)
library(modelsummary)
gm <- causaldata::gapminder
# Simply include a factor variable in the model to get it turned to binary variables. You can use factor() to ensure it's done.
m2 <- lm(lifeExp ~ factor(country) + log(gdpPercap), data = gm)
msummary(m2, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| (Intercept) | -27.773*** |
| (2.501) | |
| factor(country)Albania | 17.783*** |
| (2.195) | |
| factor(country)Algeria | 5.241** |
| (2.214) | |
| factor(country)Angola | -13.907*** |
| (2.202) | |
| factor(country)Argentina | 8.132*** |
| (2.273) | |
| factor(country)Australia | 6.400*** |
| (2.352) | |
| factor(country)Austria | 5.156** |
| (2.348) | |
| factor(country)Bahrain | -1.968 |
| (2.344) | |
| factor(country)Bangladesh | 12.401*** |
| (2.158) | |
| factor(country)Belgium | 5.617** |
| (2.349) | |
| factor(country)Benin | 7.745*** |
| (2.161) | |
| factor(country)Bolivia | 2.338 |
| (2.193) | |
| factor(country)Bosnia and Herzegovina | 17.353*** |
| (2.194) | |
| factor(country)Botswana | 3.152 |
| (2.200) | |
| factor(country)Brazil | 6.318*** |
| (2.230) | |
| factor(country)Bulgaria | 12.670*** |
| (2.239) | |
| factor(country)Burkina Faso | 6.858*** |
| (2.158) | |
| factor(country)Burundi | 12.646*** |
| (2.164) | |
| factor(country)Cambodia | 12.917*** |
| (2.160) | |
| factor(country)Cameroon | 3.062 |
| (2.171) | |
| factor(country)Canada | 5.465** |
| (2.366) | |
| factor(country)Central African Republic | 4.758** |
| (2.159) | |
| factor(country)Chad | 5.737*** |
| (2.161) | |
| factor(country)Chile | 9.915*** |
| (2.242) | |
| factor(country)China | 21.271*** |
| (2.160) | |
| factor(country)Colombia | 10.920*** |
| (2.209) | |
| factor(country)Comoros | 10.205*** |
| (2.163) | |
| factor(country)Congo, Dem. Rep. | 10.025*** |
| (2.160) | |
| factor(country)Congo, Rep. | 1.375 |
| (2.198) | |
| factor(country)Costa Rica | 14.552*** |
| (2.228) | |
| factor(country)Cote d'Ivoire | 2.598 |
| (2.173) | |
| factor(country)Croatia | 9.463*** |
| (2.270) | |
| factor(country)Cuba | 13.513*** |
| (2.243) | |
| factor(country)Czech Republic | 6.555*** |
| (2.314) | |
| factor(country)Denmark | 5.345** |
| (2.361) | |
| factor(country)Djibouti | -2.834 |
| (2.188) | |
| factor(country)Dominican Republic | 12.583*** |
| (2.186) | |
| factor(country)Ecuador | 6.361*** |
| (2.234) | |
| factor(country)Egypt | 6.577*** |
| (2.190) | |
| factor(country)El Salvador | 5.535** |
| (2.217) | |
| factor(country)Equatorial Guinea | 1.062 |
| (2.163) | |
| factor(country)Eritrea | 12.711*** |
| (2.162) | |
| factor(country)Ethiopia | 11.529*** |
| (2.163) | |
| factor(country)Finland | 6.477*** |
| (2.331) | |
| factor(country)France | 6.910*** |
| (2.342) | |
| factor(country)Gabon | -11.429*** |
| (2.290) | |
| factor(country)Gambia | 8.599*** |
| (2.159) | |
| factor(country)Germany | 5.064** |
| (2.353) | |
| factor(country)Ghana | 12.292*** |
| (2.160) | |
| factor(country)Greece | 9.832*** |
| (2.303) | |
| factor(country)Guatemala | 3.754* |
| (2.209) | |
| factor(country)Guinea | 6.163*** |
| (2.158) | |
| factor(country)Guinea-Bissau | 4.059* |
| (2.160) | |
| factor(country)Haiti | 5.867*** |
| (2.168) | |
| factor(country)Honduras | 8.154*** |
| (2.190) | |
| factor(country)Hong Kong, China | 9.664*** |
| (2.302) | |
| factor(country)Hungary | 6.901*** |
| (2.288) | |
| factor(country)Iceland | 8.393*** |
| (2.350) | |
| factor(country)India | 13.971*** |
| (2.159) | |
| factor(country)Indonesia | 10.675*** |
| (2.167) | |
| factor(country)Iran | 0.256 |
| (2.250) | |
| factor(country)Iraq | -1.896 |
| (2.251) | |
| factor(country)Ireland | 8.536*** |
| (2.309) | |
| factor(country)Israel | 9.395*** |
| (2.306) | |
| factor(country)Italy | 8.460*** |
| (2.321) | |
| factor(country)Jamaica | 11.500*** |
| (2.240) | |
| factor(country)Japan | 9.201*** |
| (2.321) | |
| factor(country)Jordan | 9.472*** |
| (2.193) | |
| factor(country)Kenya | 11.335*** |
| (2.162) | |
| factor(country)Korea, Dem. Rep. | 15.650*** |
| (2.182) | |
| factor(country)Korea, Rep. | 9.455*** |
| (2.227) | |
| factor(country)Kuwait | -10.344*** |
| (2.504) | |
| factor(country)Lebanon | 7.039*** |
| (2.254) | |
| factor(country)Lesotho | 13.976*** |
| (2.159) | |
| factor(country)Liberia | 7.804*** |
| (2.160) | |
| factor(country)Libya | -2.924 |
| (2.285) | |
| factor(country)Madagascar | 5.475** |
| (2.163) | |
| factor(country)Malawi | 9.266*** |
| (2.161) | |
| factor(country)Malaysia | 10.309*** |
| (2.216) | |
| factor(country)Mali | 7.889*** |
| (2.159) | |
| factor(country)Mauritania | 9.907*** |
| (2.164) | |
| factor(country)Mauritius | 11.610*** |
| (2.212) | |
| factor(country)Mexico | 6.405*** |
| (2.255) | |
| factor(country)Mongolia | 11.758*** |
| (2.168) | |
| factor(country)Montenegro | 12.110*** |
| (2.248) | |
| factor(country)Morocco | 9.573*** |
| (2.182) | |
| factor(country)Mozambique | 6.881*** |
| (2.162) | |
| factor(country)Myanmar | 22.143*** |
| (2.167) | |
| factor(country)Namibia | 1.226 |
| (2.205) | |
| factor(country)Nepal | 11.912*** |
| (2.158) | |
| factor(country)Netherlands | 6.690*** |
| (2.360) | |
| factor(country)New Zealand | 6.752*** |
| (2.340) | |
| factor(country)Nicaragua | 7.015*** |
| (2.199) | |
| factor(country)Niger | 7.482*** |
| (2.158) | |
| factor(country)Nigeria | 0.236 |
| (2.166) | |
| factor(country)Norway | 5.266** |
| (2.381) | |
| factor(country)Oman | -2.695 |
| (2.275) | |
| factor(country)Pakistan | 12.524*** |
| (2.164) | |
| factor(country)Panama | 11.706*** |
| (2.231) | |
| factor(country)Paraguay | 16.062*** |
| (2.196) | |
| factor(country)Peru | 2.468 |
| (2.233) | |
| factor(country)Philippines | 13.968*** |
| (2.178) | |
| factor(country)Poland | 10.312*** |
| (2.263) | |
| factor(country)Portugal | 8.670*** |
| (2.281) | |
| factor(country)Puerto Rico | 11.211*** |
| (2.279) | |
| factor(country)Reunion | 11.900*** |
| (2.221) | |
| factor(country)Romania | 9.763*** |
| (2.251) | |
| factor(country)Rwanda | 5.802*** |
| (2.159) | |
| factor(country)Sao Tome and Principe | 15.282*** |
| (2.164) | |
| factor(country)Saudi Arabia | -9.374*** |
| (2.350) | |
| factor(country)Senegal | 6.765*** |
| (2.167) | |
| factor(country)Serbia | 7.871*** |
| (2.270) | |
| factor(country)Sierra Leone | -3.283 |
| (2.160) | |
| factor(country)Singapore | 7.762*** |
| (2.298) | |
| factor(country)Slovak Republic | 8.651*** |
| (2.284) | |
| factor(country)Slovenia | 7.244*** |
| (2.307) | |
| factor(country)Somalia | 0.122 |
| (2.161) | |
| factor(country)South Africa | -4.901** |
| (2.254) | |
| factor(country)Spain | 10.452*** |
| (2.301) | |
| factor(country)Sri Lanka | 21.770*** |
| (2.170) | |
| factor(country)Sudan | 2.888 |
| (2.172) | |
| factor(country)Swaziland | -1.101 |
| (2.192) | |
| factor(country)Sweden | 7.988*** |
| (2.351) | |
| factor(country)Switzerland | 3.963* |
| (2.394) | |
| factor(country)Syria | 11.325*** |
| (2.192) | |
| factor(country)Taiwan | 12.772*** |
| (2.243) | |
| factor(country)Tanzania | 9.881*** |
| (2.159) | |
| factor(country)Thailand | 14.436*** |
| (2.181) | |
| factor(country)Togo | 10.646*** |
| (2.161) | |
| factor(country)Trinidad and Tobago | 7.972*** |
| (2.254) | |
| factor(country)Tunisia | 10.108*** |
| (2.195) | |
| factor(country)Turkey | 6.352*** |
| (2.212) | |
| factor(country)Uganda | 10.085*** |
| (2.158) | |
| factor(country)United Kingdom | 5.909** |
| (2.349) | |
| factor(country)United States | 2.477 |
| (2.386) | |
| factor(country)Uruguay | 12.138*** |
| (2.252) | |
| factor(country)Venezuela | 4.388* |
| (2.285) | |
| factor(country)Vietnam | 18.639*** |
| (2.159) | |
| factor(country)West Bank and Gaza | 8.657*** |
| (2.201) | |
| factor(country)Yemen, Rep. | 3.485 |
| (2.166) | |
| factor(country)Zambia | 3.440 |
| (2.164) | |
| factor(country)Zimbabwe | 17.592*** |
| (2.160) | |
| log(gdpPercap) | 9.769*** |
| (0.297) | |
| Num.Obs. | 1704 |
| R2 | 0.846 |
| R2 Adj. | 0.832 |
| AIC | 10649.6 |
| BIC | 11433.1 |
| Log.Lik. | -5180.806 |
| F | 60.592 |
| RMSE | 5.06 |
import numpy as np
import statsmodels.formula.api as smf
from causaldata import gapminder
gm = gapminder.load_pandas().data
gm['logGDPpercap'] = np.log(gm['gdpPercap'])
# Use C() to include binary variables for a categorical variable
m2 = smf.ols(formula = '''lifeExp ~ logGDPpercap + C(country)''', data = gm).fit()
m2.summary()| Dep. Variable: | lifeExp | R-squared: | 0.846 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.832 |
| Method: | Least Squares | F-statistic: | 60.59 |
| Date: | 周六, 23 8月 2025 | Prob (F-statistic): | 0.00 |
| Time: | 23:50:24 | Log-Likelihood: | -5180.8 |
| No. Observations: | 1704 | AIC: | 1.065e+04 |
| Df Residuals: | 1561 | BIC: | 1.143e+04 |
| Df Model: | 142 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -27.7735 | 2.501 | -11.107 | 0.000 | -32.678 | -22.869 |
| C(country)[T.Albania] | 17.7826 | 2.195 | 8.101 | 0.000 | 13.477 | 22.088 |
| C(country)[T.Algeria] | 5.2411 | 2.214 | 2.367 | 0.018 | 0.897 | 9.585 |
| C(country)[T.Angola] | -13.9071 | 2.202 | -6.316 | 0.000 | -18.226 | -9.588 |
| C(country)[T.Argentina] | 8.1322 | 2.273 | 3.578 | 0.000 | 3.674 | 12.590 |
| C(country)[T.Australia] | 6.4003 | 2.352 | 2.721 | 0.007 | 1.787 | 11.014 |
| C(country)[T.Austria] | 5.1559 | 2.348 | 2.196 | 0.028 | 0.550 | 9.762 |
| C(country)[T.Bahrain] | -1.9682 | 2.344 | -0.840 | 0.401 | -6.566 | 2.629 |
| C(country)[T.Bangladesh] | 12.4011 | 2.158 | 5.745 | 0.000 | 8.167 | 16.635 |
| C(country)[T.Belgium] | 5.6169 | 2.349 | 2.391 | 0.017 | 1.009 | 10.225 |
| C(country)[T.Benin] | 7.7450 | 2.161 | 3.584 | 0.000 | 3.506 | 11.984 |
| C(country)[T.Bolivia] | 2.3380 | 2.193 | 1.066 | 0.286 | -1.963 | 6.639 |
| C(country)[T.Bosnia and Herzegovina] | 17.3528 | 2.194 | 7.911 | 0.000 | 13.050 | 21.655 |
| C(country)[T.Botswana] | 3.1522 | 2.200 | 1.433 | 0.152 | -1.162 | 7.467 |
| C(country)[T.Brazil] | 6.3181 | 2.230 | 2.833 | 0.005 | 1.944 | 10.692 |
| C(country)[T.Bulgaria] | 12.6702 | 2.239 | 5.659 | 0.000 | 8.279 | 17.062 |
| C(country)[T.Burkina Faso] | 6.8580 | 2.158 | 3.177 | 0.002 | 2.624 | 11.092 |
| C(country)[T.Burundi] | 12.6464 | 2.164 | 5.843 | 0.000 | 8.401 | 16.892 |
| C(country)[T.Cambodia] | 12.9167 | 2.160 | 5.981 | 0.000 | 8.680 | 17.153 |
| C(country)[T.Cameroon] | 3.0620 | 2.171 | 1.411 | 0.159 | -1.196 | 7.320 |
| C(country)[T.Canada] | 5.4647 | 2.366 | 2.309 | 0.021 | 0.823 | 10.106 |
| C(country)[T.Central African Republic] | 4.7576 | 2.159 | 2.204 | 0.028 | 0.523 | 8.992 |
| C(country)[T.Chad] | 5.7370 | 2.161 | 2.655 | 0.008 | 1.498 | 9.976 |
| C(country)[T.Chile] | 9.9150 | 2.242 | 4.421 | 0.000 | 5.516 | 14.314 |
| C(country)[T.China] | 21.2712 | 2.160 | 9.846 | 0.000 | 17.034 | 25.509 |
| C(country)[T.Colombia] | 10.9197 | 2.209 | 4.943 | 0.000 | 6.587 | 15.253 |
| C(country)[T.Comoros] | 10.2055 | 2.163 | 4.718 | 0.000 | 5.963 | 14.448 |
| C(country)[T.Congo, Dem. Rep.] | 10.0249 | 2.160 | 4.640 | 0.000 | 5.787 | 14.262 |
| C(country)[T.Congo, Rep.] | 1.3751 | 2.198 | 0.626 | 0.532 | -2.936 | 5.686 |
| C(country)[T.Costa Rica] | 14.5522 | 2.228 | 6.533 | 0.000 | 10.183 | 18.922 |
| C(country)[T.Cote d'Ivoire] | 2.5978 | 2.173 | 1.195 | 0.232 | -1.665 | 6.861 |
| C(country)[T.Croatia] | 9.4630 | 2.270 | 4.169 | 0.000 | 5.011 | 13.915 |
| C(country)[T.Cuba] | 13.5129 | 2.243 | 6.025 | 0.000 | 9.114 | 17.912 |
| C(country)[T.Czech Republic] | 6.5554 | 2.314 | 2.833 | 0.005 | 2.017 | 11.094 |
| C(country)[T.Denmark] | 5.3450 | 2.361 | 2.264 | 0.024 | 0.713 | 9.977 |
| C(country)[T.Djibouti] | -2.8341 | 2.188 | -1.296 | 0.195 | -7.125 | 1.457 |
| C(country)[T.Dominican Republic] | 12.5829 | 2.186 | 5.755 | 0.000 | 8.294 | 16.872 |
| C(country)[T.Ecuador] | 6.3606 | 2.234 | 2.847 | 0.004 | 1.979 | 10.743 |
| C(country)[T.Egypt] | 6.5774 | 2.190 | 3.003 | 0.003 | 2.282 | 10.873 |
| C(country)[T.El Salvador] | 5.5348 | 2.217 | 2.497 | 0.013 | 1.187 | 9.883 |
| C(country)[T.Equatorial Guinea] | 1.0619 | 2.163 | 0.491 | 0.623 | -3.180 | 5.304 |
| C(country)[T.Eritrea] | 12.7109 | 2.162 | 5.879 | 0.000 | 8.470 | 16.952 |
| C(country)[T.Ethiopia] | 11.5290 | 2.163 | 5.331 | 0.000 | 7.287 | 15.771 |
| C(country)[T.Finland] | 6.4773 | 2.331 | 2.778 | 0.006 | 1.904 | 11.050 |
| C(country)[T.France] | 6.9102 | 2.342 | 2.950 | 0.003 | 2.316 | 11.504 |
| C(country)[T.Gabon] | -11.4293 | 2.290 | -4.992 | 0.000 | -15.920 | -6.938 |
| C(country)[T.Gambia] | 8.5989 | 2.159 | 3.983 | 0.000 | 4.364 | 12.834 |
| C(country)[T.Germany] | 5.0641 | 2.353 | 2.152 | 0.032 | 0.448 | 9.680 |
| C(country)[T.Ghana] | 12.2924 | 2.160 | 5.691 | 0.000 | 8.056 | 16.529 |
| C(country)[T.Greece] | 9.8322 | 2.303 | 4.270 | 0.000 | 5.316 | 14.349 |
| C(country)[T.Guatemala] | 3.7542 | 2.209 | 1.699 | 0.089 | -0.579 | 8.087 |
| C(country)[T.Guinea] | 6.1635 | 2.158 | 2.855 | 0.004 | 1.930 | 10.397 |
| C(country)[T.Guinea-Bissau] | 4.0586 | 2.160 | 1.879 | 0.060 | -0.177 | 8.295 |
| C(country)[T.Haiti] | 5.8666 | 2.168 | 2.706 | 0.007 | 1.613 | 10.120 |
| C(country)[T.Honduras] | 8.1544 | 2.190 | 3.723 | 0.000 | 3.858 | 12.451 |
| C(country)[T.Hong Kong, China] | 9.6640 | 2.302 | 4.198 | 0.000 | 5.149 | 14.179 |
| C(country)[T.Hungary] | 6.9008 | 2.288 | 3.016 | 0.003 | 2.413 | 11.389 |
| C(country)[T.Iceland] | 8.3929 | 2.350 | 3.571 | 0.000 | 3.783 | 13.003 |
| C(country)[T.India] | 13.9706 | 2.159 | 6.471 | 0.000 | 9.736 | 18.206 |
| C(country)[T.Indonesia] | 10.6745 | 2.167 | 4.927 | 0.000 | 6.425 | 14.924 |
| C(country)[T.Iran] | 0.2559 | 2.250 | 0.114 | 0.909 | -4.157 | 4.669 |
| C(country)[T.Iraq] | -1.8960 | 2.251 | -0.842 | 0.400 | -6.310 | 2.519 |
| C(country)[T.Ireland] | 8.5358 | 2.309 | 3.697 | 0.000 | 4.007 | 13.065 |
| C(country)[T.Israel] | 9.3947 | 2.306 | 4.073 | 0.000 | 4.871 | 13.919 |
| C(country)[T.Italy] | 8.4601 | 2.321 | 3.646 | 0.000 | 3.908 | 13.012 |
| C(country)[T.Jamaica] | 11.5001 | 2.240 | 5.133 | 0.000 | 7.106 | 15.894 |
| C(country)[T.Japan] | 9.2015 | 2.321 | 3.964 | 0.000 | 4.648 | 13.755 |
| C(country)[T.Jordan] | 9.4722 | 2.193 | 4.319 | 0.000 | 5.170 | 13.774 |
| C(country)[T.Kenya] | 11.3351 | 2.162 | 5.244 | 0.000 | 7.095 | 15.575 |
| C(country)[T.Korea, Dem. Rep.] | 15.6495 | 2.182 | 7.173 | 0.000 | 11.370 | 19.929 |
| C(country)[T.Korea, Rep.] | 9.4547 | 2.227 | 4.245 | 0.000 | 5.086 | 13.823 |
| C(country)[T.Kuwait] | -10.3442 | 2.504 | -4.132 | 0.000 | -15.255 | -5.433 |
| C(country)[T.Lebanon] | 7.0391 | 2.254 | 3.123 | 0.002 | 2.619 | 11.460 |
| C(country)[T.Lesotho] | 13.9758 | 2.159 | 6.474 | 0.000 | 9.741 | 18.210 |
| C(country)[T.Liberia] | 7.8041 | 2.160 | 3.613 | 0.000 | 3.567 | 12.041 |
| C(country)[T.Libya] | -2.9242 | 2.285 | -1.279 | 0.201 | -7.407 | 1.559 |
| C(country)[T.Madagascar] | 5.4748 | 2.163 | 2.531 | 0.011 | 1.231 | 9.718 |
| C(country)[T.Malawi] | 9.2660 | 2.161 | 4.288 | 0.000 | 5.027 | 13.505 |
| C(country)[T.Malaysia] | 10.3093 | 2.216 | 4.653 | 0.000 | 5.963 | 14.655 |
| C(country)[T.Mali] | 7.8890 | 2.159 | 3.654 | 0.000 | 3.654 | 12.124 |
| C(country)[T.Mauritania] | 9.9067 | 2.164 | 4.579 | 0.000 | 5.663 | 14.151 |
| C(country)[T.Mauritius] | 11.6102 | 2.212 | 5.250 | 0.000 | 7.272 | 15.948 |
| C(country)[T.Mexico] | 6.4046 | 2.255 | 2.840 | 0.005 | 1.981 | 10.828 |
| C(country)[T.Mongolia] | 11.7578 | 2.168 | 5.424 | 0.000 | 7.506 | 16.010 |
| C(country)[T.Montenegro] | 12.1100 | 2.248 | 5.387 | 0.000 | 7.700 | 16.520 |
| C(country)[T.Morocco] | 9.5733 | 2.182 | 4.387 | 0.000 | 5.293 | 13.853 |
| C(country)[T.Mozambique] | 6.8808 | 2.162 | 3.183 | 0.001 | 2.640 | 11.121 |
| C(country)[T.Myanmar] | 22.1428 | 2.167 | 10.219 | 0.000 | 17.892 | 26.393 |
| C(country)[T.Namibia] | 1.2264 | 2.205 | 0.556 | 0.578 | -3.098 | 5.551 |
| C(country)[T.Nepal] | 11.9117 | 2.158 | 5.519 | 0.000 | 7.678 | 16.146 |
| C(country)[T.Netherlands] | 6.6905 | 2.360 | 2.834 | 0.005 | 2.060 | 11.321 |
| C(country)[T.New Zealand] | 6.7521 | 2.340 | 2.886 | 0.004 | 2.163 | 11.342 |
| C(country)[T.Nicaragua] | 7.0148 | 2.199 | 3.190 | 0.001 | 2.701 | 11.328 |
| C(country)[T.Niger] | 7.4817 | 2.158 | 3.466 | 0.001 | 3.248 | 11.715 |
| C(country)[T.Nigeria] | 0.2362 | 2.166 | 0.109 | 0.913 | -4.012 | 4.484 |
| C(country)[T.Norway] | 5.2662 | 2.381 | 2.212 | 0.027 | 0.596 | 9.936 |
| C(country)[T.Oman] | -2.6950 | 2.275 | -1.185 | 0.236 | -7.157 | 1.767 |
| C(country)[T.Pakistan] | 12.5242 | 2.164 | 5.789 | 0.000 | 8.280 | 16.768 |
| C(country)[T.Panama] | 11.7063 | 2.231 | 5.247 | 0.000 | 7.330 | 16.083 |
| C(country)[T.Paraguay] | 16.0615 | 2.196 | 7.315 | 0.000 | 11.755 | 20.368 |
| C(country)[T.Peru] | 2.4684 | 2.233 | 1.105 | 0.269 | -1.913 | 6.849 |
| C(country)[T.Philippines] | 13.9677 | 2.178 | 6.414 | 0.000 | 9.696 | 18.239 |
| C(country)[T.Poland] | 10.3117 | 2.263 | 4.557 | 0.000 | 5.873 | 14.750 |
| C(country)[T.Portugal] | 8.6700 | 2.281 | 3.801 | 0.000 | 4.196 | 13.144 |
| C(country)[T.Puerto Rico] | 11.2107 | 2.279 | 4.920 | 0.000 | 6.741 | 15.680 |
| C(country)[T.Reunion] | 11.9001 | 2.221 | 5.358 | 0.000 | 7.543 | 16.257 |
| C(country)[T.Romania] | 9.7631 | 2.251 | 4.337 | 0.000 | 5.348 | 14.178 |
| C(country)[T.Rwanda] | 5.8015 | 2.159 | 2.687 | 0.007 | 1.566 | 10.037 |
| C(country)[T.Sao Tome and Principe] | 15.2824 | 2.164 | 7.062 | 0.000 | 11.038 | 19.527 |
| C(country)[T.Saudi Arabia] | -9.3738 | 2.350 | -3.990 | 0.000 | -13.982 | -4.765 |
| C(country)[T.Senegal] | 6.7647 | 2.167 | 3.122 | 0.002 | 2.514 | 11.016 |
| C(country)[T.Serbia] | 7.8708 | 2.270 | 3.467 | 0.001 | 3.417 | 12.324 |
| C(country)[T.Sierra Leone] | -3.2832 | 2.160 | -1.520 | 0.129 | -7.520 | 0.953 |
| C(country)[T.Singapore] | 7.7617 | 2.298 | 3.378 | 0.001 | 3.254 | 12.269 |
| C(country)[T.Slovak Republic] | 8.6514 | 2.284 | 3.788 | 0.000 | 4.172 | 13.131 |
| C(country)[T.Slovenia] | 7.2443 | 2.307 | 3.139 | 0.002 | 2.718 | 11.770 |
| C(country)[T.Somalia] | 0.1222 | 2.161 | 0.057 | 0.955 | -4.116 | 4.361 |
| C(country)[T.South Africa] | -4.9006 | 2.254 | -2.174 | 0.030 | -9.322 | -0.479 |
| C(country)[T.Spain] | 10.4524 | 2.301 | 4.542 | 0.000 | 5.939 | 14.966 |
| C(country)[T.Sri Lanka] | 21.7701 | 2.170 | 10.034 | 0.000 | 17.514 | 26.026 |
| C(country)[T.Sudan] | 2.8876 | 2.172 | 1.329 | 0.184 | -1.373 | 7.148 |
| C(country)[T.Swaziland] | -1.1013 | 2.192 | -0.502 | 0.615 | -5.401 | 3.199 |
| C(country)[T.Sweden] | 7.9883 | 2.351 | 3.398 | 0.001 | 3.377 | 12.600 |
| C(country)[T.Switzerland] | 3.9627 | 2.394 | 1.655 | 0.098 | -0.733 | 8.659 |
| C(country)[T.Syria] | 11.3248 | 2.192 | 5.167 | 0.000 | 7.026 | 15.624 |
| C(country)[T.Taiwan] | 12.7718 | 2.243 | 5.694 | 0.000 | 8.372 | 17.171 |
| C(country)[T.Tanzania] | 9.8807 | 2.159 | 4.578 | 0.000 | 5.647 | 14.115 |
| C(country)[T.Thailand] | 14.4358 | 2.181 | 6.619 | 0.000 | 10.158 | 18.714 |
| C(country)[T.Togo] | 10.6459 | 2.161 | 4.927 | 0.000 | 6.407 | 14.884 |
| C(country)[T.Trinidad and Tobago] | 7.9716 | 2.254 | 3.537 | 0.000 | 3.551 | 12.393 |
| C(country)[T.Tunisia] | 10.1077 | 2.195 | 4.605 | 0.000 | 5.802 | 14.413 |
| C(country)[T.Turkey] | 6.3521 | 2.212 | 2.872 | 0.004 | 2.014 | 10.690 |
| C(country)[T.Uganda] | 10.0855 | 2.158 | 4.673 | 0.000 | 5.852 | 14.319 |
| C(country)[T.United Kingdom] | 5.9090 | 2.349 | 2.516 | 0.012 | 1.301 | 10.517 |
| C(country)[T.United States] | 2.4768 | 2.386 | 1.038 | 0.299 | -2.204 | 7.157 |
| C(country)[T.Uruguay] | 12.1379 | 2.252 | 5.390 | 0.000 | 7.721 | 16.555 |
| C(country)[T.Venezuela] | 4.3883 | 2.285 | 1.920 | 0.055 | -0.094 | 8.870 |
| C(country)[T.Vietnam] | 18.6392 | 2.159 | 8.634 | 0.000 | 14.405 | 22.874 |
| C(country)[T.West Bank and Gaza] | 8.6575 | 2.201 | 3.933 | 0.000 | 4.340 | 12.975 |
| C(country)[T.Yemen, Rep.] | 3.4847 | 2.166 | 1.609 | 0.108 | -0.763 | 7.733 |
| C(country)[T.Zambia] | 3.4400 | 2.164 | 1.590 | 0.112 | -0.805 | 7.685 |
| C(country)[T.Zimbabwe] | 17.5924 | 2.160 | 8.146 | 0.000 | 13.356 | 21.829 |
| logGDPpercap | 9.7690 | 0.297 | 32.944 | 0.000 | 9.187 | 10.351 |
| Omnibus: | 23.105 | Durbin-Watson: | 0.626 |
|---|---|---|---|
| Prob(Omnibus): | 0.000 | Jarque-Bera (JB): | 33.240 |
| Skew: | -0.143 | Prob(JB): | 6.05e-08 |
| Kurtosis: | 3.622 | Cond. No. | 1.23e+03 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.23e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
让我们看看两种方法的结果，使用\(表16.3\)中的 \(R\) 输出。用第二种方法我们还估计了一堆国家系数，但这些没有列在表上（因为数量太多了）。
# Gapminder FE
library(tidyverse)
library(readr)
gapminder <- read_csv('FixedEffects/gapminder.csv', show_col_types = FALSE)
gapminder <- gapminder %>%
# Put GDP per capita in log format since it's very skewed
mutate(log_GDPperCap = log(gdpPercap), show_col_types=FALSE) %>%
# Perform each calculation by group
group_by(country) %>%
# Get within variation by subtracting out the mean
mutate(lifeExp_within = lifeExp - mean(lifeExp), log_GDPperCap_within = log_GDPperCap - mean(log_GDPperCap))
# We no longer need the grouping ungroup()
# Analyze the within variation
m1 <- lm(lifeExp_within ~ log_GDPperCap_within, data = gapminder)
m2 <- lm(lifeExp ~ factor(country) + log_GDPperCap, data = gapminder)
coef(m2)['factor(country)India'] ## factor(country)India
## 13.97062coef(m2)['factor(country)Brazil']## factor(country)Brazil
## 6.31808# With dummies
knitr::opts_current$set(label = 'fixedeffects-basic-regs') ## Warning in set2(resolve(...)): The object is read-only and cannot be modified.
## If you have to modify it for a legitimate reason, call the method $lock(FALSE)
## on the object before $set(). Using $lock(FALSE) to modify the object will be
## enforced in future versions of knitr and this warning will become an error.summary(list('Life Expectancy (within)' = m1, 'Life Expectancy' = m2),
output = 'FixedEffects/fixedeffects_basic_regs.html',
coef_map = c('log_GDPperCap' = 'Log GDP per Capita', 'log_GDPperCap_within' = 'Log GDP per Capita (within)'),
gof_omit = 'Adj|AIC|BIC|Lik|F',
stars = TRUE) ## Length Class Mode
## Life Expectancy (within) 12 lm list
## Life Expectancy 13 lm listm1##
## Call:
## lm(formula = lifeExp_within ~ log_GDPperCap_within, data = gapminder)
##
## Coefficients:
## (Intercept) log_GDPperCap_within
## -4.453e-17 9.769e+00m2##
## Call:
## lm(formula = lifeExp ~ factor(country) + log_GDPperCap, data = gapminder)
##
## Coefficients:
## (Intercept)
## -27.7735
## factor(country)Albania
## 17.7826
## factor(country)Algeria
## 5.2411
## factor(country)Angola
## -13.9071
## factor(country)Argentina
## 8.1322
## factor(country)Australia
## 6.4003
## factor(country)Austria
## 5.1559
## factor(country)Bahrain
## -1.9682
## factor(country)Bangladesh
## 12.4011
## factor(country)Belgium
## 5.6169
## factor(country)Benin
## 7.7450
## factor(country)Bolivia
## 2.3380
## factor(country)Bosnia and Herzegovina
## 17.3528
## factor(country)Botswana
## 3.1522
## factor(country)Brazil
## 6.3181
## factor(country)Bulgaria
## 12.6702
## factor(country)Burkina Faso
## 6.8580
## factor(country)Burundi
## 12.6464
## factor(country)Cambodia
## 12.9167
## factor(country)Cameroon
## 3.0620
## factor(country)Canada
## 5.4647
## factor(country)Central African Republic
## 4.7576
## factor(country)Chad
## 5.7370
## factor(country)Chile
## 9.9150
## factor(country)China
## 21.2712
## factor(country)Colombia
## 10.9197
## factor(country)Comoros
## 10.2055
## factor(country)Congo, Dem. Rep.
## 10.0249
## factor(country)Congo, Rep.
## 1.3751
## factor(country)Costa Rica
## 14.5522
## factor(country)Cote d'Ivoire
## 2.5978
## factor(country)Croatia
## 9.4630
## factor(country)Cuba
## 13.5129
## factor(country)Czech Republic
## 6.5554
## factor(country)Denmark
## 5.3450
## factor(country)Djibouti
## -2.8341
## factor(country)Dominican Republic
## 12.5829
## factor(country)Ecuador
## 6.3606
## factor(country)Egypt
## 6.5774
## factor(country)El Salvador
## 5.5348
## factor(country)Equatorial Guinea
## 1.0619
## factor(country)Eritrea
## 12.7109
## factor(country)Ethiopia
## 11.5290
## factor(country)Finland
## 6.4773
## factor(country)France
## 6.9102
## factor(country)Gabon
## -11.4293
## factor(country)Gambia
## 8.5989
## factor(country)Germany
## 5.0641
## factor(country)Ghana
## 12.2924
## factor(country)Greece
## 9.8322
## factor(country)Guatemala
## 3.7542
## factor(country)Guinea
## 6.1635
## factor(country)Guinea-Bissau
## 4.0586
## factor(country)Haiti
## 5.8666
## factor(country)Honduras
## 8.1544
## factor(country)Hong Kong, China
## 9.6640
## factor(country)Hungary
## 6.9008
## factor(country)Iceland
## 8.3929
## factor(country)India
## 13.9706
## factor(country)Indonesia
## 10.6745
## factor(country)Iran
## 0.2559
## factor(country)Iraq
## -1.8960
## factor(country)Ireland
## 8.5358
## factor(country)Israel
## 9.3947
## factor(country)Italy
## 8.4601
## factor(country)Jamaica
## 11.5001
## factor(country)Japan
## 9.2015
## factor(country)Jordan
## 9.4722
## factor(country)Kenya
## 11.3351
## factor(country)Korea, Dem. Rep.
## 15.6495
## factor(country)Korea, Rep.
## 9.4547
## factor(country)Kuwait
## -10.3442
## factor(country)Lebanon
## 7.0391
## factor(country)Lesotho
## 13.9758
## factor(country)Liberia
## 7.8041
## factor(country)Libya
## -2.9242
## factor(country)Madagascar
## 5.4748
## factor(country)Malawi
## 9.2660
## factor(country)Malaysia
## 10.3093
## factor(country)Mali
## 7.8890
## factor(country)Mauritania
## 9.9067
## factor(country)Mauritius
## 11.6102
## factor(country)Mexico
## 6.4046
## factor(country)Mongolia
## 11.7578
## factor(country)Montenegro
## 12.1100
## factor(country)Morocco
## 9.5733
## factor(country)Mozambique
## 6.8808
## factor(country)Myanmar
## 22.1428
## factor(country)Namibia
## 1.2264
## factor(country)Nepal
## 11.9117
## factor(country)Netherlands
## 6.6905
## factor(country)New Zealand
## 6.7521
## factor(country)Nicaragua
## 7.0148
## factor(country)Niger
## 7.4817
## factor(country)Nigeria
## 0.2362
## factor(country)Norway
## 5.2662
## factor(country)Oman
## -2.6950
## factor(country)Pakistan
## 12.5242
## factor(country)Panama
## 11.7063
## factor(country)Paraguay
## 16.0615
## factor(country)Peru
## 2.4684
## factor(country)Philippines
## 13.9677
## factor(country)Poland
## 10.3117
## factor(country)Portugal
## 8.6700
## factor(country)Puerto Rico
## 11.2107
## factor(country)Reunion
## 11.9001
## factor(country)Romania
## 9.7631
## factor(country)Rwanda
## 5.8015
## factor(country)Sao Tome and Principe
## 15.2824
## factor(country)Saudi Arabia
## -9.3738
## factor(country)Senegal
## 6.7647
## factor(country)Serbia
## 7.8708
## factor(country)Sierra Leone
## -3.2832
## factor(country)Singapore
## 7.7617
## factor(country)Slovak Republic
## 8.6514
## factor(country)Slovenia
## 7.2443
## factor(country)Somalia
## 0.1222
## factor(country)South Africa
## -4.9006
## factor(country)Spain
## 10.4524
## factor(country)Sri Lanka
## 21.7701
## factor(country)Sudan
## 2.8876
## factor(country)Swaziland
## -1.1013
## factor(country)Sweden
## 7.9883
## factor(country)Switzerland
## 3.9627
## factor(country)Syria
## 11.3248
## factor(country)Taiwan
## 12.7718
## factor(country)Tanzania
## 9.8807
## factor(country)Thailand
## 14.4358
## factor(country)Togo
## 10.6459
## factor(country)Trinidad and Tobago
## 7.9716
## factor(country)Tunisia
## 10.1077
## factor(country)Turkey
## 6.3521
## factor(country)Uganda
## 10.0855
## factor(country)United Kingdom
## 5.9090
## factor(country)United States
## 2.4768
## factor(country)Uruguay
## 12.1379
## factor(country)Venezuela
## 4.3883
## factor(country)Vietnam
## 18.6392
## factor(country)West Bank and Gaza
## 8.6575
## factor(country)Yemen, Rep.
## 3.4847
## factor(country)Zambia
## 3.4400
## factor(country)Zimbabwe
## 17.5924
## log_GDPperCap
## 9.7690两种方法在标准误差方面存在一些细微差异，\(R^2\)方面存在一些较大差异，但系数是相同的17。既然我们有了系数，我们该如何解读它们呢？
问题3：我们如何解释回归估计结果？
固定效应模型中斜率系数的解释方式与控制其他变量时的回归分析相同 ——在同一国家内部 ，人均GDP对数 的变化如何与 预期寿命 的变化相关联？
换句话说，由于我们的系数是 \(9.769\)，可以说：对于一个特定的国家，在某一年中，如果其 \(人均GDP对数值\) 比该国典型水平 高出\(1\)个单位，那么该国的 \(预期寿命\)预计会比其典型水平长 \(9.769年\)18。
表上还有什么？我们可以看到，两列中的\(R^2\)值差异很大。为什么会这样？
记住，\(R^2\) 是一种基于残差变异程度相对于因变量变异程度的度量指标19。在第一列中，我们的因变量是 \(lifeExp\_within\) ，而非 \(lifeExp\) 。所以，这里的 \(R^2\) 是基于残差变异程度相对于组内变异（而非 \(lifeExp\) 的整体变异 ）的情况，这被称为 “组内 \(R^2\)” 。另一方面，第二列告诉我们，（完全相同的 ）这些残差的变异程度相对于 \(lifeExp\) 的所有变异而言是多少。它统计了由组间变异（我们的二元变量 ）解释的 \(lifeExp\) 的部分。
最后，让我们讨论一个表格中没有呈现的内容——固定效应本身。即我们在第二次回归中获得的那些国家二元变量的系数。例如，印度的系数是\(13.971\)，巴西的是\(6.318\)。这些数值相当于各国特定拟合线中的截距项，如 \(图16.9\) 所示。
对于这些固定效应，我们能得出什么结论？从绝对值来看意义不大20，但它们的相对比较是有意义的。因此，印度的截距比巴西高出 \(13.971-6.318=7.653\)。
我们可以这样解释：如果印度和巴西的\(人均GDP\)相同，我们预计印度的预期寿命将比巴西高出\(7.653\)岁。这里我们仍然是在观察一个变量对另一个变量的控制效果——只不过现在，我们观察的是国家变量在控制\(人均GDP\)后的影响，而非相反情况。有时研究者会利用这些固定效应估计值来对评估对象（国家）作出论断。例如，他们可能会说：考虑到\(人均GDP水平\)，印度的预期寿命显得特别高。
这种解读个体差异的方式或许颇具启发性。但需谨记：固定效应法的优势在于，它能有效控制一系列不随时间变化的未观测变量，从而聚焦少数时变变量的影响。而这些个体效应本身却无法同样有效地控制那些未观测的时变变量。因此，我们很难将这些个体固定效应估计值视为因果关系的体现21。
关于固定效应，最后需要牢记的是：它聚焦于组内变异。因此，所估计的处理效应（第\(10\)章）会更侧重于那些具有显著组内变异的个体。以预期寿命与\(GDP\)的例子来说，我们得到的\(GDP\)系数更接近于估计那些\(GDP\)随时间大幅波动国家中\(GDP\)对预期寿命的影响。如果一个国家的\(GDP\)保持相对稳定，它对估计结果的贡献就会较小！这个问题可以通过加权方法来解决22。
Charles E. Gibbons, Serrato Juan Carlos Suárez, and Michael B. Urbancic. Broken or fixed effects? Journal of Econometric Methods, 8(1):1–12, 2019.
16.2.2 多重固定效应组合
前文讨论中将固定效应视为控制分类变量的一种方法。这种方法最终能够提取出该变量内部的变异量。具体而言：当我们设置个体固定效应时，实际上是在对该个体进行跨时间维度的自我比较；而当我们设置城市固定效应时，则仅限于对同一城市内不同个体之间的比较分析。
为何不采用多重固定效应组合呢？设想我们对同一批个体进行多年追踪观测时，可以：
引入个体固定效应——通过比较个体在不同时间点的自身变化来捕捉时间维度变异；
设置年份固定效应——通过比较同一时间截面上不同个体的差异来分析截面特征。
这种同时包含个体和时间双重固定效应的方法，是多重固定效应组合的典型应用，在计量经济学中被称为”双向固定效应模型23“。其回归方程可表示为：
\[ Y_{it}=\beta_i+\beta_t+\beta_1 X_{it}+\epsilon_{it} \quad\quad(16.3) \]
这一设定能带来什么优势？此时，我们既能捕捉个体内部随时间的变化，也能控制年份层面的共同趋势。因此，模型最终识别的变异量可以理解为：在控制个体固定效应和年份固定效应的条件下，相对于预期值的偏离程度24。
例如，假设我们掌握以下面板数据：每位受访者的年度收入\(Y_{it}\)及其当年接受的培训时长\(X_{it}\)。本研究旨在解析在职培训经验对薪酬水平的因果效应。
我们确实知道，\(2009\) 年是收入情况格外糟糕的一年，这与大衰退等因素有关。假设 \(2009\) 年的收入比正常年份低\(\$10,000\)。再假设我们以安东尼为例，他的年收入为\(\$120,000\)，远高于普通人的年收入\(\$40,000\)。
\(2009\) 年，安东尼的收入仅为\(\$116,000\)。这个数字依然远高于平均水平，但这可以用 “他是安东尼” 这一事实来解释 —— 安东尼本就收入丰厚。因此，单从 “他是安东尼” 这一角度看，其收入比自身常规水平低了\(\$4,000\)。但从 “\(2009年\)” 这一背景来看，大多数人的收入都减少了\(\$10,000\)，而安东尼的收入仅减少了\(\$4,000\)。所以，综合 “他是安东尼” 和 “\(2009年\)” 这两个因素，\(2009\) 年安东尼的收入比预期水平高出了\(\$6,000\)。
需要说明的是，上述解释仅为简化说明。实际上，双向固定效应模型存在更复杂的交互机制：个体固定效应与年份固定效应会相互影响（反之亦然）。与单一固定效应类似，模型最终用于估计处理效应（第\(10\)章）的变异量，主要来自那些随时间波动显著的个体。例如：若安东尼的在职培训时长保持稳定，而卡玛拉的培训量逐年大幅波动，那么双向固定效应估计结果将更多反映卡玛拉的处理效应，而非安东尼的。
Clément De Chaisemartin and Xavier d’Haultfoeuille. Two-way fixed effects estimators with heterogeneous treatment effects. American Economic Review, 110(9):2964–96, 2020.当然，双向固定效应（个体与时间维度）并非多重固定效应组合的唯一形式。如上文所述，研究者亦可采用个体与城市维度的双重固定效应——需注意的是，这两个维度均不涉及时间因素。
此种设定下的经济学解释，更接近于在模型中同时控制个体特征变量与城市特征变量的传统做法。
但有一点需要记住：我们要分离的是 “个体内变异” 与 “城市内变异”。因此，若要让数据中存在可纳入模型的变异，你需要确保（研究对象）出现在多个城市中。相较于双向固定效应（个体与时间），这种问题（两组固定效应均与时间无关）更为常见 —— 因为通常而言，每个人在每个时间段内只会出现一次。
要是我们在收入与职业培训研究中纳入个体和城市固定效应，处理效应就只会基于那些在不同时间点，处理（职业培训等干预）情况有变化的个体和城市。某个个体或城市一直处于被处理状态（比如一直参与培训 ）？这种情况不会被纳入计算！
使用与常规”单向”固定效应估计量相同的 “二元控制变量”方法，即可估计含多组固定效应的模型 —— 只需一组固定效应。毫无问题！
然而，当需要估计的固定效应类型过多时，这种方法可能面临挑战。计算复杂度会呈指数级增长——除非各类固定效应仅包含少量分类（例如按左右利手或眼睛颜色等分类），否则通常建议使用专门处理多重固定效应的计量软件命令。
这类方法通常采用称为”交替投影法”的算法技术，其核心原理类似于基础版的组内变异计算法，但关键改进在于：能同步考量第一组固定效应对第二组组内变异的影响（反之亦然）。
让我们通过编码实现这一方法，继续使用 \(Gapminder\) 数据集进行演示。
library(tidyverse)
library(modelsummary)
library(fixest)##
## Attaching package: 'fixest'## The following object is masked from 'package:scales':
##
## pvalue## The following object is masked from 'package:lfe':
##
## fepoisgm <- causaldata::gapminder
# Run our two-way fixed effects model (TWFE).
# First the non-fixed effects part of the model, Then a |, then the fixed effects we want
twfe <- feols(lifeExp ~ log(gdpPercap) | country + year,data = gm)
# Note that standard errors will be clustered by the first fixed effect by default. Set se = 'standard' to not do this
msummary(twfe, stars = c('*' = .1, '**' = .05, '***' = .01), se = 'standrad')| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| log(gdpPercap) | 1.450** |
| (0.679) | |
| Num.Obs. | 1704 |
| R2 | 0.936 |
| R2 Adj. | 0.930 |
| R2 Within | 0.019 |
| R2 Within Adj. | 0.018 |
| AIC | 9178.1 |
| BIC | 10016.0 |
| RMSE | 3.27 |
| Std.Errors | by: country |
| FE: country | X |
| FE: year | X |
# Multiple FEs
library(tidyverse)
library(jtools)
library(lfe)
gm <- read_csv('FixedEffects/gapminder.csv', show_col_types = FALSE)
gm <- gm %>%
# Put GDP per capita in log format since it's very skewed
mutate(log_GDPperCap = log(gdpPercap))
# Run our two-way fixed effects model (TWFE)
# First the non-fixed effects part of the model Then a |, then the fixed effects we want
twfe <- felm(lifeExp ~ log_GDPperCap | country + year, data = gm)
export_summs(twfe, se = 'standrad')| Model 1 | |
|---|---|
| log_GDPperCap | 1.45 *** |
| (0.27) | |
| nobs | 1704 |
| r.squared | 0.94 |
| adj.r.squared | 0.93 |
| sigma | 3.43 |
| statistic | 148.17 |
| p.value | 0.00 |
| df | 1550.00 |
| df.residual | 1550.00 |
| nobs.1 | 1704.00 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | |
# Multiple FEs
library(tidyverse)
library(jtools)
library(lfe)
gm <- read_csv('FixedEffects/gapminder.csv', show_col_types = FALSE)
gm <- gm %>%
# Put GDP per capita in log format since it's very skewed
mutate(log_GDPperCap = log(gdpPercap))
# separates each "section" of the formula.
# Clusters go in the fourth section, and we don't need the third (instrumental variables) So we skip it with 0.
# y ~ x | FE | IV | Cluster
clfe <- felm(lifeExp ~ log_GDPperCap | country | 0 | country, data = gm)
export_summs(clfe)| Model 1 | |
|---|---|
| log_GDPperCap | 9.77 *** |
| (0.70) | |
| nobs | 1704 |
| r.squared | 0.85 |
| adj.r.squared | 0.83 |
| sigma | 5.29 |
| statistic | 60.59 |
| p.value | 0.00 |
| df | 1561.00 |
| df.residual | 1561.00 |
| nobs.1 | 1704.00 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | |
import linearmodels as lm
from causaldata import gapminder
gm = gapminder.load_pandas().data
gm['logGDPpercap'] = np.log(gm['gdpPercap'])
# Set our individual and time (index) for our data
gm = gm.set_index(['country','year'])
# Specify the regression modelAnd estimate with both sets of fixed effects
# EntityEffects and TimeEffects
# (this function can't handle more than two)
mod = lm.PanelOLS.from_formula('''lifeExp ~ logGDPpercap + EntityEffects + TimeEffects''', gm)
twfe = mod.fit()
print(twfe)## PanelOLS Estimation Summary
## ================================================================================
## Dep. Variable: lifeExp R-squared: 0.0186
## Estimator: PanelOLS R-squared (Between): 0.3542
## No. Observations: 1704 R-squared (Within): 0.1127
## Date: 周六, 8月 23 2025 R-squared (Overall): 0.3514
## Time: 23:50:27 Log-likelihood -4435.1
## Cov. Estimator: Unadjusted
## F-statistic: 29.366
## Entities: 142 P-value 0.0000
## Avg Obs: 12.000 Distribution: F(1,1550)
## Min Obs: 12.000
## Max Obs: 12.000 F-statistic (robust): 29.366
## P-value 0.0000
## Time periods: 12 Distribution: F(1,1550)
## Avg Obs: 142.00
## Min Obs: 142.00
## Max Obs: 142.00
##
## Parameter Estimates
## ================================================================================
## Parameter Std. Err. T-stat P-value Lower CI Upper CI
## --------------------------------------------------------------------------------
## logGDPpercap 1.4499 0.2676 5.4190 0.0000 0.9251 1.9747
## ================================================================================
##
## F-test for Poolability: 45.214
## P-value: 0.0000
## Distribution: F(152,1550)
##
## Included effects: Entity, Time16.2.3 随机效应模型
从研究设计维度看，固定效应的核心在于提取组内变异——通过剔除可能存在混杂因素的组间变异，从而聚焦组内变异。但就统计建模本质而言，固定效应可视为截距项\(\beta_i\)在个体间自由变动的特殊模型，其基本形式如 \(式（16.4）\)所示：
\[ Y_{it}=\beta_i+\beta_1 X_{it}+\epsilon_{it}\quad\quad(16.4) \]
事实证明，从统计学角度看，固定效应是一种让截距项变化的相对薄弱的方式。毕竟，我们只允许组内变异 。而且，每个个体实际有多少个观测值呢？我们可能只用少量观测值来估计这些\(\beta_i\)，这会使它们产生噪声（波动大、不准确 ），进而可能也会让我们对\(\beta_1\)的估计产生噪声 。
随机效应模型采用不同的建模思路：通过对个体截距项\(\beta_i\)施加概率分布约束（而非允许其完全自由变动），具体而言，该模型假定所有\(\beta_i\)均来自某个预设的随机分布——例如假设\(\beta_i \sim N(\mu, \sigma^2)\)服从正态分布。
随机效应模型的优势主要体现在：
估计精度提升：降低参数估计的标准误
个体效应改进：提高各 \(\beta_i\) 截距项的估计效率25
信息综合利用：采用组内变异与组间变异的加权组合，而非仅依赖组内变异
内生性控制：当 \(\beta_i\) 与解释变量 \(X_{it}\) 不相关时，其解决混杂因素的能力与固定效应相当
最后一点尤为关键：虽然统计精度的提升颇具吸引力，但我们最初采用固定效应模型的核心动机是解决研究设计中的内生性问题。需特别注意，随机效应模型仅在个体效应与右侧变量（含处理变量）不相关时，才具备同等的内生性控制能力。
这种假设往往难以成立。以 \(Gapminder\) 数据集为例：若要使随机效应模型适用，必须满足”国家效应”与\(GDP\)无关这一强假设——即所有影响预期寿命的、不随时间变化的国别特征（除\(人均GDP\)外），都必须与\(人均GDP\)完全无关。而这类特征显然包含工业发展水平、地理条件等重要因素，它们与\(人均GDP\)大概率存在关联。
基于上述分析，固定效应模型在绝大多数情况下都优于随机效应模型26——除非能确证右侧变量与个体效应不存在相关性。典型例外情况是随机对照实验：当\(X_{it}\)为严格随机分配的实验变量时（如临床试验中的药物剂量），其与个体效应的独立性假设自然成立。
处理这一问题的一个常见方法是使用 \(Durbin-Wu-Hausman\) 检验。\(Durbin-Wu-hausman\) 检验是一类广泛的检验方法，它通过比较不同模型的估计值来判定它们是否存在显著差异。
在固定效应与随机效应的比较中：
如果检验结果显示估计值无显著差异（即无法拒绝”估计值相同”的原假设）
则表明个体效应\(\beta_i\)与解释变量\(X_{it}\)之间的相关性较弱
此时可放心采用具有优良统计特性的随机效应模型
然而，我并不推荐使用 \(Durbin-Wu-Hausman\) 检验作为模型选择的依据27。主要原因在于：若缺乏强有力的理论支撑来确信”\(\beta_i\)与\(X_{it}\)不相关”这一假设成立，仅凭”未能拒绝原假设”的统计结果，实难令人信服该假设的合理性。
其次，\(Durbin-Wu-Hausman\) 检验仅对比了固定效应与最基础版本的随机效应模型，这未能充分发挥随机效应的优势。当前主流研究采用的随机效应模型已发展出更精妙的变体（如多层模型、部分 \(pooling\) 技术等），能有效规避前述强假设问题。欲了解这些前沿方法，请参阅”专业实践指南”章节中的”高级随机效应技术”部分。
16.3 专业实践方法
16.3.1 聚类标准误
回归模型的基本假设之一是误差项相互独立，即结果变量的数据生成过程中未被模型捕捉的部分本质上是随机的。该假设是正确计算回归标准误的必要前提。
不过，我们不难想象这一假设在固定效应模型中往往难以成立。毕竟，我们获取的是来自同一主体（个体或群体）的多次观测数据。可以预见的是，数据生成过程中那些未被模型捕捉的遗漏因素，很可能在该主体所有观测值间存在共性，从而导致误差项相互关联，并造成标准误估计的偏差。
因此，经济学家们普遍建议：当使用固定效应模型时，几乎总是应该采用聚类标准误28——且聚类层级需与固定效应层级保持一致。例如：
个体固定效应模型 → 使用个体层级的聚类标准误
城市固定效应模型 → 使用城市层级的聚类标准误
聚类标准误的计算方法允许误差项在特定层级内存在相关性，从而更准确地反映数据真实结构。
然而这一常规认知存在过度推广之嫌。要使”固定效应+聚类标准误”的组合成为必要（且有效）选择，需满足若干前提条件：首要条件是存在处理效应异质性——即不同个体的处理效应必须存在显著差异。
若上述条件成立，则需满足第二项条件：要么您数据中的固定效应组别/个体属于总体的非随机样本29，要么在固定效应组别/个体内部，处理变量的分配呈现聚类特征30。
因此，在采用聚类标准误前，请审慎评估这两个条件是否可能成立。若条件满足，则放心使用聚类调整；反之则无需勉强——因为不当的聚类操作反而会导致标准误被人为放大。
Alberto Abadie, Susan Athey, Guido W Imbens, and Jeffrey M Wooldridge. When should you adjust standard errors for clustering? Technical report, NBER, 2017. NBER Working Paper No. 24003.聚类标准误通常已内置在固定效应模型的执行命令中——鉴于二者在实际研究中的高频联用。让我们再次以 \(Gapminder\) 数据集为例进行演示：
library(tidyverse)
library(modelsummary)
library(fixest)
gm <- causaldata::gapminder
# feols clusters by the first fixed effect by default
clfe <- feols(lifeExp ~ log(gdpPercap) | country, data = gm)
msummary(clfe, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| log(gdpPercap) | 9.769*** |
| (0.702) | |
| Num.Obs. | 1704 |
| R2 | 0.846 |
| R2 Adj. | 0.832 |
| R2 Within | 0.410 |
| R2 Within Adj. | 0.410 |
| AIC | 10647.6 |
| BIC | 11425.6 |
| RMSE | 5.06 |
| Std.Errors | by: country |
| FE: country | X |
import pandas as pd
import numpy as np
import linearmodels as lm
from causaldata import gapminder
gm = gapminder.load_pandas().data
gm['logGDPpercap'] = np.log(gm['gdpPercap'])
# Set our individual and time (index) for our data
gm = gm.set_index(['country','year'])
mod = lm.PanelOLS.from_formula('''lifeExp ~ logGDPpercap + EntityEffects''',gm)
# Specify clustering when we fit the model
clfe = mod.fit(cov_type = 'clustered', cluster_entity = True)
print(clfe)## PanelOLS Estimation Summary
## ================================================================================
## Dep. Variable: lifeExp R-squared: 0.4101
## Estimator: PanelOLS R-squared (Between): 0.8786
## No. Observations: 1704 R-squared (Within): 0.4101
## Date: 周六, 8月 23 2025 R-squared (Overall): 0.8731
## Time: 23:50:27 Log-likelihood -5180.8
## Cov. Estimator: Clustered
## F-statistic: 1085.3
## Entities: 142 P-value 0.0000
## Avg Obs: 12.000 Distribution: F(1,1561)
## Min Obs: 12.000
## Max Obs: 12.000 F-statistic (robust): 195.30
## P-value 0.0000
## Time periods: 12 Distribution: F(1,1561)
## Avg Obs: 142.00
## Min Obs: 142.00
## Max Obs: 142.00
##
## Parameter Estimates
## ================================================================================
## Parameter Std. Err. T-stat P-value Lower CI Upper CI
## --------------------------------------------------------------------------------
## logGDPpercap 9.7690 0.6990 13.975 0.0000 8.3978 11.140
## ================================================================================
##
## F-test for Poolability: 14.000
## P-value: 0.0000
## Distribution: F(141,1561)
##
## Included effects: Entity16.3.2 非线性模型中的固定效应
固定效应研究设计背后的理论直觉具有普适性，可推广至各类统计模型。其核心思路是：通过消除组间变异\(（between \space variation）\)来阻断潜在的后门路径 \(（back \space door）\)，从而达成因果推断的可靠性。
然而，现有消除组间变异的方法——无论是通过添加二元控制变量表征固定效应，还是通过减去个体均值——均无法普适于所有统计模型。这两种方法均基于线性模型假设，当面对非线性模型（如\(probit\)、\(logit\)、\(有序logit\)、\(possion\)、\(tobit\)等）时便会产生技术障碍。
这究竟会引发什么问题？
让我们首先审视引入二元控制变量的方法：当固定效应分类数量有限时（例如区分左右利手或眼睛颜色的固定效应），添加对应的二元控制变量尚属可行，因为这些特征的变种数量较少。然而，当涉及大量组别时——比如在我们的\(Gapminder\)数据分析中包含\(142\)个不同国家——问题便随之产生：这需要引入\(142\)个控制变量并估计\(142\)个对应系数。虽然线性模型尚能应对这种需求，但在非线性模型中，研究者将面临”伴随参数问题”\(（incidental \space parameters \space problem）\)——模型根本无法有效处理如此大量的参数估计，最终导致估计结果变得极不可靠且偏差显著。
其次，考察通过减去组内均值的方法：在非线性模型中，由于结果变量已非连续变量，均值不再能有效表征”需要消除的组间变异量”。虽然理论逻辑依然成立，但均值消减法已不再适用。
当因变量存在非线性特征时，若仍希望采用固定效应设计，可考虑以下解决方案：
实践中主要存在两种解决方案：其一是直接忽略非线性特征，对二元或有序等非线性因变量进行线性建模（配合异方差稳健标准误）。这种做法在经济学界尤为普遍，研究者认为：虽然线性模型存在设定偏误，但相较于在固定效应框架下估计真实非线性模型所面临的计量难题，前者产生的偏差反而更具可控性31。
另一种常见方法是使用专为妥善处理固定效应而设计的非线性模型变体。到这一步，事情可能就不由你掌控，而落入计量经济学家的范畴了。去寻找一种旨在实现你所需功能的方法（或软件命令）。这种方法可能并不存在！但要是存在的话，一定要阅读文档，好让自己明白在做什么，然后再加以运用 。
要深入探讨的话，相关方法数量太多，这里无法一一涉及。但当你遇到这种情况时，或许就得上网搜索，查阅与你的具体应用场景相关的资料。同时要记住，通常有多种不同的方法可达成目的，且每种方法的含义都有所不同 。
不过，此处最常见的应用场景是针对二元因变量，且希望通过概率单位模型\(（probit）\)或逻辑回归模型\(（logit）\)处理固定效应。在 \(R\) 语言中，可查看 \(fixest\) 包中的广义线性模型 \(（generalized \space linear \space model）\)相关选项；在 \(Stata\) 中，可使用 \(xtprobit\)、\(xtlogit\) 或 \(xtPoisson\) 命令；在 \(Python\)中，有 \(pylogit\) 包可用于执行条件逻辑回归\(（conditional \space logit \space regression）\)—— 这种方法与之类似，但并非完全相同。
16.3.3 高级随机效应方法
在前文关于随机效应的讨论中已明确指出：除非能确证解释变量与个体效应无关（例如通过随机化实验设计），否则固定效应模型通常应作为首选方法。然而，这一论断实则基于对随机效应的”稻草人谬误”式简化——大量研究者在明知个体效应几乎必然与解释变量相关的情况下，仍坚持采用随机效应模型。这种现象显然无法仅用”研究者忽视该问题”来解释！
事实上，现代随机效应方法已能有效处理这一问题。其核心思路是将模型构建为多层次结构——这与固定效应模型的数据层级特性天然契合：任何固定效应分析场景中，数据本就具备层次性（至少包含个体内的时间维度观测）。既然数据存在层级结构，我们就应当通过建模显性捕捉这一特征。
多层次模型：指允许系数在样本中发生变异，并可通过其他变量解释其变异模式的建模框架
例如，我们可以构建如下模型：
\[Y_{it}=\beta_i+\beta_1 X_{it}+\epsilon_{it} \quad\quad (16.5)\]
该式构成了我们的第一层模型。随后可设定\(\beta_i\)自身的生成方程：
\[\beta_i=\beta_0+\gamma_1 Z_{i}+\mu_i
\quad\quad(16.6)\]
其中：
\(Z_i\)：个体层面的特征变量集，用于解释个体效应差异
\(\mu_i\)：个体层级的随机误差项
这使我们能够处理\(\beta_i\)和\(X_{it}\)之间的相关性，因为我们在直接对其进行建模。\(X_{it}\)中可能与\(\beta_i\)相关的、个体层面上具有时间一致性的那部分 —— 希望我们能用\(Z_i\)解释那部分内容。
这种方法为我们带来了一些好处。
首先，它赋予我们随机效应相对于固定效应所具备的良好统计性质，同时又无需施加过多额外假设。
其次，它让我们能够探究个体内部恒定变量的影响。假设我们想要研究家乡对一个人是否成为发明家的影响。毕竟，每个人的家乡在其一生中都不会改变，所以固定效应永远无法让我们研究这个问题。但多层随机效应\(（multi-level \space random \space effects）\)可以。
第三，它让我们能够明确且分别地探究组间效应\(（between \space effects）\)和组内效应 \(（within \space effects）\)。固定效应的意思是”只关注组内效应”。基本的随机效应是”组间效应和组内效应的某种混合” 。而多层随机效应，如果运用得当，会表示”这是组间效应，这是组内效应，你想怎么用就怎么用” 。
牢记这种多层\(（multi-level）\)直觉，我们有两个研究方向。
其中一种方法有助于明确区分组间效应和组内效应，那就是在回归中同时纳入这两种效应。就这么简单32！
Mundlak, Yair. 1978. “On the Pooling of Time Series and Cross Section Data.” Econometrica 46 (1): 69–85.基于多层次建模思想，我们可朝两个方向拓展：
方向一：显性分离组间与组内效应，方法极为简明——直接在回归模型中同时纳入两类效应即可33。
你要做的步骤如下：
计算每个预测变量的个体内均值（\(\bar{X}_i\) ）。
计算每个预测变量的组内变异（\(X_{it} - \bar{X}_i\) ） 。
将因变量对个体内均值（\(\bar{X}_i\) ）、组内变异（\(X_{it} - \bar{X}_i\) ），以及你想用的其他个体层面控制变量（\(Z_i\) ）进行回归，并对截距项使用随机效应。
这被称为 “相关随机效应\(（correlated \space random \space effects）\)”。它的便利之处在于，能为你区分出 “组间变异效应（个体内均值\(\bar{X}_i\) 对应的系数 ）” 和 “组内变异效应（\(X_{it} - \bar{X}_i\) 对应的系数 ）” ，还能让你看到个体层面变量\(Z_i\) 的效应。反正你本来就想要固定效应中的组内效应，现在还能额外获得一些信息 。
这种方法易于实施且便于解释。要是你能计算出个体内均值和组内变异（代码可参考 “如何实施” 部分 ），那么只需把所有内容纳入带有随机效应的模型即可（在 \(R\) 语言中，可查看 \(lme4\) 包的 \(lmer\) 函数；在 \(Stata\) 中，可查看带 \(re\) 选项的 \(xtreg\) 命令；在 \(Python\) 中，可查看 \(linearmodels\) 包的 \(RandomEffecs\) 类 ）。
基于多层直觉，我们的另一个研究方向是全面投入混合模型\(（mixed \space models）\)，也叫分层线性模型\(（hierarchical \space linear \space models）\) 。混合效应模型的核心突破在于：不仅允许截距项服从随机分布，更可进一步释放其他变量斜率系数的随机变异特性——为何止步于截距？何不让所有关键参数都获得这种灵活性？
\[Y_{it}=\beta_{0i}+\beta_{1i}X_{it}+\epsilon_{it}
\quad\quad(16.7)\]
\[\beta_{0i}=\beta_0+\gamma_1Z_i+\mu_i
\quad\quad(16.8)\]
\[\beta_{1_i}=\beta_1+\delta_iZ_i+\eta_i
\quad\quad(16.9)\]
虽然混合效应模型可能变得相当复杂，但其优势在于能够刻画变量间错综复杂的交互作用——毕竟在社会科学研究中，所有现象本质上都是多层次的！我们完全有理由假设：存在诸多变量不仅直接影响\(Y_{it}\)，更会调节\(X_{it}\)对\(Y_{it}\)的作用效应，而这类复杂机制正是模型需要捕捉的重点。
混合效应模型在诸多学科领域应用广泛，若您来自这些领域，可能会惊讶于本书仅用一个小节来讨论该方法。然而，鉴于本书聚焦研究设计而非参数估计，先探讨固定效应再延伸至混合模型的编排逻辑是合理的——毕竟前者构成了后者的方法论基础34。
不过，它们值得在本书范畴之外深入探究。我推荐\(Gelman\) 和 \(Hill\) \(（2006）\)的教材以供进一步阅读。你也可以自己动手尝试混合模型：在 \(R\) 中用 \(lme4\)包的\(lmer\)；在 \(Stata\) 中用\(xtmixed\)；或在 \(Python\) 的 \(statsmodels\) 包中用 \(mixedlm\) 。更广义地说，还有 \(Stan\)，它算是一种专门用于快速、灵活估计分层模型的独立语言35。\(Stan\)可与其他语言配合使用，命名直观的 \(RStan\)、\(StataStan\)和 \(PyStan\)包均有提供 。
Gelman, Andrew, and Jennifer Hill. 2006. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.若用更技术性的表述，我们实际上是通过控制层级数据\(（hierarchical \space data）\)中更高阶的分类变量来实现因果识别。 这种方法的计量经济学本质是：当我们在模型中纳入上层分组变量（如行业代码、区域编号等）作为固定效应时，等同于自动控制了所有在该层级内保持恒定的潜在混杂因素——无论这些因素是否被实际观测或测量。↩︎
不同学科对”固定效应”的术语理解确实存在显著差异——这完全不是您的错，纯粹是学科惯例使然。 在某些领域（如计量经济学），固定效应特指允许截距项存在更大变异的模型设定（即本章讨论的情形）；而在其他学科（如混合效应建模），该术语仅作为”非随机效应”的代称（即截距无个体变异的基本模型）。这种术语分歧虽令人困扰，但只需明确当前语境即可规避误解。↩︎
我们可以像控制任何分类变量那样“控制城镇因素”。在回归分析中，这意味着添加一组二元指标，每个城镇对应一个指标。在匹配法中，这意味着对城镇进行精确匹配。↩︎
几乎保持不变的变量也会被近乎完全地控制住。例如，对大多数人来说，性别在个体层面是恒定不变的。有些人生理性别未变但更改了心理认同的性别，所以性别在个体层面并非绝对恒定。但只要这种情况非常罕见，我们仍然可以说个体固定效应控制了性别因素。↩︎
“你所拥有的数据涵盖的时间段”很关键！你的数据所覆盖的时间跨度越长，某些看似固定的变量（如地理位置）实际上并非固定不变的可能性就越大。例如，德国的地理边界在很长一段时间内是恒定的，但在20世纪30年代和40年代确实有很大的变化。在考虑一个变量是否随时间保持恒定时，实际上更重要的是要考虑它在你数据所涵盖的时间段内是否恒定。↩︎
这种逻辑让一些人得出这样的结论：如果没有某种面板数据，即随时间进行的重复测量，以便利用组内变异，就无法识别出任何效应。但事实并非如此！本书中有很多方法在没有面板数据的情况下也能很好地发挥作用。这就是我们如此努力地绘制因果图的原因——这样我们就能知道是否需要像固定效应这类方法了。↩︎
从技术上讲，这些变量的效应也必须随时间保持固定，固定效应模型才能成立。↩︎
向加利福尼亚的阿克塔\(Arcata\)问好——不是阿卡迪亚\(（Arcadia）\)，这两个名字确实容易搞混。↩︎
这个问题可以通过随机效应模型来解决——具体请参见本章后续内容。↩︎
下次初次约会时，不妨试试用这个开场白。↩︎
不过，单从图形本身来看，这个（正相关）发现的说服力相当有限。↩︎
变截距模型的定义如下：该回归模型要求所有自变量的斜率系数（即解释变量的回归系数）在各个体间必须保持一致，但允许截距项（即常数项）随不同个体而变化。这种设定既能捕捉个体间的异质性，又能保证核心变量影响效应的同质性。↩︎
其中某些方法甚至无法完全分离出纯个体内变异——这简直令人难以置信！↩︎
当模型涉及多组固定效应时（例如同时控制个体效应和时间效应），这种方法会变得异常复杂。↩︎
使用这种方法时，\(\beta_0\)通常会被省略，因为它必然等于零。想想为什么：↩︎
为什么呢？想象一下，只有两个人 —— 我和你，要估计身高。假设我随时间的平均身高是 \(66\) 英寸，你的是 \(68\) 英寸，而且我们不排除任何人，所以我们有模型\(\text{Height}_{it} = \beta_0 + \beta_1\text{Me}_i + \beta_2\text{You}_i\) 。当\(\beta_0 = 0\)、\(\beta_1 = 66\)、\(\beta_2 = 68\) 时，与\(\beta_0 = 1\)、\(\beta_1 = 65\)、\(\beta_2 = 67\) 时拟合效果完全相同。有无数种方式能得到完全一样的拟合结果，模型没办法做出选择！所以这个模型无法被估计 。↩︎
这些代码片段旨在说明两种方法是等效的。大多数情况下，即使研究人员只有一组固定效应，他们实际上也会使用下面“多重固定效应集”部分展示的命令之一。↩︎
正如你在第13章中会回想起来的那样，人均GDP对数增加1个单位意味着人均GDP出现了一定的百分比增长。具体来说，是增长了171%。但对于较小幅度的增长而言，增长的数值和对应的百分比十分接近，所以我们可能更常说，人均GDP对数比典型水平高出0.1，即人均GDP大致比典型水平高出10%，这与预期寿命比典型水平高出0.9769年相关联。↩︎
\(R^2 = 1 - \text{var}(\text{residuals}) / \text{var}(\text{dependent variable})\)，其中 \(\text{var}\) 是方差的缩写 ↩︎
我们不能孤立地看待这些固定效应，因为它们都是相对于那个未被赋予独立系数的参照国来定义的——还记得这点吗？因此巴西的6.318并不意味着它的截距就是6.318，而是表示其截距比被省略的参照国高出6.318。↩︎
此外，个体固定效应估计值是基于相对较少的观测数据得出的——仅取决于每个个体拥有的数据量。因此这些估计值往往波动较大且不够稳定。本章后文将讨论的”随机效应”方法，正是通过”收缩”个体效应来解决这一问题。↩︎
吉本斯、塞拉托与乌尔班西克（2019）的研究提出了两种解决方案，其中较简单的方法是：先计算个体内部处理变量的方差，再以其倒数作为权重。这类似于第13章讨论的逆方差加权法。↩︎
包含个体与年份双重固定效应的模型↩︎
关键区别在于：这不同于”特定个体在特定年份”的情况，因为我们仅能观测到该个体在该年份的单一数据点（不存在变异）。每个”相对预期”的计算都是独立进行的。↩︎
其原理在于：通过假定所有\(\beta_i\)源自同一概率分布，模型能够整合全样本数据来估计该分布参数。这意味着每个\(\beta_i\)的估计都能利用更丰富的信息量——具体而言：
1）全局信息共享：所有个体截距共享分布参数(μ,σ²)的先验信息
2）估计效率提升：通过收缩估计减少小样本个体的估计方差
3）数据利用优化：同时吸收组内与组间的有效信息↩︎至少就这种基础形式的随机效应模型而言（更深入的分析请继续阅读后续内容）↩︎
我不认同这一检验方案，实在难以认同——杰瑞·豪斯曼↩︎
有关聚类标准误的基础概念与应用方法，请参阅第13章与第15章的专题介绍。↩︎
换言之，某些组别被纳入样本的概率显著高于其他组别。↩︎
以城市固定效应为例：该城市中某些个体是否比其他个体更有可能接受处理？↩︎
遗憾的是，这一方法选择本质上属于学术偏好问题——在此并不存在唯一正确的标准答案。↩︎
这一方法源于（但并非与…… 完全相同）蒙德拉克\(（Mundlak）\)在 \(1978\) 年提出的理论 。↩︎
该方法源自Mundlak (1978)的研究框架，虽与之并不完全相同。↩︎
我几乎能感受到Andrew Gelman教授正穿越时空对此表述提出异议——或许将混合模型置于此处仅在我这个”经济学外行”眼中才显得合理。谨此向Andrew致歉。↩︎