公告

[公告]
2014/01/17
由於已經是faculty的關係,不太有足夠時間寫部落格。因此更新的速度會相當緩慢。再加上近幾年來SAS GLOBAL FORUM沒有出現讓我覺得驚艷的技術文件,所以能分享的文章相對也減少許多。若有人推薦值得分享的SAS技術文件,請利用『問題討論區』告知。

2013/07/19
臉書留言板的功能因為有不明原因故障,因此特此移除。而intensedebate的留言板因管理不易,也一併移除。目前已經開啟內建的 G+ 留言系統,所以請有需要留言的朋友,可直接至『問題討論區』裡面留言。


2007年3月10日 星期六

A Macro to Calculate Kappa Statistics for Categorizations by Multiple Raters

原文載點:http://www2.sas.com/proceedings/sugi30/155-30.pdf

在類別資料分析中,我們經常會看到這樣的設計:訪問一組人兩次,然後檢定兩次的答案有沒有一致性。這樣的統計分析,可以用 Kappa 統計量來檢定。可是,當重複訪問的次數超過兩次,那 Kappa 統計量還可以用嗎?Bin Chen,Dennis Zaebst 和 Lynn Seel 等人在 2005 年的 SUGI 30 發表了利用 Kappa 統計量做多重一致性比較的分析,並分享一個新的 macro 在文章裡面。


KAPPA STATISTICS


Kappa 統計量首先是在 1960 年由 Cohen 所發明出來,之後一直被廣泛應用到現在。原始的理論是架構在 2X2 table 上,現在就先稍微介紹 nXn table (n≧2) 時的理論基礎。

假設有個評量(以下用 rater 來表示)裡面有 j=1~k 個選項,則在 n 個人重複測了 m 次之後,第 j 個選項被勾選的總比例是:

sugi155_30_01

其中分子的 Xij 表示第 i 個人勾選第 j 個選項的次數。則第 j 個選項的 Kappa 統計量可以用下式表示:

sugi155_30_02

因此,總 Kappa 統計量就是上式的加權平均:

sugi155_30_03

其標準誤是:

sugi155_30_04

MKAPPA MACRO

作者們開發了一個視窗化的 macro 程式,名為 MKAPPA。使用者可以直接在視窗內輸入一些設定即可完成多重一致性的檢定,而不需要寫上任何一個程式碼。視窗畫面如下所示:

sugi155_30_05

APPLICATION EXAMPLE


作者簡單地提供了一個小小的範例來示範這個程式的用法。假設有下面這樣的資料:

sugi155_30_06

每個人一份評量做三次,所以總共有三個 rater,每個 rater 內有五個選項。經過簡單的統計之後,每個 rater 內各選項被選到的次數分配表如下所示:

sugi155_30_07

經過 MKAPPA 的計算,結果如下:

sugi155_30_08

由於 p-value 小於顯著水準 0.05,拒絕虛無假設,因此可以下一個結論,那就是受測群三次做出來的結果並沒有一致性存在。

這個 MKAPPA 還是有一些使用上的限制,主要是處理 missing data 上。作者有提到一旦出現 missing data 的話程式便會強迫終止。但是在範例中的 rater1 有一個 missing data,程式卻仍舊可以執行。初步的研判應該是作者將 missing data 獨立成另一個選項來避免程式執行失敗。

最後,MKAPPA 的 macro 如下所示:

%window setmacro rows = 18
#2 "Data location : " dloc 80
#4 "Dataset name : " dset 40
#6 "Number of raters : " Nrater 8
#8 "Number of rating categories : " Ncat 8
#10 "Output file name : " outfile
;
%let dloc = c:\documents and settings\bin chen\my documents\projects\KappaMacro\ ;
%let dset = testno1 ;
%let Nrater = 3 ;
%let Ncat = 5 ;
%let outfile = output.rtf ;
%display setmacro ;
%macro mkappa;
libname kap "&dloc";
data temp; set kap.&dset;
subj=_n_;
data temp; set temp;

%do i=1 %to &nrater;
rate=rater&i; rater=&i; output; %end;
run;
proc freq data=temp noprint; tables subj*rate / out=xij ;
run;
data subj; set xij; keep subj;
proc sort data=subj nodupkey out=subj; by subj; run;
proc sql noprint; select count(*) into :tsub from subj;
data temp2; do subj=1 to &tsub; do rate=1 to &ncat; output; end; end; run;
data xij; set xij; drop percent;
proc sort data=xij; by subj rate;
proc sort data=temp2; by subj rate;
data xij; merge xij temp2(in=a); by subj rate; if a;
data xij; set xij; if count=. then count=0; run ;
data xij; set xij; xmx=count*(&nrater-count); run;
proc sort data=xij; by rate;
data xij; set xij; by rate;
retain x xmx2;
if first.rate then do; x=0; xmx2=0; end;
x=x+count; xmx2=xmx2+xmx;
if last.rate then output;
run;
data xij; set xij; keep rate x xmx2 ; run;
data xij; set xij;
p=x/(&tsub*&nrater);
pq=p*(1-P);
pqp=p*(1-p)*(1-2*p);
kj=abs(1 - xmx2/(&tsub*&nrater*(&nrater-1)*pq));
numj=kj*pq;
run;
data xij; set xij; if p=0 then do;
kj=0; numj=0; end; run;
data final; set xij end=end; retain num den pqpsum;
if _n_=1 then do; num=0; den=0; pqpsum=0; end;
num=num+numj;
den=den+pq;
pqpsum=pqpsum+pqp;
if end then output;
run;
data kap.koutput; set final; keep Kappa SE pvalue;
Kappa=num/den;
SE=sqrt(2*(den*den-pqpsum))/(den*sqrt(&tsub*&nrater*(&nrater-1)));
pvalue=1-probnorm(kappa/se);
run;
ods rtf file="&dloc&outfile" ;
options nonumber ;
goptions dev=png target=png htitle=1.5 htext=1.2 ftext=swissb ftitle=swissb
xmax=6 ymax=5 ;
proc print data=kap.koutput noobs label;
label kappa="Kappa" SE="Standard Error" pvalue="Prob>Z";
title "Kappa statistics for &Ncat-category ratings by &Nrater raters";
run;
ods rtf close;
%mend mkappa;


