cover
yiyuezhuo 2018年5月10日 14:09

篡改数据操纵p值的正确姿势

为某些为追求显著性(* ,**, ***)而操纵数据的同学献上此文

对回归并不十分了解,但却要硬凑量化方法的人发现不显著后不禁抓耳挠腮,较为低端的人士直接修改表格中的p值与*,然后被发现与表中其他统计量不一致(如t统计量),较为高端的人士直接修改数据让其显著。本文考虑回归的显著性对数据修改之脆弱性与直接“优化”p值的正确姿势。

回归涉及的优化有时是对残差的最小化(最小二乘法)或对似然度的最大化(极大似然法),但如果我们把优化目标改成某个我们想让它最小化的p值会怎么样呢?换而言之,优化变量变为数据,通过一个回归计算过程估计出系数,此系数结合一些其他运算算出的单个系数p值,或模型显著性成为优化目标。

首先回顾OLS p值的计算方式:

P值的计算

经典回归对数据的模型为

\begin{align} \mathbf{y} &= \mathbf{X} \beta + \epsilon \\ \epsilon &\sim N(0,\sigma^2\mathbf{I}) \end{align}

最小二乘法通过正规方程估计系数,这可以表示为:

\hat{\mathbf{\beta}} = (\mathbf{X}^{\prime} \mathbf{X})^{-1} \mathbf{X}^{\prime} \mathbf{y}

容易计算出这个估计量(它的随机性完全来自随机误差向量 \epsilon )的协方差矩阵为

\textrm{Var}(\hat{\mathbf{\beta}}) = \sigma^2 (\mathbf{X}^{\prime} \mathbf{X})^{-1}

尽管分析缺失变量影响时有用,但统计表通常抛弃除了对角线上的各估计分量的方差以外的项:

\sigma(\hat{\beta}) = \sqrt{\mathrm{diag}(\mathrm{Var}(\hat{\beta}))}

这里 \sigma(\hat{\beta}) 是个向量, \mathrm{Var}(\hat{\beta}) 是个矩阵, \mathrm{diag}(\cdot) 指取矩阵的对角元素构成一个向量, \sqrt{\cdot} 指对向量里的每个元素求平方根。后两个符号来自编程中通常对这两个函数的揭示方式。

注意矩阵 \mathrm{Var}(\hat{\beta}) 与向量 \sigma(\hat{\beta}) 为参数,稍后讨论它们的估计方法。

另一条路上,估计通常被忽视的模型参数 \sigma ,首先估计隐变量 \epsilon :

e = \hat{\epsilon} = Y - X \hat{\beta}

估计 \sigma

