现在大家都知道我强烈推荐在DID中进行“平行趋势敏感性分析”。这是因为Roth（2022，AER：Insights）指出，传统的处理前趋势检验效力在50%-80%之间。这个效力是什么？这就是今天想给大家介绍的，也是我的第三篇方法论模拟分析文章《因果推断的效力分析的最新趋势与实践指导》（前两篇分别是《听Rubin和Imbens的话：处理配置机制讨论的实践指导》和《如何控制协变量？DID中协变量调整的实践指导》）。

除了Roth的文章外，最近几年，效力分析也在Top期刊的应用研究中越来越常见，例如，Black et al.（2022，JPubE）、Dench et al.（2024，JPubE）。关于效力分析更详细的技术讲解，请参见Sylvain Chabé-Ferret（2024）的计量课本《Statistical Tools for Causal Inference》第七章。

下面的内容主要来自于Nick Huntington-Klein（2020）：《Simulation for Power Analysis》。

我们有了idea后，就要去做许多统计分析。在学习计量之初，老师就告诉我们，要尽可能的增加样本量。效力分析（Power Analysis）就是一种我们至少需要多少样本量得到给定处理效应，或者给定样本量可以得到最低处理效应是多少。

无论你进行哪种研究，效力分析都是一个好主意。然而，它在两种情况下特别有用：

• 如果你正在寻找一个你认为可能不是那么核心或重要的影响，即它是一个小的影响，或是一个系统的一部分，在这个系统中，有很多其他的事情正在发生（粗略地，你可以把它想象成“一个小的 2 ”），进行效力分析可能是一个好主意——了解小效应所需的样本量通常比你预期的要大得多，最好现在就了解这一点，而不是在你已经完成所有工作之后再做
• 如果你正在进行随机实验，你实际上可以控制样本量 - 你可以选择收集多少数据以及如何随机化处理。在进行实验之前，效力分析至关重要，这样你就不会在实验结束时才意识到“哦，该死，我应该再多几百人做一次......现在太晚了！”

效力分析的作用

让X表示处理，Y表示结果，假设我们正在研究X和Y关系，效力分析设计以下五件事：

• ① 处理效应的大小
• ② 处理的变动
• ③ 其它因素引起Y的变动
• ④ 统计精度
• ⑤ 样本量

效力分析就是保持其它四项不变，第五项是多少。例如，如果我们认为效应可能是 A，并且X的变动为B，并且Y中有 C 的变动与X无关，并且你希望至少有 D% 的机会发现一个效应（如果确实存在），那么你需要至少 E 的样本量。这告诉我们获得足够统计效力所需的最小样本量是E。

或者如果你有一个样本量为 A，并且X的变动为 B ，并且Y有 C 的变动与X无关，并且你希望至少有 D% 的机会发现一个效应（如果它确实存在），那么这个效应必须至少与 E 一样大。这告诉我们可检测到的效应的下界，即在给定样本大小的情况下，你希望有机会合理测量的最小效应。

那么，“统计精度”是什么呢？通常，你有一个目标统计效力水平（因此称为“效力分析”）。统计效力是真实率。也就是说，如果确实存在影响，并且抽样变动意味着你有 80% 的机会拒绝给定样本中无效的零假设，那么你就有 80% 的统计效力。统计效力还取决于你正在运行的测试类型 - 如果你在 95% 的置信水平下进行假设检验，则与在 99% 的置信水平下进行假设检验相比，你更有可能拒绝零假设（因此统计效力更高）。

效力分析并不一定要考虑统计效力。我们可以从以上五项的任何一个角度来进行效力分析，例如标准误。

这五项数值的来源

为了进行效力分析，我们需要知道五个数值中的四个，这样效力分析才能告诉我们第五个。

另一个问题是效力标准是多大？统计效力越高，我们错过真实处理效应的可能性越小，但这也意味着我们需要更大的样本。一般的经验来看，80%的效力是一个标准，但现在越来越多的研究者使用90%的效力。

效力分析的模拟

第一步：构造模拟数据集

我们模拟1000个样本，处理变量X是0,1上的均匀分布，处理效应是0.2，标准误是0-3的正态分布。stata代码如下：

* Set the number of observations we want in our simulated data
clear

* Create our data
* Note 0 and 1 are the default start and end of runiform; I want that, so I don't put anything inside runiform()
g X = runiform()

* Now create Y based on X
* Don't forget to add additional noise!
* The *true effect* of X on Y is .2
g Y = .2*X + rnormal(0, 3)