CONTACT INFORMATION
Your comments and questions are valued and encouraged. Please contact the authors at:
Bin Chen
WESTAT Inc.
5555 Ridge Ave
MS-R44
Cincinnati, OH 45213
Phone: 513-841-4316
Fax: 513-841-4470
Email: bchen2@cdc.gov

Text Utility Macros for Manipulating Lists of Variable Names

原文載點:http://www2.sas.com/proceedings/sugi30/029-30.pdf

SAS 裡面有相當多語法,後面是要接一連串的變數,如 KEEP, RENAME,ARRAY 或 PROC FREQ 程序中的 TABLES 指令等等。在遇到具有數十個甚至數百個變數的資料庫時,有極大的機會需要輸入大量變數在這些語法中。如果這些變數在一開始設定時已經有一些規則可循,比方說有共同的開頭字母(prefix),如 VAR001,VAR002,VAR003....,那就可以用一些連字符號 (-) 來節省輸入時間。但如果這些變數名稱沒有規則可循,這就會變成一個相當頭痛的問題。Robert J. Morris 在 2005 年的 SUGI 30 上發表了一篇專門處理龐大變數的文章,並寫了幾個text utility macro 供大家使用。

在本文中,只用一個含有六個變數的資料來當範例,但實際上可以想像中上百個甚至上千個變數。 Robert 用幾個例子來介紹 text utility macro 的用法。在這三個例子所出現的程式中有一個共同點,都是要先設定下面這行程式來表示:
%let orig_vars = a01a a01b d01 d02 t1_r t2_r;

EXAMPLE 1 – RECODING VARIABLES

假設這六個變數是一些數值變數,而我們想要製造另外六個新的變數,當舊變數是正時,新變數就是 1 ,反之則是 0。此外,新變數的名稱,是在舊變數名稱前加上「r_」。一般直覺上會寫出下面這樣的程式:

data newdata;
set mydata;
array orig[*] &orig_vars;
array recode[*] r_a01a r_a01b r_d01 r_d02 r_t1_r r_t2_r;
do i=1 to dim(orig);
if (orig[i] > 0) then recode[i] = 1; else recode[i] = 0;
end;
run;


這個程式有很大的缺點,那就是當你有上百個變數時,第二個陣列 recode 後面不就要一個一個變數加上 r_ 的開頭。這樣貼上的動作要重複上百次,非常消耗時間。如果使用 Robert 的 text utility macro,程式會變成:
data newdata;
set mydata;
array orig[*] &orig_vars;
array recode[*] %add_string(&orig_vars, r_, location=prefix);
do i=1 to dim(orig);
if (orig[i] > 0) then recode[i] = 1; else recode[i] = 0;
end;
run;

唯一和原始程式不同的地方只是把 recode 陣列後面所引入的參數換成 %add_string 這個 macro。這個 macro 會把 orig_vars 所預設的所有變數一次丟入,然後便會在每個變數前面加上 r_ 的檔頭。這樣 rename 的功能就完成了。

EXAMPLE 2 – PROC MEANS OUTPUT DATA SET

在使用 PROC MEANS 來計算每個變數的一些基本統計量時,有時候會需要將得到的結果另存成一個新的資料集,讓後續的程式能夠繼續使用。在另存的過程中,不免需要將每個變數所得到的各種統計量重新命名,於是含有大量變數的資料庫又會出現如例一的問題。當然,PROC MEANS 的 OUTPUT 語法中有一個叫做 autoname 的 option,可以自動幫你命名,可是當你不滿意 SAS 幫你命的名稱時,就是 %add_string 上場的時候了。

原始程式如下:
proc means data=mydata noprint;
output out=newdata
mean(&orig_vars) = a01a_mean a01b_mean d01_mean
d02_mean t1_r_mean t2_r_mean
max(&orig_vars) = a01a_max a01b_max d01_max
d02_max t1_r_max t2_r_max;
run;

上述程式只是將每個變數的平均數和最大值另存,這樣就要多輸入 2X6=12 個新名字。如果有一百個變數,要另存每個變數的平均數、標準差、最大值、最小值、中位數,那就等於要輸入 5X100=500 個新變數名稱,的確相當煩人。如果使用 %add_string,就會大幅縮減輸入時間:
proc means data=mydata noprint;
output out=newdata
mean(&orig_vars) = %add_string(&orig_vars, _mean)
max(&orig_vars) = %add_string(&orig_vars, _max);
run;

這個 macro 讓每個舊變數的平均數後面都加上 _mean,每個最大值後面都加上 _max。事實上,如果使用 autoname 就可以產生同樣的結果,但難保在實際應用中,不會出現其他想要設定的名稱格式。因此,%add_string 的確有他存在的必要。

EXAMPLE 3 – RENAMING VARIABLES TO HAVE A COMMON PREFIX

當要在 Data step 中大量重新命名變數時,所造成的困擾就相當清楚了。Robert 提供了另一個叫做 %rename_string 的 macro 來解決這個問題。原始程式如下:
data newdata;
set mydata;
rename a01a=r_a01a a01b=r_a01b d01=r_d01 d02=r_d02 t1_r=r_t1_r t2_r=r_t2_r;
run;

如果變數變成一百個,新的變數名稱是 r_ 加上舊的變數名稱,大概 ctrl+v 會按不完。用 %rename_string 輕鬆搞定:
data newdata;
set mydata;
rename %rename_string(&orig_vars, r_, location=prefix);
run;

原本可能要花十幾分鐘的重複貼上時間,現在只要二十秒就可完成批次改變數名稱的動作。

EXAMPLE 4 – PROC FREQ TABLE STATEMENT

