Stata：多重假设修正-rwolf

发布时间：2022-07-04 阅读 1150

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下载 - 推文合集

作者：文海铭 (广西大学)
邮箱：hming_wen@sina.com

编者按：本文主要摘译自下文，特此致谢！
Source：Clarke D, Romano J P, Wolf M. The Romano–Wolf multiple-hypothesis correction in Stata[J]. The Stata Journal, 2020, 20(4): 812-843. -PDF-

1. 问题背景

当同时考虑多个假设检验时，除非明确检验框架的多重性，否则标准的统计方法会对无效假设过度拒绝。在本文中，我们讨论了 Romano-Wolf 多重假设修正，并记录了它在 Stata 中的实现。

Romano-Wolf 校正 (渐进地) 控制了族内误差率，即在被测假设族中至少拒绝一个真实无效假设的概率。这种校正比早期的多重检验程序，如 Bonferroni 和 Holm 校正，要强大得多，它通过对原始数据的重新取样，考虑到了检验统计的依赖结构。

本推文的目的在于介绍实现这种校正的命令 rwolf，并结合案例加以说明。

2. Romano-Wolf 的多重假设修正

假设我们想要用数据 $X$ 来检验 $S$ 个假设，每个假设 $H_{s}$ 都有一个感兴趣的参数 $θ_{s}$ ，这个参数的估计值 ${\hat{θ}}_{s}$ 和对应的标准误 ${\hat{δ}}_{s}$ 关联。对于 $s = 1, \dots, S$ ，我们通常假设 $θ_{s}^{0} = 0$ 。进一步假设替代假设都是 $H_{s}^{'} : θ_{s} > 0$ 的单边类型，或者都是 $H_{s}^{'} : θ_{s} \neq 0$ 的双边类型。 $H_{s}$ 的 $S t u d e n t i z e d$ 检验统计量由以下公式给出：

t_{s} := \frac{{\hat{θ}}_{s}}{{\hat{σ}}_{s}}

接下来考虑 $X$ 的 $M$ 个重采样数据矩阵，用 $X_{*}^{1} = 0 \dots X_{*}^{M} = 0$ 表示。它们会产生估计值，用 $θ_{*}^{m s}$ 表示，以及相应的标准误差，用 ${\hat{σ}}_{*}^{m s}$ 表示，且 $m = 1, \dots, M$ 。对于每个样本 $m$ 和每个假设 $H_{s}$ ，一个 $S t u d e n t i z e d (n u l l)$ 统计量由以下公式给出：

t_{s}^{*, m} := \frac{{\hat{θ}}_{s}^{*, m} - {\hat{θ}}_{s}}{{\hat{σ}}_{s}^{*}, m}

这些统计量 $t_{s}^{*, m}$ 以 0 为中心，因为分子是由再抽样估计值减去原始估计值 (而不是无效参数) 组成的，因此 $t_{s}^{*, m}$ 的分布将构成程序的 $(n u l l)$ 分布。

