內舒拿(Nasonex)與深睡時間（更新為2017整年記錄）

當醫師發現我的鼻子常會塞住，早上容易咳嗽，建議使用類固醇鼻噴劑「內舒拿水溶性鼻用噴液劑 （NASONEX® Aqueous Nasal Spray）」。

我有戴「小米手環2」的習慣，小米手環在睡覺的時候，會紀錄睡眠狀況，包括睡眠時間、深睡時間、淺眠時間等，所謂深睡時間，應該就是睡覺時候沒有太大的動作。
我也很好奇鼻噴劑的效果，如果真的可以改善鼻塞的效果，應也可改善睡眠中斷的狀況，進而增加深睡的時間，因此本次的重點於比較使用Nasonex前後深睡時間的改變。

【於2017/7/6所撰寫的結果觀察時間共150天，而本次更新至2017整年的睡眠記錄，並重新分析結果】

---

KAMERA 急診資料挑戰賽

臉友轉貼上這個比賽，開啟了人生第一場數據分析比賽，KAMERA 急診資料挑戰賽。

原本只是抱著看熱鬧的心情，把資料下載下來後，稍微研究一下資料的格式，開始思考要如何分析這份資料，從不同月份、星期、時段、醫院別等思考，原本還打算整理一些外部的資料，例如溫度、濕度、雨量，但取得方式比預期的困難許多，也可能要花額外的費用，也就做罷了。

回到主辦單為提供資料的原本架構，從中思考有那些變項可以分析，可以套用哪些模型，首先當然先把資料分布做一個基本分析，最後我的第一個模型就是『以前一年同月份資料為基礎，計算各變項組合下的平均值』，這嚴格來說不算模型，算算平均值而已。

第一次上傳後得到0分！！！這個狀況持續了幾天（一天只能上傳兩次），最後發現上傳的檔案格式必須完全符合主辦單位所提供的資料架構才能計算分析，在成功得到分數後，竟然就得到97分（滿分100），過來幾天稍微調整了一下成長率的問題，分數就非常漂亮了。因為原本沒有打算很認真，就把這個比賽放了幾天沒用管他。

等到某一天回來發現有一堆高手分數都高於98分，也激起了我的鬥志，因為不想要再嘗試統計模型了，就將結果依據估計值的大小慢慢調整每一個變項組合的加權（註1），最後得到分數為98.842，最後結算為第5名（佳作）。

發表會那天，前三名得獎者在台上說著他們的分析方法和策略，還真的都滿有趣的。當天行政院長張善政還專程參加（張揆出席高屏地區急診資料挑戰賽成果分享會）。

發表會後的晚宴也相當精緻，得獎者及DSP知識長互相間的討論，討論到資料分析技巧的問題和產業的現況，整體來說我覺得是一個滿有意義的活動，不過比賽真的會要花心力，需要休息一下，過一段時間會考慮挑戰國外的比賽！！

PS：
註1：我調整加權的方法不是有效率的方法，得獎者在其方法中均有不同的加權方式，有人使用K-MEAN或類神經網路或許都是更有效率的方法。

嬰兒猝死症候群（SIDS）與寶寶的年齡

很多新手爸媽應該都會看過類似的文章，嬰兒猝死症候群的預防建議、7重點‧降低嬰兒猝死症發生等文章，文章中都強調如何預防嬰兒猝死症的風險，相關文章都寫得很專業，可以自行參考。

SIDS有許多未知的原因所造成，但對我來說寶寶在何時容易有SIDS的發生，才是我想瞭解的。

彙整了政府的2008-2013年的統計資料，這六年間有231名不幸的寶寶發生SIDS，死亡率為每萬人口數1.95人，也就是平均來說每1萬個寶寶有1.95人會死於SIDS，這個數值並不算高（註1）。

依據這231位寶寶的資料，未滿月寶寶佔SIDS的18.2％，在滿1個月但未滿2個月寶寶佔SIDS的23.4％，占所比率為最高。

試管嬰兒的成功率（以政府資料分享）

所謂試管嬰兒的技術是人工生殖技術的其中之一（如果有誤請告知），在國民健康署的「100年台灣地區人工生殖施行結果分析報告」(註1)中提到政府所針對人工生殖技術包括了四種方法，最常用的是所謂「試管嬰兒 (In Vitro Fertilization, IVF/ET )」，在資料中有超過99%接受治療的都是使用這個技術。
而人工生殖技術（因為有99%為試管嬰兒，所以後面就只接稱試管嬰兒）成功率在報告中其實分為兩種，「懷孕率」及「活產率」(註2)。
先不論治療者何種狀況平均來說懷孕率為36.9%，活產率為27.4%。（因為有99%均採用試管嬰兒技術，幾乎可以直接反應為試管嬰兒的成功率）
但其實這個成功率與接受治療者的年齡有很大的關係，下面的圖就表示治療者年齡和兩種成功率的關係，在過30歲後成功率就一直下滑。

(原)台中市房價的對話－買房屋要多少錢？

希望…