在 PROC FREQ 程序中,如果要產生不同變數排列組合的次數分配表,那就是如下所示:
proc freq data=newdata;
tables a01a*r_a01a a01b*r_a01b d01*r_d01 d02*r_d02 t1_r*r_t1_r t2_r*r_t2_r;
run;

上述程式只是簡單的做出每個舊的變數和其對應的新變數(開頭加上 r_ 的那幾個)的次數分配表。我們可以一次使用 %parallel_join 和 %add_string 這兩個 macro 來完成這個動作:
proc freq data=newdata;
tables %parallel_join(
&orig_vars,
%add_string(&orig_vars, r_, location=prefix),
*
);
run;

%parallel_join 可以將 orig_vars 中每一個變數對應到相對位置的另一群變數,而這另一群變數就是由 %add_string 所製造出來(讓每個舊變數前面加上 r_),最後,用「*」把他們分別黏起來,就完成了原始程式碼中 tables 後面那一串文字的輸入了。

EXAMPLE 5 – RENAMING VARIABLES TO HAVE A COUNTER SUFFIX

最後一個例子算是終極大絕招。如果變數名稱太長,用 %rename_string 只會讓變數更長。理想的新變數是像 qvar1、qvar2、qvar3.... 這樣。這次 Robert 一次用三個 text utility macro 來完成這個動作。原始沒有用 macro 的程式碼如下:
data newdata;
set mydata;
rename a01a=qvar1 a01b=qvar2 d01=qvar3 d02=qvar4 t1_r=qvar5 t2_r=qvar6;
run;

新的程式碼看起來比舊的程式碼要來的長,不過一旦變數變成超級多時,新的程式碼看起來就相當短了:
data newdata;
set mydata;
rename %parallel_join(
&orig_vars,
%suffix_counter(qvar, %num_tokens(&orig_vars)),
=
);
run;

在新的程式碼中,一開始也是先引入 %parallel_join 把新舊變數用等號連結起來。在設定新變數名稱時,使用 %suffix_counter 來把新變數設定成 qvar 開頭,而 qvar 後面接著的數字,就由 %num_tokens 來計算 &orig_vars 所設定的每個舊變數的次序。第一個舊變數就會在 qvar 後面加上 1,第二個就加上 2,以此類推。如此一來又大功告成了!!!

原文中還有對每個 text utility macro 做更進一步的參數解說,可從原文載點下載來看。在此僅列出所有 macro 的程式碼,有興趣的人可以自由下載回去使用。

[num_tokens]
%macro num_tokens(words, delim=%str( ));
%local counter;

%* Loop through the words list, incrementing a counter for each word found. ;
%let counter = 1;
%do %while (%length(%scan(&words, &counter, &delim)) > 0);
%let counter = %eval(&counter + 1);
%end;

%* Our loop above pushes the counter past the number of words by 1. ;
%let counter = %eval(&counter - 1);

%* Output the count of the number of words. ;
&counter
%mend num_tokens;

[add_string]
%macro add_string(words, str, delim=%str( ), location=suffix);
%local outstr i word num_words;

%* Verify macro arguments. ;
%if (%length(&words) eq 0) %then %do;
%put ***ERROR(add_string): Required argument 'words' is missing.;
%goto exit;
%end;
%if (%length(&str) eq 0) %then %do;
%put ***ERROR(add_string): Required argument 'str' is missing.;
%goto exit;
%end;
%if (%upcase(&location) ne SUFFIX and %upcase(&location) ne PREFIX) %then %do;
%put ***ERROR(add_string): Optional argument 'location' must be;
%put *** set to SUFFIX or PREFIX.;
%goto exit;
%end;

%* Build the outstr by looping through the words list and adding the
* requested string onto each word. ;
%let outstr = ;
%let num_words = %num_tokens(&words, delim=&delim);
%do i=1 %to &num_words;
%let word = %scan(&words, &i, &delim);
%if (&i eq 1) %then %do;
%if (%upcase(&location) eq PREFIX) %then %do;
%let outstr = &str&word;
%end;
%else %do;
%let outstr = &word&str;
%end;
%end;
%else %do;
%if (%upcase(&location) eq PREFIX) %then %do;
%let outstr = &outstr&delim&str&word;
%end;
%else %do;
%let outstr = &outstr&delim&word&str;
%end;
%end;
%end;

%* Output the new list of words. ;
&outstr

%exit:
%mend add_string;

[rename_string]
%macro rename_string(words, str, delim=%str( ), location=suffix);

%* Verify macro arguments. ;
%if (%length(&words) eq 0) %then %do;
%put ***ERROR(rename_string): Required argument 'words' is missing.;
%goto exit;
%end;
%if (%length(&str) eq 0) %then %do;
%put ***ERROR(rename_string): Required argument 'str' is missing.;
%goto exit;
%end;
%if (%upcase(&location) ne SUFFIX and %upcase(&location) ne PREFIX) %then %do;
%put ***ERROR(rename_string): Optional argument 'location' must be;
%put *** set to SUFFIX or PREFIX.;
%goto exit;
%end;

%* Since rename_string is just a special case of parallel_join,
* simply pass the appropriate arguments on to that macro. ;
%parallel_join(
&words,
%add_string(&words, &str, delim=&delim, location=&location),
=,
delim1 = &delim,
delim2 = &delim
)

%exit:
%mend rename_string;

[suffix_counter]
%macro suffix_counter(base, end, start=1, zpad=0);
%local outstr i counter;

%* Verify macro arguments. ;
%if (%length(&base) eq 0) %then %do;
%put ***ERROR(suffix_counter): Required argument 'base' is missing.;
%goto exit;
%end;
%if (%length(&end) eq 0) %then %do;
%put ***ERROR(suffix_counter): Required argument 'end' is missing.;
%goto exit;
%end;
%if (&end < &start) %then %do; %put ***ERROR(suffix_counter): The 'end' argument must not be less; %put *** than the 'start' argument.; %goto exit; %end; %* Construct the outstr by looping from &start to &end, adding the counter * value to &base in each iteration. To handle the zero-padding, use the * putn function to format the counter variable with the Z. format. ; %let outstr=; %do i=&start %to &end; %if (&zpad > 0) %then %do;
%let counter = %sysfunc(putn(&i, z&zpad..));
%end;
%else %do;
%let counter = &i;
%end;
%let outstr=&outstr &base&counter;
%end;

%* Output the new list. ;
&outstr

%exit:
%mend suffix_counter;

[parallel_join]
%macro parallel_join(words1, words2, joinstr, delim1=%str( ), delim2=%str( ));
%local i num_words1 num_words2 word outstr;

%* Verify macro arguments. ;
%if (%length(&words1) eq 0) %then %do;
%put ***ERROR(parallel_join): Required argument 'words1' is missing.;
%goto exit;
%end;
%if (%length(&words2) eq 0) %then %do;
%put ***ERROR(parallel_join): Required argument 'words2' is missing.;
%goto exit;
%end;
%if (%length(&joinstr) eq 0) %then %do;
%put ***ERROR(parallel_join): Required argument 'joinstr' is missing.;
%goto exit;
%end;

%* Find the number of words in each list. ;
%let num_words1 = %num_tokens(&words1, delim=&delim1);
%let num_words2 = %num_tokens(&words2, delim=&delim2);

%* Check the number of words. ;
%if (&num_words1 ne &num_words2) %then %do;
%put ***ERROR(parallel_join): The number of words in 'words1' and;
%put *** 'words2' must be equal.;
%goto exit;
%end;

%* Build the outstr by looping through the corresponding words and joining
* them by the joinstr. ;
%let outstr=;
%do i = 1 %to &num_words1;
%let word = %scan(&words1, &i, &delim1);
%let outstr = &outstr &word&joinstr%scan(&words2, &i, &delim2);
%end;

%* Output the list of joined words. ;
&outstr

%exit:
%mend parallel_join;


CONTACT INFORMATION
The author welcomes any comments, questions, or suggestions for improvements. Contact the author at:
Robert J. (Joey) Morris
RTI International
3040 Cornwallis Rd.
Research Triangle Park, NC 27709
Email: rjmorris@rti.org

2007年3月6日 星期二

Identifying Plant Species: A Botanical Analysis Using PROC DISCRIM

原文載點:http://www2.sas.com/proceedings/sugi27/p199-27.pdf

鑑別分析是多變量分析中一個相當有趣的統計分析,他決定哪些變數可將一個群體完全切割成一些小群體,並將那些變數線性化,以提供後續的資料進行分類的動作。Robert G. Downer 和 Philip E. Hyatt 利用他做一些植物鑑別的分析,並將結果發表在 2002 年的 SUGI 27 上。兩人並沒有提出什麼特別的程式,但這的確是個不錯的範例介紹,可供初學者做簡易的入門。

話說 C. texensis 和 C. retroflexa 這兩類植物在以往相當難分別,於是兩人收集了六十株這類植物的一些資料,包含 perigynium length, spongy layer length, width of the widest leaf, length of the longest bract 和 height of the tallest culm 來讓鑑別分析做出分類的準則。

文內有一小段理論說明,但還是請有興趣的人直接去看教科書比較好。

在 SAS 中,鑑別分析的語法如下:

PROC DISCRIM DATA = calibration CROSSVALIDATE OUTCROSS = resultset ;
CLASS group;
BY ….;
PRIORS….;
VAR ; ….
RUN;


CORSSVALIDATE 可以提供鑑別分析裡面所能算出的估計值和誤差。分析結果可用 OUTCROSS 儲存起來。 CLASS 則是宣告哪些變數是離散型的。BY 則提供不同組採個別的分析處理。PRIORS 可以提前給定一些東西,比方說強迫某些觀測值進入某個群體。VAR 就放入提供鑑別分析的變數,但如果不寫 VAR 的話,SAS 會將全部的變數都丟進去作分析。

以下她們提供兩個分析,一個是指定某些變數拿來做鑑別的準則,另一個是將全部的變數拿來做鑑別的準則。

[程式一:放入全部變數]
proc discrim data = avgdat crossvalidate outcross = avgcvout;
class specid;
run;


[結果]
sugi199_27_01

[程式二:放入 avgpwid 變數]

proc discrim data = avgdat crossvalidate outcross = avgcvout2;
class specid;
var avgpwid;
run;


[結果]
sugi199_27_02

其實上面的結果不是太重要,只是在講經過鑑別分析後歸類正確的比例。當然正確比例越高表示結果越好,這是無庸置疑的。另一個更重要的是如何用圖像來清楚表示兩群是分開來的。Robert G. Downer 和 Philip E. Hyatt 只有秀出第二個程式所跑出來的鑑別分析圖。從[圖一]可以清楚看見被分開的兩群。這邊有兩個問題,一是兩人沒有提供做圖的程式。由於這篇文章是在 2002 年發表,而當時的 ODS/GRAPHICS 可能還不能做出這種圖,但我確信目前的 ODS/GRAPHICS 已經可以輕鬆地產生這個圖檔。有興趣的人可以往前找那一篇關於 ODS/GRAPHICS 的文章。另一個問題是這張圖並沒有標示分隔線。我記得以前用 Statistica 這套統計軟體做過鑑別分析,他會將鑑別線畫出來,大概就像[圖二]那樣。不過這是我自己模擬加上去的,不代表真正的鑑別線就是如此。總而言之,那條線的線性方程式可以提供後續資料來鑑別種類,相當實用。

[圖一]
sugi199_27_03

[圖二]
sugi199_27_04

CONTACT INFORMATION
Contact the authors at:
Robert G. Downer
Dept. Experimental Statistics
161 Agricultural Administration Building
Baton Rouge, LA 70803-5606
(225) 578-8373
rdowner@lsu.edu
Philip E. Hyatt
U.S. Forest Service
2500 Shreveport Hwy
Pineville, LA 71360
(318) 473-7262
phyatt@fs.fed.us

2007年3月4日 星期日

The Power of CALL SYMPUT – DATA Step Interface by Examples

原文載點:http://www2.sas.com/proceedings/sugi29/052-29.pdf

CALL SYMPUT 是一個使用在 Data step 的指令。其主要功能是提供一個良好的橋樑讓變數可以自由使用在 Data step 和 Macro 之間。Yunchao Tian 在 2004 年的 SUGI 29 上利用三個例子簡單地來介紹這個語法的強大功能。

EXAMPLE 1: CREATE A SERIES OF VARIABLE NAMES FROM ANOTHER VARIABLE'S VALUES

在模式配適中,我們經常需要將一些離散變數轉成 dummy variable。但有時候會遇到離散變數內的分項太多,這樣要製造 dummy variable 便會花上很多時間。舉例來說,有個離散變數 Con 含有 506 層,傳統方法要轉成 dummy variable 要用到以下程式碼:

IF CON = 1 THEN CON1 = 1; ELSE CON1 = 0;
IF CON = 2 THEN CON2 = 1; ELSE CON2 = 0;
. . . . . .
IF CON = 506 THEN CON506 = 1; ELSE CON506 = 0;


總共 506 行,在撰寫上極為消耗時間。而 SYMPUT 可以大幅簡化上述的程式碼,自動生成所有的 dummy variable。

以下有個簡短的例子。假設有個資料夾裡面有一個變數,12 個觀測值。如下所示:
DATA TESTDATA;
INPUT CON @@;
CARDS;
1 7 34 115 7 1 487 34 506 57 7 43
;
RUN;


欲新增數個變數,其名稱是 CON 後面接上觀測值的數字(CON1, CON7, ..., CON43)。當新的變數等於 CON 後面的數字時,其值等於 1,反之為 0。要完成這個目的,首先先將資料排序,並刪除重複數據(使用 NODUPKEY),然後把結果存成新資料集 UNIQUE:
PROC SORT DATA=TESTDATA OUT=UNIQUE NODUPKEY;
BY CON;
RUN;


建立一個新的暫存資料庫 _NULL_,將 UNIQUE 資料集導入,並將 CON 的最大值(也是最後一個數值)設定成一個 macro 變數,名為 N。
DATA _NULL_;
SET UNIQUE END=LAST;
IF LAST THEN CALL SYMPUT('N', PUT(CON, 3.));
RUN;


再設一個新的暫存資料庫,並將 0 設定給所有的 macro 變數 M1~M8。
DATA _NULL_;
DO I = 1 TO &N;
CALL SYMPUT('M'||LEFT(PUT(I, 3.)), '0');
END;
RUN;


再設一個新的暫存資料庫,讓 UNIQUE 資料庫裡面的 CON 值設定到對應的 macro 變數。
DATA _NULL_;
SET UNIQUE;
CALL SYMPUT('M'||LEFT(PUT(CON, 3.)), PUT(CON, 3.));
RUN;


接下來用一個 macro 程式來生成 dummy variable。
%MACRO GETCON;
%DO I = 1 %TO &N;
%IF &M&I = 0 %THEN %GOTO OUT;
IF CON = &M&I THEN CON&I = 1;
ELSE CON&I = 0;
%OUT: %END;
%MEND GETCON;


在原始的資料裡面導入剛剛做好的 macro 程式,這樣就大功告成了。
DATA TESTDATA;
SET TESTDATA;
%GETCON
RUN;


把資料列印出來看看。
PROC PRINT DATA=TESTDATA;
TITLE 'Table 1. List of CON with dummy variables';
RUN;


sugi052_29_01

EXAMPLE 2: GENERATE LABELS FOR A SERIES OF VARIABLES USING EXISTING FORMATS

這個例子可以讓人輕鬆的一次貼上大量的變數標籤。作者僅用一個含有六個變數的例子簡化程式的說明。

有一個資料內有六個變數(FLAG1~FLAG6)。
DATA FLAGS;
FLAG1 = 1; FLAG2 = 1; FLAG3 = 0;
FLAG4 = 1; FLAG5 = 0; FLAG6 = 1;
RUN;


先把要用到的標籤用 PROC FORMAT 設定好。
PROC FORMAT;
VALUE FMTFLAG
1='Red'
2='Purple'
3='Blue'
4='Yellow'
5='Orange'
6='Green';
RUN;


設定一永久變數 N 等於六,再製造一個新的暫存資料集,讓 SYMPUT 將剛剛 PROC FORMAT 弄好的標籤依序分配到新的 macro 變數 FMT1~FMT6。
%LET N = 6;
DATA _NULL_;
DO I = 1 TO &N;
CALL SYMPUT('FMT'||LEFT(PUT(I, 3.)), PUT(I, FMTFLAG.));
END;
RUN;


接著寫一小段 macro 程式,讓剛剛生好的 macro 變數再依序設定到原始變數 FLAG1~FLAG6。
%MACRO LABELS;
%DO I = 1 %TO &N;
FLAG&I = "I"
%END;
%MEND LABELS;


最後,再原始資料集裡面執行上一段 macro 程式,就大功告成了。
DATA FLAGS;
SET FLAGS;
LABEL %LABELS;
RUN;


將資料內容列印出來看看。
PROC CONTENTS DATA=FLAGS;
TITLE 'Table 2. Contents of SAS data set FLAGS showing variable labels';
RUN;


sugi052_29_02

最後一個例子有點是在講如何做批次改檔名,但我覺得沒有相當實用,所以不多加說明。有興趣的人去看本文。

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author at:
Yunchao (Susan) Tian
Social & Scientific Systems, Inc.
8757 Georgia Avenue, 12th Floor
Silver Spring, MD 20910
Work Phone: (301) 628-3285
Fax: (301) 628-3201
Email: Stian@s-3.com

