IDNote: 2010

2010年11月28日星期日

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. 我們的研究團隊跟老闆。

2010年8月2日星期一

可轉債自動報價查詢 II

可轉債, CB, 自動化查詢, 股票, 金融

改進了之前的 script ，全部改成 python 來寫。可由標準輸入檔案讀入想查詢的台灣股市上市上櫃公司代號，完畢後輸出成 html table 並發信至指定的電郵信箱。

usage:

python thisscript.py < stocklist.txt

這支 script 基本上做了下面的事情:

1. 依代號查網頁 (tw.stock.yahoo.com)

2. 解析網頁資料

3. 輸出成 html table

4. 寄信 with attachment

#! /usr/bin/env python
# -*- coding: big5 -*-
import urllib2, re, sys, csv
import smtplib
from datetime import datetime
from email.mime.multipart import MIMEMultipart
from email.mime.text import MIMEText
from email import encoders
from email.utils import COMMASPACE

html = '<html><table border="1" style="border-collapse:collapse;" borderColor=black>\n'
html = html + '<tr><th>代碼</th><th>HTTP</th><th>名稱</th><th>成交價</th><th>漲跌</th><th>買進</th><th>賣出</th>'
html = html + '<th>張數</th><th>昨收</th><th>開盤</th><th>最高</th><th>最低</th></tr>\n'
data = []
for line in sys.stdin:
F1, F2, F3, F4, F5 = '', '', '', '', ''
F1 = re.match('\d+', line).group(0)
url = 'http://tw.stock.yahoo.com/q/q?s=' + F1
myrequest = urllib2.Request(url)
myrequest.add_header('User-Agent', 'Mozilla 5.0')
try:
content = urllib2.urlopen(myrequest)
F2 = 'OK'
except:
F2 = 'Failed'
c = content.read()
m = re.search(r'href="/q/bc\?s=\d+">(\d+.*?)</a>', c)
F3 = m.group(1) if m else '－'
m = re.search(r'nowrap><b>(\d+.?\d+|\d+|－)</b></td>', c)
F4 = m.group(1) if m else '－'
m = re.search(r'nowrap><font color=#......>(▽\d+.\d+|△\d+.\d+|\d+.\d+|－)', c)
F5 = m.group(1) if m else '－'
F6 = re.findall(r'wrap>(\d+.?\d+|\d+|－)</td>', c)
data.append([F1, F2, F3, F4, F5]+F6[1:])
for line in sorted(data, key=lambda x: float(x[3]) if (x[3] <> '－') else x[3]):
result = ''
for x in line:
result = result + x + ','
result = re.sub(',' , '</td><td>', result)
result = '<tr><th>'+result+'</tr>\n'
result = re.sub(r'<th>(\d+)</td>', r'<th>\1</th>', result)
html = html + result

html = html + '</table></html>\n'

tt = datetime.now().timetuple()

me = "sender@email.address"
you = ['recipient1@email.address', 'recipient2@email.address']

# Create message container - the correct MIME type is multipart/alternative.
msg = MIMEMultipart('alternative')
msg['Subject'] = "CB Report " + str(tt[0]) + "-" + str(tt[1]) + "-"  + str(tt[2]) + " "  + str(tt[3]) + ":"  + str(tt[4]) + ":"  + str(tt[5])
msg['From'] = me
msg['To'] = COMMASPACE.join(you)

# Create the body of the message (a plain-text and an HTML version).
# text = """ <html><b>test</b></html> """
# html = """ this is html data """
# Record the MIME types of both parts - text/plain and text/html.
# part1 = MIMEText(text, 'plain')
part2 = MIMEText(html, 'html', _charset='Big5') # _charset is important to avoid bad characters

# Attach parts into message container.
# According to RFC 2046, the last part of a multipart message, in this case
# the HTML message, is best and preferred.
#msg.attach(part1)
msg.attach(part2)

# Send the message via local SMTP server.
s = smtplib.SMTP('localhost')
# sendmail function takes 3 arguments: sender's address, recipient's address
# and message to send - here it is sent as one string.
s.sendmail(me, you, msg.as_string())
s.quit()

2010年7月21日星期三

可轉債自動報價查詢

關鍵字: perl, python, 可轉債 (convertible bond)

上市公司及上櫃公司的可轉債約有200多卷，可由公開資訊觀測站查詢相關的訊息。以可轉債的特性而言，比較有用的是跌到票面價格以下的債卷，但必須要常常追蹤價格以及市場消息，對於非專業人士實在力有未逮。如果要一筆一筆輸入 Yahoo 的投資組合來 monitor，也是耗時又費力。因此在一台 linux 上面寫了個 perl 的查詢器，再放到 crontab 來定時自動查詢，方便性增加不少。在這邊做個紀錄。

下面是一開始用的 perl script，可以用標準輸入讀入可轉債代號文字檔案 (一個代號一列)，查詢後輸出到標準輸出 (csv tab 分隔文字格式)。以下的方法也可以用來查詢一般的股票。

