(完整版)應(yīng)用Beta-二項分布模型及MCMC算法評估醫(yī)院醫(yī)療事故數(shù)據(jù)

1、統(tǒng)計計算案例2，呂曉玲應(yīng)用 Beta-二項分布模型及MCMC 算法評估醫(yī)院醫(yī)療事故數(shù)據(jù)1. 問題提出與數(shù)據(jù)外科手術(shù)直接影響著人民群眾的身體健康和生命安全。因此所有醫(yī)生都須盡職盡責(zé)，每家醫(yī)院都要嚴格控制事故發(fā)生率，醫(yī)療管理部門也需要做好嚴格的監(jiān)管工作。表一給出了九家醫(yī)院的醫(yī)療事故數(shù)據(jù)，哪家醫(yī)院的手術(shù)更安全，事故發(fā)生率更低呢？當(dāng)然一個最最簡單的辦法就是按照事故比例對醫(yī)院進行排序。不過我們都知道觀測數(shù)據(jù)帶有隨機誤差，簡單的按照比例排序不能反映真實的情況。為此我們考慮應(yīng)用 Beta-二項分布模型及MCMC 算法分析此數(shù)據(jù)，希望可以考慮到數(shù)據(jù)的觀測誤差這個因素，從而對醫(yī)院的醫(yī)療事故進行準確的評估。表一、