An Alternative to PROC MI for Large Samples

原文載點:http://www2.sas.com/proceedings/sugi28/271-28.pdf

許多統計方法在進行時,多半會把 missing data 完全移除。好處是處理方便,壞處是,有時missing data 仍有一些有用但未知的資訊,我們永遠也不知道刪除他們會產生什麼重大的影響。另一方面,有時 missing data 太多,導致樣本變很小,power 一整個降低。在 SAS 中,PROC MI 和 PROC MIANALYZE 這兩個程序可以用一些內插法來補 missing data,以降低刪除 missing data 對資料本身的傷害性。不過 David Lanning 和 Doug Berry 利用更簡單的方法來做到內插的效果,並將結果發表於 2003 年的 SUGI 28 上。

簡單來說,David Lanning 和 Doug Berry 是用 PROC MEANS 程序來抓出每個變數的最小值和全距(同時間最大值也被決定了),並進行資料插補的動作,使的新增資料會在這個全距中隨機生成,並讓新的變數平均數和舊的變數平均數不會差異太多,減低填補的資料可能會造成資料巨大變異的可能性。

David 和 Doug 只有在文內發表他們所寫的程式,而沒有詳細解釋。我試著來解譯每個部分的功能。

資料如下:

data one;
input abbrid $ cnt amt_rnd amt_dec age @@;
datalines;
DTL 1 110 30.67 21 KRL 2 320 24.56 22
ARL 3 330 . 23 GLM . . . .
CRM 5 150 45.32 . MGM 6 . 28.12 26
PRG . 270 . 27 RDS . . 34.00 .
JKS 9 . . . DTL 2 160 54.98 20
KRL 4 180 67.54 22 ARL 6 140 34.13 24
GLM . . . . CRM 8 170 23.87 .
MGM 2 . 43.76 26 PRG . 130 12.54 28
RDS . . . . JKS 5 . . .
DTL 3 190 43.67 23 KRL 7 110 32.87 26
ARL 3 120 34.34 29 GLM . . . .
CRM 5 130 54.12 . MGM 2 . 76.43 21
PRG . 140 19.92 25 RDS . . . .
JKS 8 . 53.13 . JTS 6 . 73.93 30
;
run;


利用 PROC MEANS 來計算每個數值變數的最小值和全距,並存成新的資料集 three。

proc means data=one noprint;
var _numeric_;
output out=three (drop=_freq_ _type_) min= range= / autoname;
run;


製造一個暫時的資料集 _null_(不會顯示在library內),並把資料集 three 引入。裡面有個 symput 指令,可以將 data step 中的變數轉成 macro 可用的變數。我會在之後的文章引用另一篇 SUGI 的文章來討論 symput。但簡單來講,這個資料集將剛剛用 PROC MEANS 生出的資料集 three 中的所有結果存成新的變數以讓之後的程序可以直接用到。

data _null_;
set three;
array ass{*} _numeric_;
do i=1 to dim(ass);
call symput(vname(ass(i)),ass(i));
end;
drop i;
run;


新增一個資料集 four,並引入舊資料集 one。這一個步驟最主要就是在計算插補到 missing data 的數據。增補的根據是用 Uniform 分配生出亂數,並乘上之前算出來的全距,最後再加上最小值(粗體字部分)。

data four;
set one;
array wit{*} _numeric_;
do i=1 to dim(wit);
if wit(i)=. then
wit(i)=round((ranuni(0)*(symget(vname(wit(i))||'_range')))+(symget(vname(wit(i))||'_min')));
end;
drop i;
run;


現在我們來比較新舊資料集的結果。

[舊資料]
proc means data=one min max range mean;
title 'Original dataset with missing values';
run;




[新資料]
proc means data=four min max range mean;
title 'New dataset using random generated values for missing observations';
run;




由上面兩張圖表可知,新舊資料集的 Min, Max 和 Range 完全一樣,只有 Mean 的部分有些微不同。這種方法並沒有牽涉到太複雜的統計模擬運算,純粹是要解決 PROC MI 在插補龐大資料時所可能會消耗的冗長時間。由於兩人是在知名的保險公司 State Farm 工作,可想而知龐大的客戶群資料也許沒有辦法讓他們慢慢等待 PROC MI 模擬生成的效果。

CONTACT INFORMATION
Your comments and questions are valued and encouraged.
Contact the author at:
David Lanning
State Farm Insurance Companies
One State Farm Plaza – SC-3
Bloomington, IL 61710-001
Work Phone: (309)735-2723
Email: david.lanning.lyus@statefarm.com

ODS Sttisticl Graphics for Clinical Research

原文載點:http://www2.sas.com/proceedings/sugi31/095-31.pdf

Wei Cheng 於 2006 年的 SUGI 31 發表了一篇關於 ODS (Output Delivery System)繪製高解析度圖形的報告,相當具有建設性!!SAS 在以前最被人詬病的就是複雜的做圖程式還有與所花時間不成比例的低劣品質。這一直到現在還是一些人拒絕使用 SAS 的理由,尤其當不需要使用複雜資料整理和模式配適,而只要簡單生一個 box plot 或 histogram 的時候,讓許多人轉而用 SPSS 或其他統套軟體。現在經由 ODS 就可以在彈指之間,一口氣地畫出很複雜的圖形。現在,就來介紹 ODS/GRAPHICS 最基本的語法:

ODS GRAPHICS ON [/ MAGEFMT = image-file-type | STATIC | STATICMAP
IMAGENAME = filename RESET ];

procedures or data steps
ODS GRAPHICS OFF;


一開始,需用 ODS GRAPHICS ON 來啟動 ODS/GRAPHICS,後面可以加上一些圖檔格式和名稱的設定,但完全不加也沒有關係。中間就可以放任何的 Data procedure 和 Proc procedure,最後再用 ODS GRAPHICS OFF 來關閉 ODS。

