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

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 个假设,每个假设 Hs 都有一个感兴趣的参数 θs,这个参数的估计值 θ^s 和对应的标准误 δ^s 关联。对于 s=1,,S,我们通常假设 θs0=0。进一步假设替代假设都是 Hs:θs>0 的单边类型,或者都是 Hs:θs0 的双边类型。Hs 的 Studentized 检验统计量由以下公式给出:

接下来考虑 X 的 M 个重采样数据矩阵,用 X1=0XM=0 表示。它们会产生估计值,用 θms 表示,以及相应的标准误差,用 σ^ms 表示,且 m=1,,M。对于每个样本 m 和每个假设 Hs,一个 Studentized(null) 统计量由以下公式给出:

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

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

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

对于一个给定的 jC^(1α,j) 表示统计量 maxt,j,mmM=1 的经验 1α 四分位数。例如,在测试 S=4 个假设的情况下,maxt,1,m 表示 t1,t(2),t(3),t(4) 的最大值,maxt,2,m 表示子向量 t(2),t(3),t(4) 的最大值 (即只对应于三个最不显著假设的检验统计数据向量),以此类推。

最后,我们简单地得到 max(4),mt(4),m。一个重要的结果是 C^(1α,j) 在 j 中是弱递减的,也就是说 C^(1α,j)C^(1α,j+1)j=1,,S1

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

  1. 对于 s=1,,S,拒绝 H(s)ift(s)>C^(1α,1)
  2. 用 R1 表示拒绝的假设数量。如果 R1=0,则停止;否则,让 j=2
  3. 对于 s=Rj1+1,,S,拒绝 H(s)ift(s)>C^(1α,Rj1+1)
  4. 如果没有进一步的假设被拒绝,则停止。否则,用 Rj 表示到目前为止被拒绝的假设数量,之后让 j=j+1。然后返回到第 3 步。

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

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

对 s=2,,S,令

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

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 时被拒绝。特别是,变量 y4 和 y7 在两个未校正程序中的 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 直方图显示

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

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

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

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

5. 相关推文

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

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

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

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

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

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

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