#!/usr/bin/perl -w
use LWP::UserAgent;
use HTTP::Request;
use HTTP::Response;
$ua = new LWP::UserAgent; # 產生 UserAgent 物件
my $year = (localtime)[5] + 1900; # year
my $month = (localtime)[4] + 1; # month
my $day = (localtime)[3]; # day
my $date = $year."-".$month."-".$day;
print "代碼"."\t"."HTTP"."\t"."名稱"."\t"."成交價"."\n";
foreach (<>) {
chomp($_);
$request = new HTTP::Request('GET', 'http://tw.stock.yahoo.com/q/q?s='.$_); # 產生 Request 物件
$response = $ua->request($request); # 開始抓取網頁，並將結果傳會 $response
if ($response->is_success) { # 若抓取網頁成功，則印出 HTML 原始碼
   my $html = $response->content;
   if ($html =~ m{href="/q/bc\?s=\d+">(\d+.*?|－)}) {
   print $_."\tOK\t".$1."\t";}
   else {print $_."\tOK\t"."抓不到\t";}
   if ($html =~ m{nowrap>(\d+.?\d+)}) {
   print $1."\n";}
   else {print "抓不到\n";}
   } else { # 若抓取網頁不成功，則印出錯誤訊息
   #print $response->error_as_HTML;
   print $_."\t"."網頁讀取有問題"."\n";
}
}

配合上 crontab 執行以下的 perl script 可以固定每天定時查詢後發信至指定的信箱，包括日期跟時間

#!/usr/bin/perl -w
my $year = (localtime)[5] + 1900; # year
my $month = (localtime)[4] + 1; # month
my $day = (localtime)[3]; # day
my $hour = (localtime)[2]; # hour
my $minute = (localtime)[1]; # minute
my $second = (localtime)[0]; # second
my $cbdate = $year."-".$month."-".$day." ".$hour.":".$minute.":".$second;
my $output = "/tmp/"."cblist-".$year."-".$month."-".$day."-".$hour."-".$minute."-".$second.".txt";

`/path-to-script/queryscript.pl < /convertible_bond_list.txt > $output`;
sleep(25);
open IN, ("< $output");

$to='your@email.address';
$from= 'your@email.address';
$subject='CB Report '.$cbdate;

open(MAIL, "|/usr/sbin/sendmail -t");

## Mail Header
print MAIL "To: $to\n";
print MAIL "From: $from\n";
print MAIL "Subject: $subject\n\n";
## Mail Body
foreach () {
print MAIL $_;
}

close(MAIL);

查完後發現還是不夠完美，我希望可以針對成交價由小至大做排序，這樣只要看前面的幾個可轉債就好了，所以把查詢的 script 用 python 來改寫，輸出 csv 逗號分隔檔至標準輸出裝置，另外增加了 yahoo 上可查詢到的其他欄位。

#! /usr/bin/env python
# -*- coding: big5 -*-
import urllib2, re, sys, csv
print '代碼,HTTP,名稱,成交價,漲跌,買進,賣出,張數,昨收,開盤,最高,最低'
data = []
for line in sys.stdin:
   F1, F2, F3, F4, F5 = '', '', '', '', ''
   F1 = re.match('\d+', line).group(0)
   url = 'http://tw.stock.yahoo.com/q/q?s=' + F1
   myrequest = urllib2.Request(url)
   myrequest.add_header('User-Agent', 'Mozilla 5.0')
   try:
   content = urllib2.urlopen(myrequest)
   F2 = 'OK'
   except:
   F2 = 'Failed'
   c = content.read()
   m = re.search(r'href="/q/bc\?s=\d+">(\d+.*?)', c)
   F3 = m.group(1) if m else '－'
   m = re.search(r'nowrap>(\d+.?\d+|\d+|－)', c)
   F4 = m.group(1) if m else '－'
   m = re.search(r'nowrap>(▽\d+.\d+|△\d+.\d+|\d+.\d+|－)', c)
   F5 = m.group(1) if m else '－'
   F6 = re.findall(r'wrap>(\d+.?\d+|\d+|－)', c)
   data.append([F1, F2, F3, F4, F5]+F6[1:])
for line in sorted(data, key=lambda x: float(x[3]) if (x[3] <> '－') else x[3]):
   result = ''
   for x in line:
   result = result + x + ','
   print result[:-1]

這個 script 有幾個問題
1. 在比大小的時候，python 內建的排序遇到浮點數會被當成字串，因為一開始我們查詢回來的資料就是字串，因此比大小的時候要 type casting。
2. 有用到 key，所以必須要用 python2.4 以上。

2010年5月27日星期四

單一比例的信賴區間估計 (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:

Newcombe RG. Two-sided confidence intervals for the single proportion: comparison of seven methods. Stat Med. 1998;17(8):857-872.
Clopper-Pearson Interval. Available at: http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Clopper-Pearson_interval [Accessed May 26, 2010]
Agresti A, Coull BA. Approximate Is Better than "Exact" for Interval Estimation of Binomial Proportions. The American Statistician. 1998;52(2):119-126.
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.

2010年5月24日星期一

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.

2010年1月17日星期日

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]

2010年1月1日星期五

費雪精密檢定（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

訂閱：文章 (Atom)

2010年11月28日 星期日

2010年8月2日 星期一

2010年7月21日 星期三

2010年5月27日 星期四

2010年5月24日 星期一

2010年1月17日 星期日

2010年1月1日 星期五