以下就來介紹一些例子。

Regression

ods html;
ods graphics on / imagename = 'regression';

proc reg data = angina;
model y_impr = x_dur;
run;quit;
ods graphics off;
ods html close;


上述只是一個很簡單的迴歸程式,但只要在 PROC REG 前後粗體黑字部分的指令,就可以一口氣生出下列圖形:

sugi095_31_01sugi095_31_02sugi095_31_03

如果想要將這些圖分別儲存的話,只要在 PROC REG 後面加上 plots(unpack) 即可。

ods html;
ods graphics on;
proc reg data = angina plots(unpack);
model y_impr = x_dur;
run; quit;
ods graphics off;
ods html close;


附帶一提的是,所有的輸出報表和圖形,因為 ODS HTML 被呼叫的關係,都會被儲存為 html 檔。

另外,可以用 ods select 的指令來指定只生成某一圖檔。以上述圖形而言,想只生成最左邊那張圖,可用 ods select fit 來限定:

ods html;
ods graphics on;
ods select fit;
proc reg data = angina;
model y_impr = x_dur;
run;
quit;
ods graphics off;
ods html close;


如果描繪的點太多,沒有辦法把每一點的 ID 都標示在圖片上,這時可以用 imagefmt = staticmap 的指令,讓觀測值相關訊息可透過滑鼠來顯示。程式如下:

ods html;
ods graphics on / imagefmt = staticmap;
ods select fit;
proc reg data = angina;
model y_impr = x_dur;
run;
quit;
ods graphics off;
ods html close;


從下圖可知,當滑鼠移到某一點時,就會彈出新的訊息框顯示相關訊息。

sugi095_31_04

總而言之,在 PROC REG 下可以用 ODS/GRAPHICS 產生八種圖形:

• Residuals versus the predicted values
• Studentized residuals versus the predicted values
• Studentized residuals versus the leverage
• Normal quantile plot of the residuals
• Dependent variable values versus the predicted values
• Cook's D versus observation number
• Histogram of the residuals
• A "Residual-Fit" (or RF) plot consisting of side-by-side quantile plots of the centered fit and the residuals.

GLM

GLM 的分析中並沒有用到太多的圖形,但利用 ODS/GRAPHICS 仍可在 PROC GLM 中生出簡單的 box plot。同樣也是只要呼叫 ods graphics 即可。

ods html;
ods graphics on;
proc glm data = angina;
class x_dur;
model y_impr = x_dur;
run; quit;
ods graphics off;
ods html close;


sugi095_31_05

@ANCOVA @

共變數分析同樣利用 PROC GLM 完成,不同的地方在於模式裡面有個連續變數。因此當 SAS 發現有連續變數放入 PROC GLM 時,就會啟動共變數分析。此時若同時啟動 ODS/GRAPHICS 系統,則會產生 covariance plot。程式和圖形如下:

ods html;
ods graphics on / imagename = 'ancova';
proc glm data = tri;
class trt;
model trichg = trt hgba1c / solution;
run; quit;
ods graphics off;
ods html close;


sugi095_31_06

@ Log Rank Test @

在倖存分析中,可用 PROC LIFETEST 來進行 Log Rank 檢定。同樣地,也可利用 ODS/GRAPHICS 將最後的倖存機率畫出來。

ods html;
ods graphics on / imagename = 'lifetest';
proc lifetest data = hsv;
time wks * cens (1);
strata vac;
run;
ods graphics off;
ods html close;


sugi095_31_07

如果想看各細部的圖形,可利用下面的程式來增生:

ods html;
ods graphics on;
proc lifetest data = hsv;
time wks * cens (1);
strata vac;
survival plots = (survival,density,epb,hazard,loglogs,logsurv,hwb,cl,stratum);
run;
ods graphics off;
ods html close;


sugi095_31_08sugi095_31_09
sugi095_31_10sugi095_31_11
sugi095_31_12sugi095_31_13
sugi095_31_14sugi095_31_15

以上所有的 ODS/GRAPHICS 範例,都是用預設的設定來繪圖,而這些設定都放在一個叫做 Stat.Reg.Graphics 的模版裡面。當然,這個模版也是可以做調整的,只要利用 PROC TEMPLATE 的程序,就可以做進一步的更動。PROC TEMPLATE 的基本語法如下:

PROC TEMPLATE;
DEFINE STATGRAPH name-of-graph-definition;
[declaration-statements;]
[layout-statements;]
[plot-statements;]
[text-statements;]
END;
RUN;


要嵌入自訂的模版,可在 Data procedure 中加上下面黑色粗體字的那兩行。

ODS GRAPHICS ON;
DATA _NULL_;
SET my-data;
FILE PRINT
ODS = (TEMPLATE = “my-graphics-template”);
PUT _ODS_;
RUN;


文內並沒有詳細說明 PROC TEMPLATE 的所有語法,僅列了兩個例子,因此我只寫其中一個範例在這,詳細情況請參見原文。

@ TWO-SAMPLE T-TEST @

程式:

proc format;
value $trt 'A' = 'Active' 'P' = 'Placebo';
run;
data fev;
length trtgrp $ 7;
input patno trt $ fev0 fev6 @@;
chg = fev6 - fev0; if chg = . then delete; trtgrp = put(trt, $trt.);
datalines;
101 A 1.35 . 103 A 3.22 3.55 106 A 2.78 3.15
……
;
run;
proc template;
define statgraph mygraphs.meanchg;
layout gridded;
entrytitle 'Bar Chart of Mean Change by Treatment Group' ;
entrytitle 'with Upper and Lower CLM' ;
barchartparm x = trtgrp y = chg_mean / yerrorupper=uclm yerrorlower=lclm;
endlayout;
end;
run;
ods html;
ods graphics on / imagename='meanchange';
proc summary data = fev nway;
class trtgrp;
var chg;
output out = meanchg (drop = _type_ _freq_) mean = chg_mean lclm = lclm uclm = uclm;
run;
data _null_;
set meanchg;
label chg_mean = "Mean Change" trtgrp = "Treatment Group";
uclm = uclm - chg_mean; lclm = chg_mean - lclm;
file print ods = (template='mygraphs.meanchg');
put _ods_;
run;
ods graphics off;
ods html close;


