多元回归系数：我们都解释错了？

发布时间：2020-06-04 阅读 8098

Stata 连享会主页 || 视频 || 推文

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

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

作者：金钊 (中山大学)
E-Mail: 980510243@qq.com

编者按： 在 Stata Journal (各期 SJ) 2016 年第 1 期中连续登载了 4 篇「吵架」论文。其中，首篇为 David Hoaglin 撰写的长文 (p.5-22)，题为「Regressions are Commonly Misinterpreted」。从标题上来看，这无异于「挑战」我们的常识。三位知名的计量经济学家 (James Hardin, Scott Long, David Drukker) 撰写了两篇短文对此进行评论/批评。而同期第四篇论文刊登的就是 David Hoaglin 的「反驳 (Rejoinder)」。显然，这四篇文章是 Stata Journal 的编辑们蓄意之作，同时，也凸显出此问题的重要性。本文对其中的一些核心观点进行梳理，感兴趣的读者可以阅读原文以便品尝原味红茶。

Hoaglin David C., 2016, Regressions are Commonly Misinterpreted, Stata Journal, 16(1): 5–22. [PDF]
Hardin James W. , 2016, Regressions are Commonly Misinterpreted: Comments on the Article, Stata Journal, 16(1): 23–24. [PDF]
Long J. Scott, David M. Drukker, 2016, Regressions are Commonly Misinterpreted: Comments on the Article, Stata Journal, 16(1): 25–29. [PDF]
Hoaglin David C., 2016, Regressions are Commonly Misinterpreted: A Rejoinder, Stata Journal, 16(1): 30–36. [PDF]

同主题阅读：

1. 引言

多元回归模型一直被广泛运用，也是最常见和最基础的计量模型。多元回归模型中各个变量间关系相对复杂，其回归系数惯常解释为：

「当其他变量保持不变或控制其他变量不变时， $X_{i}$ 每改变一个单位时因变量 $Y$ 的平均变化量」。

然而，Hoaglin (2016) 指出，这种常见的解读存在错误。这些问题常常出现在 OLS 回归、logistic 回归和其他广义线性模型以及生存分析、纵向分析和层次分析回归中。

Hoaglin (2016) 认为，这些解释既没有体现多元回归的基本原理，也不符合现实情况。他从图形、多元正态分布和最小二乘几何特征等角度解释「控制」和「保持不变」的不适性。为此，他们提出一直全新的解释「调整 $X_{i}$ 和其他预测因子对 $Y$ 的共同线性影响后， $X_{i}$ 每改变一个单位时因变量 $Y$ 的平均变化量」。

2. 多元线性回归系数的常见误解

2.1 多元线性回归计量模型

我们常见的多元回归的总体 (population) 模型为:

Y_{i} = β_{1} X_{i 1} + . . . + β_{p} X_{i p} + ε_{i} (1)

其中， $X_{i k}$ 为解释变量或预测因子； $Y_{i}$ 为被解释变量或响应变量。我们通常令第一个解释变量为单位向量，即 $X_{i 1} \equiv 1$ 。 $β_{1}, . . ., β_{p}$ 为待估参数，称为「回归系数」(regression coefficients)。

在多元回归中，预测因子间不可能是完全独立的，每个回归系数的下标应该包含方程中的其他预测因子。为此，预测因子 $X_{i 2}$ 的回归系数可以写为： $β_{y 2 \cdot 13. . . p}$ ，第一个下标表示响应变量，第二个下标表示系数所附的预测因子，而在「·」之后的下标表示其他预测因子。则回归模型可以变为：

Y_{i} = β_{y 1 \cdot 2. . . p} X_{i 1} + . . . + β_{y p \cdot 1. . . p - 1} X_{i p} + ε_{i} (2)

运用数据可以对回归系数 $β_{y k \cdot 12. . . k - 1, k + 1. . . p}$ 进行估计，可以得到估计值 $b_{y k \cdot 12. . . k - 1, k + 1. . . p}$ 。则 (2) 的拟合方程为：

y_{i} = b_{y 1 \cdot 2. . . p} X_{i 1} + . . . + b_{y p \cdot 1. . . p - 1} X_{i p} + e (3)

