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 系列推文:

PDF下载 - 推文合集

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

编者按:本文整理自 Types Of Transformations For Better Normal Distribution,特此致谢!


目录


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

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

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

2. 判断指标:偏度(Skewness)和峰度(Kurtosis)

偏度衡量随机变量概率分布的不对称性,是相对于平均值不对称程度的度量。通过对偏度系数的测量,我们能够判定数据分布的不对称程度以及方向,其计算公式为 S=1ni=1n(XiX¯)31ni=1n(XiX¯)23。标准正态分布的偏度为 0,而非标准的数据分布分别被称为 “左偏” (偏度小于零,有极小值) 和 “右偏” (偏度大于零,有极大值)。

峰度是研究数据分布陡峭或平滑的统计量,通过对峰度系数的测量,我们能够判定数据相对于正态分布而言是更陡峭 / 平缓,其计算公式为 K=1ni=1n(XiX¯)41ni=1n(XiX¯)243 。标准正态分布的峰度为 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 为底数的 Log(X) 对数值作为新的分布数据。当原始数据中有负值及零值时,亦可做 Y=Log(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 变换是一种广义幂变换方法,可以明显地改善数据的正态性和方差齐性,其计算公式为:

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

  • θ=0 时相当于对数变换;
  • θ=2 时等同于平方变换;
  • θ=1 时相当于没有进行变换;
  • θ=12 时相当于平方根变换;
  • θ=1 时等同于倒数变换。

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

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

depvarindepvars:被解释变量和解释变量。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 变换在实数域的推广。当随机变量非负时,其表达式为:

当随机变量为负时,其表达式为:

由于 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小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh