shapleyx 命令如何使用 bootstrap 方法得到置信区间

最近有个小伙伴遇到了这样一个问题，他想使用 shapleyx 命令，但是这个命令的结果里面没有置信区间：

sysuse auto, clear
replace price = price/1000
shapleyx weight i.foreign i.rep78 mpg length, result(e(r2)):regress price @

*>     Shapley decomposition
*>
*>  Factors | 1st round |  Shapley
*>          |  effects  |   value
*> ---------+-----------+-----------
*> weight   |    .2901  |    .2473
*> foreign  |  .002374  |  .092834
*> rep78    |  .014495  |  .025609
*> mpg      |   .21958  |  .086631
*> length   |   .18648  |   .11039
*> ---------+-----------+-----------
*> Residual |  -.15027  |
*> ---------+-----------+-----------
*>    Total |   .56276  |   .56276

不过在 shapleyx 命令的帮助文档里面有提到可以使用 bootstrap 计算置信区间：

No standard errors for factor contributions have been provided so far. A reasonable thing to do would be to bootstrap your data; see bootstrap. Doing so, however, would result in even more computationally intensive calculations…

但是 shapleyx 命令又没有提供 scalar 类型的返回值：

ret list

*> macros:
*>             r(factors) : "(weight) (i.foreign) (i.rep78) (mpg) (length)"
*>               r(names) : "weight foreign rep78 mpg length"

*> matrices:
*>            r(decompos) :  6 x 2

mat list r(decompos)

*> r(decompos)[6,2]
*>           OneStage    Shapley
*>  weight   .2901023  .24730103
*> foreign  .00237359  .09283407
*>   rep78  .01449477  .02560857
*>     mpg  .21958286  .08663138
*>  length  .18647823  .11038516
*>   total  .71303175  .56276021

对于这种情况，我们可以考虑自编一个带 r 类返回值的命令：

cap prog drop myprog
prog def myprog, rclass
 shapleyx weight i.foreign i.rep78 mpg length, result(e(r2)):regress price @
 ret scalar weight = r(decompos)[1,2]
 ret scalar foreign = r(decompos)[2,2]
 ret scalar rep78 = r(decompos)[3,2]
 ret scalar mpg = r(decompos)[4,2]
 ret scalar length = r(decompos)[5,2]
 ret scalar total = r(decompos)[6,2]
end

然后对这个命令应用 bootstrap：

bootstrap weight = r(weight) foreign = r(foreign) ///
 rep78 = r(rep78) mpg = r(mpg) length = r(length) ///
 total = r(total): myprog

*> (running myprog on estimation sample)
*>
*> warning: myprog does not set e(sample), so no observations will be excluded from the resampling
*>          because of missing values or other reasons. To exclude observations, press Break, save
*>          the data, drop any observations that are to be excluded, and rerun bootstrap.
*>
*> Bootstrap replications (50)
*> ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
*> ..................................................    50
*>
*> Bootstrap results                                           Number of obs = 74
*>                                                             Replications  = 50
*>
*>       Command: myprog
*>        weight: r(weight)
*>       foreign: r(foreign)
*>         rep78: r(rep78)
*>           mpg: r(mpg)
*>        length: r(length)
*>         total: r(total)
*>
*> ------------------------------------------------------------------------------
*>              |   Observed   Bootstrap                         Normal-based
*>              | coefficient  std. err.      z    P>|z|     [95% conf. interval]
*> -------------+----------------------------------------------------------------
*>       weight |    .247301   .0690901     3.58   0.000     .1118869    .3827152
*>      foreign |   .0928341   .0337649     2.75   0.006      .026656    .1590121
*>        rep78 |   .0256086   .0253655     1.01   0.313     -.024107    .0753241
*>          mpg |   .0866314   .0302501     2.86   0.004     .0273423    .1459204
*>       length |   .1103852   .0297727     3.71   0.000     .0520318    .1687385
*>        total |   .5627602   .0941743     5.98   0.000     .3781819    .7473385
*> ------------------------------------------------------------------------------

这样我们就解决了这个问题。

点击这里跳转到 RStata 短书平台获取附件：shapleyx 命令如何使用 bootstrap 方法得到置信区间