Stata：正态转换的五种方法

发布时间：2023-06-14 阅读 1064

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

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

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

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

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

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：王硕 (山东大学)
邮箱：shweung@163.com

编者按：本文整理自 Types Of Transformations For Better Normal Distribution，特此致谢！

1. 数据需要正态转换的原因

以概率分布为核心的研究大都聚焦于正态分布。例如一般线性模型假定 $Y = X β + ϵ$ ，其中误差项 $ϵ$ 服从或近似服从正态分布，但由于 $ϵ$ 可能和 $Y$ 有关，而 $Y$ 不服从正态分布，这会给线性回归的最小二乘估计系数结果带来误差。此时，就需要对数据正态化处理以有效降低模型的均方根误差 (RMSE) 或均方根对数误差 (RMSLE)。

在本推文中，我们将给出 2 个判断指标以及 5 种类型的数据转换方案，以更好地辨别和适应正态分布。

2. 判断指标：偏度（Skewness）和峰度（Kurtosis）

偏度衡量随机变量概率分布的不对称性，是相对于平均值不对称程度的度量。通过对偏度系数的测量，我们能够判定数据分布的不对称程度以及方向，其计算公式为 $S = \frac{\frac{1}{n} \sum_{i = 1}^{n} (X_{i} - \bar{X})^{3}}{{\sqrt{\frac{1}{n} \sum_{i = 1}^{n} (X_{i} - \bar{X})^{2}}}^{3}}$ 。标准正态分布的偏度为 0，而非标准的数据分布分别被称为 “左偏” (偏度小于零，有极小值) 和 “右偏” (偏度大于零，有极大值)。

峰度是研究数据分布陡峭或平滑的统计量，通过对峰度系数的测量，我们能够判定数据相对于正态分布而言是更陡峭 / 平缓，其计算公式为 $K = \frac{\frac{1}{n} \sum_{i = 1}^{n} (X_{i} - \bar{X})^{4}}{{\sqrt{\frac{1}{n} \sum_{i = 1}^{n} (X_{i} - \bar{X})^{2}}}^{4}} - 3$ 。标准正态分布的峰度为 3，而非标准的数据分布分别被称为 “尖峰” (峰度大于 3，数据点集中于均值) 和 “厚尾” (峰度小于 3，数据点较分散) 。

2.1 方法一：计算统计量

在查看数据偏度和峰度统计量时，可以借助 Stata 中的 tabstat 命令，并加入 count mean p50 skewness kurtosis 等选项呈现出观测值、均值、中位数、峰度和偏度等指标。以 Stata 软件自带的数据文件 nlsw88.dta 中的 wage 工资变量为例进行演示，该数据包含了 1988 年 2246 个美国妇女的资料，其命令的语法结构如下：

. sysuse nlsw88, clear
. tabstat wage, stat (count mean p50 skewness kurtosis)

 Variable |    N      Mean      p50  Skewness  Kurtosis
----------+---------------------------------------------
     wage | 2246  7.766949  6.27227  3.096199  15.85446
--------------------------------------------------------

从呈现结果可以看出，偏度为 3.096，均值 7.767 高于中位数 6.272，因此很明显工资变量是一个右偏的变量。峰度为 15.854，远大于标准正态分布的峰值 3，表明工作变量较为集中地分布在均值附近。

2.2 方法二：图示法

图形可以便于我们理解连续变量分布，从而直观地评估偏度和峰度，以建立对数据的初步了解。Stata 中的直方图可以较好地展示数据的偏度和峰度特征，其命令的语法结构如下：