圖形:
sugi095_31_16

CONTACT INFORMATION
I welcome and appreciate your comments and questions. Contact the author at:
Wei Cheng,
Isis Pharmaceuticals, Inc.,
1896 Rutherford Rd., Carlsbad, CA 92008
(760) 603-3807
Email: wcheng@isisph.com

Identifying Continuity in Longitudinal Data

原文載點:http://www2.sas.com/proceedings/sugi28/089-28.pdf

這一篇由 Merle E. Hamburger 和 Thomas Sukalac 二人一同於 2003 年的 SUGI 28 所發表的文章,主要是針對在分析長期趨勢資料時所經常會碰到的一個問題,即觀測值是否在實驗期間不間斷地都有觀測到數據或是長期間連續地暴露在某一種狀態下。也就是說,我們所關心的,是倒底這段不中斷的時間到底持續了多久。兩人寫了一段 SAS 程式,來說明如何從原始資料得到這種數據。

兩人使用一個有關 HIV 感染者的世代研究資料來當作範例。此資料包含下列變數:

HRNID = study ID
VLOAD = HIV viral load
VBDATE = date of viral load assessment

首先,利用 PROC SORT 將資料 依照 HRHID 變數排序。

proc sort data=VLDATA ;
by hrnid ;
run ;


接下來就是使用一個 Data procedure 把每個病人從一開始到 viral load 成功的時間算出來。我將這個程式拆成數段分別解釋。

data INTERVAL ;
set VLDATA ;
by hrnid ;
retain inint begdate enddate ;
keep hrnid begdate enddate interval undetect ;


建立一個新的資料集名為 INTERVAL,除了引用原始資料 VLDATA 外,retain 可以使變數在經過任何迴圈或條件式之後,在下一次進入迴圈或條件式前仍能繼續保有數值。若沒有 retain 的話,所有的變數在一進入 Data procedure 時都會被視為 missing。By hrnid 是希望舊的資料能夠依照 hrnid 的順序讀進,並順便讓 SAS 產生兩個隱藏的變數:first.hrnid 和 last.hrnid。當該筆資料是某個病人的第一個觀測值,則其 first.hrnid = 1,反之為 0。 同道理,若該筆資料是某個病人的最後一個觀測值,則 last.hrnid = 1,反之為 0。這個程序會生出一些新的變數,如 begdate, enddate, interval (最重要的)和 undetect。連同 hrnid ,將會被 keep 指令給繼續留在新的資料集當中。

接著,新增一個 inint 的變數,且讓他在 first.hrnid =1 時為零。

if first.hrnid then inint = 0 ;

再來,設定兩個條件子 inint = 0 和 vload = 1,產生一個新變數 begdate,使的滿足這兩個條件的病人的 begdate 值等於 vbdate,並讓其 inint 值改變為 1。

if (inint = 0) and (vload = 1) then do ;
begdate = vbdate ;
inint = 1 ;
end ;


完成上述動作後,新增一個 enddate 變數,使其在滿足 inint = 1 時,其值會等於 vbdate。

if inint = 1 then enddate = vbdate ;

接著,當滿足兩個條件子 inint = 1 和 vload > 1 時,先讓 inint 歸零,再來就是要算時間區間 interval 了(除以 365.25 以轉換其單位成年)。但若 interval 大於二,則新增一個變數 undetect 並使其為 1。

if inint = 1 and vload > 1 then do ;
inint = 0;
interval = (enddate - begdate) / 365.25 ;
if interval ge 2 then do ;
undetect = 1 ;
output;
end ;
end ;


下面的程式功能和上面那個一樣,唯一不同的是考量當 last.hrnid = 1 時該怎樣去算 interval。

if last.hrnid and inint = 1 then do ;
interval = (enddate - begdate) / 365.25 ;
if interval ge 2 then do ;
undetect = 1 ;
output;
end ;
end ;


最後利用 format 語法將算出來的 begdate 和 enddate 轉換成「月/日/年」的時間型態,而 interval 則轉成 4.2 格式型態。

format begdate mmddyy10. enddate mmddyy10. interval 4.2 ;

如此一來就大功告成了!!

完整程式:

data INTERVAL ;
set VLDATA ;
by hrnid ;
retain inint begdate enddate ;
keep hrnid begdate enddate interval undetect ;
if first.hrnid then inint = 0 ;
if (inint = 0) and (vload = 1) then do ;
begdate = vbdate ;
inint = 1 ;
end ;
if inint = 1 then enddate = vbdate ;
if inint = 1 and vload > 1 then do ;
inint = 0;
interval = (enddate - begdate) / 365.25 ;
if interval ge 2 then do ;
undetect = 1 ;
output;
end ;
end ;
if last.hrnid and inint = 1 then do ;
interval = (enddate - begdate) / 365.25 ;
if interval ge 2 then do ;
undetect = 1 ;
output;
end ;
end ;
format begdate mmddyy10. enddate mmddyy10. interval 4.2 ;
run ;


CONTACT INFORMATION
Your comments and questions are valued and encouraged.
Contact the author at:
Merle E. Hamburger
CDC-P
1600 Clifton Rd., MS E-45
Atlanta, GA 30333
(404) 639-5263
(404) 639-6127 (fax)
mhamburger@cdc.gov
CODE { display: block; /* fixes a strange ie margin bug */ font-family: Courier New; font-size: 8pt; overflow:auto; background: #f0f0f0 url(http://klcintw.images.googlepages.com/Code_BG.gif) left top repeat-y; border: 1px solid #ccc; padding: 10px 10px 10px 21px; max-height:200px; height:200px; // for IE6 line-height: 1.2em; }