$e$ 为残差， $y_{i}$ 为 $Y_{i}$ 的观测值。在多元回归中，每个预测因子的系数都说明了其他预测因子的贡献，也就是说，它反映了对这些预测因子的调整。

2.2 多元线性回归系数的图形解释

我们运用 Stata 自带的 1978 汽车数据集 auto.dta 中的进口汽车 (foreign) 数据来解释多元回归系数估计值的含义。

首先，我们把汽车的百英里油耗 (100/mpg) 当作被响应变量，汽车的重量 (weight) 和排量 (displacement) 为预测因子。通过散点图 (图 1)，我们可以看到汽车油耗与重量和排量的相关性很高，汽车的重量和排量的相关性也很强。

. sysuse auto, clear
(1978 Automobile Data)

. generate gp100m = 100/mpg

. label var gp100m "Gallons per 100 miles"

*-相关系数
. pwcorr gp100m weight displacement if foreign==1

             |   gp100m   weight displa~t
-------------+---------------------------
      gp100m |   1.0000 
      weight |   0.8172   1.0000 
displacement |   0.8444   0.9507   1.0000 

*-散点图矩阵
. graph matrix gp100m weight displacement if foreign==1

图 1 汽车油耗、重量和排量的散点图矩阵

首先，我们看二元回归的估计结果如下，可以发现汽车重量回归系数的估计值为 0.396，而汽车排量回归系数的估计值为 0.032。

. regress gp100m weight displacement if foreign == 1

      Source |       SS       df       MS      Number of obs  =      22
   ----------+------------------------------   F(2, 19)       =   23.86
       Model |  19.6704568     2  9.83522842   Prob > F       =  0.0000
    Residual |  7.83165119    19  .412192168   R-squared      =  0.7152
   ----------+------------------------------   Adj R-squared  =  0.6853
       Total |   27.502108    21  1.30962419   Root MSE       =  .64202

-----------------------------------------------------------------------
      gp100m |  Coef.   Std. Err.      t    P>|t|  [95% Conf. Interval]
-------------+---------------------------------------------------------
      weight |  0.396      1.044     0.38   0.708    -1.788       2.580
displacement |  0.032      0.018     1.78   0.091    -0.006       0.070
       _cons | -0.196      0.811    -0.24   0.812    -1.893       1.501
-----------------------------------------------------------------------

其次，我们单独对汽车重量进行回归，可以发现，在第 (1) 列中，汽车重量回归系数的估计值为 2.160，比第二列所呈现的二元回归的系数估计值高 (0.396)：

 . regress gp100m weight              if foreign == 1
 . est store m1
 . regress gp100m weight displacement if foreign == 1
 . est store m2

 . esttab m1 m2, nogap b(%6.3f) s(N r2_a) 
   --------------------------------------------
                         (1)             (2)   
                      gp100m          gp100m   
   --------------------------------------------
   weight              2.160***        0.396   
                      (6.34)          (0.38)   
   displacement                        0.032   
                                      (1.78)   
   _cons              -0.689          -0.196   
                     (-0.86)         (-0.24)   
   --------------------------------------------
   N                  22.000          22.000   
   r2_a                0.651           0.685   
   --------------------------------------------
   t statistics in parentheses
   * p<0.05, ** p<0.01, *** p<0.001

接着，我们用图示显示在一元回归后加入汽车重量的估计结果。

图 2 显示了，汽车油耗 (gp100m) 对排量 (displacement) 回归后的残差，与汽车重量 (weight) 对排量 (displacement) 回归后的残差之间的关系。可以发现，两类残差的相关性与二元回归中汽车重量 (weight) 的系数估计值和 $t$ 值完全一致。命令如下：

*-部分回归图
 . ssc install reganat, replace  //下载外部命令
 . reganat gp100m weight displacement if foreign == 1, dis(weight) biline
 . graph export "图2-reganat部分回归图.png", replace  //保存图片

我们也可以手动计算上述残差，并进而用这两个残差做回归，得到与二元回归中一样的系数估计值 ${\hat{β}}_{1} = 0.396$ 。在这段命令中，第 [1]-[2] 行的目的在于从 gp100m 中「滤掉」displacement 的影响 (也就是前文反复提及的「调整」)；第 [3]-[4] 行的目的相似：从 weight 中「调整掉」displacement 的影响。调整后的的两个残差 e_y_x2 和 e_x1_x2 中已经不再包含 displacement 变量的信息了，因此，我们执行 reg e_y_x2 e_x1_x2 时，产生的效果与 regress gp100m weight displacement 是完全一致的。

