Stata：多重假设修正-rwolf

👇 连享会 · 推文导航 | www.lianxh.cn

• 🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令
• 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享
• 💹 计量：回归分析 | 交乘项-调节  | IV-GMM | 时间序列 | 面板数据  | 空间计量 | Probit-Logit | 分位数回归
• ⛳ 专题：SFA-DEA | 生存分析  | 爬虫 | 机器学习 | 文本分析
• 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching
• 🔨 工具：工具软件 | Markdown | Python-R-Stata
• 🎧 课程：公开课-直播 |  计量专题 | 关于连享会

连享会 · 2022暑期班

作者：文海铭 (广西大学)
邮箱：[e-mail]

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：

 

编者按：本文主要摘译自下文，特此致谢！
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. 问题背景

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

• 3. 命令介绍

• 4. 具体应用

• 4.1 案例 1

• 4.2 案例 2

• 5. 相关推文

---

1. 问题背景

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

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

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

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

假设我们想要用数据 来检验 个假设，每个假设 都有一个感兴趣的参数 ，这个参数的估计值 和对应的标准误 关联。对于 ，我们通常假设 。进一步假设替代假设都是 的单边类型，或者都是 的双边类型。 的 检验统计量由以下公式给出：

接下来考虑 的 个重采样数据矩阵，用 表示。它们会产生估计值，用 表示，以及相应的标准误差，用 表示，且 。对于每个样本 和每个假设 ，一个 统计量由以下公式给出：

这些统计量 以 0 为中心，因为分子是由再抽样估计值减去原始估计值 (而不是无效参数) 组成的，因此 的分布将构成程序的 分布。

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

和

如上所述，我们将按照显著性顺序重新标记假设，但现在是基于它们的检验统计量 ，而不是它们的 值 ，就像 Holm 程序所做的那样。此时， 指的是具有最大的相应检验统计量的假设 (标记为 )，而 指的是具有最小的假设 (标记为 )。在下文中，我们用 表示向量 。

对于一个给定的 ， 表示统计量 的经验 四分位数。例如，在测试 个假设的情况下， 表示 的最大值， 表示子向量 的最大值 (即只对应于三个最不显著假设的检验统计数据向量)，以此类推。

最后，我们简单地得到 = 。一个重要的结果是 在 中是弱递减的，也就是说 ，。

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

1. 对于 ，拒绝 。
2. 用 表示拒绝的假设数量。如果 ，则停止；否则，让 。
3. 对于 ，拒绝 。
4. 如果没有进一步的假设被拒绝，则停止。否则，用 Rj 表示到目前为止被拒绝的假设数量，之后让 。然后返回到第 3 步。

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

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

对 ，令

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

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)：指定计算基于单边检验的 值；
• bl(string)：允许在每个模型中加入因变量的基线测量值作为控制；
• iv(varlist)：只有在指定 method(ivregress) 时才是必要的；
• otherendog(varlist)：如果在 ivregress 模型中需要一个以上的内生变量，那么 otherendog(varlist) 将指定额外的内生变量；
• indepexog：如果指定了 method(ivregress)，但 indepvar() 是一个外生变量，则应该指定 indepexog；
• holm：指定与标准输出一起提供 Holm(1979) 校正对应的 值。

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
------------------------------------------------------------------------------

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

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

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 直方图显示

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

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

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

在最不显著变量 的情况下，这两个分布是相似的，因为 Null 分布只基于单一回归的 回归估计值，而这是结果变量。

5. 相关推文

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

• 专题：Stata教程
• Stata：自己动手做组间系数差异检验-bootstrap-bdiff
• Stata检验：AIC-BIC-MSE-MAE-等信息准则的计算
• 专题：Stata命令
• Stata：堆叠回归及组间差异联合检验
• Stata：线性趋势检验-lintrend
• ebalance：基于熵平衡法的协变量平衡性检验
• 品头论足-distcomp：组间分布差异检验
• reldist-相对分布：分布差异分析和检验
• 专题：数据处理
• Stata：如何正确检验U型关系的存在
• Stata：多变量均值和中位数差异检验lianxh.cn/news/051e3a01cdb19.html)
• 专题：面板数据
• Stata：xtbunitroot-允许结构突变的面板单位根检验
• xtgranger：面板格兰杰因果检验-C103
• Stata：异方差和自相关稳健F检验和t检验-T316
• Stata面板：suest支持面板数据的似无相关检验
• Stata面板：Granger-因果检验
• Sargan+Hansen：过度识别检验及Stata实现
• 专题：IV-GMM
• IV专题- 内生性检验与过度识别检验
• Stata：工具变量的秩检验-bootrantest
• IV专题: 内生性检验与过度识别检验
• 专题：倍分法DID
• 多期DID之安慰剂检验、平行趋势检验
• 多期DID：平行趋势检验图示

课程推荐：因果推断实用计量方法
主讲老师：丘嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp

 

New！ Stata 搜索神器：lianxh 和 songbl  GIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
  . ssc install lianxh
  . ssc install songbl
👉  使用：
  . lianxh DID 倍分法
  . songbl all

🍏 关于我们

• 连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。
• 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。