【应用计量系列103】线性回归、标准误与stata案例

原文信息：Fernando Frios-Avila，2023："Linear Regressions, OLS and Standard Errors"

一、引言

线性回归（LR）是经验研究中最常用的工具。我们有许多方法可以估计感兴趣的参数，最常用的是普通最小二乘（OLS）。在一些假设之下，当保持其它因素不变时，OLS可以得到因变量X对结果变量Y的无偏、一致估计量。OLS的优势在于易于实施。

考察下列线性回归模型（矩阵形式），且所有经典假设均成立：
$y =X \beta + e$
那么，OLS估计量为：
$\hat \beta_{OLS}=(X'X)^{-1}X'y$

我们也可以推导出估计量 $\hat \beta_{OLS}$ 与真实 $\beta$ 值之间的关系：
$\begin{aligned} \hat \beta_{OLS}&=(X'X)^{-1}X'(X\beta + e) \\ &=(X'X)^{-1}X'X\beta + (X'X)^{-1}X'e \\ &= \beta + (X'X)^{-1}X'e \end{aligned}$

也就是说，当 $E((X'X)^{-1}X'e)=0$ ，即感兴趣的因变量与扰动项无关时，OLS估计量是无偏的。但是，当使用有限样本时，OLS估计量与真实值还是存在一些差异，这主要是由于抽样误差导致。

二、估计量有多精确

利用OLS估计量与真实值之间的关系，我们还可以推导出估计量的精度。在上式两边同时取方差Var：
$\begin{aligned} Var(\hat\beta_{OLS})=Var(\beta + (X'X)^{-1}X'e) \\ Var(\hat\beta_{OLS})=Var(\beta) +Var((X'X)^{-1}X'e) \end{aligned}$

假设 $X$ 是固定的，我们可以重写上述方程为方差-协方差矩阵的形式：
$Var(\hat\beta_{OLS})=Var(\beta) +(X'X)^{-1}X' Var(e) X (X'X)^{-1}$
其中，
$Var(e) = \left( \begin{matrix} \sigma^2_1 & \sigma_{1,2} & \sigma_{1,3} & ... & \sigma_{1,n} \\ \sigma_{1,2} & \sigma_2^2 & \sigma_{2,3} & ... & \sigma_{2,n} \\ \sigma_{1,3} & \sigma_{2,3} & \sigma^2_{3} & ... & \sigma_{3,n} \\ ... & ... & ... & ...& ... \\ \sigma_{1,n} & \sigma_{n,2} & \sigma_{n,3} & ... & \sigma^2_{n} \\ \end{matrix}\right)$

在截面数据中，我们仅仅只能观测到一个个体的单次样本。但理论上来说，我们所观测到的样本仅仅时候个体可能得所有状态的一种实现值。这意味着，在每一个可能得状态，个体都可能有不同的 $e_i$ 。因此，如果我们能观测到所有的状态，就有可能获得个体不可观测的方差 $\sigma_i^2$ ，以及与其它个体的协方差 $\sigma_{ij}$ 。

我们可以将上式转换成我们常见的方差协方差公式：
$Var(\hat\beta_{OLS}) = \color{brown}{(X'X)^{-1}} \color{green}{X'} \color{red} \Omega \color{green}{X}\color{brown}{(X'X)^{-1}}$
上述公式也成为“三明治”公式：

$\color{brown}{(X'X)^{-1}}$ ：是面包；
$\color{green}{X'}$ ：生菜
$\color{red} \Omega$ ：它是三明治最好吃的部分——肉类。

下面，我们来看看不同类型的标准误。为了简化阐述，忽略自由度的矫正。

（一）同方差标准误

我们所学的第一种类型的标准误就是同方差假设，且个体间不相关。这意味着：

影响个体i的不可观测因素 $e$ 完全与影响个体j的不可观测因素无关。
不可观测因素的方差在所有个体间是相同的。

数学上，这意味着：
$\begin{aligned} \sigma_i^2 &=\sigma_j^2 = \sigma^2 \\ \sigma_{ij} &=0 \ \forall \ i\neq j \end{aligned}$
这极大地简化了矩阵 $\Omega$ :
$\Omega_0 = \left( \begin{matrix} \sigma^2 & 0 & 0 & ... & 0 \\ 0 & \sigma^2 & 0 & ... & 0 \\ 0 & 0 & \sigma^2 & ... & 0 \\ ... & ... & ... & ...& ... \\ 0 & 0 & 0 & ... & \sigma^2 \\ \end{matrix}\right)= \sigma^2 I(N)$
也极大地简化了估计量的方差-协方差：
$Var(\hat\beta_{OLS})_0 = (X'X)^{-1} X'\sigma^2 I(N) X (X'X)^{-1}=\sigma^2 (X'X)^{-1}$

（二）稳健标准误

稳健标准误（Huber-White稳健标准误）放松了一个假设：不同个体的方差不同。但个体间仍不相关。
$\begin{aligned} \sigma_{ij} &=0 \ \forall \ i\neq j \\ \exists i,j &:\sigma^2_{i}\neq \sigma^2_{j} \end{aligned}$
这并非意味着所有两个个体总是有不同的方差，仅仅是可能不同。在这个条件下，我们能将 $\Omega$ 矩阵表示成：
$\Omega_1 = \left( \begin{matrix} \sigma_1^2 & 0 & 0 & ... & 0 \\ 0 & \sigma_2^2 & 0 & ... & 0 \\ 0 & 0 & \sigma_3^2 & ... & 0 \\ ... & ... & ... & ...& ... \\ 0 & 0 & 0 & ... & \sigma_n^2 \\ \end{matrix}\right)$

（三）一维聚类标准误

聚类标准误放松了另一个假设——不可观测因素在个体间并不必然独立。例如，同一家庭内的人可能有相同的经历，因此，他们的不可观测因素是相关的。同一教室的学生也是一样，他们有相同的学习经历。同理，相同企业的职工有相同的工作环境等等。

当然，我们仍然需要对数据施加一些限制：

首先，个体彼此相关仅仅在一个维度上；
第二，如果个体不属于相同的群体，他们的标准误仍然是独立的。

为了用数学表示，定义函数 $g()$ 表示个体i的群体成员。换言之，如果个体i属于组群 $k$ ，那么， $g(i)=k$ 。利用这个函数，我们可以顶级个体之间的相关性：
$\begin{aligned} \sigma_{i,j} \neq 0 \text{ if } g(i)=g(j) \end{aligned}$
在这种情形下， $\Omega$ 矩阵看起来像稳健标准误的形式。我们可以将其分块表示，块的元素是主对角线非零元素，其它元素为0：
$\Omega_2 = \left( \begin{matrix} \sigma_1^2 & \sigma_{1,2} & 0 & ... & 0 \\ \sigma_{1,2} & \sigma_2^2 & 0 & ... & 0 \\ 0 & 0 & \sigma_3^2 & ... & \sigma_{3,n} \\ ... & ... & ... & ... & ... \\ 0 & 0 & \sigma_{3,n} & ... & \sigma_n^2 \\ \end{matrix}\right) = \left( \begin{matrix} \Sigma_1 & 0 \\ 0 & \Sigma_2 \\ \end{matrix}\right)$
stata可视化 $\Omega$ 矩阵：

clear
set scheme white2
color_style tableau
set seed 1
set obs 50
gen r1=runiformint(1,4)
gen r2=runiformint(1,4)
gen id=_n
sort r1  r2
qui:mata:
r1=st_data(.,"r1")
r2=st_data(.,"r2")
rr1=J(rows(r1)*rows(r2),4,0)
k=0
for(i=1;i<=50;i++){
    for(j=1;j<=50;j++){
        if ((r1[i]==r1[j]) | (r2[i]==r2[j])) {
            k++
            rr1[k,]=(51-i,j,(r1[i]==r1[j]),(r2[i]==r2[j]) )         
        }
    }   
}
rr1=rr1[1..k,]
end
getmata rr1*=rr1, replace force

two (scatter rr11 rr12 if rr13==1,  ms(s) msize(2.1))  ///
    (scatter rr11 rr12 if 51-rr11 == rr12, ms(s) msize(2.1) color(gs1)  ) ///
    , aspect(1) legend(off)  xtitle("") ytitle("") yscale(off) xscale(off)

image.png

在上图中，蓝色的块显示，我们允许一些相同组群（聚类）的观测样本之间非零的相关性。隐含地，假设所有主对角线元素（红色）也是非零的，这就是为什么聚类标准误对异方差也是稳健的原因。

（四）多维聚类标准误

正如上文提到的，一维聚类标准误允许不同个体不可观测误差是相关的，但仅仅只在他们属于同一组群。如果他们来自不同组群，那么，仍然不相关。一维聚类标准误也假设同一组群的个体仅仅通过一维渠道相互联系。

后面这个假设是非常强的。例如，人们可能相互之间产生联系。在学术圈，我可能不认识你，但我们可能有一些共同的朋友或老师或学生。

好消息是，多维聚类标准误就是应对这种情形。尤其是，两个个体至少有一个共同的组群，那么，我们会假设个体间不可观测因素的相关性不等于0。

为了从数学上展示，我们定义函数 $g_h()$ 表示基于聚类变量h，个体所属的组群。如果i和j没有任何联系， $gg(i,j)$ 为0，否则为1：
$\begin{aligned} gg(i,j)&=0 \text{ if } \forall h:g_h(i)\neq g_h(j) \\ &\text{ and 1 otherwise} \end{aligned}$

这就意味着， $\Omega$ 矩阵不再是对角块，因为主对角线外的元素也不等于0：
$\Omega_3 = \left( \begin{matrix} \sigma_1^2 & \sigma_{1,2} & 0 & ... & \sigma_{1,n} \\ \sigma_{1,2} & \sigma_2^2 & \sigma_{2,3} & ... & 0 \\ 0 & \sigma_{2,3} & \sigma_3^2 & ... & \sigma_{3,n} \\ ... & ... & ... & ... & ... \\ \sigma_{1,n} & 0 & \sigma_{3,n} & ... & \sigma_n^2 \\ \end{matrix}\right)$
stata可视化：

* 一维聚类
two (scatter rr11 rr12 if rr13==1,  ms(s) msize(2.1))  ///
    , aspect(1) legend(off)  xtitle("") ytitle("") yscale(off) xscale(off) name(m1, replace) 


* 一维聚类
two (scatter rr11 rr12 if rr14==1,  ms(s) msize(2.1))  ///
    , aspect(1) legend(off)  xtitle("") ytitle("") yscale(off) xscale(off) name(m2, replace)    


* 二维聚类
two (scatter rr11 rr12 if rr14==1 | rr13==1,  ms(s) msize(2.1))  ///
    , aspect(1) legend(off)  xtitle("") ytitle("") yscale(off) xscale(off) name(m3, replace)

第一个一维聚类

第二个一维聚类

二维聚类

上图显示了二维聚类（两个聚类变量）。前两个一维聚类是分别基于单一聚类变量的聚类标准误。而第三个图是同时考虑两个聚类变量：蓝色的区域是个体间相关性不等于0。

四、如何估计方差协方差？

那么我们如何估计 $\sigma_i^2$ 和 $\sigma_{i,j}$ ？

因为我们只能在截面上观测到一次个体，最佳估计就是
$\hat \sigma_{ij} = \varepsilon_i*\varepsilon_j$

其中， $\varepsilon_i$ 是个体的不可观测部分。只要我们有了 $\hat \sigma_{ij}$ ，我们就只需要定义我们想要使用的标准误假设（同方差、异方差、聚类）。但是在实践应用中，大多数软件仍然有两个限制：

① 上述阐述完全忽略了自由度的作用和对应的矫正。
② 上述内容 $\Omega$ 矩阵计算适合于小样本集。

五、Stata 18更新了线性回归中的标准误命令

在Stata的线性回归命令中，选项vce(vcetype)声明标准误类型，它包括渐进理论得到的标准误（ols），异方差稳健标准误（robust），聚类标准误（cluster clustervallist），bootstrap或jackknife的标准误（bootstrap，jackknife）。

vce(ols)是默认选项，ols回归的标准方差估计量；
vce(robust)，无聚类情形，它用 $\hat{\sigma}_j^2=\frac{n}{n-k} u_j^2$ 作为第j个观测值方差估计量。其中，n表示观测量，k表示解释变量数量， $u_j$ 表示计算的余值，包括 $\frac{n}{n-k}$ 是为了改善总估计量的小样本性质；
vce(cluster clustvarlist)声明聚类变量(clustvarlist)层面的聚类标准误；
vce(hc2 [clustvar] [,dfadjust])用 $\frac{u_j^2}{1-h_{jj}}$ 作为观测值的方差估计量，其中 $h_{jj}$ 是矩阵的对角元素。如果模型是同方差的，那么，这个估计量是无偏的。vce(hc2)趋向于得到更加保守的置信区间。clustvar允许产生多个聚类变量的聚类标准误。dfadjust计算Bell and McCaffrey(2002)的调整自由度。
vce(hc3)用 $\frac{u_j^2}{(1-h_{jj})^2}$ 作为方差估计量，这是由Davidson and MacKinnon (1993)提出的标准误，作者认为当存在异方差时，这个标准误更好。且它产生更加保守的置信区间。

下面，我们来看看stata代码：

*加载数据
sysuse auto, clear

generate gpmw = ((1/mpg)/weight)*100*1000
summarize gpmw

*(1)同方差标准误
regress gpmw foreign

同方差标准误

下面，用vce（robust）来展示一下异方差稳健标准误：

regress gpmw foreign, vce(robust)

异方差稳健标准误

我们能看到，同方差和异方差稳健标准误下，点估计量都是相同的（外国汽车的耗油量是0.25），但是标准误相差20%（0.055 vs 0.068）。传统回归报告的95%置信区间为[0.14, 0.36]，而在稳健标准误的置信区间[0.11, 0.38]。

那我们应该相信哪个？注意，耗油量可能存在很大的异方差：

tabulate foreign, summarize(gpmw)

结果变量的方差

所以稳健标准误的结果可能更加可信。

下面，我们来看看其他的一些稳健标准误，用vce（hc2）和vce（hc3）：

regress gpmw foreign, vce(hc2)

regress gpmw foreign, vce(hc3)

hc2稳健标准误

hc3稳健标准误

从上述结果可以看出，hc2和hc3稳健标准误得到更加保守的置信区间。

下面，我们用年龄-工资收入的数据来看看聚类标准误：

use https://www.stata-press.com/data/r18/regsmpl, clear

regress ln_wage age c.age#c.age tenure

无聚类标准误

我们有理由相信无聚类标准误的结果没有意义。因为更高工资的女性可能在其它年份也有更高的工资，因此，余值并不是独立的。这个时候用聚类稳健标准误更好。

regress ln_wage age c.age#c.age tenure, vce(cluster id)

聚类标准误

对比来看，关注tenure的系数——继续工作的回报率。在无聚类下，置信区间为[0.038, 0.041]。而聚类稳健标准误下置信区间更大[0.036, 0.042]。我们还可以用vce(hc2 id):

regress ln_wage age c.age#c.age tenure, vce(hc2 id)

hc2聚类标准误

下面，我们来看看多维聚类标准误——个体和教育水平：

regress ln_wage age c.age#c.age tenure, vce(cluster idcode grade) clustertable

多维聚类标准误

选项clustertable可以让结果包含描述联合聚类的表，也报告每个聚类水平的数量。上述结果聚类两个变量。t统计量的自由度是从联合聚类中来的，选择最小的聚类水平数。在上述例子中个，自由度是19-1=18。wald检验F统计量的分母包含自由度18。现在tenure系数的95%的置信区间变得更大[0.035, 0.043]。

我在我的课程中详细地讲过：