*-手动计算二元回归的系数
 . keep if foreign==1
 . reg gp100m displacement      // [1]
 . predict e_y_x2, res          // [2] 
 . reg weight displacement      // [3]
 . predict e_x1_x2, res         // [4]
 
 . reg e_y_x2  e_x1_x2
   Source |       SS           df       MS      Number of obs   =        22
----------+----------------------------------   F(1, 20)        =      0.15
    Model |  .059470475         1  .059470475   Prob > F        =    0.7009
 Residual |  7.83165211        20  .391582605   R-squared       =    0.0075
----------+----------------------------------   Adj R-squared   =   -0.0421
    Total |  7.89112258        21  .375767742   Root MSE        =    .62577
---------------------------------------------------------------------------
   e_y_x2 |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
----------+----------------------------------------------------------------
  e_x1_x2 |      0.396      1.017     0.39   0.701       -1.725       2.518
    _cons |      0.000      0.133     0.00   1.000       -0.278       0.278
---------------------------------------------------------------------------

编者注：有关这部分内容的详情，参见「图示线性回归系数：Frisch-Waugh定理与部分回归图」。这里使用 reganat 命令绘制的图形是编者自行添加的。

由此，我们可以清晰的知道，二元回归中汽车重量 ( $x_{1}$ ) 的系数估计值 ${\hat{β}}_{1}$ 是指汽车油耗 ( $y$ ) 和重量 ( $x_{1}$ ) 同时调整了汽车排量 ( $x_{2}$ ) 对其的线性变化后的关系。从一元和二元回归结果中汽车重量的系数来看(分别为 2.160 和 0.396)，这个调整的影响是很大的。

图 2 增加变量后的散点图

2.3 多元线性回归系数的代数和矩阵解释

我们从简单的二元回归模型来看系数的代数解释。首先，我们假设回归计量模型为:

Y_{i} = β_{0} + β_{y 1.2} X_{i 1} + β_{y 2.1} X_{i 2} + ε_{i} (4)

通过最小二乘估计，我们可以得到系数的代数表达式如下。可以发现，系数 $β_{y 1.2}$ 正是 $Y - ρ_{y 2} X_{2}$ 对 $X_{1} - ρ_{12} X_{2}$ 回归的斜率。也就是说， $Y$ 对 $X_{1}$ 的回归系数 $β_{y 1.2}$ 是 $Y$ 和 $X_{1}$ 对 $X_{2}$ 回归后的残差再回归的结果。可以理解为是调整了 $Y$ 和 $X_{1}$ 对 $X_{2}$ 的共同线性变化后， $X_{1}$ 变动一单位对 $Y$ 的平均变化。_

\begin{aligned} b_{y 1.2} & = \frac{ρ_{1 y} - ρ_{12} ρ_{2 y}}{1 - ρ_{12}^{2}} = \frac{cov (X_{1} - ρ_{12} X_{2}, Y - ρ_{y 2} X_{2})}{var (X_{1} - ρ_{12} X_{2})} \\ b_{y 2.1} & = \frac{ρ_{2 y} - ρ_{12} ρ_{1 y}}{1 - ρ_{12}^{2}} = \frac{cov (X_{2} - ρ_{12} X_{1}, Y - ρ_{y 1} X_{1})}{var (X_{1} - ρ_{12} X_{2})} \end{aligned}