當我準備在(原)台中市購屋後，也看一些待售的房屋後，一些打算購屋朋友（但還沒行動），就開始出現許多有趣的問題，例如

台中房價到底如何？

什麼那個地方的房屋不值這個價錢啦！

這個價錢都可以買獨棟房屋了！

房價有那麼貴嗎？應該是亂開的啦！

聽到這些問題後，我心中有許多的疑問，房價真有如朋友口中得那麼高呢，還是朋友們的房價記憶能停留在先前的年代？

在早期的年代，所有房價的訊息只能透過仲介，我自己看房子也愛問，這樣狀況的房屋應該會賣多少錢？多少錢買到才會合理？不管仲介給怎樣的數據我都很懷疑。

在內政部推不動產實價登記前，有內政部曾邀請房仲業者自行將售屋資訊提供出來（http://ehouse.land.moi.gov.tw/A0/A10/A1001/A1001A.aspx），這是在實價登錄前相對完整的房屋成交資訊了。

因此~我試圖對台中市房價的瞭解了。

-----------故事開始-------------

英國孩童過多？？

這是一篇網路上看到的新聞「學童太多英國小學教室不夠」，主要是說英國這幾年出生率提高，有些地方政府，在數年前為了因應出生率過低的情形，將部分的教室廢除了，而使得今年數千名學童沒有教室可用。

這一篇新聞打破我的一個固有的印象，英國(或是說英格蘭)的出生率不是非常低嗎，甚至看過一篇報到提到「英國人睡覺寧可抱著熱水袋，也不願意抱著另一半睡，英國人如何繁殖後代真是一個謎」，但怎麼突然之間……他們的孩童變成過多了！！那台灣呢？

使用ODS GRAPHICS和PROC LIFETEST繪製存活曲線圖 (Survival Curve)

太久沒寫有的沒有的SAS程式了，這次來寫寫關於Kaplan-Meier Estimate、Log-Rank Test及KM Plot的心得吧。
這篇主要的目的為兩個，如何使用PROC LIFETEST進行存活資料的的分析（跳過統計理論的部分）， 另一個為使用SAS 9.2新的ODS GRAPHICS，快速產生勘用的存活機率圖，並將圖檔自動放入RTF(Rich Text Format)檔中。

在SAS中使用已經整理好的統計量

實務上常遇到研究者只提供已算好的數字(統計量)，例如列聯表中的個數和百分比或兩樣本的人數、平均數和標準差等，如果要利用這些統計量做某些檢定，例如卡方檢定或雙樣本T檢定等，可使用PROC FREQ和PROC TTEST處理僅有統計量的問題。

PROC FREQ

有兩群病患分別使用A藥和B藥，再觀察兩組實驗組症狀有無改善，結果如下表，如研究者想問A藥和B藥病患症狀改善的比率有無不同，可採用卡方檢定(Chi-square test)

	A藥	B藥
症狀無改善	9	12
症狀改善	21	18

當然第一步要先把資料輸入至SAS中，程式可以這樣寫：

DATA TRIAL_1;
INPUT DRUG $ EFFECT COUNT;
DATALINES;
A 0 9
A 1 21
B 0 12
B 1 18
;
RUN;
PROC SQL;
SELECT *
FROM TRIAL_1;
QUIT;

接下來是卡方檢定的部分

PROC FREQ DATA=TRIAL_1;
WEIGHT COUNT;
TABLE EFFECT*DRUG/CHISQ;
RUN;

PROC TTEST

同樣兩群病患使用A藥和B藥，使用後測量某個臨床指標（連續型變項），研究者欲使用T檢定，探討兩組病患的臨床指標平均值有無差異。資料如下表

	A藥(N=30)	B藥(N=30)
臨床指標	89(12)	95(16)

mean(std)

因為這裡我們只有拿到兩組病患的人數、平均值和標準差，所以在輸入資料時要稍微複雜一點，必須符合SAS程式的規定，一定要有一個變項名稱為『_TYPE_』或『_STAT_』，在這個變項中，必須要放入‘N’、‘MEAN’、 ‘STD’這三筆資料，資料輸入程式如下

DATA TRIAL_2;
INPUT DRUG $ _TYPE_ $ EFFECT;
DATALINES;
A N 30
A MEAN 89
A STD 12
B N 30
B MEAN 95
B STD 18
;
RUN;

接下來執行PROC TTEST程序就簡單了，和一般寫法沒有太大差異

PROC TTEST DATA=TRIAL_2;
CLASS DRUG;
VAR EFFECT;
RUN;

以上程式可於SAS 9.1.3和SAS 9.2執行。

繪製ROC曲線和比較兩條ROC曲線下面積

SAS實用巨集簡介：繪製ROC曲線和比較兩條ROC曲線下面積

目的為介紹SAS官方網頁中關於ROC曲線的巨集(Macro)程式，程式分別為「繪製ROC曲線」和「以無母數方式比較兩條ROC曲線下面積」。本文中所提及之程式，僅為原作者的部分語法，詳細語法請見原文說明，兩篇原文出處連結如下：

