aStat 實證醫學統計小工具介紹

(部份內容已先行發表於 mobile01.com)


前言

因為工作/研究上的需求，常常得算些二維列表的統計量，比如說: 卡方檢定，費雪精密檢定等。雖然已經有太多的軟體可以計算這些東西，但是很多資料我都是整理成乾淨的 excel 列表再使用樞紐分析表的功能去檢視變數之間的相關性，網路上工具一堆，可是上網有時也不是頂方便的。後來更因為要計算敏感度，特異性的信賴區間，以及 positive likelihood ratio (PLR) 等等，發現常用的 SPSS 並不提供比例或是 PLR 的信賴區間估計功能，按計算機或是用 excel 表又有點麻煩，因而萌生自己寫個小工具的念頭。

aStat 統計小工具 (A Statistical Calculator for Android) 就是這樣的狀況下的產物。

特點

1. 二維列連表 (2x2 contingency table) 的相關統計量計算，包括: 

• 統計意義檢驗 (Test of significance): 卡方檢定 (Chi-square test), 連續性校正的卡方檢定 (Yate's corrected Chi-square test), 費雪精密檢定 (Fisher's exact test)
• 相關性強度檢驗 (Test of association): 勝算比 (odds ratio), 相對危險性 (risk ratio or relative risk), 危險差 (risk difference), 包含 95% 信賴區間的數值。
• 檢驗的正確性評價 (Evaluation of diagnostic power): 計算敏感度 (sensitivity), 特異性 (specificity), 陽性預測值 (PPV: positive predictive value), 陰性預測值 (negative predictive value), 陽性相似比 (positive likelihood ratio), 陰性相似比 (negative likelihood ratio), 亦包含 95% 信賴區間的數值。

2. Person-time 統計推論。

• 比較不同的組別 (暴露 vs. 非暴露) 的 incidence rate 是不是有差? 
• 計算 incidence rate ratio 以及 incidence rate difference

3. NNT/NNH (number needed to treat/harm)