2、九家醫(yī)院的醫(yī)療事故數(shù)據(jù)醫(yī)院醫(yī)療事故次數(shù) (y i )外科病例數(shù) (T i )11102010311245605156103207301528505109273202. 模型及估計方法2.1 Beta-二項分布模型s,N) 家醫(yī)院外科手術(shù)的事故率。如果對于第i 家醫(yī)院進行重復(fù)觀假設(shè)i 是第 i(i=1,測的次數(shù) Ti是已知的，那么第 i家醫(yī)院發(fā)生外科手術(shù)事故的數(shù)量Yi 服從二項分布BinNom (T ,s) :iisTisy i(1sT i -y i(1)p(y i | i) (yi) ( i)- i )如果概率is 之間存在差異，處理這個問題的一種方法是假設(shè)其服從一定的概率分布：is B( ,

3、 )(2)其中( ,) 是未知參數(shù)。那么is 的均值和方差可以表示為：E( is )(3)1Qvar(is )() 2 (1)(1 )1(4)(1)2.2先驗分布在貝葉斯估計中，我們要有一個, 的先驗分布p( , ) 。我們假設(shè),是相互獨立的。對每個參數(shù)我們選取尺度參數(shù)為1 但自由度不同的伽瑪分布為,的先驗分布： (e,1), (e,1)(5)這說明 E()e , E()e。我們可以運用以下結(jié)果選擇e ,e：可以得出，e , e 越大，會使得(4)式的方差的系數(shù)越小。而 e/( ee ) 這一比例又直接對總體均值造成很大影響。如果選取e , e 使得e (e e ) ，則總體均值將會趨于 0。

4、2.3 MCMC算法進行統(tǒng)計推斷在 給 定 觀 測 數(shù) 據(jù) y( y1,., yN ) 的 情 況 下 ， 統(tǒng) 計 推 斷 是 根 據(jù) 后 驗 分 布p(,1s ,.,Ns | y) 進行 ( ,1s ,.,Ns ) 的估計。在此對后驗分布應(yīng)用MCMC方法進行抽樣模擬。從 (0 ) , ( 0) 開始，運行一個三階段抽樣，主要步驟如下：（a）從分布密度p(s |,( m1) , ( m 1) )中獨立抽取s；iyi（）從(|, (s )( m),., (s) (m ) ,( m 1) )中抽取( m)bpy1N；（c）從(|, (s ) (m),., (s) (m) ,( m) )中抽取( m

5、)py1N。2.3.1聯(lián)合后驗分布根據(jù)貝葉斯估計方法可以知道，聯(lián)合后驗分布p(, 1s,.,Ns | y) 由下式給出：ssp( , , |y)1N(6)ssssP( y|, , )p( , , |,)p( )p( )1N1N從式 (2)中可以得到先驗分布p( 1s ,.,Ns |, )如下：ssN1s( -1)(1 -s( -1)p( 1, , N | ,) i=1( i)i)B( , )2(7)由式 (1)可以得到似然函數(shù) p( y |1s ,., Ns ) 如下：ssNss(T i -y i )p(y| , ) ( )y i (1 - )1Ni=1ii(8)將（ 7）和（ 8）代入（ 6

6、），得ssp( , , |y)1N(9)1p( )p(Nss(T i -y i )s( -1)(1 -s( -1)N)( ) y i (1 -)( ) )B( i=1iiii, )()() ( )其中 B , = ( +)2.3.2 步驟 a：對獨立的概率值進行抽樣：首先需要從分布p(s,.,s| ,) 中分別抽取概率值s, 已知。根1Nyi ，其中假設(shè)據(jù)式 (9)，我們可以得到如下關(guān)系：ssN(s(y i + -1)(1s(T i -y i+ -1)p( , , |,y) )- )1Ni=1ii(10)p( is | yi ,)為給定,后的后驗分布，其中is 相互獨立。對任意i （ i=1

7、， ， N），后驗分布密度p(is | yi ,)為 B(yi ,Tiyi ) 。該結(jié)論在0,0 的情況下對任意 yi適用。2.3.3 步驟 b：對參數(shù)抽樣根據(jù)式 (9)，我們得到一個不是很熟悉的后驗分布：ss()NNs( +)（ 11）p( |, , , y)p( )( ( )1N( )i=1i采用 Metropolis-Hasting算法中的對數(shù)正態(tài)隨機游走Metropolis-Hasting算法對問題進行處理。假設(shè)lognew N (logold , c ) ，然后以下述概率接受新值：min(1, (newold)(old)N (Nis )newold p(new)new(12)(new

8、 )i 1p(old )old接受新值與否可通過如下過程判斷：抽取服從均勻分布的隨機數(shù)U m U 0,1 ，當(dāng) U m 小于式(12) 給出的概率值時，接受新值，否則不接受新值。為了簡化運算，提高運算速度，可以將數(shù)據(jù)進行對數(shù)處理后再進行比較，即當(dāng)：3LogU mN (log(new) log(old)log( old )log (new )N( newold )logislog p(new )log p(old )log newlogoldi 1（ 13）此時，選擇接受新值new 。2.3.4 步驟 c：對參數(shù)抽樣根據(jù)式 (9)，我們得到一個不是很熟悉的的后驗分布：ss()NNs, y)( +

9、p( |, , )p( )( ( 1 - ) )1N()i=1i(14)采用 Metropolis-Hasting算法中的對數(shù)正態(tài)隨機游走Metropolis-Hasting算法對問題進行處理。假設(shè) lognew N (logold ,c) ，抽取服從均勻分布的隨機數(shù)Vm U 0,1 ，當(dāng)LogVmN (log(new )log (old )log( old ) log (new )N( newold)log(1 is )log p(new ) log p(old )log newlogoldi1(15)此時，選擇接受新值new 。3. 數(shù)據(jù)分析3.1 第一步：為,選擇合適的先驗分布我們通過數(shù)

10、據(jù) y i/T i 的均值 m和方差 s 確定,先驗分布的兩個參數(shù)e ， e ：令 m =E()e 1可以得到： e = m (m1-m = E( )+E( ) =e+e， s = Q = （ 1 - ）e+e+1s - 1) = 1-m1.87 ， e= e m= 17.28 。3.2 第二步： MCMC抽樣在使用 Metropolis-Hasting算法對,進行抽樣時，使用不同的變量c , c 可以求得不同的接受概率，從而確定使接受概率最佳（介于0.40.7 之間）的 c ,c。在這里 c,c 的取值從 0.1 到 1，可以得到,的接受率分別如下：的接受概率如下：4c.beta=0.10.

11、20.30.40.50.60.70.80.91.0c.alpha=0.10.5936 0.5868 0.5884 0.5866 0.5932 0.5932 0.5880 0.5898 0.5850 0.59120.20.4762 0.4778 0.4822 0.4834 0.4830 0.4834 0.4824 0.4818 0.4778 0.48680.30.4202 0.4132 0.4148 0.4110 0.4178 0.4246 0.4050 0.4016 0.4058 0.41620.40.3826 0.3672 0.3876 0.3746 0.3772 0.3688 0.3780

12、 0.3512 0.3734 0.38600.50.3334 0.3460 0.3420 0.3658 0.3374 0.3292 0.3436 0.3468 0.3436 0.35680.60.3100 0.3248 0.3176 0.3132 0.3218 0.3066 0.3096 0.3198 0.3044 0.32480.70.2934 0.2966 0.2998 0.3030 0.2864 0.2922 0.2904 0.2962 0.2874 0.30000.80.2900 0.2734 0.2790 0.2828 0.2800 0.2834 0.2830 0.2776 0.27

13、26 0.28480.90.2764 0.2582 0.2580 0.2792 0.2674 0.2630 0.2622 0.2694 0.2722 0.26261.00.2618 0.2532 0.2538 0.2626 0.2538 0.2436 0.2456 0.2508 0.2534 0.2616的接受概率如下：c.beta=0.10.20.30.40.50.60.70.80.91.0c.alpha=0.10.5354 0.4186 0.3448 0.3200 0.2706 0.2684 0.2468 0.2374 0.2288 0.20800.20.5326 0.4246 0.354

14、6 0.3164 0.2934 0.2668 0.2554 0.2252 0.2208 0.21000.30.5182 0.4294 0.3576 0.3176 0.2926 0.2706 0.2490 0.2374 0.2084 0.21040.40.5274 0.4202 0.3462 0.3308 0.2808 0.2656 0.2538 0.2384 0.2184 0.21040.50.5278 0.4192 0.3636 0.3286 0.2932 0.2776 0.2588 0.2306 0.2212 0.20560.60.5214 0.4080 0.3512 0.3116 0.2

15、912 0.2680 0.2428 0.2342 0.2292 0.20820.70.5156 0.4184 0.3522 0.3168 0.2778 0.2678 0.2460 0.2376 0.2120 0.20880.80.5338 0.4256 0.3634 0.3156 0.2792 0.2538 0.2604 0.2312 0.2216 0.20520.90.5130 0.4134 0.3576 0.3232 0.2968 0.2588 0.2488 0.2418 0.2106 0.21081.00.5234 0.4174 0.3524 0.3172 0.2710 0.2614 0

16、.2402 0.2352 0.2312 0.2088選定 c0.3,c0.2 ， 這樣 可以使得 和 的 接受概率 在 40-43%。 并且選取初始值050, 0100 ，使用 MCMC 抽樣得到(m) ,(m) 的序列如圖 1。圖 1：(m) ,(m) 的序列 ( 初始值 ? = 50,0= 100)05從圖 1 中可以看出，序列很快就趨于收斂了，不妨認為從第501 個抽樣開始趨于收斂。去除數(shù)據(jù)趨于收斂前的樣本，然后計算,的估計值，得到的結(jié)果如下：的均值 = 2.04的方差 = 0.36的中位數(shù) = 1.97的百分之九十五的置信區(qū)間= 1.09 ， 3.39的均值= 17.93的方差 var

17、iance = 15.34的中位數(shù) median = 17.49的百分之九十五的置信區(qū)間= 11.29，26.51畫出,的后驗分布直方圖，如圖2。圖 2： 和 的分布3.3 第三步：分析不可觀測的異質(zhì)性根據(jù)(3) 和 (4) 式，通過 MCMC 抽樣得到的(m) ,(m ) 的值計算均值以及方差 Q 的后驗分布，得到的估計值如下：的均值= 0.1031的方差= 0.0006的中位數(shù) Median = 0.1009的百分之九十五的置信區(qū)間= 0.0614， 0.15556Q 的均值 = 0.0046Q 的方差 = 0.0000026Q 的中位數(shù) = 0.0044Q 的百分之九十五的置信區(qū)間= 0

18、.0024， 0.0080畫出和 Q 的后驗分布直方圖，如圖3。圖 3：u 和 Q 的分布3.4第四步：根據(jù)各個醫(yī)院的事故率si 對醫(yī)院進行排序根據(jù) MCMC 抽樣中各個醫(yī)院的概率值s(m)s( )，計算 的中位數(shù)，并根據(jù)此值對醫(yī)院進sii行排序，結(jié)果如表二。畫出(m)的盒型圖，如圖4。( )i表二、九家醫(yī)院的醫(yī)療事故數(shù)據(jù)醫(yī)院醫(yī)療事故次數(shù)外科病例數(shù)均值按均值排名后驗分布中按后驗分布(y i )(T i )（ yi /T i）位數(shù)中位數(shù)排名11100.170.091662010010.0579231120.08330.0868545600.08340.085445150.290.110586103200.03120.034717301520.19780.1