对于多元回归系数，可以写成矩阵的形式， $b_{y 1 \cdot 2. . . p} = (X^{'} X)^{- 1} X^{'} Y$ 。令 $Y$ 的预测值为 $\hat{y}$ ，则 $\hat{y} = X b = X (X^{'} X)^{- 1} X^{'} Y$ ,其中 $X (X^{'} X)^{- 1} X^{'}$ 被称为「投影矩阵」，也就是说 $\hat{y}$ 是 $Y$ 在 $X$ 超平面上的投影。

2.4 常见的错误解读

在拟合方程中，预测因子的系数估计值 $b_{k}$ 不仅仅代表斜率，其包含更复杂的关系。

y_{i} = b_{1} X_{i 1} + b_{2} X_{i 2} + . . . + b_{p} X_{i p} + e (5)

常见的对 $b_{2}$ 的解释为：

「控制其他解释变量不变, $X_{2}$ 变化一单位， $Y_{i}$ 的平均变化」。

然而，这样的解释只有当 $X_{2}$ 为虚拟变量，即其值由 0 变 1 时才成立。这样运用「控制」来描述 $X_{2}$ 与其他预测因子之间的关系显然存在问题

(1) 「控制」的描述可能意味着在数据收集中对预测因子运用了随机化原则 (randomization rules)，即总体或样本中每个个体发生的概率均等。

(2) 「控制」的描述并不能反映多元回归的工作原理，即忽视了其它预测因子对回归系数的影响。「控制」的解释通常是从偏导数的定义来的。简单说， $β_{2}$ 为 $Y$ 对 $X_{2}$ 的偏导数，即 $\partial Y_{i} / \partial X_{2} = β_{2}$ 。然而，运用偏导数的概念来解释 $β_{2}$ 存在两方面缺陷，一方面，实际数据是无法观测的， $Y$ 对 $X_{2}$ 的偏导只是形式上的；另一方面，在微积分中的「保持不变」是偏导的假设，而这里是把假设当作结论。偏导数并不能解释系数在多大程度上反映其他预测因子的贡献。

(3) 在很多模型中，我们没法保证 $X_{2}$ 的变动不会引起其它预测因子的变动。可以从两个例子来进一步理解「保持不变」的说法是不合理的。第一个例子，在模型 (6) 中，我们不可能要求在 $X^{2}$ 不变的情况下仅让 $X$ 变动。第二个例子，在模型 (7) 中，在改变 $X_{1}$ 时，若想保持 $X_{1} X_{2}$ 不变，就必须让 $X_{2}$ 能够变动。

Y_{i} = β_{0} + β_{1} X + β_{2} X^{2} + ε_{i} (6)

Y_{i} = β_{0} + β_{1} X_{1} + β_{2} X_{2} + β_{3} X_{1} X_{2} + ε_{i} (7)

(4) 允许一个预测因子变化而其他预测因子固定在其平均值所获得的预测值可能没有意义。一方面，某些预测因子的均值可能缺乏实际经济意义；另一方面，在数据中预测时使用的任何特定预测因子模式都可能没法显示。

2.5 正确的解释

为此，Hoaglin (2016) 提出对于多元回归系数的「正确解释」应该为：

$b_{y 2 \cdot 13. . . p}$ 表示，调整了 $X_{i 2}$ 和其他预测因子对 $Y_{i}$ 的共同线性影响后， $X_{i 2}$ 变化一单位， $Y_{i}$ 的平均变化。

他们运用「调整」来代替「控制」，同时突出预测因子间对响应变量的共同线性影响。 Long and Drukker(2016)

3. Stata 命令：margins 运用问题

背景知识：参见「Stata：因子变量全攻略」，「Stata：边际效应分析」。

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。

对多元回归系数估计值的重新解读会对 margins 命令的运用带来重要影响。下面是 margins 命令官方说明文档中的一个例子，所使用的数据为虚构的网络数据 margex.dta：

. webuse margex, clear
(Artificial data for margins)

. tab group sex,column

+-------------------+
| Key               |
| ----------------- |
| frequency         |
| column percentage |
+-------------------+

           |          sex
     group |      male     female |     Total
-----------+----------------------+----------
         1 |       215        984 |     1,199
           |     14.35      65.51 |     39.97
-----------+----------------------+----------
         2 |       666        452 |     1,118
           |     44.46      30.09 |     37.27
-----------+----------------------+----------
         3 |       617         66 |       683
           |     41.19       4.39 |     22.77
-----------+----------------------+----------
     Total |     1,498      1,502 |     3,000
           |    100.00     100.00 |    100.00

可以看到样本在不同组别中的性别分布是截然不同的，接下来我们做因变量 $y$ 对性别 $s e x$ 和组别 $g r o u p$ 的简单回归，数据中并未说明 $y$ 的具体含义，不妨假设其代表小时工资.

. reg y i.sex i.group

      Source |       SS           df       MS      Number of obs   =     3,000
-------------+----------------------------------   F(3, 2996)      =    152.06
       Model |  183866.077         3  61288.6923   Prob > F        =    0.0000
    Residual |  1207566.93     2,996  403.059723   R-squared       =    0.1321
-------------+----------------------------------   Adj R-squared   =    0.1313
       Total |  1391433.01     2,999  463.965657   Root MSE        =    20.076

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         sex |
     female  |     18.322      0.893    20.52   0.000       16.571      20.073
             |
       group |
          2  |      8.038      0.914     8.80   0.000        6.246       9.829
          3  |     18.639      1.160    16.08   0.000       16.366      20.913
             |
       _cons |     53.321      0.935    57.06   0.000       51.489      55.154
------------------------------------------------------------------------------

在不加任何选项的情况下， margins 命令默认计算的是平均调整过的预测值 (Average adjusted predictions, AAPs) ，将样本视为每个人都是男性 (或是女性) ， margins 命令的结果如下：

. margins sex

Predictive margins                              Number of obs     =      3,000
Model VCE    : OLS

Expression   : Linear prediction, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        sex  |
       male  |     60.560      0.578   104.74   0.000       59.427      61.694
     female  |     78.882      0.577   136.65   0.000       77.751      80.014
------------------------------------------------------------------------------

从结果可以看出，两个 AAPs 之差刚好等于性别的估计系数 (78.88−60.56 = 18.32) 。然而，如果让 AAPs 有意义，就必须假定每个个体刚好以 39.97%，37.27% 和 22.77% 的概率分布在 group1、group2 和 group3。

并且，在这三个群体中，男性和女性都应该具有相同的分布。然而，样本的「预测空间」为六个点，分别对应为 ${m a l e f e m a l} \times {g r o u p 1 g r o u p 2 g r o u p 3}$ 。在交叉表中可以看出，三个组中的男性和女性的分布有显著差异。

以上分析属于线性回归，接着，我们运用 nhanes2f.dta 数据集，讨论 logistic 回归的情况。

. webuse nhanes2f, clear

. logit diabetes black female age

Iteration 0:   log likelihood = -1999.0668
Iteration 1:   log likelihood = -1841.3525
Iteration 2:   log likelihood = -1812.3671
Iteration 3:   log likelihood = -1811.9834
Iteration 4:   log likelihood = -1811.9828
Iteration 5:   log likelihood = -1811.9828

Logistic regression                             Number of obs     =     10,335
                                                LR chi2(3)        =     374.17
                                                Prob > chi2       =     0.0000
Log likelihood = -1811.9828                     Pseudo R2         =     0.0936

------------------------------------------------------------------------------
    diabetes |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       black |      0.718      0.127     5.66   0.000        0.469       0.966
      female |      0.155      0.094     1.64   0.101       -0.030       0.339
         age |      0.059      0.004    15.93   0.000        0.052       0.067
       _cons |     -6.405      0.237   -27.00   0.000       -6.870      -5.940
------------------------------------------------------------------------------

我们知道，logit 回归的系数并不代表边际效应。对于非线性模型，边际效应不是常数，而是随着解释变量而变化。

我们接下来看年龄 (age) 在 20、30、40、50、60 和 70 处，性别 (female) 和种族 (black) 在样本均值处的边际效应。

可以发现，当性别和种族处于均值时 (0.525 和 0.105) ，70 年龄组的边际效应是 20 年龄组的 18 倍 (11% 和 0.6%) 。样本的「预测空间」为四个点： $(b l c a k, f e m a l e) = (0, 0), (0, 1), (1, 0), (1, 1)$ 。要使这个结果有解释意义，我们必须假定 20 岁年龄组和 70 岁年龄组的女性占比为 0.525，黑人占比为 0.105。然而，实际上，在nhanes2f.dta 数据集中，20 岁年龄组的女性占比为 0.578，而黑人占比 0.123；70 岁年龄组的女性占比 0.5，而黑人占比 0.064。

. margins, at(age=(20 30 40 50 60 70)) atmeans

Adjusted predictions                            Number of obs     =     10,335
Model VCE    : OIM

Expression   : Pr(diabetes), predict()

1._at        : black           =    .1050798 (mean)
               female          =    .5250121 (mean)
               age             =          20

2._at        : black           =    .1050798 (mean)
               female          =    .5250121 (mean)
               age             =          30

3._at        : black           =    .1050798 (mean)
               female          =    .5250121 (mean)
               age             =          40

4._at        : black           =    .1050798 (mean)
               female          =    .5250121 (mean)
               age             =          50

5._at        : black           =    .1050798 (mean)
               female          =    .5250121 (mean)
               age             =          60

6._at        : black           =    .1050798 (mean)
               female          =    .5250121 (mean)
               age             =          70

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         _at |
          1  |      0.006      0.001     6.38   0.000        0.004       0.008
          2  |      0.011      0.001     8.25   0.000        0.009       0.014
          3  |      0.020      0.002    11.42   0.000        0.017       0.024
          4  |      0.036      0.002    16.99   0.000        0.032       0.041
          5  |      0.064      0.003    22.50   0.000        0.059       0.070
          6  |      0.110      0.006    18.82   0.000        0.099       0.122
------------------------------------------------------------------------------

样本取不同值的边际效应差别会很大。然而，一般文献里常计算样本均值处的边际效应。从下图我们可以发现，当性别和人种取不同值时，20 岁年龄组和 70 岁年龄组的边际效应差别非常大。尽管 margins 可以为很多模型的预测提供更大的灵活性，然而，在分析时不能一味的选择「样本均值处边际效应」，而是应该详细分析样本的预测区间，谨慎选取样本代表值。

. margins, at(age=(20 70) black=(0 1) female=(0 1))

Adjusted predictions                            Number of obs     =     10,335
Model VCE    : OIM

Expression   : Pr(diabetes), predict()

1._at        : black           =           0
               female          =           0
               age             =          20

2._at        : black           =           0
               female          =           0
               age             =          70

3._at        : black           =           0
               female          =           1
               age             =          20

4._at        : black           =           0
               female          =           1
               age             =          70

5._at        : black           =           1
               female          =           0
               age             =          20

6._at        : black           =           1
               female          =           0
               age             =          70

7._at        : black           =           1
               female          =           1
               age             =          20

8._at        : black           =           1
               female          =           1
               age             =          70

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         _at |
          1  |      0.005      0.001     5.99   0.000        0.004       0.007
          2  |      0.096      0.007    13.51   0.000        0.082       0.110
          3  |      0.006      0.001     6.10   0.000        0.004       0.008
          4  |      0.110      0.007    15.05   0.000        0.096       0.125
          5  |      0.011      0.002     5.24   0.000        0.007       0.015
          6  |      0.179      0.020     9.08   0.000        0.140       0.217
          7  |      0.013      0.002     5.32   0.000        0.008       0.018
          8  |      0.203      0.021     9.66   0.000        0.161       0.244
------------------------------------------------------------------------------

4. 小结

总的来说，Hoaglin (2016) 认为在解释多元回归系数时，文献和教科书中常用的「保持不变」，「控制」和「固定」等词语均不能很好的反应多元回归的基本原理，也不符合现实情况。

他认为，「调整共同线性影响」的表述会更适合。此外，在运用 margins 命令时应该谨慎选择「样本均值处」的边际效应。需要根据实际情况，弄清楚「样本预测空间」，并据此在合适的「点位」上求取边际效应，才能保证我们所着力解释的系数有真正的经济含义。

参考文献

Hoaglin David C., 2016, Regressions are Commonly Misinterpreted, Stata Journal, 16(1): 5–22. [PDF]
Hardin James W. , 2016, Regressions are Commonly Misinterpreted: Comments on the Article, Stata Journal, 16(1): 23–24. [PDF]
Long J. Scott, David M. Drukker, 2016, Regressions are Commonly Misinterpreted: Comments on the Article, Stata Journal, 16(1): 25–29. [PDF]
Hoaglin David C., 2016, Regressions are Commonly Misinterpreted: A Rejoinder, Stata Journal, 16(1): 30–36. [PDF]
胡雨霄, 连玉君. 2019. 「图示线性回归系数：Frisch-Waugh定理与部分回归图」，连享会 ( https://www.lianxh.cn ) 推文。
杨柳, 连玉君. 2019. 「Stata：因子变量全攻略」, 连享会 ( https://www.lianxh.cn ) 推文。