应用Beta-二项分布模型及MCMC算法评估医院医疗事故数据

1、统计计算案例 2，吕晓玲应用 Beta-二项分布模型及 MCMC 算法评估医院医疗事故数据1. 问题提出与数据 外科手术直接影响着人民群众的身体健康和生命安全。 因此所有医生都须尽职尽责， 每 家医院都要严格控制事故发生率， 医疗管理部门也需要做好严格的监管工作。 表一给出了九 家医院的医疗事故数据， 哪家医院的手术更安全， 事故发生率更低呢？当然一个最最简单的 办法就是按照事故比例对医院进行排序。 不过我们都知道观测数据带有随机误差， 简单的按 照比例排序不能反映真实的情况。 为此我们考虑应用 Beta-二项分布模型及 MCMC 算法分析 此数据，希望可以考虑到数据的观测误差这个因素， 从而

2、对医院的医疗事故进行准确的评估。表一、九家医院的医疗事故数据医院医疗事故次数 (y i )外科病例数 (T i )11102010311245605156103207301528505109273202. 模型及估计方法2.1 Beta- 二项分布模型si 家医院进行重复观假设 is是第 i(i=1, ,N) 家医院外科手术的事故率。如果对于第测的次数 Ti 是已知的，那么 第 i 家医院发生外科手术事故的数量Yi 服从二项分布BinNom (Ti , is) :p(y i| is) ( yi) ( is)yi(1 - is)Ti-yi(1)如果概率 is 之间存在差异，处理这个问题的一种方法

3、是假设其服从一定的概率分布： is B( , ) (2)其中 ( , )是未知参数。那么 is 的均值和方差可以表示为：(3)E( is)Q var( is)( ) 2(1)(1 ) (1)(4)2.2 先验分布在贝叶斯估计中，我们要有一个, 的先验分布p( , ) 。我们假设 , 是相互独立的。对每个参数我们选取尺度参数为1但自由度不同的伽玛分布为, 的先验分布： (e ,1), (e ,1)(5)这说明 E( ) e ,E( ) e 。我们可以运用以下结果选择 e,e：可以得出，e ,e 越大，会使得 (4) 式的方差的系数越小。而 e /(e e ) 这一比例又直接对总体均值造成很大影响

4、。如果选取 e ,e 使得 e (e e ) ，则总体均值将会趋于 0 。2.3 MCMC算法进行统计推断在 给 定 观 测 数 据 y ( y1,., yN ) 的 情 况 下 ， 统 计 推 断是 根 据 后 验 分 布 s s s sp( , , 1s ,., Ns |y)进行 ( , , 1s,., Ns )的估计。在此对后验分布应用 MCMC方法进行抽 样模拟。从 (0), (0) 开始，运行一个三阶段抽样，主要步骤如下：a)从分布密度p( is |(m 1)y, ,(m1) 中独立抽取 is；b)从 p(| y,(s )(m)1,., ( s )(m) N, (m 1)中抽取 (m

5、)；c)从 p(| y,(s)(m)1),., ( s )(m),., ( N ), (m)中抽取 (m) 。2.3.1 联合后验分布根据贝叶斯估计方法可以知道，联合后验分布ssp( , , 1 ,., N |y) 由下式给出：p( ,s1, , sN|y) P(y|s1,sN)p(s1,sN|,)p( )p( ) (6)从式 (2) 中可以得到先验分布 p( 1s ,., Ns | , ) 如下：p( s1 , , sN | ,) iN=1 B(, ()is ) ( -1) (1 - is) (-1)(7)由式(1) 可以得到似然函数 p(y| 1s,., Ns )如下： p(y| s1,s

6、N) iN=1 ( is )y i (1 - is )(T i -y i)(8)将( 7)和( 8)代入( 6)，得p( ,s1, sN |y) B(1,N)p( )p( )iN=1(is)yi(1- is) (T i -y i)(is)(-1)(1 - is)(-1)(9)其中 B(,) =()( )( +)2.3.2 步骤 a：对独立的概率值进行抽样已知。根首先需要从分布 p( 1s ,., Ns | y, , )中分别抽取概率值 is ，其中假设据式(9) ，我们可以得到如下关系：p( 1s,sN|(10)p( is | yi ,y,) iN=1(is)(yi+-1)(1 - is )(

7、T i -y i+-1)(i=1 ， N)，) 为给定, 后的后验分布，其中is 相互独立。对任意 i后验分布密度p( is |yi , ,)为 B( yi, Tiyi ) 。该结论在0,0的情况下对任意 yi 适用。2.3.3 步骤 b：对参数 抽样 根据式 (9) ，我们得到 一个不是很熟悉的后验分布：( +) Np( | ,s1 , , sN, y) ( () ) p ( )( iN=1 (is)(11)采用 Metropolis-Hasting 算法中的对数正态随机游走 Metropolis-Hasting 算法对问题进 行处理。假设 log new N(log old,c ) ，然后

8、以下述概率接受新值：min(1,( (newold)()(old N new)N( new) i 1is)new old p(p(new) old )newold(12)接受新值与否可通过如下过程判断： 抽取服从均匀分布的随机数 Um U 0,1 ，当U m小于式(12) 给出的概率值时，接受新值，否则不接受新值。为了简化运算，提高运算速度，可以 将数据进行对数处理后再进行比较，即当：LogUmN(lognew) log ( old ) log ( old ) log ( new)newNold )i1lognew oldlog p( ) log p( )lognewlogold13)此时，选