. sysuse nlsw88,clear
. sum wage, d 
. local p50  = r(p50)
. local mean = r(mean)
. twoway (hist wage, vertical color("136 187 203")                      ///
>     xline(`mean', lcolor("220 146 20")) xline(`p50', lcolor(black)))  ///
>     (kdensity wage, lpattern(dash)),                                  ///
>     ylabel(0(0.035)0.14) graphregion(color(white))                    ///
>     xtitle(Wage,place(right)) ytitle(Density, place(top))             ///
>     legend(ring(0) col(1) position(2.5))                              
. gr save wage0.gph, replace

直方图横轴表示数据类型，纵轴表示分布情况。从上图可以看出，均值 (橙色虚线) 位于中位数 (黑色虚线) 右侧，且右侧分布较多明显的极大值，这同样可以得出工资变量是一个右偏的变量的结论。从峰度来看，数据集中于 3-6 范围内，且高出同均值和方差的正态分布曲线，呈现出 “尖峰” 样式。

2.2 方法三：外部命令

首先，采用由开发者提供的命令链接进行下载安装。

. ssc install jb6, replace
. sysuse nlsw88, clear
. jb6 wage
Jarque-Bera normality test:  1.9e+04 Chi(2)      0
Jarque-Bera test for Ho: normality: (wage)

从上述结果可以看出，该命令的原假设 H0 为数据服从正态分布，而结果中的 $p$ 值为 0 表明拒绝原假设，即数据并非服从正态分布。

3. 数据转换方案

3.1 对数转换 (Log Transformation)

对数变换即将原始数据 $X$ 转换为以 10 为基数、以 2 为基数或以自然对数 $e$ 为底数的 $L o g (X)$ 对数值作为新的分布数据。当原始数据中有负值及零值时，亦可做 $Y = L o g (X + K)$ 变换。以 10 为底数的对数转换为例，其 Stata 命令语法结构如下：

. sysuse nlsw88,clear
. gen Logwage=log(wage)
. twoway (hist wage, vertical color("136 187 203"))        ///
>     (kdensity wage, lpattern(dash)),                     ///
>     ylabel(0(0.035)0.14) graphregion(color(white))       ///
>     xtitle(Wage,place(right)) ytitle(Density,place(top)) ///
>     legend(ring(0) col(1) position(1))
. gr save wage.gph, replace

. twoway (hist Logwage,vertical color("96 133 149"))          ///
>     (kdensity Logwage, lpattern(dash)),                     ///
>     graphregion(color(white))                               ///
>     xtitle(Logwage,place(right)) ytitle(Density,place(top)) ///
>     legend(ring(0) col(1) position(1))
. gr save Logwage.gph, replace
. gr combine wage.gph Logwage.gph, ysize(3)

上图是原始数据 (左) 和对数转换数据 (右) 的比较。可以看出转换后的数据的偏度明显降低且峰度与正态分布接近。

3.2 平方根转换 (Square-Root Transformation)

平方根变换即将原始数据 $X$ 的平方根作为新的分布数据， $Y = sqrt (X)$ 。相较于对数转换，平方根转换的正态效果相对较弱，但其优点是适用于存在零值的数据。其命令语法结构如下：

. sysuse nlsw88,clear
. gen Sqrtwage = sqrt(wage)
. twoway (hist wage, vertical color("136 187 203"))        ///
>     (kdensity wage, lpattern(dash)),                     ///
>     ylabel(0(0.035)0.14) graphregion(color(white))       ///
>     xtitle(Wage,place(right)) ytitle(Density,place(top)) ///
>     legend(ring(0) col(1) position(1))
. gr save wage.gph, replace
. twoway (hist Sqrtwage, vertical color("96 133 149"))         ///
>     (kdensity Sqrtwage, lpattern(dash)),                     ///
>     graphregion(color(white))                                ///
>     xtitle(Sqrtwage,place(right)) ytitle(Density,place(top)) ///
>     legend(ring(0) col(1) position(1))
. gr save Sqrtwage.gph, replace
. gr combine wage.gph Sqrtwage.gph, ysize(3)

上图是原始数据 (左) 和平方根转换数据 (右) 的比较。从上图可以看出，无论是从偏度还是峰度角度，平方根转换效果均低于对数转换效果。

3.3 倒数转换 (Reciprocal Transformation)

在倒数转换中， $X$ 将转化为自身的倒数。此转换对分布的形状影响不大，且只能用于非零值数据。其命令语法结构如下：

. sysuse nlsw88,clear
. gen Reciwage=1/wage
. twoway (hist wage, vertical color("136 187 203"))         ///
>     (kdensity wage, lpattern(dash)),                      ///
>     ylabel(0(0.035)0.14) graphregion(color(white))        ///
>     xtitle(Wage,place(right)) ytitle(Density,place(top))  ///
>     legend(ring(0) col(1) position(1))
. gr save wage.gph, replace
. twoway (hist Reciwage, vertical color("96 133 149"))         ///
>     (kdensity Reciwage, lpattern(dash)),                     ///
>     graphregion(color(white))                                ///
>     xtitle(Reciwage,place(right)) ytitle(Density,place(top)) ///
>     legend(ring(0) col(1) position(1))
. gr save Reciwage.gph, replace
. gr combine wage.gph Reciwage.gph, ysize(3)

上图是原始数据 (左) 和倒数转换数据 (右) 的比较。从上图可以看出，无论是从偏度还是峰度角度，倒数转换未发生明显变动，其正态效果皆弱于对数转换和倒数转换。

3.4 伯克斯-考克斯转换 (Box-Cox Transformation)

Box-Cox 变换是一种广义幂变换方法，可以明显地改善数据的正态性和方差齐性，其计算公式为：

Y (θ) = {\begin{cases} \frac{X^{θ} - 1}{θ} ， & θ \neq 0 \\ \ln (X) ， & θ = 0 \end{cases}

其中， $X$ 为连续变量，且要求取值为正 (当原始数据中有负值及零值时，亦可取 $X^{'} = X + K$ ，然后对 $X^{'}$ 做 Box-Cox 变换)。 $θ$ 为变换参数，不同的 $θ$ 对应不同的变换方式。

$θ = 0$ 时相当于对数变换；
$θ = 2$ 时等同于平方变换；
$θ = 1$ 时相当于没有进行变换；
$θ = \frac{1}{2}$ 时相当于平方根变换；
$θ = - 1$ 时等同于倒数变换。

通过求解 $θ$ 值即可确定具体的变换方式， $θ$ 值的估计方法可采用最大似然估计。Stata 里面直接提供的是 Box-Cox 回归，可以利用 boxcox 命令获得参数 $θ$ 。boxcox 命令的语法结构如下：

boxcox depvar [indepvars] [if] [in] [weight] [, options]

depvar 和 indepvars：被解释变量和解释变量。options 选择项如下：

model(lhsonly)：指定只对因变量进行 Box-Cox 变换，而不对自变量进行变换。这种情况下，Box-Cox 转换的指数 $θ$ 是仅基于被解释变量 $Y$ 的最大似然估计 (MLE) 得到的。使用此选项可以在保持解释变量不变的情况下对被解释变量进行变换。
model(rhsonly)：指定只对解释变量进行 Box-Cox 变换，而不对被解释变量进行变换。这种情况下，Box-Cox 转换的指数 $θ$ 是仅基于自变量 $x$ 的最大似然估计 (MLE) 得到的。使用此选项可以在保持被解释变量不变的情况下对解释变量进行变换。
model(lambda)：指定 Box-Cox 转换的指数 $θ$ 是通过基于被解释变量和解释变量的联合最大似然估计 (MLE) 得到的。这是 boxcox 命令的默认选项，它同时对因变量和自变量进行 Box-Cox 变换。
model(theta)：指定 Box-Cox 转换的指数 $θ$ 是通过估计数据的 Mills 比例，基于被解释变量和解释变量的联合最大似然估计 (MLE) 得到的。
notrans：指定不对数据进行任何变换，仅计算出基于因变量和自变量的联合最大似然估计 (MLE) 得到的 Box-Cox 变换的指数 $θ$ 。使用此选项可以仅获得指数而不进行实际的变换。

下面仅对 wage 工资变量为例进行转换说明：

. sysuse nlsw88, clear
. boxcox wage, model(lhsonly)

                                       Number of obs =  2,246
                                       LR chi2(0)    =   0.00
Log likelihood = -6108.0084            Prob > chi2   =      .
----------------------------------------------------------------
wage | Coefficient Std. err.     z   P>|z|  [95% conf. interval]
-----+-----------------------------------------------------------
/theta| -.2241236  .0288636   -7.76  0.000  -.2806952    -.167552
-----------------------------------------------------------------
 
Estimates of scale-variant parameters
----------------------------
             | Coefficient
-------------+--------------
Notrans      |
       _cons |   1.502707
-------------+--------------
      /sigma |   .3727356
----------------------------

-----------------------------------------------------
   Test       Restricted    LR statistic
    H0:     log likelihood      chi2      Prob > chi2
-----------------------------------------------------
theta = -1    -6449.0767      682.14          0.000
theta =  0     -6138.868       61.72          0.000
theta =  1    -7117.2949     2018.57          0.000
-----------------------------------------------------

第一个表是 Box-Cox 变换参数 $θ$ 的结果， $θ = - 0.2241$ ，即 wage 变换值为 (wage^(-0.2241)-1)/(-0.2241)；

第二个表包含尺度变量参数的估计值，会显示各个自变量的参数估计值和模型标准差。本例只是对因变量进行了变化，该表实际意义不大；

第三个表包含了 Box-Cox 变换中三个常用函数的限制性似然比检验输出结果，即倒数变换 ( $θ = - 1$ )、对数变换 ( $θ = 0$ ) 和线性回归 ( $θ = 1$ )。

接下来，对 Box-Cox 变换参数 $θ = - 0.2241$ 进行验证，相应命令为：

. sysuse nlsw88, clear
. gen Boxwage=(wage^(-0.2241)-1)/(-0.2241)
. twoway (hist wage, vertical color("136 187 203"))        ///
>     (kdensity wage, lpattern(dash)),                     ///
>     ylabel(0(0.035)0.14) graphregion(color(white))       ///
>     xtitle(Wage,place(right)) ytitle(Density,place(top)) ///
>     legend(ring(0) col(1) position(1))
. gr save wage.gph, replace
. twoway (hist Boxwage, vertical color("96 133 149"))         ///
>     (kdensity Boxwage, lpattern(dash)),                     ///
>     graphregion(color(white))                               ///
>     xtitle(Boxwage,place(right)) ytitle(Density,place(top)) ///
>     legend(ring(0) col(1) position(10))
. gr save Boxwage.gph, replace
. gr combine wage.gph Boxwage.gph, ysize(3)

上图是原始数据 (左) 和 Box-Cox 转换数据 (右) 的比较。从上图可以看出，Box-Cox 转换的效果较好且与对数转化结果相似，实际上，对数转化模型是 Box-Cox 转换在模型选择时的一种情况 ( $θ = 0$ )。理论上，通过对 $θ$ 的选择，Box-Cox 转换能够选出不差于对数转化的模型。

3.5 杨-约翰逊转换 (Yeo-Johnson Transformation)

Yeo-Johnson 转换也是一种广义幂变换方法，通过构建一组单调函数对随机变量进行数据变换，其特点在于其可被应用于包含 0 值和负值的样本中，因此其也被认为是 Box-Cox 变换在实数域的推广。当随机变量非负时，其表达式为：

Y (θ) = {\begin{cases} \frac{(X + 1)^{θ} - 1}{θ} ， & θ \neq 0 ， X \geq 0 \\ \ln (X + 1) ， & θ = 0 ， X \geq 0 \end{cases}

当随机变量为负时，其表达式为：

Y (θ) = {\begin{cases} \frac{- [(- X + 1)^{2 - θ} - 1]}{2 - θ} ， & θ \neq 2 ， X < 0 \\ - \ln (- X + 1) ， & θ = 2 ， X < 0 \end{cases}

由于 Stata 中暂时没有专门处理 Yeo-Johnson 转换的命令，我们借助 Python 中的 scipy.stats 模块中的 yeojohnson 函数进行处理，相应命令为：

. sysuse nlsw88, clear
. python
---- python (type end to exit) ---
>>> from sfi import Data
>>> from scipy.stats import yeojohnson
>>> Wage=Data.get(va = 'wage')
>>> YJwage,lam = yeojohnson(Wage)
>>> print(YJwage,lam)
>>> Data.addVarFloat("YJwage")
>>> Data.store("YJwage",None,val=YJwage)
>>> end
-----------------------------------
. twoway (hist wage, vertical color("136 187 203"))        ///
>     (kdensity wage, lpattern(dash)),                     ///
>     ylabel(0(0.035)0.14) graphregion(color(white))       ///
>     xtitle(Wage,place(right)) ytitle(Density,place(top)) ///
>     legend(ring(0) col(1) position(1))
. gr save wage.gph, replace
. twoway (hist YJwage, vertical color("96 133 149"))         ///
>     (kdensity YJwage, lpattern(dash)),                     ///
>     graphregion(color(white))                              ///
>     xtitle(YJwage,place(right)) ytitle(Density,place(top)) ///
>     legend(ring(0) col(1) position(10))
. gr save YJwage.gph, replace
. gr combine wage.gph YJwage.gph, ysize(3)

上图是原始数据 (左) 和 Yeo-Johnso 转换数据 (右) 的比较。从上图可以看出，Yeo-Johnso 转换的效果与 Box-Cox 转换非常相似，也是一种较好的转化方式。

4. 总结

通过前述讨论，我们有以下几方面的发现：

通过偏度 (Skewness) 和峰度 (Kurtosis) 判断数据是否接近正态分布以及是否需要进行进行正态转化。
进行判断方法有计算统计量、图示法以及使用外部命令等。
能够较好地进行正态转化的方法有对数转换、Box-Cox 转换以及 Yeo-Johnso 转换，并且对数转化是后两者的一种特例。
平方根转换能够对数据进行正态转化但效果不十分明显，倒数转化的作用十分有限。

5. 参考资料

Tamil Selvan S，博客，Types Of Transformations For Better Normal Distribution
Yeo, I.K. and Johnson, R.A., 2000. A new family of power transformations to improve normality or symmetry. Biometrika, 87(4), pp.954-959. -PDF- , -Link-
Box, G.E. and Cox, D.R., 1964. An analysis of transformations. Journal of the Royal Statistical Society: Series B (Methodological), 26(2), pp.211-243. -PDF- , -Link-
Joanes, D.N. and Gill, C.A., 1998. Comparing measures of sample skewness and kurtosis. The Statistician, 47, pp.183–189. -PDF- , -Link-

6. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 对数 python交互, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]