Stata：二维估计的可视化-confcomptwo

发布时间：2023-09-15 阅读 359

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

作者：雷诺 (新加坡国立大学大学)
邮箱：leinuo@u.nus.edu

编者按：本文主要摘译自下文，特此致谢！
Source：Eckert, M., & Vach, W. (2023). Visualizing uncertainty in a two-dimensional estimate using confidence and comparison regions. The Stata Journal, 23(2), 455–490. -Link- -PDF-

1. 引言

1.1 二维置信区域的含义

在统计分析中，我们已经习惯使用置信区间 (confidence intervals) 来描述单参数估计的不确定性。即使在需要同时考虑多个参数估计时，通常也会将不确定性分别附加在每个单参数估计中。然而，通过使用置信区域 (confidence regions)，我们能够分析多元回归中两个回归系数估计的联合不确定性。

以图 1 为例，首先，我们发现握力 (grip strength) 与年龄 (age) 两个回归系数之间存在负相关关系。在解释这两个估计的大小时，需要考虑到对一个参数过高估计可能意味着对另一个参数的低估。

其次，值得注意的是， $(0, 0)$ 点位于置信区域之外，这意味着我们可以拒绝两个协变量均无效的零假设。当两个回归参数都与 0 的显著性存在差异时，置信区域提醒我们不应错误地将结果解释为这两个协变量之间没有任何关联。这种方法有助于更全面地理解参数估计的不确定性以及参数之间的关系。

1.2 比较区域的必要性

置信区域可以用于将参数固定在某个值的零假设的事后检验。如对于某一 $θ_{0} \in R^{2}$ ，其零假设为 $H_{0} : θ = θ_{0}$ 。如果点 $θ_{0}$ 没有被置信区域 $(1 - α)$ 覆盖，我们则在显著性水平 $α$ 下拒绝零假设。

然而，在某些统计应用中，我们关注的是两个参数的相关假设，如它们的平均值是否超过了某个阈值。此时，我们关注的是零假设 $H_{0} : θ \notin R$ 。当然，在事后检验中我们可以检查 $(1 - α)$ 置信区域是否完全覆盖了我们关心的区域。若成立，则可以在显著性水平 $α$ 下拒绝零假设。然而，这种测试方法非常保守，因为实际检验的显著性水平远低于 $α$ 。

为了解决这个问题，Eckert & Vach (2019) 提出了比较区域 (comparison region) 的概念。显著性水平 $α$ 下的比较区域是一个依赖数据的区域 $C \subseteq R^{2}$ ，并具有如下特征：

$I (C) \subseteq R$ 定义了一个在显著性水平 $α$ 下的零假设 $H_{0} : θ \notin R$ , 并对于任何凸集 $R \in R^{2}$ 成立

因此，比较区域可以用于对我们所关心的零假设进行事后检验 (图 2)。同时，比较区域的构建方法与置信区域非常相似。

2. 命令介绍

confcomptwo 命令是 Eckert & Vach (2023) 基于上述理论所编写的 Stata 的新命令。

命令安装：

net install st0716.pkg, replace // 安装命令
net get st0716.pkg, replace     // 下载文件

命令语法：

使用 confcomptwo 作为估计后命令并应用 Wald test 原理。

confcomptwo parname1 parname2 [, confcomp_options twoway_options]

使用 concomptwo 的直接形式以应用 Wald test 原理。其中 \#1 和 \#2 表示两个参数估计值，se1() 和 se2() 提供参数估计的标准误，默认相关性为 corr(0)。

confcomptwo #1 #2, se1(real) se2(real) [corr(real) confcomp_options twoway_options]

应用 likelihood-ratio (LR)

confcomptwo #1 #2, call(expr) [confcomp_options twoway_options linesearch_options]

call() 中的表达式描述了调用一个用户编写的程序，该程序会计算某种统计模型的对数似然值 (log-likelihood)，并将这个对数似然值存储在变量 r(ll) 中。

其中，confcomp_options 包括：

parname1 parname2：先前拟合的模型的参数名；
alpha1 (real)：定义比较区域的显著性水平或 1 减去置信区域的显著性水平；
points(numlist)：定义要绘制的点；
rescale(expr)：允许修改方向的选择；
loptsconf(line_options)：指定用于绘制置信区域的线的外观；
loptscomp(line_options)：指定用于绘制比较区域的线条的外观；
mopts(marker_options)：指定用于绘制点估计值的标记的外观；
only(conf|comp)：抑制绘制比较区域或置信区域；
nograph：禁止绘制图形；
savepoints(filename[, saveoptions])：将计算的点数保存为 Stata 文件；
reverse：交换 $y$ 轴和 $x$ 轴；
addplot(plot ... [|| plot ... [...]][, below])：向图中添加双向图。