• 給定實驗組與對照組的 events rate (i.e., bad outcomes) 計算 absolute risk reduction/increase
• 提供 NNT/NNH 信賴區間 (Newcombe's method)

4. 單一比例的信賴區間估計

• 使用三個方法計算比例的信賴區間: exact method, Wilson's score method, Agresti-Coull method

5. 自由設定您所要的信賴區間: 90%, 95%, 99%，能符合大多數人的需求
6. 設定顯示非零小數位數: 2-8 位數。
7. 除了提供一般教科書讓所使用的信賴區間估計法，另外還提供諸如 Agresti-Coull, Wilson, Newcombe 等人所提出的信賴區間方法，特別適用於小樣本數的狀況 (如: N < 40-50)。
8. 免費 (自由使用)

對象
研究人員，公衛/統計相關人員，實證醫學領域相關人士

當前版本
0.6

系統需求
Android 1.6 以上

近期計畫
1. Chi-square test for trend。之後應該暫時不會再更新。
中文介面的部份考量到詞彙的統一性，以及還在增加功能的過程，暫時還沒有做出來，還請見諒。

下載
您可以在 android market 搜尋關鍵字 aStat 就可以了，或是掃描以下的 QR code。


開發網頁
http://aStat.twbbs.org

展示圖片



結語

aStat並不同於一般 android market 上的統計工具，其目的並非要取代 PC 上的統計軟體，而是期待能提供一個快速，方便而準確的評估數據的方法。雖然是小眾市場，但是秉持著好東西要分享的精神，希望對您能有所幫忙! 若您覺得好用，請在 market 上給個好評吧! 此外，若您對區間估計的 algorithms 有興趣，程式碼亦可提供給您。目前所使用的機率分佈函式庫為 JSci，如果有人找到更小的歡迎再告訴我，謝謝!

致謝
1. Google。提供了這麼一個開放的平台，讓想像與時間是唯一的限制。
2. gasolin, 就是因為您的教材: 深入淺出 android ，讓非資訊相關產業的我有辦法完成這個小品!
3. 我們的研究團隊跟老闆。

單一比例的信賴區間估計 (Confidence Interval Estimation for Proportions)

整理一下幾種單一比例的信賴區間估計方法。


令

X 代表成功次數

N 代表試驗次數

p = X/N 代表比例

z (alpha/2) = 1.960 at alpha level = 0.05 (i.e., 95% confidence interval)

sqrt = 開平方根

I. Wald method (Wald interval)

教科書裏面最常提到的為 Wald method，又稱為 Wald Interval。

下限 (LL) = p + z * sqrt (p(1-p)/N)

上限 (UL) = p - z * sqrt (p(1-p)/N)

Wald Interval 雖便於計算，但還是有使用上的限制。由於是透過常態分佈近似而來的區間，在樣本數不足，或是 p 接近 0 或 1 時，計算得出的並不是真正的 95% 信賴區間，有 overshooting (大於一) 或 degeneracy (小於零) 的問題。甚至有人證明即使樣本數很大，也不能保證其區間機率 (coverage probability) 近似 95%。但是好算好教好懂是它的最大優點，大部分教科書也是以 Wald Interval 為例子。

II. Exact method (Clopper-Pearson Interval)

由 Clopper 與 Pearson 提出的方法，計算方法較為複雜，可由查表或數值方法或跟 F distribution 的關係作計算。現在電腦很方便，excel 運算功能強大，利用 excel 的 F 分佈反函數即可算出特定單一比例的上下限。以 excel 的內建函數 FINV 則 95% 上下限可用下面的式子表示:

LL = IF(X=0,0,X/(X+(1+N-X)*FINV(0.025,2*(1+N-X),2*X)))
UL = IF(X=N,1,(X+1)*FINV(0.025,2*(X+1),2*(N-X))/(N-X+(X+1)*FINV(0.025,2*(X+1),2*(N-X))))

對於一般的醫學研究，可以用這個方法快速求出某個比例的 95% 信賴區間，例如二維列聯表的 sensitivity 或 specificity， 這也是 Graphpad 統計軟體計算所使用的方法。缺點是可能 over conservative。

III. Wilson's score method (Wilson interval)

LL = (2*n*p + z*z - z*sqrt(z*z+4*n*p*(1-p)))/(2*(n+z*z))
UL = (2*n*p + z*z + z*sqrt(z*z+4*n*p*(1-p)))/(2*(n+z*z))

這個 interval 被近來許多的統計學家用模擬的方法證實，coverage probability 較為接近 95%，即使在樣本數少的狀況下。另外，計算上也還不算太繁複，相較於 exact method 所估計出來的區間比較窄一些。

IV. Modified Wald method (Agresti-Coull Interval)

在 Newcombe 比較七種方法和推崇 Wilson's score method 的差不多時間 (1998)，UF 的一個統計學教授 Alan Agresti 發表文章，原本的 Wald method 經由分子加二，分母加四 (作者稱為 pseudo-observations) 可以大大地改進原本的 Wald method 的 coverage probability，使之接近 95% 名目機率。作者還進一步闡述，經過變換，其實加上 pseudo-observations 使得原本的 Wald interval 可以作為 Wilson's score interval 的近似。也就是說，大部分的狀況，都可以用這個改變後的 p' 數值代入原本的 Wald method 得到接近 95% coverage probability 的信賴區間。這個方法有可能產生 overshooting，也就是區間大於一或小於零的狀況，這時候要把他裁掉。

令 p'= (X+2)/(N+4)

下限 (LL) = p' + z * sqrt (p'(1-p')/N)

上限 (UL) = p' - z * sqrt (p'(1-p')/N)

其實加上的 pseudo-observations:
2 是 z*z/2 的近似值
4 是 z*z 的近似值
因此也可以拿來套用在不同的 alpha level。

在 Agresti 與 Caffo 合著的文章 (2000) 中提到，即使沒有資料也可以計算信賴區間，因為分子為二分母為四可以帶入運算，這是 Agresti 提出的方法一個有趣的地方。

CONCLUSIONS

以應用上來說，若要好算又嚴謹，可以直接使用 Agresti-Coull interval。若要名字好聽，那就用所謂的 exact confidence interval，這是多數統計軟體提供的，也是所謂的金標準。雖然它已經被證實不比起 Agresti-Coull interval 跟 Wilson's interval 精確，而且它的信賴區間過度保守 (conservative ~ 比較寬)，不過一般的接受度應該是很高的。

個人是傾向 Newcombe 提倡的 Wilson's interval 或  Agresti-Coull interval，除了用 excel 方便計算，區間在小樣本的時候看起來窄一點之外，更重要的是它們都可以應用在建立兩個比率的差異的信賴區間。最後，忘了在哪一篇文章中曾提到，樣本數若偏小 (N<40)，Wilson's interval 可能是最好的選擇。

REFERENCES:

1. Newcombe RG. Two-sided confidence intervals for the single proportion: comparison of seven methods. Stat Med. 1998;17(8):857-872.
2. Clopper-Pearson Interval. Available at: http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Clopper-Pearson_interval [Accessed May 26, 2010]
3. Agresti A, Coull BA. Approximate Is Better than "Exact" for Interval Estimation of Binomial Proportions. The American Statistician. 1998;52(2):119-126.
4. Agresti A, Caffo B. Simple and Effective Confidence Intervals for Proportions and Differences of Proportions Result from Adding Two Successes and Two Failures. The American Statistician. 2000;54(4):280-288. 
   
    

A 2x2 contigency table calculator

http://sites.google.com/a/idman.twbbs.org/linkoucgmhid/files/2x2.zip

Updated. 2010/05/24

1. Confidence interval estimation for two proportions (unpaired samples) using Newcombe's procedure (1998), which was suitable for both equal n or unequal n.

2. Confidence interval estimation for a single proportion with three methods: Exact method (reverse F distribution), modified Wald method (Agresti, 1998), and Wilson's score method.

Cox 比例風險回歸 (Cox Proportional Hazards Model) 的應用

The prototype of hazard function：h(t,X) = h0(t) exp (b1X1 + b2X2 + b3X3 + ... + bkXk)

h0(t): baseline hazard function
t: time to event
X: covariates
b: coefficients
exp: exponential function

例子：加護病房病人產生黴菌尿後，有些接受治療有些不治療。經過一段時間後，他們的預後有沒有差別？換句話說，他們的死亡跟有沒有用藥到底有沒有關聯？

此例並非隨機雙盲試驗，我們無法只比較 Kaplan-Meier curves 來得到結論。用不用藥物受到病人的本身身體狀況，以及疾病嚴重度的影響。有沒有用藥又進一步受到時間因素的干擾：有些人來不及用藥，有些剛開始沒有用藥，後來因為狀況不好又把藥物加上去。 在這裡我們可以用 Cox proportional hazards model 來考慮這些干擾因素的影響，以評估藥物與病人預後（死亡）的相關性。


Cox porportianl hazards model 的特性 [1-3]：
1. 不需事先知道 ho(t) 這個 baseline hazard function，也不要求服從特定的機率分佈。因此它有 non-parametric model 的特性。

2. 解釋變數或共變數（covariates）的部份有 parametric model 的特性。所以有些人說 Cox 回歸屬於 semi-parametric models。

3. 解釋變數的風險（hazards）是一個定值，不隨著時間變化。亦即，符合比例風險假定（proportinal hazards assumption）。

4. 解釋變數可以是連續變項，類別變項或是有次序的變項。 


應用 Cox proportional hazards model 來做回歸我們必須注意[1]：

1. 研究樣本數要足夠，一般需要自變數個數的十倍以上。若研究得出的結果違背常理，則必須考慮到模型設計上面的缺陷，不一定是科學上的新發現。

2. 要符合比例風險假定 (proportional hazards assumption)。任意兩個個體在同一個變數的風險比 (risk ratio) 是一個定值，不隨著時間變化。不符合比例風險假定則可分層或是作 Cox regression with time-depedent covariates。以前面黴菌尿的例子，治療會受到時間因素的影響，因此必須作時間協變的 Cox regression。
3. covariates 若是連續變數，可考慮轉成類別變數再分析。比如說血球數目，回歸的時候單位是用顆來計算，所得到的結果意義就代表，血球數每增加或減少一個所帶來的風險變化。
4. covariates 若是沒有次序性的，例如居住地等等，則需加入 dull variables，以免造成解釋上的困難。
5. 變數的篩選需考慮：臨床上有意義的重要的變數要放入模型，無關的不擺進去，或是先作單變數分析再篩選等等。

判斷是否符合風險比例假定有幾個方法 [1-3]：
1. 圖形法（graphic methods）。

• 直接比較 survival curves 是最直觀的方法之一。就在 SPSS 裡面直接作 Cox regression 把懷疑受時間影響的變數放進 Strata，covariates 不放，然後畫圖。如果圖形有明顯交錯就不符合風險比例假定。
• Log-minus-log plots
• Partial residual plots
  

2. Goodness-of-fit test。查不到怎麼實做。
3. Computational methods, like time-interaction test. 一個符合比例風險假定的 covariates，和時間交互作用後的變項，將其加入 Cox 回歸模型，其回歸係數應該很接近零（和 0 沒有顯著差異）[3]。若 covariate 對時間的交互作用項在 Cox 回歸模型中有統計意義，covariate 就不符合比例風險假定 [1]。

REFERENCES:
1. 生存分析。醫學統計學，第二十三章。人民衛生出版社 2005 年第一版。
2. Proportional hazards model. Available at: http://en.wikipedia.org/wiki/Proportional_hazards_models [Accessed January 16, 2010]
3. Cox regression. Available at: http://faculty.chass.ncsu.edu/garson/PA765/cox.htm [Accessed January 17, 2010]

費雪精密檢定（Fisher's Exact Test）

好文章：
1. http://www.math.nsysu.edu.tw/outstanding/use/report/h.pdf
2. http://en.wikipedia.org/wiki/Fisher%27s_exact_test
計算工具：

1. http://www.langsrud.com/fisher.htm
2. http://cgmh.idman.twbbs.org/Home/tools

布瓦松分佈與區間估計 (Poisson Distribution and Estimation of Confidence Interval)

Poisson Distribution 可視為二項式分佈 (Binomial Distribution) 的極限情形。亦即，當試驗次數 (n) 趨近無窮大，p 很小，但一段時間裏面 np 為一不大不小的定值時。從品質管理的角度，如：每日不良產品的數目，不良事件的發生數目等等，每週醫療照顧相關感染的發生次數等等，皆適用於布瓦松分佈。[1]

依從 Poisson 分布的群體，最佳的點估計即為平均值 u。當我們只做了一個單次的觀察，可利用卡方分佈 (Chi square Distribution) 與 Poisson 分布的關係，很快的求出對應的 95% 信賴區間估計 (exact confidence interval)。[2-4]

令 x 為單次的觀察結果 (發生次數)，InvChiSqare (v,y) 函數代表自由度 v 的 Chi Sqare 分布下面積為 y 的機率分布對應的 Chi Sqare 值。則依從 Poisson 分布的 95 % 信賴區間上下限 (exact confidence interval; UL; LL) 與卡方分布有如下的關係：

LL = 1/2 * InvChiSqare (2x, 0.025)
UL = 1/2 * InvChiSqare (2x + 2, 0.975)

當觀察次數不止一次時 (observations = N)，所有的觀察次數合 (n) 依從 Poisson 分布。其 95% 信賴區間上下限可如此表示：[5]

LL = 1/2 * InvChiSqare (2n, 0.025) / N
UL = 1/2 * InvChiSqare (2n + 2, 0.975) / N

References:

1. Poisson 分布. Available at: http://episte.math.ntu.edu.tw/articles/sm/sm_16_07_1/index.html [Accessed Sep 30, 2009]
2. Confidence intervals for the mean of a Poisson distribution. Available at: http://www.math.mcmaster.ca/peter/s743/poissonalpha.html [Accessed Sep 30, 2009]
3. Poisson confidence interval. Available at: http://www.statsdirect.com/help/parametric_methods/pest.htm [Accessed Sep 30, 2009]
4. An exact method for calculating a confidence interval of a poisson parameter. Mulder et al. American Journal of Epidemiology 117 (3): 337. Available at: http://aje.oxfordjournals.org/cgi/pdf_extract/117/3/377 [Accessed Sept 30, 2009]
   
5. Estimating the mean of a Poisson population. Chapter 6.3. Available at Google Books [Accessed Sep 30, 2009]
   

在 Excel 內增加色票 (Add new colors to the default excel palette)

Refer to http://www.bloggpro.com/add-new-colors-to-the-default-excel-palette/

選擇顏色的工具: http://colorschemedesigner.com/

APACHE II Score and Predicted Mortality Rate

The following result was based on original APACHE II (Acute Physiology and Chronic Health Evaluation II) score, published by Knaus et al. in 1985.
----
score = predicted mortality reate
0 = 0.0288323804597415
1 = 0.0332141827839264
2 = 0.0382356916987282
3 = 0.0439818435946882
4 = 0.0505461568401609
5 = 0.0580307274067857
6 = 0.0665458910150112
7 = 0.0762094426987653
8 = 0.0871452904871902
9 = 0.0994814113171375
10 = 0.113346978732479
11 = 0.12886854890866
12 = 0.146165230272512
13 = 0.165342828461315
14 = 0.186487056810361
15 = 0.209656033245229
16 = 0.234872441095319
17 = 0.262115898694886
18 = 0.291316235325043
19 = 0.322348474912558
20 = 0.355030345948782
21 = 0.389123033352462
22 = 0.424335648915215
23 = 0.460333532118856
24 = 0.496750045768218
25 = 0.533201073292596
26 = 0.569301043624338
27 = 0.604679084771601
28 = 0.638993887971074
29 = 0.671946054939952
30 = 0.70328705911182
31 = 0.732824402144069
32 = 0.760423002108014
33 = 0.78600323342325
34 = 0.809536301764071
35 = 0.831037763630525
36 = 0.85056000246354
37 = 0.868184382009093
38 = 0.884013650945879
39 = 0.898165005907147
40 = 0.910764060077794
41 = 0.921939828549029
42 = 0.93182073729233
43 = 0.940531590674707
44 = 0.948191389177323
45 = 0.954911868587507
46 = 0.960796628155722
47 = 0.96594072236451
48 = 0.970430604427108
49 = 0.974344325979221
50 = 0.977751914291392
51 = 0.980715864241253
52 = 0.983291696442564
53 = 0.985528544995485
54 = 0.987469748266567
55 = 0.989153424071641
56 = 0.990613016854848
57 = 0.991877809199794
58 = 0.992973393538575
59 = 0.993922102483462
60 = 0.994743398004825
61 = 0.995454220894235
62 = 0.996069302727867
63 = 0.99660144299736
64 = 0.9970617542933
65 = 0.997459878479777
66 = 0.997804176739593
67 = 0.998101896237682
68 = 0.998359315973963
69 = 0.998581874196749
70 = 0.99877427953855
71 = 0.998940607827617

An Opensource Equation Editor

http://mathcast.sourceforge.net/home.html

Binary Logistirc Regression (SPSS)

基本算式: Ln (Odd) = Ln (Y'/1-Y') = a + bX
Y' represents the predicted probability of the event which is coded with 1 (one).

For dichomotous variables, the calculated odds ration (OR) of a logistic model is

• (the odd in favor of the event coded with 1) / (the odd in favor the event coded with 0)

筆記：

• Continuous variables in a logistic regression model - the importance of the "scale" - description 對應的某單位變動的勝算比，使用的單位不同，勝算比也不同。如：早產兒的體重使用 1 gram vs. 1 kilogram，odds ratio 有不同的大小。

參考資料：

1. 生物谷的教學
2. Multivariate Analysis with SPSS - 一個教學資訊豐富的網頁，包括如何使用 spss 做 binary logistic regression analysis。