如果替代假设是 $H_{s}^{'} : θ_{s} \neq 0$ 类型的双侧假设，则必须使用检验统计量的绝对值。为了在下面的算法中保持符号的统一，我们将在双侧情况下使用以下计算惯例 (但在单侧情况下不使用)：

t_{s} := | t_{s} |

和

t_{s}^{*, m} := | t_{s}^{*}, m |

如上所述，我们将按照显著性顺序重新标记假设，但现在是基于它们的检验统计量 $t_{s}$ ，而不是它们的 $p$ 值 $p_{s}$ ，就像 Holm 程序所做的那样。此时， $H (1)$ 指的是具有最大的相应检验统计量的假设 (标记为 $t (1)$ )，而 $H (S)$ 指的是具有最小的假设 (标记为 $t (S)$ )。在下文中，我们用 ${max}_{t, j}^{*, m}$ 表示向量 ${t_{(j)}^{*, m}, \dots, t_{(S)}^{*, m}}$ 。

{max}_{t, j}^{*, m} := max {t_{(j)}^{*, m}, \dots, t_{(S)}^{*, m}}; f o r j = 1, \dots, S 和 m = 1, \dots, M

{max}_{t, j}^{*, m} := max {t_{(j)}^{*, m}, \dots, t_{(S)}^{*, m}}; f o r j = 1, \dots, S 和 m = 1, \dots, M

对于一个给定的 $j$ ， $\hat{C} (1 - α, j)$ 表示统计量 ${max}_{t, j}^{*, m}_{m}^{M} = 1$ 的经验 $1 - α$ 四分位数。例如，在测试 $S = 4$ 个假设的情况下， ${max}_{t, 1}^{*, m}$ 表示 $t_{1}^{*}, t_{(2)}^{*}, t_{(3)}^{*}, t_{(4)}^{*}$ 的最大值， ${max}_{t, 2}^{*, m}$ 表示子向量 $t_{(2)}^{*}, t_{(3)}^{*}, t_{(4)}^{*}$ 的最大值 (即只对应于三个最不显著假设的检验统计数据向量)，以此类推。

最后，我们简单地得到 ${max}_{(4)}^{*, m}$ = $t_{(4)}^{*, m}$ 。一个重要的结果是 $\hat{C} (1 - α, j)$ 在 $j$ 中是弱递减的，也就是说 $\hat{C} (1 - α, j) \geq \hat{C} (1 - α, j + 1)$ ， $j = 1, \dots, S - 1$ 。

基于上述情况， $α$ 级的主要递减多重检验程序 (基于 Romano 和 Wolf (2016) 的算法) 可以总结为以下几点:

对于 $s = 1, \dots, S$ ，拒绝 $H_{(s)} if t_{(s)} > \hat{C} (1 - α, 1)$ 。
用 $R 1$ 表示拒绝的假设数量。如果 $R 1 = 0$ ，则停止；否则，让 $j = 2$ 。
对于 $s = R_{j - 1} + 1, \dots, S$ ，拒绝 $H_{(s)} if t_{(s)} > \hat{C} (1 - α, R_{j - 1} + 1)$ 。
如果没有进一步的假设被拒绝，则停止。否则，用 Rj 表示到目前为止被拒绝的假设数量，之后让 $j = j + 1$ 。然后返回到第 3 步。

如同 Holm (1979) 的程序一样，这个修正是一个递减程序， $\hat{C} (1 - α, j)$ 在 $j$ 中是弱递减的。因此，拒绝的标准对较显著的假设要求更高，对较不显著的假设要求更低。基于 ${max}_{t, j}^{*, m}$ 的 Null 分布是基于原始数据的再抽样，它们隐含地考虑了检验统计量的基本依赖结构，导致了比 Holm 程序的潜在的显著增加，该程序假定了最坏情况下的依赖结构。

上述算法导致了在给定的显著性水平 $α$ 下，对每个无效假设 $H_{s}$ 的拒绝或不拒绝的决定。然而，也许更方便的是，可以直接计算每个 $H_{s}$ 的多重检验调整的 $p$ 值，如下面的算法所述。这个算法是对单一无效假设的基于重抽样的 $p$ 值的概括。

p_{(1)}^{adj} := \frac{# {{max}_{t, 1}^{*, m} \geq t_{(1)}} + 1}{M + 1}

对 $s = 2, \dots, S$ ，令

p_{(s)}^{initial} := \frac{# {{max}_{t, m}^{*, m} \geq t_{(s)}} + 1}{M + 1}

然后通过定义来增加单调性，即

p_{(s)}^{a d j} := max {p_{(s)}^{i n i t i a l}, p_{(s - 1)}^{a d j}}

3. 命令介绍

命令安装：

ssc install rwolf, replace

命令语法：

rwolf depvars [if] [in] [weight], [options]

其中，depvars 为需要分析的多个结果变量。主要选项如下：

indepvar(varlist)：表示包含在多假设检验中的独立 (处理) 变量；
method(string)：向 Stata 指出每个多假设检验 (即基线模型) 是如何进行的；
controls(varlist)：表示在每个 depvar 对感兴趣的 indepvar 进行回归时应该包括的额外控制；
strata(varlist)：指定了识别分层的变量；
cluster(varlist)：指定了识别重采样聚类的变量；
onesided(string)：指定计算基于单边检验的 $p$ 值；
bl(string)：允许在每个模型中加入因变量的基线测量值作为控制；
iv(varlist)：只有在指定 method(ivregress) 时才是必要的；
otherendog(varlist)：如果在 ivregress 模型中需要一个以上的内生变量，那么 otherendog(varlist) 将指定额外的内生变量；
indepexog：如果指定了 method(ivregress)，但 indepvar() 是一个外生变量，则应该指定 indepexog；
holm：指定与标准输出一起提供 Holm(1979) 校正对应的 $p$ 值。

4. 具体应用

4.1 案例 1

. set seed 130319     // 设置种子数
. set obs 100         // 设置观测值数
. generate treat = runiform()>0.5 
. foreach num of numlist 1(1)10 { 
  2.     generate c`num'= rnormal()
  3. }                // 生成 treat 值，均匀分布
. local r=0.25        // 存储 r 值
. #delimit ; 
. matrix corr =  
> (1,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' \ 
> `r' ,1 ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' \ 
> `r' ,`r' ,1 ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' \ 
> `r' ,`r' ,`r' ,1 ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' \ 
> `r' ,`r' ,`r' ,`r' ,1 ,`r' ,`r' ,`r' ,`r' ,`r' \ 
> `r' ,`r' ,`r' ,`r' ,`r' ,1 ,`r' ,`r' ,`r' ,`r' \ 
> `r' ,`r' ,`r' ,`r' ,`r' ,`r' ,1 ,`r' ,`r' ,`r' \ 
> `r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,1 ,`r' ,`r' \ 
> `r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,1 ,`r' \ 
> `r',`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,`r' ,1 ) ;
.  #delimit cr

先看看我们定义的 corr 矩阵：

. mat list corr

symmetric corr[10,10]
      c1   c2   c3   c4   c5   c6   c7   c8   c9  c10
 r1    1
 r2  .25    1
 r3  .25  .25    1
 r4  .25  .25  .25    1
 r5  .25  .25  .25  .25    1
 r6  .25  .25  .25  .25  .25    1
 r7  .25  .25  .25  .25  .25  .25    1
 r8  .25  .25  .25  .25  .25  .25  .25    1
 r9  .25  .25  .25  .25  .25  .25  .25  .25    1
r10  .25  .25  .25  .25  .25  .25  .25  .25  .25    1

继续进行后续分析：

. matrix corsim = cholesky(corr)
. foreach num of numlist 1(1)10 { 
  2.     matrix eps`num'= corsim[`num',1..10] 
  3.     matrix score epsilon`num'= eps`num'
  4.     generate y`num'= 1 + 0*treat + epsilon`num'
  5. }
. rwolf y1 y2 y3 y4 y5 y6 y7 y8 y9 y10, indepvar(treat) reps(1000) nodots holm graph

Independent variable:  treat
Outcome variables:   y1 y2 y3 y4 y5 y6 y7 y8 y9 y10
Number of resamples: 1000
------------------------------------------------------------------------------
                    |     Model       Resample     Romano-Wolf        Holm
   Outcome Variable |    p-value      p-value        p-value         p-value
--------------------+---------------------------------------------------------
                 y1 |    0.1142        0.1009         0.5534         0.8072
                 y2 |    0.8906        0.8931         0.9940         0.8931
                 y3 |    0.2750        0.2797         0.7872         1.0000
                 y4 |    0.0292        0.0280         0.1938         0.2517
                 y5 |    0.1914        0.1818         0.6883         1.0000
                 y6 |    0.8683        0.8741         0.9940         1.0000
                 y7 |    0.0137        0.0100         0.1009         0.0999
                 y8 |    0.8337        0.8372         0.9940         1.0000
                 y9 |    0.3849        0.3966         0.8382         1.0000
                y10 |    0.1199        0.1149         0.5534         0.8042
------------------------------------------------------------------------------

该命令返回一个与每个多重结果相关的 $p$ 值列表。模型 $p$ 值一栏列出了直接来自估计回归模型的分析性 $p$ 值，在每种情况下都是基于 $t$ 统计量和具有适当自由度的 $t$ 分布的 (反) 累积分布函数。再抽样 $p$ 值一栏列出了基于再抽样的 $p$ 值，它也没有对多重检验进行校正。

在这两个未校正程序的情况下，尽管所有真实的 $β_{1}$ 值都是 0，但许多 $β_{1} = 0$ 的假设都在 $α = 0.05$ 时被拒绝。特别是，变量 $y 4$ 和 $y 7$ 在两个未校正程序中的 $p$ 值都低于 0.05。第三列显示的是 rwolf 调整后的 $p$ 值，我们注意到现在没有任何一个空假设被拒绝 (即使是在 $α = 0.10$ 的时候)。

4.2 案例 2

之前展示的基于 rwolf 的模拟例子提供了一个标准的实现，其中多个结果变量被回归到一个独立 (处理) 变量上，每个变量都使用普通最小二乘法回归。然而，rwolf 的标准命令语法允许对这个基线实现进行许多扩展，包括使用替代估计方法 (例如，IV 回归、probit 回归和其他 Stata 估计命令)，实施单边检验，或使用 bootstrap 程序。

在这里，我们扩展到一个有多个结果的案例，并进行单边和双边的测试。

4.2.1 双边假设检验

. set seed 120113031 
. use http://www.stata-press.com/data/r13/hsng, clear 
. rwolf rent popden popgrow hsng, indepvar(hsngval) method(ivregress) ///
>     iv(faminc i.region) reps(10000) graph nodots controls(pcturban) 

Independent variable:  hsngval
Outcome variables:   rent popden popgrow hsng
Number of resamples: 10000
------------------------------------------------------------------------------
   Outcome Variable | Model p-value    Resample p-value    Romano-Wolf p-value
--------------------+---------------------------------------------------------
               rent |    0.0000             0.0157              0.0157
             popden |    0.0654             0.0848              0.0848
            popgrow |    0.0019             0.0119              0.0200
               hsng |    0.0236             0.0120              0.0515
------------------------------------------------------------------------------

4.2.2 单边假设检验

. use http://www.stata-press.com/data/r13/hsng, clear 
. replace popden=-popden 
. replace hsng =-hsng 
. rwolf rent popden popgrow hsng, indepvar(hsngval) method(ivregress) ///
>     iv(faminc i.region) reps(10000) graph onesided(negative) nodots ///
>     controls(pcturban) 

Independent variable:  hsngval
Outcome variables:   rent popden popgrow hsng
Number of resamples: 10000
------------------------------------------------------------------------------
   Outcome Variable | Model p-value    Resample p-value    Romano-Wolf p-value
--------------------+---------------------------------------------------------
               rent |    0.0000             0.0003              0.0003
             popden |    0.0327             0.0193              0.0193
            popgrow |    0.0010             0.0012              0.0017
               hsng |    0.0118             0.0060              0.0110
------------------------------------------------------------------------------

4.2.3 直方图显示

下图显示了双边情况下的这些分布。这里的直方图显示的是每一个 $(S t u d e n t i z e d)$ 估计值的绝对值，这些估计值来自于强加空值的 bootstrap 复制，黑色虚线显示的是一个精确的半正态分布。回归系数的实际 $S t u d e n t i z e d$ 值显示为一条垂直实线。

在左上角的面板上，第一个 Null 分布比理论分布要求高得多，因为它累积了每个结果的最大估计系数。在右上角和左下角的面板中，这些 Null 分布的要求越来越低，因为先前测试的变量被从池中移除以形成 Null 分布。最后，在右下角 (对于最不重要的变量)，Null 分布是基于仅来自该变量的 bootstrap，因此，Null 分布非常接近于理论上的半正态分布。

下图展示了同样的 Null 分布，但现在是基于单边测试。这里的直方图记录了每个自举复制中多个变量的最大值，而黑色虚线呈现了理论上的正态分布。再一次，当我们考虑到观察到更多显著关系的结果时，用于计算校正 $p$ 值的经验分布比黑点分布要求更高，黑点分布将在没有校正和正态性假设的情况下使用。

在最不显著变量 $(p o p d e n)$ 的情况下，这两个分布是相似的，因为 Null 分布只基于单一回归的 $s t u d e n t i z e d$ 回归估计值，而这是结果变量。

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 检验, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

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

Stata：多重假设修正-rwolf

1. 问题背景

2. Romano-Wolf 的多重假设修正

3. 命令介绍

4. 具体应用

4.1 案例 1

4.2 案例 2

4.2.1 双边假设检验

4.2.2 单边假设检验

4.2.3 直方图显示

5. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