其中，twoway_options 包括：

added_line_options：在指定的 $y$ 或 $x$ 值处绘制线条；
added_line_options: 以指定的 $(y, x)$ 值显示文本；
axis_options：标签，刻度，网格，对数刻度；
title_options: 标题、副标题、注释、说明；
legend_options：解释图例；
scale(#)：调整文本、标记和行宽的大小；
region_options：轮廓，阴影，图形大小；
aspect_option：限制纵横比；
scheme(schemename)：整体外观；
play(recordingname)：从 recordingname 展示编辑的内容；
by(varlist, ...)：重复子组；
nodraw：抑制图形显示；
name(name, ...)：为图形指定名称；
saving(filename, ...)：将图形保存到文件中。

3. Stata 案例

3.1 案例1：两个回归系数的联合置信区域

首先，我们运用 Stata 提供的 auto.dta 来研究里程等级与汽车重量和原产地 (国外或国内) 的关系。

. sysuse auto, clear
. regress mpg weight foreign

  Source |       SS          df       MS      Number of obs =     74
---------+---------------------------------   F(2, 71)      =  69.75
   Model |   1619.2877        2  809.643849   Prob > F      = 0.0000
Residual |  824.171761       71   11.608053   R-squared     = 0.6627
---------+---------------------------------   Adj R-squared = 0.6532
   Total |  2443.45946       73  33.4720474   Root MSE      = 3.4071
--------------------------------------------------------------------
     mpg | Coefficient Std. err.      t    P>|t| [95% conf. interval]
---------+----------------------------------------------------------
  weight |     -0.007     0.001   -10.34   0.000   -0.008     -0.005
 foreign |     -1.650     1.076    -1.53   0.130   -3.796      0.495
   _cons |     41.680     2.166    19.25   0.000   37.362     45.998
---------------------------------------------------------------------

. confcomptwo weight foreign, only(conf) rescale(1666) legend(cols(1))  ///
>     loptsconf(lpattern(solid)) ytitle({&beta}(weight)) xtitle({&beta} ///
>     (foreigh)) aspect(1) xsize(4) ysize(4.5)

上图中两个回归系数呈正相关关系，这是因为两个协变量之间负相关：平均而言，外国汽车的重量低于国内汽车。

3.2 案例2：配对诊断准确性研究分析

Ng et al., (2008) 对两种诊断工具进行了比较，这两种工具分别是 18F-Fluor-2- 脱氧葡萄糖正电子发射断层扫描和远场多探测器计算机断层扫描，用于检测口咽癌或下咽鳞状细胞癌患者的远处恶性肿瘤。研究共招募了 160 名患者，这些患者接受了至少 12 个月的随访或直到死亡。

该研究采用了配对设计，即每位患者都分别使用两种诊断工具进行了测试。研究结果呈现在下面的 2×2×2 列表中，该列表包含了要比较的两种指数检验的每种结果组合的频率和参考标准。

. use studyng, clear
. list

     +---------------------------------+
     | test1   test2   refere~e   freq |
     |---------------------------------|
  1. |     1       1          1     12 |
  2. |     1       1          0      2 |
  3. |     1       0          1      1 |
  4. |     1       0          0      1 |
  5. |     0       1          1      8 |
     |---------------------------------|
  6. |     0       1          0      6 |
  7. |     0       0          1      5 |
  8. |     0       0          0    125 |
     +---------------------------------+

首先，我们探究将 test 1 替换为 test 2 对灵敏度和特异性的影响。为了基于 Wald 检验计算区域，我们使用了 confcomptwo 的更新版本进行估计。随后，绘制了比较区域和置信区域，以图形化方式展示了估计值和感兴趣差异的不确定性。

. generate correct1 = test1 == reference
. generate correct2 = test2 == reference
. /* 根据测试变量 test1 和 test2 是否等于参考值 reference 
> 来生成两个新的变量 correct1 和 correct2，如果相等则为 1，否则为 0。*/

. quietly glm correct1 if reference == 1 [fw=freq], family(bernoulli) link(id)
. /* 执行广义线性模型 (GLM) 对 correct1 和 correct2 进行回归分析，
> 筛选条件是 reference 等于 1。其中 quietly 表示屏蔽输出，
> glm 表示拟合广义线性模型，family(bernoulli) 表示使用伯努利分布，
> link(id) 表示使用恒等链接函数。*/
. estimates store sens1 

. quietly glm correct2 if reference == 1 [fw=freq], family(bernoulli) link(id)
. estimates store sens2
. quietly glm correct1 if reference == 0 [fw=freq], family(bernoulli) link(id)
. estimates store spec1
. quietly glm correct2 if reference == 0 [fw=freq], family(bernoulli) link(id)
. estimates store spec2
. suest sens1 sens2 spec1 spec2

Simultaneous results for sens1, sens2, spec1, spec2
                                                             Number of obs = 8
--------------------------------------------------------------------------------
               |               Robust
               | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
sens1_correct1 |
         _cons |      0.500      0.098     5.08   0.000        0.307       0.693
---------------+----------------------------------------------------------------
sens2_correct2 |
         _cons |      0.769      0.083     9.28   0.000        0.607       0.932
---------------+----------------------------------------------------------------
spec1_correct1 |
         _cons |      0.978      0.013    76.25   0.000        0.952       1.003
---------------+----------------------------------------------------------------
spec2_correct2 |
         _cons |      0.940      0.021    45.80   0.000        0.900       0.981
------------------------------------------------------------------------------

. nlcom (deltasens:[sens2_correct2]_cons-[sens1_correct1]_cons) ///
>     (deltaspec:[spec2_correct2]_cons-[spec1_correct1]_cons), post
. /* 利用 nlcom 计算非线性组合 (nonlinear combination) 的估计和置信区间。
> 其中前半部分计算两个敏感性 (sensitivity) 估计值之间的差异，
> 后半部分计算两个特异性 (specificity) 估计值之间的差异。
> post 表示将计算的结果存储为后续可以使用的变量。*/

   deltasens: [sens2_correct2]_cons-[sens1_correct1]_cons
   deltaspec: [spec2_correct2]_cons-[spec1_correct1]_cons
---------------------------------------------------------------------
          | Coefficient Std. err.    z    P>|z|  [95% conf. interval]
----------+----------------------------------------------------------
deltasens |      0.269     0.103   2.62   0.009     0.068       0.471
deltaspec |     -0.037     0.020  -1.91   0.056    -0.076       0.001
----------------------------------------------------------------------

. confcomptwo deltasens deltaspec,                                 ///
>     loptsconf(points, n(100) yxratio(2) msize(*0.2) mcol(black)) ///
>     xtitle({&Delta} specificity) ytitle({&Delta} sensitivity)    ///
>     addplot((scatteri 0.15 -0.15 -0.1 0.1, connect(1) msymbol(i) ///
>     lpattern(dot) lcolor(gs7))                                   ///
>     (scatteri 0.225 -0.15 -0.1 0.066666,                         ///
>     connect(1) msymbol(i) lpattern(shortdash) lcolor(gs7))       ///
>     (scatteri 0.3 -0.15 -0.1 0.05,                               ///
>     connect(1) msymbol(i) lpattern(longdash) lcolor(gs7)))       ///
>     legend(cols(1) order(1 2 3))                                 ///
>     xline(0, lcolor(gs14)) yline(0, lcolor(gs14))                ///
>     xlabel (-.15(.1).15) ylabel(-.1(.1).5) aspect(2) xsize(2.5) ysize(4)

在上图中，我们可以观察到灵敏度有所提高，而特异性略微下降。由于所研究疾病的患病率较低，因此对于灵敏度变化的估计精度相对较低。我们还可以看到三条额外的线，它们分别代表了 Wald test 的 5% 比较，以及敏感性和特异性变化的 95% 置信区域。这些线条用于表示灵敏度和特异性变化的不同加权平均值。研究结果表明，灵敏度的提高程度超过了特异性的下降程度，这意味着灵敏度得到了显著的改善。

接下来，计算二元分类模型中的真阳性 (True Positive)、假阳性 (False Positive) 等统计量，并进行估计和比较。

. use studyng, clear
. generate TP1 = test1 & reference == 1
. generate TP2 = test2 & reference == 1
. /* 根据测试变量 test1 和 test2 是否为真阳性生成两个新的变量 TP1 和 TP2。
> 即在 test 为 1 且 reference 为 1 的情况下，将对应的变量值设置为 1，否则设置为 0。*/

. generate FP1 = test1 & reference == 0
. generate FP2 = test2 & reference == 0
. /* 同理，生成两个新的变量 FP1 和 FP2，用来表示假阳性情况。*/

. quietly glm TP1 [fw=freq], family(bernoulli) link(id)
. estimates store TP1
. quietly glm TP2 [fw=freq], family(bernoulli) link(id)
. estimates store TP2
. quietly glm FP1 [fw=freq], family(bernoulli) link(id)
. estimates store FP1
. quietly glm FP2 [fw=freq], family(bernoulli) link(id)
. estimates store FP2
. suest TP1 TP2 FP1 FP2

Simultaneous results for TP1, TP2, FP1, FP2            Number of obs = 8
------------------------------------------------------------------------
             |               Robust
             | Coefficient  std. err.   z    P>|z|  [95% conf. interval]
-------------+----------------------------------------------------------
TP1_TP1      |
       _cons |      0.081      0.022  3.75   0.000     0.039       0.124
-------------+------------------------------------------------------------
TP2_TP2      |
       _cons |      0.125      0.026  4.77   0.000     0.074       0.176
-------------+----------------------------------------------------------
FP1_FP1      |
       _cons |      0.019      0.011  1.74   0.081    -0.002       0.040
-------------+----------------------------------------------------------
FP2_FP2      |
       _cons |      0.050      0.017  2.89   0.004     0.016       0.084
------------------------------------------------------------------------

. nlcom (deltaTP:[TP2_TP2]_cons-[TP1_TP1]_cons) ///
>     (deltaFP:[FP2_FP2]_cons-[FP1_FP1]_cons), post
. /* 使用 nlcom 命令，计算两个真阳性估计值之间的差异和两个假阳性估计值之间的差异，
> 然后将这些差异结果存储为后续可以使用的变量。*/

     deltaTP: [TP2_TP2]_cons-[TP1_TP1]_cons
     deltaFP: [FP2_FP2]_cons-[FP1_FP1]_cons
--------------------------------------------------------------------
        | Coefficient  Std. err.     z    P>|z| [95% conf. interval]
--------+-----------------------------------------------------------
deltaTP |      0.044      0.018    2.37   0.018    0.008       0.080
deltaFP |      0.031      0.016    1.91   0.057   -0.001       0.063
--------------------------------------------------------------------

. confcomptwo deltaTP deltaFP, xtitle({&Delta} FP) ytitle({&Delta} TP) ///
>    loptsconf(points, n(250) yxratio(2) msize(*0.1) mcol(black))      ///
>    addplot((scatteri -0.01 -0.01 0.08 0.08, connect(line) msymbol(i) ///
>     lpattern(dot) lcolor(gs7))                                       ///
>     (scatteri -0.01 -0.02 0.04 0.08, connect(line) msymbol(i)        ///
>     lpattern(shortdash) lcolor(gs7))                                 ///
>     (scatteri -0.0067 -0.02 0.0266 0.08, connect(line) msymbol(i)    ///
>     lpattern(longdash) lcolor(gs7)))                                 ///
>    legend(cols(1) order(1 2 3))                                      ///
>    xline(0, lcolor(gs14)) yline(0, lcolor(gs14))                     ///
>    legend(col(1) order(1 2 3))                                       ///
>    xlabel(-.02(.02).08) ylabel(0(.02).1)                             ///
>    yscale(range(-0.005 0.1))                                         ///
>    xsize(4) ysize(5)

从上图中可知，真阳性 (TP) 决策的频率上升了 $4.4 ％$ ，相对于假阳性 (FP) 决策的增加幅度更高。这种情况在疾病状态的患病率较低时，灵敏度和特异性的变化显得尤为显著。在决策制定中，不同的利益相关者对于在获取 TP 决策的同时接受 FP 决策的容忍程度存在差异。然而，新的测试方法仅在需要在获得一个 TP 决策之前接受三个 FP 决策的情况下才显示出优势。这强调了在低患病率情况下，灵敏度和特异性变化对于权衡决策的重要性。

4. 参考资料

Eckert, M., & Vach, W. (2023). Visualizing uncertainty in a two-dimensional estimate using confidence and comparison regions. The Stata Journal: Promoting Communications on Statistics and Stata, 23(2), 455–490. doi:10.1177/1536867x231175314. -PDF-
Eckert, M., & Vach, W. (2019). On the use of comparison regions in visualizing stochastic uncertainty in some two‐parameter estimation problems. Biometrical Journal, 62(3), 598–609. doi:10.1002/bimj.201800232. -PDF-
Ng, S.-H. et al. (2008) ‘Distant metastases and synchronous second primary tumors in patients with newly diagnosed oropharyngeal and hypopharyngeal carcinomas: Evaluation of 18F-FDG PET and extended-field multi-detector row CT’, Neuroradiology, 50(11), pp. 969–979. doi:10.1007/s00234-008-0426-2. -PDF-

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 置信区间系数可视化, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

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

Stata：二维估计的可视化-confcomptwo

1. 引言

1.1 二维置信区域的含义

1.2 比较区域的必要性

2. 命令介绍

3. Stata 案例

3.1 案例1：两个回归系数的联合置信区域

3.2 案例2：配对诊断准确性研究分析

4. 参考资料

5. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