\hat{\sigma^2} = \frac{e'e}{n-k}

其中 n 为样本规模,或者说 \epsilon 的长度,或者说自由度。 k 为模型用掉的自由度,即自变量个数,估计中对 \epsilon 施加的线性约束个数,据此,在一元回归中为 k=2 (这 k 包不包括常数项不同的地方有不同的假定)

有了 \sigma^2 的估计量就可以估计 参数 \mathrm{Var}(\hat{\beta}) 了(Emmm...估计只是随机变量以某种方式计算的函数,至于这么计算为什么有用,性质哪来的,如这种常见的替换模式为什么是有效的是另一回事,事实证明也很少有人在乎它们)

\hat{\textrm{Var}}(\hat{\mathbf{\beta}}) = \hat{\sigma}^2 (\mathbf{X}^{\prime} \mathbf{X})^{-1}

\mathrm{se(\hat{\beta})} = \hat{\sigma}(\hat{\beta}) = \sqrt{\mathrm{diag}(\mathrm{\hat{Var}}(\hat{\beta}))}

这个 \mathrm{se}(\hat{\beta}) 就是统计表里显示的那个数,当然既然它也是个随机的估计量,我们还能讨论它对应的参数 \sigma(\hat{\sigma}(\hat{\beta})) 与估计量 \hat{\sigma}(\hat{\sigma}(\hat{\beta})) 甚至 \sigma(\hat{\sigma}(\hat{\sigma}(\hat{\beta}))) 之类的,不过看上去这并不是很有用。

假设 \sigma(\hat{\beta}) 并不用估计的话,可以造出 z 统计量

z(\hat{\beta}) = \frac{\hat{\beta}-\beta^H}{\sigma(\hat{\beta})}

这里除是对应元素除, \beta^H 为假设的 \beta 值。假设 E(\hat{\beta}) = \beta = \beta^H 成立的话有

z(\hat{\beta}) \sim N(0,I)

对应的双尾p值有

p_z(\hat{\beta}) = 2(1-\Phi(|(z(\hat{\beta})|))

其中 \Phi(\cdot) 为标准正态分布累积密度分布函数(cdf)

假设成立则有

p_z(\hat{\beta}) \sim U(0,1)

拿单个出来看就有 z(\hat{\beta}_i) \sim N(0,1) ,于是可以看 z 是不是在 N(0,1) 的较极端的位置——不太可能出现的情况——较小的双尾p值——某一显著性水平下的最优假设检验拒绝域来判断这个假设 (常见的是 \beta^H_i = 0 这种形式)是不是不太合理。

然而参数 \sigma(\hat{\beta}) 并不是已知的(知道 参数 \sigma 的值就可以知道这个值,这并不是没有意义,统计在讨论“一致性”和某些性质时是假设这些值已知的。虽然使用者可能根本不知道这些性质,以及不知道自己的使用方法是错误的),这时用估计量替换它得到其“估计量”, t 统计量:

t(\hat{\beta}) = \hat{z}(\hat{\beta})= \frac{\hat{\beta}-\beta^H} {\mathrm{se}(\hat{\beta})}

假设成立时有(常见的假设是 \beta^H = \mathbf{0}

t(\hat{\beta}) \sim t(n-k)

这里指统计量向量 t(\hat{\beta}) 每个都服从自由度为 n-k 的分布

对应的 p 值统计量为

p_t(\hat{\beta}) = 2(1-t_{n-k}(|(t(\hat{\beta})|))

其中 t_{n-k}(\cdot) n-k 自由度的 t 分布的cdf。

假设成立时有

p_t(\hat{\beta}) \sim U(0,1)

p 值就是一个(假设成立下)服从 (0,1) 均匀分布的统计量,就像 z 值(统计量)就是服从 N(0,1) 分布的统计量一样,这些统计量通常通过某种“标准化”导出,将原来统计量的极端性以标准的方式显出,较小的p值就是个常见的“极端”的标准)

直接优化p值

这里利用梯度下降优化感兴趣的系数的p值。因为最小二乘法估计是由矩阵计算出来的而不是似然度梯度上升之类的方法,所以系数仅仅可以看做计算中的可导的一步。利用梯度下降的另一个困难是t分布的分布函数( t_{n-k}(\cdot) )不好计算,我这直接用方差与对应t分布一致的正态分布函数替换了(其实直接优化t值好像也行,因为p值这个操作是单调的)。正态分布函数虽然的计算虽然也不是解析的,但有高效的近似算法(这个算法不应该是迭代的,否则梯度计算会很麻烦)。我这直接用pytorch自动求导了。

单个点优化

给定10个二维数据点如图所示,x的p值为0.6969,极其不显著

删除样本点的效果是显然的就不展示了,无约束增加样本点会出现最优值不收敛的情况,如图所示,一个坐标对应的数表示在该点插入样本所计算出的p值

这很容易理解,因为这实际上相当于在回归两个点,一个点是新插入的点,另一个点是原先的“聚在一团”的等价点。类似下图所示:

一个比较有意思的是不修改自变量,而只修改因变量,有时候也许也只能这么改,如时间序列数据。

下面用贪婪法修改数据,每步迭代每个还没修改的样本点,求它使得p值最小化的取值,取其中使得p值下降最多者。然后下一步再在剩余样本点中这么找,以此类推,结果如下图所示:

可以看出,只需修改3个数就可以使得模型在0.05显著性水平下显著了。另外之所以修改单个因变量值为何不会产生不收敛情况?比如可以考虑修改第一个点一样把这个点拉到更高?因为这也会拉高其他点的预测值,增加其余点的残差,从而增加模型方差的估计量的估计值,从而使得t统计量变小,尽管偏斜率的估计值上升了。

全局优化

施工中。。。

代码: jupyter notebook: yiyuezhuo/Mathematical-statistics

相关阅读
  • 推荐阅读
  • 文章导航