第二步：执行分析

假设我们计划获得回归中的稳健标准误：

* reg is short for regress
reg Y X, robust

* And if we're just interested in significance, say at the 95% level,
* we can compare the p-value to .05 and store the result as a local variable  (i.e. just store the single number, not as a regular Stata variable)
local sig = 2*ttail(e(df_r),abs(_b[X]/_se[X])) <= .05

* di is short for display
di `sig'

第三步：重复上述过程

forvalues i = 1(1)500 {
  * Set the number of observations we want in our simulated data
  clear
  set obs 1000
  
  * Create our data
  g X = runiform()
  g Y = .2*X + rnormal(0, 3)
  
  * Run analysis quietly
  qui reg Y X, robust
  
  * Pull out results
  local sig = 2*ttail(e(df_r),abs(_b[X]/_se[X])) <= .05
  di `sig'
  local coef = _b[X]
  di `coef'
}

这个重复过程将会产生500次数据，然后他将获取X的系数以及95%置信水平上是否显著。

第四步：储存结果

我们首先创建一个数据集来存储结果。然后，随着我们反复进行分析，我们将该preserve结果数据集清除，进行数据生成，并将结果存储在本地变量中。然后我们restore返回结果数据集并将这些本地变量放入实际变量中。最后，我们将拥有一个充满结果的数据集！

clear
* We want a number of observations equal to the number of times we will simulate the data
set obs 500

* Blank variables
g coef_results = .
g sig_results = .

* Let's wrap the whole thing in quietly{} - we don't need the output!
quietly {
  forvalues i = 1(1)500 {
    
    * Preserve our results data set so we can instantly come back to it
    preserve
    
    * Set the number of observations we want in our simulated data
    * Now we're no longer in the results data set; we're in the simulated-data data set
    clear
    set obs 1000
    
    * Create our data
    g X = runiform()
    g Y = .2*X + rnormal(0, 3)
    
    * Run analysis quietly
    qui reg Y X, robust
    
    * Pull out results
    local sig = 2*ttail(e(df_r),abs(_b[X]/_se[X])) <= .05
    local coef = _b[X]
    
    * Now restore to come back to results
    restore
    
    * And store our results in the results variables, using "in" to tell it which row to store the data in
    replace coef_results = `coef' in `i'
    replace sig_results = `sig' in `i'
  }
}

* summ is short for summarize
summ sig_results

结果如下：

结果显示，6%的统计效力（均值）。

第五步：更多的结果

我们还可以考察系数的分布：

tw kdensity coef_results, xti("Coefficient on X") yti("Density") lc(red)

image.png

X系数的密度分布

label def sig 0 "Not Significant" 1 "Significant"
label values sig_results sig
graph bar, over(sig_results) yti("Percentage")

image.png

第七步：高级效力分析

效力分析还可以获得更多的信息：（1）我们想要分析不同效应下的效力是多少。首先固定1000个样本，尝试分析不同效应的大小的效力水平：

foreach effect in .4 .8 1.2 1.6 2 {
  foreach sample_size in 1000 {
    quietly {
      clear
      * We want a number of observations equal to the number of times we will simulate the data
      set obs 500
      
      * Blank variables
      g sig_results = .
      
      * Let's wrap the whole thing in quietly{} - we don't need the output!
      forvalues i = 1(1)500 {
        
        preserve
        
        clear
        * NOTICE I'm putting the "sample_size" variable from above here
        set obs `sample_size'
        
        * Create our data
        g X = runiform()
        * and the "effect" variable here
        g Y = `effect'*X + rnormal(0, 3)
        
        * Run analysis quietly
        reg Y X, robust
        
        * Pull out results
        local sig = 2*ttail(e(df_r),abs(_b[X]/_se[X])) <= .05
        local coef = _b[X]
        
        * Now restore to come back to results
        restore
        
        * And store our results in the results variables, using "in" to tell it which row to store the data in
        replace sig_results = `sig' in `i'
      }
      
      summ sig_results
      * since we're inside a quietly{} we need a noisily to see this
      noisily di "With a sample size of `sample_size' and an effect of `effect', the mean of sig_results is " r(mean)
    }
  }
}

image.png

因此，看起来我们需要一个介于 0.8 和 1.2 之间的效应，才能有 80% 的机会找到显著的结果。如果我们认为效应实际上不太可能那么大，那么我们需要想其他办法---找到更大的样本，使用更精确的估计方法，等等！否则我们可能应该放弃。