這一篇是置頂讓大家問問題的！
1.有要問語法問題的請在此留言。
2.即日起不回覆私底下email給我的問題，請善加利用此問題討論區。
3.不保證能夠回答所有問題，若不見我回覆則盼熱心網友能協助解答。
4.請勿貼出一大串程式碼說跑不出來要我幫忙debug。
5.SAS軟體本身的問題請自行聯繫SAS客服部門。
(PS)原Question and Discussion因回文過多被拆成兩頁導致使用不便，所以刪除了。
2108年9月30日星期日
問題討論區
加入書籤：
2009年11月18日星期三
SAS® Abbreviations Are Your Friends; Use a Template Method to Code!
加入書籤：
原文載點：SAS® Abbreviations Are Your Friends; Use a Template Method to Code!
SAS 造成許多初學者學習障礙的主要原因是因為語法很多，對於沒有程式撰寫基礎的人來說，門檻可說是相當的高。坊間有許多不同的 SAS 參考書，但內容程度不一，不過有個共同的特點就是，頁數都很多，內容也很雜。通常參考書一開始會從資料格式與輸入開始教，再來是敘述統計量的輸出方法，還有一些簡單的 Proc 程序，然後是作圖，最後通常會用一些函式和 macro 的介紹。印象中自己也買過兩本中文的 SAS 參考書，但必須承認的一件事是：我從來沒有從頭到尾看完過一遍，通常都是要用到某些語法時才去翻。即便自己對 SAS 已屬稍微熟悉，但仍舊沒有辦法記住所有語法，每隔一段時間，某些語法就會從腦海裡面消失，等到要用的時候再去翻書。後來網路發達了，可以用孤狗去搜尋語法，看似比翻書快了一些，但網路搜尋的結果有時候一開始並不是自己真正想要去尋找的目標，也許需要翻個幾頁，開過幾個連結看一下才能找到真正想要的內容。這種情況所消耗的時間可能跟翻書是差不多的。此外，即便是已經熟悉語法了，但可能因為分析內容的需要，必須大量重複使用一些相當簡單比方說 Proc sort 或 Proc print 等程序。Macro 雖然提供了簡化程式行數的功能，但像 Proc sort 或 Proc print 程序本來就不需要花什麼行數，先宣告 macro 並沒有特別節省掉多少程式碼。如果有一種功能，就是當寫程式寫到一半，知道某個語法但忘記裡面要放什麼參數，SAS 能適時給予提醒，或者是利用更簡短的自訂名稱來取代頻繁使用的簡單程序，便可以讓程式撰寫的過程更省時省力。這種功能在 SAS 裡面，叫做「SAS abbreviation」。
2009年6月18日星期四
Examining Mediator and Moderator effect using Rural Women HIV Study
加入書籤：
原文載點：http://support.sas.com/resources/papers/proceedings09/191-2009.pdf
這一篇技術文件是簡單地利用一個真正的女性HIV資料來教如何使用 SAS 檢定 mediator（或稱 mediation） 和 moderator。關於 mediator 和 moderator 的定義請參照：
Mediator：http://davidakenny.net/cm/mediate.htm
Moderator：http://davidakenny.net/cm/moderation.htm
首先，這個資料背景是來自一個cross-sectional的長期研究裡面所抽出來的第一次面訪資料，總計有 280 位遭到 HIV 感染的女性。
Baron & Kenny (1986) 首度發表檢測 mediator 的方法，整個流程要跑三個迴歸模式（Eq1, Eq2, Eq3），並且要符合四個準則（C1, C2, C3, C4）。三個迴歸模式分別為：
(Eq1) IV -> DV
(Eq2) IV -> M (Mediator)
(Eq3) IV + M -> DV
前兩個準則 C1 和 C2 是，如果 Eq1 和 Eq2 都出現顯著的結果，基本上就表示 Mediator 可能是存在的。此外，還需要另兩個存在於 Eq3 的準則需要符合：
(C3) M 在 Eq3 一定要顯著。
(C4) IV 的估計參數在 Eq3 要降到 0。
若這四個準則都達成了，則可以說這個 M 變數是 full mediator。如果 C4 沒有達成，則稱 M 為 partial mediator。
最後，再用 Sobel test 來檢定一個 mediator 是否顯著地將 IV -> DV 的總效應完全取代。
因此，在這個 HIV Women 的資料庫中，定義 Available Social Support (tssqav) 為 IV，Reason for Missing Medication (treas) 為 DV，而 mediator 變數則為 Spiritual Activity (tcopesa)。
我們可以連續用三個 PROC REG 程序把 Eq1~Eq3 給建立起來：
ods rtf;
ods listing close;
proc reg data=two;
model treas = tssqav / stb pcorr2 scorr2;
title ' Regression model / step1 y=x' ;
run;
proc reg data=two;
model tcopesa = tssqav / stb pcorr2 scorr2;
title ' Regression model / step2 m=x' ;
run;
proc reg data=two;
model treas = tssqav tcopesa / stb pcorr2 scorr2;
title ' Regression model / step3 y=m x' ;
run;
ods rtf close;
ods listing;
quit;
run;不過本篇技術文件並沒有繼續去做 Sobel test，但前北卡大心理系教授 Dr. Preacher 和 OSU 教授 Dr. Hayes 曾經在 2004 年發表了一篇關於 Sobel test 的論文，裡面有附完整的 Sobel test SAS macro。
原文：http://www.comm.ohio-state.edu/ahayes/BRMIC2004.pdf
程式：here
語法：%sobel(data=file, y=dv, x=iv, m=med, boot=z);
其中，data表示想要呼叫進來使用的資料，y 是放 DV 變數名稱，x 是放 IV 變數名稱， m 是放 mediator 變數名稱，而 boot 則是指定要做 bootstrap resampling 的次數，從 1000 到 1000000 之間任一數字皆可。如果不想用的話就直接寫 0，這樣一來 %sobel 會自動關閉 bootstrap resampling 的功能。
網路上也有不少的 Sobel test calculator，可自行用 google 搜尋。
那麼，要進行 moderator 的檢定，則需要配適這個迴歸模型：IV+M+IV*M -> DV。如果 IV*M 的估計參數是顯著的話，則表示 M 是 IV 和 DV 的 moderator。程式如下：
ods rtf;
ods listing close
proc reg data=two;
model treas = tssqav tcopesa sscopesa/ stb pcorr2 scorr2;
title ' Regression model / testing moderator effect' ;
run;
ods rtf close;
ods listing;
quit;
run;最後我們得到 IV*M 的估計參數 β=0.00175，其 p-value=0.5172 並不顯著，因此可以測得 Spiritual Activity 並不是 Moderator。
2009年6月11日星期四
Using PROC SGPLOT for Quick High-Quality Graphs
加入書籤：
原文載點：http://support.sas.com/resources/papers/proceedings09/158-2009.pdf
雖然不曉得有多少人已經拿到 SAS V9.2，不過由於我已經拿到了，所以之後會開始陸續介紹一些新版的功能。
首先先來展示一個 V9.2 最新的繪圖程序—PROC SGPLOT。舊版的 SAS 雖然有提供繪圖程序，但是他們都分散在不同的程序裡面，反而造成使用者的不便。此外，他們的老毛病還是存在，那就是畫出來的圖品質不佳，後來雖然有 ODS 的協助，稍微改善了這方面的缺失，不過 V9.2 版把這些舊的繪圖程序都打包在 PROC SGPLOT 裡面。SGPLOT 顧名思義就是 sophisticated graphical plot 的縮寫，讓我們先來看看這個新繪圖程序的功能。
。HISTOGRAMS
舊版的長條圖要用 PROC GCHART 或 PROC NORMAL 裡面的 HISTOGRAM statement 才能畫出，而且還需要加上許多語法。現在在 PROC SGPLOT 裡面只要使用 HISTOGRAM statement，後面加上變數名稱，就可以完成一個精美的長條圖。
範例：
成果：PROC SGPLOT DATA = Freestyle;
HISTOGRAM Time;
TITLE "Olympic Men's Swimming Freestyle 100";
RUN; 
若想要加上機率密度曲線，則只要多寫一行 DENSITY option，並宣告變數名稱即可。
範例：
成果：PROC SGPLOT DATA = Freestyle;
HISTOGRAM Time;
DENSITY Time;
TITLE "Olympic Men's Swimming Freestyle 100";
RUN; 
。BAR CHARTS
要畫柱狀圖的的話，原本 PROC GCHART 裡面的 VBAR 和 HBAR 被完全移植過來。
範例：
成果：PROC SGPLOT DATA = Countries;
VBAR Region;
TITLE 'Olympic Countries by Region';
RUN; 
若想要顯示每根 bar 裡面不同群體所佔的比例，則只要在後面加上 GROUP option 即可。
範例：
成果：PROC SGPLOT DATA = Countries;
VBAR Region / GROUP = PopGroup;
TITLE 'Olympic Countries by Region and Population Group';
RUN; 
如果要計算次數的目標是另一個變數的話，可以用 RESPONSE option來另外累計。
範例：
成果：PROC SGPLOT DATA = Countries;
VBAR Region / RESPONSE = NumParticipants;
TITLE 'Olympic Participants by Region';
RUN;
。SERIES PLOTS
畫 X-Y 座標圖的功能被整合進 SERIES statement 中。而且重複呼叫 SERIES statement 的話，可以自動完成重疊圖形的功能，不用像以前一樣一定得加上 overlay option。
範例：
成果：PROC SGPLOT DATA = Weather;
SERIES X = Month Y = BRain;
SERIES X = Month Y = VRain;
SERIES X = Month Y = LRain;
TITLE 'Average Monthly Rainfall in Olympic Cities';
RUN; 
接下來看看怎樣替圖形做一些細部的調整。
。XAXIS AND YAXIS STATEMENTS
這兩個 statement 便是拿來更改 X 軸和 Y 軸設定的語法，功能就和以前的 AXISn statement 一樣。從上面的圖可以發現，X-軸代表月份，但是刻度卻是 2.5, 5.0, 7.5, 10.0, 12.5，顯然是不合理的，若加上 TYPE=DISCRETE 則會以實際資料裡面的數據來刻畫度數。GRID 則是在每個刻度上面劃上一條淡灰色的準線。如果 XAXIS 和 YAXIS 都加上 GRID option 的話就可以畫出格狀的底圖。LABEL 自然就是將軸重新命名，而 VALUES 則可以自己定義刻度的起始點和間距大小。
範例：
成果：PROC SGPLOT DATA = Weather;
SERIES X = Month Y = BRain;
SERIES X = Month Y = VRain;
SERIES X = Month Y = LRain;
XAXIS TYPE = DISCRETE GRID;
YAXIS LABEL = 'Rain in Inches' GRID VALUES = (0 TO 10 BY 1);
TITLE 'Average Monthly Rainfall in Olympic Cities';
RUN; 
。PLOT STATEMENT OPTIONS
如果想要針對座標軸內的線條圖形或是圖例說明做調整，則需要在 SEREIS statement 後面加上一些 option。LEGENDLABEL 可以更改圖例說明內的標籤，MARKERS 的功能就和以前的 SYMBOL 一樣，可以在資料點上標上符號。LINEATTRS 則是可以規範線條型態和粗細。
範例：
成果：PROC SGPLOT DATA = Weather;
SERIES X = Month Y = BRain / LEGENDLABEL = 'Beijing' MARKERS LINEATTRS = (THICKNESS = 2);
SERIES X = Month Y = VRain / LEGENDLABEL = 'Vancouver' MARKERS LINEATTRS = (THICKNESS = 2);
SERIES X = Month Y = LRain / LEGENDLABEL = 'London' MARKERS LINEATTRS = (THICKNESS = 2);
XAXIS TYPE = DISCRETE;
TITLE 'Average Monthly Rainfall in Olympic Cities';
RUN; 
。REFLINE STATEMENT
如果想要在座標圖上加上一些參考線，則可以 REFLINE statement來完成。同樣地，參考線也可以做一些細部設定，比方說透明度可以用 TRANSPARENCY option 來設定，每條參考線也可以用 LABEL option 寫上標籤。
範例：
成果：PROC SGPLOT DATA = Weather;
SERIES X = Month Y = BRain;
SERIES X = Month Y = VRain;
SERIES X = Month Y = LRain;
XAXIS TYPE = DISCRETE;
REFLINE 2.03 4.78 1.94 / TRANSPARENCY = 0.5 LABEL = ('Beijing(Mean)' 'Vancouver(Mean)' 'London(Mean)');
TITLE 'Average Monthly Rainfall in Olympic Cities';
RUN; 
。INSET STATEMENT
座標圖內可以用 INSET statement 寫上一些註釋文字，至於位置當然就得用 POSITION option 來設定。若要將註釋文字加框，則簡單地用 BORDER option 即可搞定。
範例：
成果：PROC SGPLOT DATA = Weather;
SERIES X = Month Y = BRain;
SERIES X = Month Y = VRain;
SERIES X = Month Y = LRain;
XAXIS TYPE = DISCRETE;
INSET 'Source Lonely Planet Guide'/ POSITION = TOPRIGHT BORDER;
TITLE 'Average Monthly Rainfall in Olympic Cities';
RUN; 
。THE SGPANEL PROCEDURE
以前如果要針對一個資料集裡面不同的群組或個體畫出各自的圖形並且放在同一張圖裡面，是相當大費周章的事情。現在這種苦差事利用 PROC SPGPANEL 程序即可輕鬆解決。假設我們要畫一個資料集的迴歸線，用 PROC SGPLOT 的 REG statement 就可以畫出。程式如下：
圖形如下：PROC SGPLOT DATA=sg.countries;
REG X=NumParticipants Y=TotalMedals;
TITLE 'Number of Participants by Total Medals Won for Each Country';
RUN; 
如果這個資料集裡面共包含六個區域，該如何分別製作迴歸圖，並且用 2X3 排成一張圖。程式如下：
首先呼叫 PROC SGPANEL，再用 PANELBY statement 定義 Region 是類別變數。這功能就很像 CLASS statement。其餘的程式都和之前的 PROC SGPLOT 一樣，結果如下：PROC SGPANEL DATA=sg.countries;
PANELBY Region;
REG X=NumParticipants Y=TotalMedals;
TITLE 'Number of Participants by Total Medals Won for Each Country';
RUN; 
這份技術文件裡面有製作幾個表格，分別說明每種圖形要在 PROC SGPLOT 裡面用哪一種 statement，以及選了一些比較重要的 option，如下所示：
由於裡面的敘述都很簡單，各位可以自行下載原始檔案來看。不過這些語法都只是整個 PROC SGPLOT 裡面的九牛一毛而已，當然如果都學會的話應該也夠用了。如果想要知道全部的語法，可以到 SAS 官網去看。
網址：http://support.sas.com/documentation/cdl/en/grstatproc/61948/HTML/default/sgplot-stmt.htm
ABOUT THE AUTHORS
Lora Delwiche and Susan Slaughter are the authors of The Little SAS Book: A Primer, and The Little SAS Book for
Enterprise Guide which are published by SAS Institute. The authors may be contacted at:
Lora D. Delwiche
(530) 752-9321
llddelwiche@ucdavis.edu
Susan J. Slaughter
(530)756-8434
susan@avocetsolutions.com
2009年4月27日星期一
Updates to SAS® Power and Sample Size Software in SAS/STAT® 9.2
加入書籤：
原文載點：http://www2.sas.com/proceedings/forum2008/368-2008.pdf
SAS V9.2 釋出已經有一段時間，雖然仍舊有很多學校機關沒有升級到最新的版本，而仍舊沿用 V9.1.3，不過我已經用新版差不多快四個月的時間，所以之後會慢慢來介紹 V9.2 的新功能。這篇技術文件首先是來介紹新版的 proc power 以及 proc glmpower 所帶來的新的功能。
1. LOGISTIC REGRESSION
新版的 proc power 已經可以計算 logistic regression 下的 power analysis。範例如下：proc power;
logistic
alpha = 0.05
vardist(’Duration’) = normal(4, 1.5)
testpredictor = ’Duration’
testoddsratio = 1.7
responseprob = 0.65
ntotal = 50 60 70
power = . ;
run;
語法的使用方法如下：
用法和一般的 proc power 一樣，如果要計算某個樣本下的 power，則 power option 後面就要打個 dot。如果已經決定了 power，想要知道多少的樣本才能滿足，則在 power option 後面填上設定值，在 ntotal option 後面打個 dot。
此範例的報表結果如下： 由於此例一口氣設定了三個樣本數（50, 60, 70），所以報表裡面也列出三個樣本數所產生的 power 值（0.664, 0.742, 0.805）。
2. CONFIDENCE INTERVAL FOR ONE PROPORTION
新版的 proc power 可以計算一個二項變數比例的信賴區間。此功能跟 proc freq 一樣。範例如下：proc power;
onesamplefreq ci = Wilson
alpha = 0.05
proportion = 0.3
halfwidth = 0.1
ntotal = 70
probwidth = .;
run;
語法解釋如下：
此例是計算一個比例為 0.3 的信賴區間寬度。由於程式中設定整個信賴區間的長度是 0.2（halfwidth*2），所以 probwidth 算出來的結果一定會小於 halfwidth*2。結果如下： 
3. WILCOXON MANN-WHITNEY TEST FOR TWO INDEPENDENT GROUPS
Wilcoxon Mann-Whitney test 底下的 power analysis 在新版中也可以計算出來了，但只限定於兩個獨立樣本的檢定。範例如下：proc power;
twosamplewilcoxon
alpha = 0.05
vard ist(’lidocaine’) =
ordinal( (-3 -2 -1 0 1 2 3): (.01 .04 .20 .50 .20 .04 .01) )
vardist(’mironel plus lidocaine’) =
ordinal( (-3 -2 -1 0 1 2 3): (.01 .03 .15 .35 .30 .10 .06) )
variables = ’lidocaine’ | ’mironel plus lidocaine’
sides = u
ntotal = .
power = 0.85;
run;
語法解釋如下：
此例兩個要比較的群組分別是 lidocaine 和 mironel plus lidocane，要計算在 power=0.85 底下要多少樣本才夠，結果如下： 
此外，SAS V9.2 中也設計了視窗介面來執行 power analysis。其實前幾個版本的 SAS 大概覺得受到 SPSS 或其他類似視窗介面型的統計軟體的威脅，在 V6.0 之後開始加入視窗分析介面，但總是讓人有種畫虎不成反類犬的情況。V9.2 版的情況大概也不可能一下跳升到 SPSS 或 Statistica 那種程度，更何況個人是覺得真要用視窗分析介面的話，SAS 的副產品「JMP」就可以辦到了，SAS 也不用特定來搞這塊。無論如何，從這份技術文件上，我們可以看到新版的視窗分析介面的確有點改進。以下是一些截圖： 
但由於不覺得真的有到很方便的地步，所以建議還是用寫程式的方法來進行 power analysis。
CONTACT INFORMATION
Wayne Watson
Building S, Room 3040
SAS Institute, Inc.
SAS Campus Drive
Cary, NC 27513
Work Phone: 919-531-6770
E-mail: wayne.watson@sas.com
2008年11月10日星期一
Combining, Combining, Combining, Splitting, Splitting, Splitting
加入書籤：
原文載點：http://www.nesug.info/Proceedings/nesug06/dm/da30.pdf
這是一篇教導如何使用 data step 和 proc sql 合併和分割資料的 SAS 技術文件，由 Emmy Pahmer 於 NESUG 2006 發表。
本文所使用的範例如下：
這筆資料總共只有五個觀測值，每個觀測值包含三個變數：年齡、姓名和性別。其中第五個觀測值的性別是打錯的 G。
。分割資料
若想要將上述資料依照性別切割成兩個分開的資料集，則有下列幾種方式。
1. 使用兩個 data step：這是最笨但也是最直接的方式。程式碼如下所示：data males;
set everyone;
if sex = 'M';
run;
data females;
set everyone;
if sex = 'F';
run;
2. 使用兩個 data step 配合 where：這個方法僅僅是縮短程式碼行數。程式碼如下：data female;
set everyone (where=(sex=’F’)) ;
run;
data male;
set everyone (where=(sex=’M’)) ;
run;
有此程式可知，這只是把 IF 的指令改成 WHERE 並放在 SET 後面。對行數來說，的確是減少了，不過打的字變多了，因為 WHERE 後面的條件式要加上括弧。
3. 使用兩個 data step 配合 where：這個方法僅僅是把 WHERE 從 SET 移到 DATA 後面，並沒有太特別的地方。程式碼如下：data female (where= (sex = ‘F’)) ;
set everyone ;
run;
data male (where= (sex = ‘M’)) ;
set everyone ;
run;
4. 使用一個 data step 並配合 IF...ELSE... 和 OUTPUT：這是比較進階的方式，可以大幅縮短程式碼的行數。如下所示：data males females;
set everyone;
if sex = 'F' then output females;
else if sex = 'M' then output males;
run;
由於本資料中某觀測值的性別是打錯的。為了確定有沒有這類的情況，當資料相當龐大時，建議使用下列程式碼：data males females;
set everyone;
if sex = 'F' then output females;
else if sex = 'M' then output males;
else put “Neither F nor M - check “ _all_; *or output to another dataset ;
run;
其中紅色那行的程式碼會讓性別變數不是 F 和 M 的觀測值通通分類到 _all_ 這個資料集中。亦或是另外定義一個資料集把他們 output 進去。當打開這個資料集時，如果裡面是空的，就可以確定沒有打錯的情況產生。
5. 使用一個 data step 並配合 WHERE：這是從方法二改良來。程式碼如下：data female (where=(sex=’F’)) male (where=(sex=’M’)) ;
set everyone ;
run;
這個程式碼比方法四 要來的更精簡精簡。如果想要擁有方法三第二個可以偵測有沒有打錯的資料，則可以改進如下：data female (where=(sex=’F’)) male (where=(sex=’M’)) checkothers (where = (sex not in (‘M’,’F’))) ;
set everyone ;
run;
道理同方法四，把 sex 不是 M 和 F 的通通丟到 checkothers 這個資料集裡面，然後再去看看該資料集是不是空的。
6. 使用 proc sql：使用 proc sql 看起來好像比較高檔，但是行數並沒有減少。proc sql;
create table males as
select *
from everyone (where=(sex='M'));
create table females as
select *
from everyone (where=(sex='F'));
quit;
7. 使用一個 data step 把 不同的變數放到不同的資料集：若想要把 who 和 age 放到新資料集 age，然後再把 who 和 sex 放到新資料集 sex，則可仿照方法五來切割。程式如下：data age (keep = who age) sex (keep = who sex);
set everyone;
run;
8. 使用 proc sql 完成方法七：程式如下：proc sql;
create table age as
select who, age
from everyone;
create table sex as
select who, sex
from everyone;
quit;
感覺行數沒有減少很多，只是寫法比較接近口語。
。合併資料
這回使用兩筆資料，如下所示：
其中 EVERYONE 這筆和之前用的一樣，而新的 ACTIVITY 資料有一點要特別注意的是它並沒有排序過。為什麼要特別強調這一點，理由是在 SAS 的合併過程中，一定要設定一個 index variable，這樣 SAS 才有辦法依照那個 index variable 來進行資料合併。而那個被設定成 index variable 的變數一定要經過排序，否則 SAS 在合併的過程中會錯亂掉。這個錯亂有時候還是會給你 output，只是結果是錯的。如果一時忽略沒有看到 log 視窗上面的警告訊息，就完蛋了。因此若要依照「who」這個變數來合併這兩組資料，則必須要先用 PROC SORT 把該變數排序：proc sort data=activity;
by who;
run;
而通常要養成一個好習慣就是所有合併的資料最好都給他排序一下，免得有漏網之魚。反正 PROC SORT 的程式碼很簡單，多寫幾行比較安心：proc sort data=everyone;
by who;
run;
然後用 merge 和 by 來合併：data combined_11;
merge everyone activity;
by who;
run;
結果如下：
從上表得知，Annie, Bill, Chandra, Igor, Jose 和 Karen 在 ACTIVITY 裡面有資料，所以合併時會顯示在 activity 這個變數底下，但 David 和 Eleanor 則沒有出現在 ACTIVITY 裡面，所以合併後他們兩人的 activity 變數就變成 missing data 了。Age 也是同樣的道理。
如果只想顯示同時出現在兩個資料的觀測值，則必須啟用 in 這個指令。程式碼如下：data combined_12;
merge everyone (in = in_a) activity (in = in_b);
by who;
if in_a and in_b;
run;
使用 in 這個指令會讓 SAS 在合併的過程中，於兩組資料裡面各加上一個隱藏的變數，分別名為 in_a 和 in_b，其數值都預設為 1。合併之後在程式裡面加上「if in_a and in_b;」來讓 SAS 挑出同時具有 in_a=1 和 in_b=1 的觀測值（要打成「if in_a=1 and in_b=1;」也可以），缺少任一個變數的觀測值則會自動被剔除。結果如下：
如果想要知道哪些觀測值被剔除，可使用下列程式碼：data combined_14;
merge everyone (in = in_a) activity (in = in_b);
by who;
if in_a and in_b then output;
else if in_a then put "In A only: " Who=; *or output to another dataset ;
else if in_b then put "In B only: " Who=;
run;
最後那兩個 else if... 會讓程式在 log 視窗印出下列字樣：
同樣地，proc sql 也可達成同樣效果：proc sql;
create table combined_15a as
select a.*, b.activity, b.sex
from everyone as a, activity as b
where a.who = b.who
order by activity
;
quit;
或者proc sql;
create table combine_15b as
select a.*, b.activity, b.sex
from everyone as a inner join activity as b
on a.who = b.who
;
quit;
但由於 proc sql 的指令比較麻煩，所以還是建議使用 data step 來完成。
CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author at:
Emmy Pahmer
MDS Pharma Services
St. Laurent, Québec
Work Phone: (514) 333-0042 ext. 4222
E-mail: emmy.pahmer@mdsinc.com
2008年11月1日星期六
Funny ^Stuff~ in My Code – Using ODS ESCAPECHAR
加入書籤：
原文載點：http://www2.sas.com/proceedings/forum2007/099-2007.pdf
SAS 的 ODS 系統通常是拿來把報表另存成新的資料或者是列印出比較精緻的格式出來。但大多數的人可能不知道 ODS 還可以像 Word 一樣編輯報表，舉凡加上顏色、字型粗細大小、甚至是上標下標、放特殊符號，通通可以在 ODS 裡面完成。Cynthia L. Zender 於 SAS Global Forum 2007 發表的技術文件教導使用者如何來完成這些動作。
在決定使用 ODS 來編輯文件時，必須先宣告要使用哪種特殊符號來命令 ODS 開始進行一些動作。這個宣告得用 ods escapechar 來完成。比方：
ods escapechar='^';
ods escapechar='~';
ods escapechar='#';
接下來的範例接使用「^」符號來命令 ODS 編輯文件。
範例一：編輯文字
一個很簡單例子。假設要列印出某個字串變數，資料如下所示：data use_esc;
length var1 $200;
var1 = "The quick brown fox";
output;
run;
ods html;
proc print data=use_esc;
run;
ods html close;
則 ODS 只會印出下面這種相當單調的結果：
下面的程式將可改變字體的粗細、字型和字的顏色：data use_esc;
length var1 $200;
var1 = "The quick brown fox";
output;
var1 = "The ^S={font_weight=bold}quick ^S={} brown fox";
output;
var1 = "The quick ^S={foreground=brown}brown ^S={}fox";
output;
var1 = "The quick brown ^S={font_face='Courier New'}fox ^S={}";
output;
var1 = "The ^S={font_weight=bold}quick ^S={}" ||
"^S={foreground=brown}brown ^S={}"||
"^S={font_face='Courier New'}fox ^S={}";
output;
run;
ods html style=sasweb;
ods escapechar='^';
proc print data=use_esc;
run;
ods html close;
首先，在原先的 data step 裡面加上紅色部分的程式碼。其中，「^」就是用來宣告之後的 ODS 要進行某些動作。S={} 即是用來編輯。以第一個紅色程式碼來看，^S={font_weight=bold} 接在 quick 前面即是要讓 quick 變成粗體的 quick。第二個紅色程式碼 ^S={foreground=brown} 是將 brown 的顏色變更為棕色的 brown。第三個紅色程式碼 ^S={font_face='Courier New'} 即是變更 fox 的字型為 Courier New 的 fox。所有的動作完成後都必須要在後面加上 ^S={}，ODS 才會停止編輯。如果沒有 ^S={} 來中斷，則之前所宣告的編輯動作將會繼續套用在後面的內容中。第四個紅色程式碼比較複雜，這是將前三個紅色程式碼合併。此處將 quick, brown 和 fox 拆成三段，編輯過後用連接符號「||」將結果串起來。但根據我自己試驗的結果，不需要那麼麻煩，直接寫成下列的程式碼也可以做出同樣的效果：
var1 = "The ^S={font_weight=bold}quick ^S={} ^S={foreground=brown}brown ^S={} ^S={font_face='Courier New'}fox ^S={}";
最後，在宣告 ods html 之後，不要忘了加上 ods escapechar='^' 來讓 ODS 知道之前使用的「^」符號是要叫 ODS 做一些編輯動作。
執行後的結果如下所示：
範例二：加上特殊符號
特殊符號的通常用在數學方程式、註解、或化學式裡面。而使用方式也不像編輯文字一樣要用 ^S={options} 來宣告，而只要用 ^={options} 即可。假設想要做出下面這種表格：
則可使用下列程式：data sup_sub;
length myvar $200;
myvar = "Pythagorean Theorem: a^{super 2} + b^{super 2} = c^{super 2}";
output;
myvar = "This is something that needs a footnote. ^{super 1}";
output;
myvar = "Macbeth: 'Is this a dagger I see before me?' ^{dagger}";
output;
myvar = "The Caffeine molecule is an alkaloid of the methylxanthine family: " ||
"C^{sub 8}H^{sub 10}N^{sub 4}O^{sub 2}";
output;
run;
ods html file='inline2.html' style=sasweb;
ods rtf file='inline2.rtf' notoc_data;
ods pdf file='inline2.pdf';
ods escapechar='^';
proc print data=sup_sub;
title j=r 'PDF & RTF: Page ^{thispage} of ^{lastpage}';
title2 j=c 'RTF only: ^{pageof}';
footnote '^{super 1}If this were a real footnote, there would be something very academic here.';
footnote2 '^{dagger} Macbeth talked to himself a lot. This quote is from Macbeth: Act 2, Scene 1, Lines 33-39.';
run;
ods _all_ close;
在第一個紅色程式碼中，畢氏定理需要在 a, b, c 各上標一個「2」，這個動作用 ^{super 2} 即可完成。換句話說，super 這個指令在 ^{} 裡面就是拿來做上標用的。第二個紅色程式碼同樣也是要將最後的 1 給上標，所以使用 ^{super 1} 來完成。第三個紅色程式碼同樣也是要加上一個特殊符號當作註解符號。而這個小小的十字架需要使用 ^{dagger} 才能呼叫出來。第四個紅色程式碼是要編輯一個化學方程式，讓一些數字在英文字旁邊坐下標，所以這個動作就得用 ^{sub n} 來完成。接著看最後 proc print 裡面的紅色程式碼。裡面用兩個 title 指令和 footnote 指令來完成表格標題和註解的動作。其中，在 title 裡面，^{thispage} 和 ^{lastpage} 只會作用於 ods rtf 或 ods pdf 裡面。這兩個指令會自動將整份文件的頁數已經該頁是屬於第幾頁列印出來。但這個功能在 ods html 顯示不出來。所以 ^{thispage} 和 ^{lastpage} 在 html 報表裡面會變成空白。同樣地，title2 是讓標題顯示這屬於整份文件的第幾頁，使用 ^{pageof} 來執行。可是這個指令只能在 rtf 文件裡面才會有作用，在 pdf 裡面將會把 ^{pageof} 原封不動地印出來，但在 html 文件裡面則會變成空白。
這是 rtf 文件的結果：
這是 pdf 文件的結果：
範例三：折行
如果一串文字太常，如上例的第四個化學方程式，則可使用折行的指令來讓字串在某處斷行。這個動作只要簡單地在宣告符號後面加上個「n」或「an」(a=任意數字)即可。程式如下：data linebr;
length myvar $200;
myvar ="The Caffeine molecule is an alkaloid of the methylxanthine family: ~n" ||
" C~{sub 8}H~{sub 10}N~{sub 4}O~{sub 2}";
output;
myvar ="The Caffeine molecule is an alkaloid of the methylxanthine family: ~2n" ||
" C~{sub 8}H~{sub 10}N~{sub 4}O~{sub 2}";
output;
myvar ="The Caffeine molecule is an alkaloid of the methylxanthine family: ~3n" ||
" C~{sub 8}H~{sub 10}N~{sub 4}O~{sub 2}";
output;
myvar ="The Caffeine molecule is an alkaloid of the methylxanthine family: ~4n" ||
" C~{sub 8}H~{sub 10}N~{sub 4}O~{sub 2}";
output;
run;
ods listing;
ods html file='inline3.html' style=sasweb;
ods rtf file='inline3.rtf' notoc_data;
ods pdf file='inline3.pdf';
ods escapechar='~';
proc print data=linebr;
title 'Using the ODS ESCAPECHAR for line break';
title2 'Title 2';
title3 'Title 3';
title4 'Title 4';
title5 'Title 5';
title6 'Title 6';
title7 'Title 7';
title8 'Title 8';
title9 'Title 9';
title10 'Title 10 ~n Title 11 ~n Title 12 ~n Title 13 ~n Title 14';
run;
ods _all_ close;
為了與前兩例做一點區別，此例將宣告符號「^」改成「~」。在 data step 裡面，分別將 The Caffeine molecule is an alkaloid of the methylxanthine family: 和 化學方程式之間各斷一、二、三、四行，因此只要在冒號後面打上「~n」,「~2n」,「~3n」,「~4n」即可。設定完後不要忘了在 proc print 前面加上 ods escapecha='~'; ，否則之前的宣告不會被執行。至於下面很多 title 指令是要展示折行的效果。如果不使用折行，折需要連續呼叫 titlen statement。但有了折行指令，則只需要打一行，裡面再連續使用數個 ~n 即可（如最後的紅色程式碼所示）。僅列出 html 的結果如下：
例四：更改欄位設定
以 ods html 來說，如果沒有特別的設定，則一般的輸出報表會像下圖所示：
如果想要將整行從第一個 Sales 處斷行，並把 Sales Total Sales 改成斜體的話，可使用下列程式碼：ods html style=sasweb;
ods escapechar='^';
proc means data=sashelp.shoes sum min mean max;
class Region;
var Sales;
label Sales='^S={font_style=italic}^nShoe Sales';
run;
ods html close;
紅色部分就是用來更改那個表格內標籤的設定。^S={font_style=italic} 是用來將後面的 Total Sales 改成斜體字。而緊接著的 ^n 就是用來斷行的。結果如下：
範例五：特殊符號
有些特殊符號需要使用特定的指令來印出。承上例，我們想要在表格後面加上註腳，並打上 trade mark、copy right mark 和美分符號，程式如下：data _null_;
hexcode = 'AE'x;
call symput ('hexreg',hexcode);
hexcode1 = 'A9'x;
call symput('hexcopy' , hexcode1);
hexcode2 = 'A2'x;
call symput('hexcent', hexcode2);
run;
** Display the hex codes in the SAS log;
%put hexreg= &hexreg;
%put hexcopy= &hexcopy;
%put hexcent= &hexcent;
ods html file='inline5.html' style=sasweb;
ods rtf file='inline5.rtf' notoc_data startpage=no;
ods pdf file='inline5.pdf' startpage=no;
ods escapechar='^';
ods noptitle;
proc means data=sashelp.shoes sum min mean max;
class Region;
var Sales;
label Sales='^S={font_style=italic}^nShoe Sales';
run;
ods pdf text="^S={just=c font_size=18pt}PDF File^nOpens with Adobe Reader^n
&hexreg &hexcopy My 2&hexcent";
ods rtf text="^S={just=c font_size=18pt}RTF File^nOpens with Microsoft Word^n
&hexreg &hexcopy My 2&hexcent";
ods html text="^S={just=c font_size=18pt}HTML File^nOpens with Microsoft Internet
Explorer^n &hexreg &hexcopy My 2&hexcent";
proc freq data=sashelp.shoes;
tables Product /nocum;
label Product='^S={font_style=italic}Shoe Types Sold';
run;
ods _all_ close;
使用者必須先用一個 _null_ 的 data step，呼叫 %symput 函式把這三個特殊符號叫進來。這三個特殊符號在 SAS 系統裡都有自己的代稱，trade mark 的代碼是 'AE'x，copy right mark 的代碼是 'A9'x，而美分的代碼是 'A2'x。然後用 %put 指令把剛剛用 %symput 叫進來的三個代碼設定為巨集參數，供之後使用。而那些註解可以直接寫在 ods 後面，使用 text 來印出，而所有諸如編輯文字、斷行以及安插特殊符號的指令全部都可以在 text 後面使用。以 ods html 來講， ^S={just=c font_size=18pt} 即是設定整個註腳的字型大小為 18pt，並且全部置中。放在 File 和 Opens 中間的 ^n 即表示從這個位置斷行，最後使用 &hexreg、&hexcopy 和 &hexcent 把 trade mark, copy right mark 和美分符號安插進去。結果如下：
範例六：安插圖片
ODS 也可以拿來安插圖片。先看看程式：proc format;
value $prdimg 'Boot' = 'c:\temp\boot.jpg'
"Women's Dress" = 'c:\temp\dressheels.jpg'
'Slipper' = 'c:\temp\slipper.jpg';
run;
ods html file='inline6.html' style=sasweb;
ods rtf file='inline6.rtf' notoc_data;
ods pdf file='inline6.pdf';
ods escapechar='^';
proc report data=sashelp.shoes nowd style(column) = {vjust=b};
title '^S={preimage="c:\temp\shoe_logo.jpg"}';
where product in ('Boot', 'Slipper', "Women's Dress");
column Product Sales N;
define Product / group
style(header)={background=white foreground=black}
style(column)={postimage=$prdimg.};
define Sales / sum;
define N / '^S={background=white foreground=black}Number of Sales';
run;
ods _all_ close;
如果只要在 title 上安插圖形的話，則可簡單地於 title statement 後面加上 ^S={preimage="路徑\檔名.jpg"} 來把圖形放上。如果要在文件內容部分的話，得需要先用 proc format 把每一張圖的路徑和檔名分別設定給某個變數的不同的 value，並用一個 $prdimg 來表示（這個字串可以自行設定）。然後在 proc report 裡面的 define 後面放上選項 style(column)={postimage=prdimg.}，則每一欄裡面都會安插不同的圖形。若只想單純的更改欄位背景或前景顏色的話，則使用 ^S={backgroun=color1 foreground=color2} 來改變。此例使用 ^S={background=white foreground=black} 把總數欄位的前景顏色改成黑色，背景顏色改成白色。
三種輸出格式的結果如下：
CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author:
Cynthia L. Zender
SAS Institute Inc.
Home Office
Las Cruces, NM 88011
Work Phone: (505) 522-3803
Fax: (505) 521-9328
E-mail: Cynthia.Zender@sas.com
2008年10月23日星期四
Performing Iterative Processes with the Macro Facility
加入書籤：
原文載點：http://www.nesug.info/Proceedings/nesug07/cc/cc21.pdf
有時候要針對一組資料裡面不同類別的觀測值進行分析，在所有的 PROC 程序裡面都可以利用 BY 這個語法來完成。但是 BY 並沒有辦法在某些特殊情況下使用，因此利用 macro 進行遞迴式的處理就變成一另一種可行的方式。Katie Joseph 和 Taylor Lewis 在 NESUG 2007 發表了一篇利用 macro + do loop + %scan 來完成 BY 所無法完成的工作，可供有需要的人參考。
假設有一組資料如下所示：
每個 agency 有數筆重複觀測資料，Q1-Q73 為 73 個類別變數，weight 是連續變數表「權重」，而 strata 則是一個獨立的分層變數，同一個 agency 可能有兩種以上不同的分層變數。
假設我們要計算抽樣母體的平均值，則可使用下列程式：proc surveymeans data=survey total=frametotals;
strata STRATA;
var Q1-Q73;
weight weight;
domain agency;
ods output domain=outstats;
run;
這個程式特殊的地方在於他有使用 domain statement 來執行 domain analysis。關於 domain analysis 可以參考下面這個網址：
http://support.sas.com/rnd/app/da/new/801ce/stat/chap13/sect9.htm#smeansdomainanalysis
Domain analysis 在變數分層多的時候會消耗大量電腦記憶體，而此例的 agency 總共有 87 層，一般電腦可能會出現「out of memory」的訊息而中斷程序。因此得轉個彎讓每個 agency 跑一次 proc surveymeans，然後一個一個去跑 domain analysis 就不會有記憶體不足的問題。假設只跑 agy1 這一層，則程式如下：data split;
set survey;
length split $1;
if agency='agy1' then split='Y';
else split='N';
run;
proc surveymeans data=split total=frametotals;
strata STRATA;
var Q1-Q73;
weight weight;
domain split;
ods output domain=outstats_agy1 (where=(split='Y'));
run;
第一個 data step 設定一個新的變數 split，當 agency=agy1 時則 split = 'Y'，反之則 split = 'N'。然後執行 proc surveymeans 時，把 split 放進 domain statement 裡面。由於 split 只是個二項變數，因此 domain 只要處理兩層就好。然後用一個 ods output 把 domain analysis 的表格輸出並存成 outstats_agy1。由於我們只關心 agy1 的結果，所以只要留下 split='Y' 的部分，至於 split='N' 的部分就不是我們所關心的了，這就是 ods output 後面要用一個 where 來抑制 split='N' 部分的輸出。
根據這個程式，我們可用一個 macro + do loop 讓他對每個 agency 都執行一次 proc surveymeans 程序。程式如下：%macro agydoloop();
%do agynum=1 %to 87;
data split;
set survey;
length split $1;
if agency="agy&agynum" then split='Y';
else split='N';
run;
proc surveymeans data=split total=frametotals;
strata STRATA;
var Q1-Q73;
weight weight;
domain split;
ods output domain=outstats_agy&agynum (where=(split='Y'));
run;
%end;
%mend agydoloop;
這個程式首先先用 %macro agydoloop(); 和 %mend agydoloop; 把整個程式包起來。由於沒有額外的巨集參數需要使用，所以 %macro agydoloop(); 裡面就流空白即可。然後裡面再用一個 %do agynum=1 %to 87; 和 %end; 把那個 data step 以及 proc surveymeans 程序給包起來。特別注意是 do, to 和 end 若使用在 macro 裡面時，前面需加上百分比符號「%」。之後再把 1 改成 &agynum 即可，這樣讓 do loop 在跑時能把 1 到 87 遞迴地代入 agynum 裡面。如此一來 87 次 proc surveymeans 就可以輕鬆完成了。
如果不想一次把 Q1 到 Q73 跑完，而也想用 do loop 分批跑的話，只要再加上一層 do loop 即可。程式如下：%macro agyandQdoloop();
%do agynum=1 %to 87;
data split;
set survey;
length split $1;
if agency="agy&agynum" then split='Y';
else split='N';
run;
%do qnum=1 %to 73;
proc surveymeans data=split total=frametotals;
strata STRATA;
var Q&qnum;
weight weight;
domain split;
ods output domain=outstats_agy&agynum._Q&qnum (where=(split='Y'));
run;
%end; /* qnum do loop */
%end; /* agynum do loop */
%mend agyandQdoloop;
由於 Qn 變數只出現在 proc surveymeans，所以新的 do loop 只需要放在 proc surveymeans 的頭尾即可。然後用 &qnum 來替換數字 1 到 73。特別一提的是，此處的 ods output 在設定輸出資料名稱時，原先的 &agynum 由於後面接了一個特殊符號「_」，因此 &agynum 後面要先打上一個句點「.」再接「_」。如果沒有打上句點，則這個 macro 會無法分辨出 &agynum 是個參數而非定值。這個規則通用於各種情況，只要巨集變數後面接上符號如「_」,「.」或「\」，都需要打上一個句點先。
可是，當類別變數全都是字串時，do loop 就無法處理了。因此需要利用 %scan 來輔助完成。%scan 函式是用來切割一長串的字串。其語法為：%scan(參數, n [,分隔符號]);
其中，n 代表要抽出切割後的第幾個字串，而[,分隔符號]是一個選擇性的參數。在沒有額外設定的情況下，預設值為「空白」,「.」,「<」,「(」,「+」,「|」,「&」,「$」,「*」,「)」,「;」,「﹁」,「-」,「/」,「,」,「%」,「￠」。
假設 agency 變數不再是 agy1-agy87，而是 ED, BO, CM 這三個字串，則程式如下：%macro agyscanloop(agylist=);
%let num=1;
%let agy=%scan(&agylist,&num);
%do %while (&agy ne );
data split;
set survey;
length split $1;
if agency="&agy" then split='Y';
else split='N';
run;
proc surveymeans data=split total=frametotals;
strata STRATA;
var Q1-Q73;
weight weight;
domain split;
ods output domain=outstats_&agy(where=(split='Y'));
run;
%let num=%eval(&num+1); *** increase by one the position pointer;
%let agy=%scan(&agylist,&num);
%end;
%mend agyscanloop;
其中，agy=%scan(&agylist,&num)是要把巨集參數 agylist 裡面所打的字串切割並取第 &num 個部分。由於一開始 &num 設定為 1，所以每跑完一次之後要把 &num 加一，這樣跑第二次時就會去取切割後的第二個字串。而遞迴的過程是由 %do %while (&agy ne); 來驅動。這個道理很簡單，只要當 &agy=%scan(&agylist, &num) 的結果不是空白（即 missing data），則這個 do while 就會一直跑，直到 &agy 是 missing data 為止。要怎樣讓這個 do while 結束，就是靠之後得 %let num=%eval(&num+1); 和 %let agy=%scan(&agylist, &num); 來控制。當跑完第一次時，num=%eval(&num+1) 會讓 &num = 2。之所以要用 %eval 的原因是因為任何四則運算在 macro 裡面都會被視為字串。只有包在 %eval 函式裡面的四則運算過程才會真正被計算出來並輸出結果。所以如果沒有用 %eval 函式的話，num 會變成字串型的「1+1」，而非數值型的「2」。當 num=2 後代入下面那個 %let agy=%scan(&agylist, &num); 時，這行就會把 &agylist 所代表的字串再次切割，然後取出第二個部分。所以這個 do while 就會繼續去執行，直到 num=3 時才會跳出。
因此當執行下面這行時：%agyscanloop(agylist=ED BO CM);
此巨集會把「ED BO CM」從中間的空白處切成「ED」,「BO」以及「CM」，然後再一個一個去執行 proc surveymeans。
雖然 %scan 可以幫忙切割字串，但是如果 agency 這個變數有數十個字串型的類別，而非上例的只有三個類別，則在輸入巨集參數 &agylist 時便會顯得相當相當耗時。作者提出了一個 proc sql 程序來解決這項煩人的工作。程式改寫如下：proc sql;
select distinct agy into :agylist separated by ' '
from survey;
%put &agylist;
quit;
%macro agySQLscanloop();
%let num=1;
%let agy=%scan(&agylist,&num);
%do %while (&agy ne );
*** PROC SURVEYMEANS Code - uses &agy as the agency code;
%let num=%eval(&num+1);
%let agy=%scan((&agylist,&num);
%end;
%mend agySQLscanloop;
這個 proc sql 會把 agy 所有的類別黏成一串，並用空白符號分隔開來，然後把黏成一串的結果用 %put 函式放進一個巨集變數 &agylist 裡面。由於這個 &agylist 變數已經在外部設定好了，所以之後寫 macro 時就不用在 %macro agySQLscanloop(); 裡面放 agylist 了。至於 macro 裡面的內容和設定則完全跟前例一樣。
CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the authors at:
Katie Joseph
US Office of Personnel Management
1900 E St, NW, Room 7439
Washington, DC 20415
Work Phone: (202) 606-1817
Fax: (202) 606-1719
E-mail: Katie.Joseph@opm.gov
Taylor Lewis
U.S. Office of Personnel Management (OPM)
1900 E St., NW, Room 7439
Washington, DC 20415
Work Phone: (202) 606-1309
Fax: (202) 606-1719
E-mail: Taylor.Lewis@opm.gov
2008年10月21日星期二
A macro for nearest neighbor imputation
加入書籤：
原文載點：http://www.mediafire.com/?sharekey=00ab0205b2465593dd8b33b5aa27078d
長期以來 SAS 內部只能使用 PROC MI 來執行多重插補法的資料插補動作，而當然世界上不可能只有這麼一種方法。由於多重插補法的使用時機需要配合他強大的假定，所以一旦資料不適用 PROC MI 時，SAS 使用者往往沒有別的替代方案。另外，有時缺失值的情況不嚴重，但使用 Complete Case Analysis 直接刪除有缺失值的觀測值，往往對其他沒有缺失值的變數造成一些資料上的浪費。使用 PROC MI 又需要驗證許多假定，還得在事後進行 PROC MIANALYSIS，實在相當浪費時間。如果有個簡單但又不失可靠的替代方案，將可大大增加處理缺失值的效率。Nearest neighbor imputation（以下簡稱 NNI）是個形之多年的資料插補法，但 SAS 現存的程序並不支援。這篇於 SESUG 2008 由 Lung-Chang Chien (我本人)和 Mark Weaver (我的committee member)所發表的技術文件，提供了一個簡易的 macro 來執行 NNI 的插補動作。
首先，要瞭解 Nearest neighbor imputation 的理論，可以參考這一篇文章：
Nearest。Neighbor。Imputation
我所設計的 %NNI 一共只有五個語法便可完成插補的動作。原始碼如下：%MACRO NNI(INDATA=, /*INPUT DATA SET(INCLUDING LIBRARY NAME) */
MISSVAR=, /*INPUT SINGLE VARIABLE WITH MISSING DATA*/
RESPVAR=, /*RESPONSE VARIABLE*/
IDVAR=, /*SUBJECT VARIABLE*/
OUTDATA= /*OUTPUT DATA SET WITH COMPLETE DATA*/
);
DATA OBSY_MISSX;
SET &INDATA;
IF &MISSVAR=.;
KEEP &IDVAR &RESPVAR;
RUN;
PROC MEANS DATA=&INDATA N NMISS NOPRINT;
VAR &MISSVAR;
OUTPUT OUT=OBSN N=OBSX NMISS=MISSX;
RUN;
DATA _NULL_;
SET OBSN;
CALL SYMPUT('OBSN',OBSX);
CALL SYMPUT('MISSN',MISSX);
RUN;
DATA SIMDATA_NOMISS;
SET &INDATA;
IF &MISSVAR NE .;
RUN;
PROC IML;
USE SIMDATA_NOMISS;
READ ALL VAR {&MISSVAR} INTO XMAT;
READ ALL VAR {&RESPVAR} INTO YMAT;
USE OBSY_MISSX;
READ ALL VAR {&IDVAR &RESPVAR} INTO MISSMAT;
TXMAT=T(XMAT);
TYMAT=T(YMAT);
DISTANCE=J(&MISSN,&OBSN,.);
MIND=J(&MISSN,1,.);
MISSVAR=J(&MISSN,&OBSN,.);
IMP=J(&MISSN,1,.);
DO I = 1 TO &MISSN;
DO J = 1 TO &OBSN;
DISTANCE[I,J]=ABS(MISSMAT[I,2]-TYMAT[,J]);
END;
MIND[I,]=MIN(DISTANCE[I,]);
END;
DO I = 1 TO &MISSN;
DO J = 1 TO &OBSN;
IF DISTANCE[I,J]=MIND[I,] THEN MISSVAR[I,J]=TXMAT[,J];
END;
IMP[I,]=MISSVAR[I,:];
END;
CNAME={"&IDVAR" "&RESPVAR" "&MISSVAR"};
IMPX=MISSMAT||IMP;
CREATE NNI FROM IMPX[C=CNAME];
APPEND FROM IMPX;
QUIT;
DATA &OUTDATA;
MERGE &INDATA NNI;
BY &IDVAR;
RUN;
%MEND;
語法使用方式如下：
這邊使用一個實際範例來說明如何使用 %NNI。我使用一個位於北卡 Raleigh 的空氣品質監測站的資料，使用裡面的 PM2.5 和二十四小時日均溫這兩個變數，時間點取在西元兩千年，所以總共有 366 個觀測值。其中氣溫變數是完整的，但是 PM2.5 有十七個 missing data。我用這個資料來配適下面這個模式：
程式如下：%NNI(INDATA=CASE.RALFORNNI,
MISSVAR=PM25TMEAN,
RESPVAR=TMPD,
IDVAR=DATE,
OUTDATA=CASE.RAL_NNI);
這張散佈圖說明了插補資料所在的位置：
由於 Nearest neighbor imputation 是一個比較保守的插補法，所以插補值很依賴既有資料的變異程度。如果既有資料的變異程度比較小，那能夠拿來當作插補值的變化也會跟著變小。不過至少不可能插補到離群值，這一點是絕對可以保證的。
再來看看 PM2.5 在插補前後的差異：
可以發現數值都很逼近。
經過 PROC GAM 程序的配適後，比較一下最重要的兩個參數估計值的結果：
從上表可知彼此的差異不大。顯著情況也沒有改變。
再來看看平滑曲線圖：
兩張圖的情況也非常近似。不過這邊要特別強調一點，那就是如果 missing data 是出現在要配適平滑曲線的變數時，Nearest neighbor imputation 的成效並不好。這一點已經經過我用另外一個模擬結果證實了，只是沒放在這篇技術文件裡面。這個實例之所以會有相似的平滑曲線是因為 missing data 只出現在線性的應變數 PM2.5 裡面，而非使用在平滑曲線的時間變數裡面。
Nearest neighbor imputation 還是有一些限制，整理如下：
1. 不太適用於插補離散變數，但如果資料夠大，能來搭配的完整變數內的數值變化也夠多的話，還是可以勉強使用。
2. 資料一定要多，太小的樣本數會大幅降低插補功效。
3. 資料裡面一定至少要有一個完整的變數，如果所有變數都有 missing data，則此方法和程式都會失效。
4. 目前這個 macro 的版本只能在 MISSVAR 和 RESPVAR 裡面各放一個變數。如果有數個變數含有 missing data，則需要重複使用這個 %NNI 來插補。如果情況是完整的變數有很多個，那到底要放哪一個變數在 RESPVAR 裡面。目前為止並沒有什麼準則來決定最好的完整變數。根據我個人的經驗，最好先用 PROC CORR 把所有變數的相關矩陣弄出來，然後挑出那個和有 missing data 變數有最高相關程度的完整變數即可。
此外，這次在 St. Pete Beach 舉辦的 SESUG 2008 會議，我也有去會場發表。整個過程可以看下面這篇文章：St。Pete。Beach
比較令人印象深刻的是我在 Tampa International Airport 等交通車時遇到了 Wendi Wright 女士。本部落格曾經發表過兩篇她寫的 SAS 技術文件：
A Legend is Not Just a Legend
Checking for Duplicates
她今年也在 SESUG 發表了兩篇技術文件，有空我再來分享她的著作。
CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author at:
Lung-Chang Chien
University of North Carolina at Chapel Hill
Chapel Hill, NC 27599
E-mail: cchien@email.unc.edu
2008年10月14日星期二
A Legend is Not Just a Legend
加入書籤：
原文載點：http://www.nesug.org/Proceedings/nesug07/po/po21.pdf
這是一篇關於介紹 SAS/GRAPH 裡面圖例（legend）的基本寫法和一些變化，由 Wendi Wright 在 NESUG 2007 所發表。
假設有個程式如下所示：TITLE font=’Times New Roman’ height=1.5 ‘Number of Hits on Websites 1, 2, and 3’;
TITLE2 font=’Times New Roman’ height=1.5 ‘For the Month of ‘ COLOT=red ‘March 2007’;
FOOTNOTE JUSTIFY=left ‘Educational Testing Service’ JUSTIFY=right ‘April 1, 2007’;
AXIS1 LABEL=(ANGLE=270 ROTATE=90 HEIGHT=1.5 ‘Number of Hits’)
ORDER=(0 to 1800 by 200) MINOR=(NUMBER=3);
AXIS2 REFLABEL=(POSITION=top JUSTIFY=center ‘Email Ad‘)
VALUE=(‘1’ ‘2’ ‘3’ ‘4’ ‘5’ ‘6’ ‘7’ ‘8’ ‘9’);
SYMBOL1 COLOR=blue INTERPOL=join LINE=1 VALUE=dot;
SYMBOL2 COLOR=red INTERPOL=join LINE=2 VALUE=star;
SYMBOL3 COLOR=green INTERPOL=join LNEe=3 VALUE=circle;
PROC GPLOT DATA=perm.hits;
PLOT Web1*day Web2*day Web3*day
/ OVERLAY HREF = ‘17’ FRAME VAXIS=axis1 HAXIS=axis2;
RUN;
能夠畫出這樣的圖：
以此圖為基準，我們來看看如果加上圖例，並且做一些變化。
最基本的圖例就是依照 SAS 內部的設定將三條折線的點和線的型態標註在整張圖的最下方。這個動作只需要在 plot statement 後面加上 legend 即可，如下所示：PROC GPLOT DATA=perm.hits;
PLOT Web1*day Web2*day Web3*day
/ OVERLAY HREF = ‘17’ FRAME VAXIS=axis1 HAXIS=axis2 LEGEND;
Run;
SAS 會輸出下面這種圖：
如果想要換掉 legend 裡面三條折線的標籤，則可以在 PROC GPLOT 程序前面加上一條 legendn 的指令來修改。其中 n=1,2,3... 表示可以寫 n 條 legend 的指令。之後再把寫好的 legend statement 掛在 legend option 的後面，如下所示：LEGEND1 VALUE=(COLOR=blue HEIGHT=1 ‘Web Site 1’ ‘Web Site 2’ ‘Web Site 3’);
PROC GPLOT DATA=perm.hits;
PLOT Web1*day Web2*day Web3*day
/ OVERLAY HREF = ‘17’ FRAME VAXIS=axis1 HAXIS=axis2 LEGEND=legend1;
Run;
在 legend1 指令中，value 可以設定標籤的顏色（color option）、大小（height option）還有標籤字樣（此例用 'Web Site 1' 'Web Site 2' 'Web Site 3'，記得要加引號，不要加逗號）。改好後得圖例會變成這樣：
如果也想一併把 legend 的名字（此例為「PLOT」）改掉，則可使在 legendn 指令裡面再加上一個 label 的選項，如下所示：LEGEND1 LABEL=(COLOR=blue HEIGHT=2 ‘Our Web Pages’)
VALUE=(COLOR=blue HEIGHT=1 ‘Web Site 1’ ‘Web Site 2’ ‘Web Site 3’);
PROC GPLOT data=perm.hits;
PLOT Web1*day Web2*day Web3*day
/ OVERLAY HREF = ‘17’ FRAME VAXIS=axis1 HAXIS=axis2 LEGEND=legend1;
Run;
同樣地，label 裡面也可以設定顏色大小和想要用的字，語法和 value 一模一樣。改出來的效果如下所示：
如果不想讓 legend 橫的排列，而想改成直的，則需使用 across 和 down 兩個選項控制。across 可以表示行數，down 可以表示列數。因此此例要讓 legend 直的排，就等於是要宣告 legend 裡面只要一行三列，程式改寫如下：LEGEND1 LABEL=(HEIGHT=1 ‘Our Web Pages’)
VALUE=(‘Web Site 1’ ‘Web Site 2’ ‘Web Site 3’)
ACROSS=1 DOWN=3;
PROC GPLOT DATA=perm.hits;
PLOT Web1*day Web2*day Web3*day
/ OVERLAY HREF = ‘17’ FRAME VAXIS=axis1 HAXIS axis2 LEGEND=legend1;
Run;
此時 legend 就會站起來了！
若覺得這樣還不夠，想要讓 legend 名稱「Our Web Pages」也放在圖例的上頭，則需要在 legend 裡面的 label option 多加兩個指令。一個是 position，用來指定 legend 名稱的位置，有上下左右（top, bottom, left, right）可選，然後再用 justify 指令來對齊，有置左(left)、置中(center)、和置右(right)可選。程式如下：LEGEND1 LABEL=(HEIGHT=1 POSITION=top JUSTIFY=center ‘Web Sites’)
VALUE=(‘Web Site 1’ ‘Web Site 2’ ‘Web Site 3’)
ACROSS=1 DOWN=3;
PROC GPLOT DATA=perm.hits;
PLOT Web1*day Web2*day Web3*day
/ OVERLAY HREF = ‘17’ FRAME VAXIS=axis1 HAXIS=axis2 LEGEND=legend1;
Run;
此程式是將他放在上頭並置中，所以用 position=top 和 justify=center 來調整。結果如下：
若覺得這樣的圖例太大很佔空間，想要把他移到圖內，只要在 position 裡面加上一個 inside 選項即可。如果要把圖例移到左上角，則可以再同時宣告 top 和 left 於 position 裡面。換句話說，position 可以一次宣告三種位置。另外，如果擔心圖例會遮到既有的點或線，則可以使用 mode 來決定要不要覆蓋。要的話設定為 protect，不要的話則設定為 share。此例設定為要覆蓋，所以使用 mode=protect。程式改寫如下：LEGEND1 LABEL=(HEIGHT=1 POSITION=top JUSTIFY=center ‘Web Sites’)
VALUE=(‘Web Site 1’ ‘Web Site 2’ ‘Web Site 3’)
ACROSS=1 DOWN=3
POSITION = (top left inside)
MODE=protect;
PROC GPLOT DATA=perm.hits;
PLOT Web1*day Web2*day Web3*day
/ OVERLAY HREF = ‘17’ FRAME VAXIS=axis1 HAXIS=axis2 LEGEND=legend1;
Run;
圖形如下：
如果覺得這個 legend 太單調，想要多加一個有顏色的框，則需使用 cframe 來設定顏色，如下所示：LEGEND1 LABEL=(HEIGHT=1 POSITION=top JUSTIFY=center ‘Web Sites’)
VALUE=(‘Web Site 1’ ‘Web Site 2’ ‘Web Site 3’)
ACROSS=1 DOWN=3
POSITION = (top left inside)
MODE=protect
CFRAME = yellow;
PROC GPLOT DATA=perm.hits;
PLOT Web1*day Web2*day Web3*day
/ OVERLYA HREF = ‘17’ FRAME VAXIS=axis1 HAXIS=axis2 LEGEND=legend1;
Run;
然後你就會得到一個黃色的框架：
當然框架也是可以做一些微調的，比方說若覺得他太靠近 Y 軸的話，可利用 offset 來進行調整。以此圖為例，如果想要將 legend 框架和 Y 軸拉開差不多整張圖的 3% 的距離，則可以下列程式：LEGEND1 LABEL=(HEIGHT=1 POSITION=top JUSTIFY=center ‘Web Sites’)
VALUE=(‘Web Site 1’ ‘Web Site 2’ ‘Web Site 3’)
ACROSS=1 DOWN=3
POSITION = (top left inside)
MODE=protect
CFRAME= yellow
OFFSET = (3 pct);
PROC GPLOT DATA=perm.hits;
PLOT Web1*day Web2*day Web3*day
/ OVERLAY HREF = ‘17’ FRAME VAXIS=axis1 HAXIS=axis2 LEGEND=legend1;
Run;
成效如下：
特別一提的是，這種調整通常很難一次到位，所以使用者需要輸入不同的數值以求達到自己最滿意的結果。
這框架也可以更改長寬大小，只要在 legend 裡面使用 shape=(w,h) 語法就可設定，其中 w 代表寬度，h 代表高度。預設值是 (w,h)=(5,1)。此例若要加長寬度，則需用 shape=(10,1)，如下所示：LEGEND1 LABEL=(HEIGHT=1 POSITION=top JUSTIFY=center ‘Web Sites’)
VALUE=(‘Web Site 1’ ‘Web Site 2’ ‘Web Site 3’)
ACROSS=1 DOWN=3
POSITION = (top left inside)
MODE=protect
CFRAME = yellow
OFFSET = (3 pct)
SHAPE=(10,1);
PROC GPLOT DATA=perm.hits;
PLOT Web1*day Web2*day Web3*day
/ OVERLAY HREF = ‘17’ FRAME VAXIS=axis1 HAXIS=axis2 LEGEND=legend1;
RUN;
當然 legend 指令還有很多，不過這些基本的設定可以應付大部分的情況，如果想要再深入研究其他指令，可以依照下面的路徑去找網路上的手冊：
AUTHOR CONTACT
Your comments and questions are valued and welcome. Contact the author at:
Wendi L. Wright
1351 Fishing Creek Valley Rd.
Harrisburg, PA 17112
Phone: (717) 513-0027
E-mail: wendi_wright@ctb.com
訂閱：
文章 (Atom)
發表文章