9、择接受新值new。2.3.4 步骤 c ：对参数 抽样根据式 (9) ，我们得到一个不是很熟悉的 的后验分布：( +) Np( | ,s1 , , sN, y) ( () ) p()(iN=1(1 - is) (14)采用 Metropolis-Hasting算法中的对数正态随机游走Metropolis-Hasting算法对问题进行处理。假设 lognewN(log old ,c） ，抽取服从均匀分布的随机数Vm U 0,1 ，当LogVm N(log(new) log (old) log( old ) log (new)N( new old )log(1is) log p(new) log

10、p(old new) loglogold(15)此时，选择接受新值 new3. 数据分析3.1 第一步：为, 选择合适的先验分布我们通过数据y i/Ti 的均值 m和方差 s确定, 先验分布的两个参数 e， e：令 m = E( )= ee+e ，s= Q= ( 1- ) e+e e+e1+1可以得到： e= m(m 1-ms - 1) = E( )+E( )1-m1.87 ， e= e m = 17.28 。3.2 第二步： MCMC 抽样在使用 Metropolis-Hasting 算法对 , 进行抽样时，使用不同的变量 c ,c 可以求得不同的接受概率，从而确定使接受概率最佳（介于 0.

11、40.7 之间）的 c ,c 。在这里 c ,c 的取值从 0.1到 1，可以得到 , 的接受率分别如下：的接受概率如下：c.beta=0.10.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0c.alpha=0.1 0.5936 0.5868 0.5884 0.5866 0.5932 0.5932 0.5880 0.5898 0.5850 0.59120.2 0.4762 0.4778 0.4822 0.4834 0.4830 0.4834 0.4824 0.4818 0.4778 0.4868 0.3 0.4202 0.4132 0.4148 0.4110 0.4178

12、0.4246 0.4050 0.4016 0.4058 0.41620.4 0.3826 0.3672 0.3876 0.3746 0.3772 0.3688 0.3780 0.3512 0.3734 0.38600.5 0.3334 0.3460 0.3420 0.3658 0.3374 0.3292 0.3436 0.3468 0.3436 0.3568 0.6 0.3100 0.3248 0.3176 0.3132 0.3218 0.3066 0.3096 0.3198 0.3044 0.32480.7 0.2934 0.2966 0.2998 0.3030 0.2864 0.2922

13、0.2904 0.2962 0.2874 0.3000 0.8 0.2900 0.2734 0.2790 0.2828 0.2800 0.2834 0.2830 0.2776 0.2726 0.28480.9 0.2764 0.2582 0.2580 0.2792 0.2674 0.2630 0.2622 0.2694 0.2722 0.26261.0 0.2618 0.2532 0.2538 0.2626 0.2538 0.2436 0.2456 0.2508 0.2534 0.2616的接受概率如下：c.beta=0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.

14、0c.alpha=0.10.5354 0.4186 0.3448 0.3200 0.2706 0.2684 0.2468 0.2374 0.2288 0.20800.2 0.5326 0.4246 0.3546 0.3164 0.2934 0.2668 0.2554 0.2252 0.2208 0.2100 0.3 0.5182 0.4294 0.3576 0.3176 0.2926 0.2706 0.2490 0.2374 0.2084 0.21040.4 0.5274 0.4202 0.3462 0.3308 0.2808 0.2656 0.2538 0.2384 0.2184 0.210

15、40.5 0.5278 0.4192 0.3636 0.3286 0.2932 0.2776 0.2588 0.2306 0.2212 0.20560.6 0.5214 0.4080 0.3512 0.3116 0.2912 0.2680 0.2428 0.2342 0.2292 0.20820.7 0.5156 0.4184 0.3522 0.3168 0.2778 0.2678 0.2460 0.2376 0.2120 0.20880.8 0.5338 0.4256 0.3634 0.3156 0.2792 0.2538 0.2604 0.2312 0.2216 0.20520.9 0.5

16、130 0.4134 0.3576 0.3232 0.2968 0.2588 0.2488 0.2418 0.2106 0.21081.0 0.5234 0.4174 0.3524 0.3172 0.2710 0.2614 0.2402 0.2352 0.2312 0.2088选定 c 0.3,c0.2，这样可以使得和的接受概率 在 40-43% 。并且选取初始值0 50, 0 100，使用 MCMC 抽样得到 (m), (m) 的序列如图 1。图 1： (m), (m)的序列(初始值?0? = 50, 0 = 100)从图 1 中可以看出，序列很快就趋于收敛了，不妨认为从第501 个抽样开始

17、趋于收敛。去除数据趋于收敛前的样本，然后计算 , 的估计值，得到的结果如下：的均值 = 2.04的方差 = 0.36的中位数 = 1.97的百分之九十五的置信区间 = 1.09 ， 3.39的均值 = 17.93的方差 variance = 15.34的中位数 median = 17.49 的百分之九十五的置信区间 = 11.29 ，26.51画出 , 的后验分布直方图，如图 2。图 2： 和 的分布3.3 第三步：分析不可观测的异质性根据 (3) 和 (4) 式，通过 MCMC 抽样得到的后验分布，得到的估计值如下：(m), (m)的值计算均值 以及方差 Q 的的均值 = 0.1031的方差

18、 = 0.0006的中位数 Median = 0.1009 的百分之九十五的置信区间 = 0.0614 ， 0.1555Q的均值 = 0.0046Q的方差 = 0.0000026Q的中位数 = 0.0044Q的百分之九十五的置信区间 = 0.0024 ， 0.0080 画出 和 Q 的后验分布直方图，如图3。图 3 ： u 和 Q 的分布s3.4 第四步：根据各个医院的事故率is对医院进行排序根据 MCMC 抽样中各个医院的概率值 ( is) (m) ，计算 is的中位数，并根据此值对医院进 行排序，结果如表二。画出 ( is) (m) 的盒型图，如图 4。表二、九家医院的医疗事故数据医院医疗事故次数 (y i )外科病例数 (Ti)均值 (yi/Ti)按均值排名后验分布中 位数按后验分布 中位数排名11100.170.091662010010.0579231120.08330.0868545600.08340.085445150.290.110586103200.03120.034717301520.1978