Sample 25017: Nonparametric comparison of areas under correlated ROC curves (http://support.sas.com/kb/25/017.html)
　roc.sas程式：http://support.sas.com/kb/25/addl/fusion_25017_6_roc.sas.txt

Sample 25018: Plot ROC curve with labelled points for a binary-response model (http://support.sas.com/kb/25/018.html)
　rocplot.sas程式：http://support.sas.com/kb/25/addl/fusion25018_4_rocplot.sas.txt


請先將程式下載至本機硬碟中，建議將Sample 25017的程式取名為roc.sas，而Sample 25018的程式取名為rocplot.sas。必須要注意的是，如果是想比較兩條ROC曲線下面積，也就是roc.sas這個巨集程式，必須安裝SAS/IML。

繪製ROC曲線

在SAS說明檔中其實已經有介紹如何繪製出簡易的ROC曲線，但在我們諮詢的個案中比較常被問到的問題為，圖上的某個點其「原始數值為何」，利用rocplot.sas 便可輕鬆完成這個工作。
假設有一組資料，研究者想以年齡和某一個BW指標來判斷罹患某疾病的機會。
問題：繪製以年齡作為判斷此是否罹患此疾病的ROC曲線。

第一步：先將資料讀入SAS中
data Age_data;
input disease age bw@@;
datalines;
0 50 65 0 39 61 0 21 70 0 61 67
0 30 55 0 35 63 0 25 72 0 41 66
0 43 52 0 36 54 0 37 76 0 25 61
0 41 53 0 62 55 0 28 70 0 33 68
1 52 45 1 49 61 1 47 42 1 62 31
1 55 67 1 70 61 1 75 55 1 77 52
1 81 58 1 64 53 1 62 41 1 39 57
1 61 51 1 61 55 1 57 49 1 79 47
;
run;

第二步：將rocplot.sas程式預先讀入SAS中
%include "c:\temp\ROC\rocplot.sas";
/*
%include “程式所放的位置”;
*/
透過%include這個巨集指令，可將rocplot.sas呼叫至SAS程式中。

第三步：執行proc logistic程式，並將部分結果輸出，以配合rocplot.sas程式
proc logistic data=age_data;
model disease (event='1')=age / outroc=roc1 roceps=0;
output out=outp p=phat;
run;
這個程式中，我們將disease變項中的「1」視為事件(罹患疾病)，並將繪製ROC曲線所需的資料輸出為「roc1」，和每個原值(年齡)的預測機率資料輸出為「outp」並將預測機率的變項命名為「phat」。基本上，使用時只需要改動依變項(disease)和自變項(age)這兩個部分即可。

第四步：使用%rocplot這個巨集程式。
%rocplot (outroc = roc1, out = outp, p = phat, id = age);
這個巨集中，outroc= 資料檔，配合proc logistic的outroc輸出。
out= 資料檔，也同樣配合proc logistic 的out輸出。
p= 預測機率值變項名稱，在輸出預測機率時的命名。
Id= 變項名稱，要依據那個變項繪製ROC曲線上所標記的原始數值，通常應與自變項相同。

第五步：圖形輸出結果



比較兩條ROC曲線下面積

問題：比較以年齡或BW指標對於判斷是否罹患疾病，此兩指標ROC曲線下面積。

第一步：先將roc.sas程式呼叫入SAS中
%include "c:\temp\ROC\roc.sas";


第二步：執行自變項為age的logistic 模型
proc logistic data=age_data;
model disease(event='1')=age / outroc=roc1_age roceps=0;
output out=outp_age p=phat_age;
run;
為了區別age和bw這兩個結果，這裡我們將age的output輸出的檔案命名為「outp_age」而預測機率命名為「phat_age」。

第三步：執行自變項為BW的logistic 模型
proc logistic data=age_data;
model disease(event='1')=bw / outroc=roc1_bw roceps=0;
output out=outp_bw p=phat_bw;
run;


第四步：執行%roc 巨集
%roc (data = outp_age outp_bw,
var = phat_age phat_bw,
response = disease);
在data=資料檔，這後面主要放前面步驟透過output所輸出的檔案，因此目前所放的為「outp_age」和「outp_bw」兩個檔案。
var=變項名稱，在output檔中，預測機率的變項名稱。
response= 依變項名稱，宣告依變項名稱為何。

第五步：結果
於「ROC Curve Areas and 95% Confidence Intervals」的表格中，age的AUC(area under curve, 曲線下面積)為0.912，而bw的AUC為0.6836。age和bw這兩個自變項AUC比較，可於「Contrast Test Results」中看到p-value為0.0378，因此這兩條ROC曲線下面積為顯著的不同。

The Rule of Defining Dummy Variables

1. If an inercept is used in the regression equation, proper difinition of the k-1 dummy variables automatically indexes an k categories.
2. If k dummy variables are used to describe a nominal variable with k categories in a model containing an intercept, all the coefficients in the model cannot be uniquely estimated because collinearity is present(Why??).