SUGI CLUB — SAS User Group International Club

問題討論區

2108-09-30T12:21:00.011-04:00

[版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規]
1.有要問語法問題的請在此留言。
2.即日起不回覆私底下email給我的問題，請善加利用此問題討論區。
3.不保證能夠回答所有問題，若不見我回覆則盼熱心網友能協助解答。
4.請勿貼出一大串程式碼說跑不出來要我幫忙debug。
5.SAS軟體本身的問題請自行聯繫SAS客服部門。
6.不幫忙代寫程式。
[版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規][版規]

Fitting Cox Model Using PROC PHREG and Beyond in SAS

2011-10-14T18:52:00.001-04:00

http://support.sas.com/resources/papers/proceedings09/236-2009.pdf

Cox PH 模型在倖存分析（survival analysis）被廣泛的使用。一篇發表在SAS GLOBAL FORUM 2009的技術文件解說了在 SAS 底下如何用 Cox PH 模型來計算一些重要估計量的方法，並且提供一個方便的巨集程式來簡化繁複的程式寫作和運行。

C-index
C-index原本是在ROC curve裡面用來解釋預測鑒別能力的測量指標，但後來在Cox PH model裡也可以應用。在倖存分析裡，C-index可以解釋為一個機率，是從某個事件(如：死亡)來的一個病人，他會比非該事件內的任何一個病人都擁有較高的機會會發生該事件。以下列程式來說明如何用 SAS 算出 C-index：

proc phreg data=sample;
  id idn;
  model combdays*combfv(0)=mihx diabhx lowef;
  output out=obs survival=surv;
run;

這個PROC PHREG程序並沒有什麼特別，使用者唯一需要額外做的就是利用一個ODS OUTPUT指令把個體的ID, 觀測到的倖存時間以及倖存函數估計值另存一個新檔(ods)。接著，用下面兩到程序去整理出一個concordance data。

proc sql;
      create table allset as
      select idn_j, y_j, x_j, idn as idn_i, surv as y_i, combdays as x_i
      from evtset, obs
      where idn_j<>idn;
quit;
data concord;
      set allset;
      if (x_iy_j) or (x_i>x_j and y_i

然後用下面這一串程式直接把C-index及其95%信賴區間該算出來：

data _null_;
    set concord end=eof;
    retain nch ndh;
    if _N_=1 then do;
      nch=0;
      ndh=0;
      end;
    if concord=1 then nch+1;
    if concord=0 then ndh+1;
    if eof=1 then do;
      call symput('ch',trim(left(nch)));
      call symput('dh',trim(left(ndh)));
      call symput('uspairs',trim(left(_n_)));
      end;
run;
data _null_;
     set sample end=eof;
     if eof=1 then call symput('totobs',trim(left(_n_)));
run;
%put &ch &dh &uspairs &totobs;
data calculat;
    ch=input("&ch",12.0);
    dh=input("&dh",12.0);
    uspairs=input("&uspairs",12.0);
    totobs=input("&totobs",10.0);
    pc=ch/(totobs*(totobs-1));
    pd=dh/(totobs*(totobs-1));
    c_hat=pc/(pc+pd);
    w=(2*1.96**2)/(totobs*(pc+pd));
    low_ci_w=((w+2*c_hat)/(2*(1+w)))-(sqrt((w**2+4*w*c_hat*(1-c_hat))/(2*(1+w))));
    upper_ci_w=((w+2*c_hat)/(2*(1+w)))+(sqrt((w**2+4*w*c_hat*(1-c_hat))/(2*(1+w))));
run;

結果如下圖所示：

Calculating adjusted survival rates using corrected group prognosis method
作者自製了一段巨集程式 %CGPM_adj 可將倖存機率以及倖存曲線圖直接輸出。其原始碼如下：

*****************************************************************************;
* Program: CGPM_adj.sas;
* Purpose: Calculate the unadjusted and adjusted survival rates and plot the
*survival curves using corrected group prognosis method;
*****************************************************************************;
*----------------------------------------------------------------------------;
* 1. All variables used in the model have to be binary variable.
*Categorical and numeric variables need to be re-coded as dummy
*variable to meet this requirement;
* 2. The maximum number of variables in the model is limited to 20;
*(not including the control variable). But this number can be
*easily expanded as needed.;
* 3. The analysis dataset for modeling was assign to be in SAS
*default library WORK. ;
*----------------------------------------------------------------------------;
%macro CGPM_adj(outputpath, figname, tbname, dsn, covars, ctrlvar, outcom, cnsrval=0, timevar, rlalpha=0.05);
*----------------------------------------------------------------------------;
* Calculate unadjusted survival rate;
*----------------------------------------------------------------------------;
*create dataset for controlled variable;
proc sql;
      create table ctrlset as
      select distinct &ctrlvar
      from &dsn
      order by &ctrlvar;
quit;
*get the unadjusted survival rates for controlled variable;
proc phreg data=&dsn noprint;
      model &timevar * &outcom (&cnsrval) = &ctrlvar;
      baseline covariates=ctrlset out=unadjset survival=survival / nomean;
run;

*----------------------------------------------------------------------------;
* Calculate adjusted survival rates using corrected group prognosis method;
*----------------------------------------------------------------------------;
*create identifier set for all possible combinations of covariates;
data _null_;
      cnt=1+count("&covars",' ');
      call symputx('varcnt', cnt);
run;
%put &varcnt;

%macro idx;
      %let i=1;
      %do %until(%scan(&covars,&i,' ')=%str());
      %local varnam;
            %let varnam=%scan(&covars,&i,' ');
            xp=&i+1;
            if &varnam=1 then substr(idx,xp,1)=1;
      %let i=%eval(&i+1);
      %end;
%mend;
%macro idxrv;
      %let i=1;
      %do %until(%scan(&covars,&i,' ')=%str());
      %local varnam;
            %let varnam=%scan(&covars,&i,' ');
            xp=&i+1;
            if substr(idx,xp,1)='1' then &varnam=1; else &varnam=0;
      %let i=%eval(&i+1);
      %end;
%mend;

data modelset;
      set &dsn;
      idx=left(put(10**(input("&varcnt",2.0)),21.));
      %idx;
      keep idx &ctrlvar &covars &timevar &outcom;
run;
proc freq data=modelset noprint;
      tables idx / out=idx;
run;
data popset (drop=percent xp);
      if _n_=1 then do until(last);
            set idx end=last;
            &ctrlvar=0;
            %idxrv;
            output;
            end;
      set idx;
            &ctrlvar=1;
            %idxrv;
            output;
run;
ods output censoredsummary=censum parameterEstimates=parasum;
proc phreg data=modelset;
      id idx;
      model &timevar * &outcom (&cnsrval) = &ctrlvar &covars / rl alpha=&rlalpha;
      baseline covariates=popset out=adjset0 survival=survival / nomean;
run;
ods output close;
proc sort data=popset out=count (keep=idx count);
by idx;
run;
proc sort data=adjset0;
by idx;
run;
data adjset1;
      merge adjset0 count;
      by idx;
run;
proc sort data=adjset1;
      by &ctrlvar &timevar;
run;
data adjset;
      set adjset1;
      by &ctrlvar &timevar;
            if first.&timevar then frqsum=0;
            frqsum+count;
            if first.&timevar then sursum=0;
            sursum+survival*count;
      if last.&timevar;
      adjsurv=sursum/frqsum;
      keep &ctrlvar &timevar adjsurv;
run;
data allsurv;
      length group $22;
      set adjset (in=a rename=(adjsurv=survival)) unadjset (in=b);
      if a=1 and &ctrlvar=0 then group='Adjusted'||' '||"&ctrlvar"||'=0';
     else if a=1 and &ctrlvar=1 then group='Adjusted'||' '||"&ctrlvar"||'=1';
      else if b=1 and &ctrlvar=0 then group='Unadjusted'||'
'||"&ctrlvar"||'=0';
      else if b=1 and &ctrlvar=1 then group='Unadjusted'||'
'||"&ctrlvar"||'=1';
      eventrate=1-survival;
run;
data summary;
      set allsurv;
      by group;
      if last.group;
run;
*----------------------------------------------------------------------------;
* You can change the code below to modify the title and graph axis scale etc.;
*----------------------------------------------------------------------------;
goption reset=all device=emf gsfname=grfa hsize=6.0 vsize=4 ftext='Arial/bo' htext=11pt display;
filename grfa "&outputpath.\&figname..emf";

symbol1 i=steplj v=none l=1 ci=green w=1;
symbol2 i=steplj v=none l=1 ci=black w=1;
symbol3 i=steplj v=none l=3 ci=green w=1;
symbol4 i=steplj v=none l=3 ci=black w=1;
axis1 offset=(0.5 cm,0.5 cm);

axis2 label=(angle=90 j=c 'Survival') order=(0.8 to 1 by 0.05) offset=(0.0cm,0.0 cm);
   legend1 across=2
   mode=share
   position=(bottom left inside)
   label=none
   value=(j=l h=10pt)
   offset=(0.5cm, 0cm);
proc gplot data=allsurv;
   plot survival * &timevar = group / legend=legend1 haxis=axis1 vaxis=axis2;
   label eventrate='Survival Rate' &timevar='Days After Randomization';
run;
quit;
title;
footnote;
ods rtf style=journal file="&outputpath.\&tbname..rtf" bodytitle
startpage=no;
proc print data=parasum noobs;
title "Analysis of Maximum Likelihood Estimates (alpha=&rlalpha for CI)";
run;
proc print data=summary noobs;
var group survival eventrate;
title 'Estimated survival and event rates at the time of last event observed';
run;
ods rtf close;
%mend;

巨集參數定義如下：

outputpath=報表與圖形輸出的路徑
figname=圖形檔名
tbname=報表檔名
dsn=存放在SAS裡面的資料檔名
covars=所有預測變數的名稱，不包含控制變數的名稱
ctrlvar=控制變數的名稱
outcom=輸出變數(事件)的名稱
cnsrval=設限的代碼，預設值為0
timevar=時間變數名稱
rlalpha=顯著水準alpha，預設值為0.05

承上面的範例，巨集程式執行的寫法如下：

%CGPM_adj(outputpath=L:\SGF2009\Examples, 
          figname=Adjdiab_fig, 
          tbname=adjdiab_summary, 
          dsn=sample, 
          covars=age65 lowef chfhx, 
          ctrlvar=diabhx, 
          outcom=deathfv, 
          cnsrval=0, 
          timevar=deathdays, 
          rlalpha=0.05);

結果如下：

Presenting the effect of a continuous covariate on estimated survival
當模型內包含連續型變數時，需要額外多兩個步驟來計算n-year倖存率。首先，用PROC PHREG程序把基底函數(baseline function)等一干相關的估計量另存成兩個新檔 yrout & sample(紅色部分)：

proc phreg data=sample;
model combdays*combfv(0)=age male;
baseline out=yrout lower=low_ci upper=up_ci survival=surv covariates=sample/alpha=0.05 nomean;
run;

再依照想要分類的變數(male)，去和倖存時間變數去排序(combdays)後再另存新檔(maxday)：

proc sort data=sample out=maxday (where=(combfv=1) keep=idn male combfv combdays);
      by male combdays;
run;

然後利用下面兩個資料步驟程序去計算男人和女人的5-year倖存率：

data _null_;
      set maxday;
      by male;
      if male=0 and last.male then call symput('lastday_f',combdays);
      if male=1 and last.male then call symput('lastday_m',combdays);
run;
%put &lastday_f &lastday_m;
data yrest;
      set yrout;
      evtrate=(1-surv)*100;
      evtlow=(1-low_ci)*100;
      evtup=(1-up_ci)*100;
      if combdays=&lastday_m and male=1 then output;
      if combdays=&lastday_f and male=0 then output;
run;

最後畫出的圖形如下所示：

(PS) 原文內並沒有附此圖的繪製程式。

Contact information
Lea Liu, Maryland Medical Research Institute
600 Wyndhurst Ave. Baltimore, MD 21210
(410) 435-4200, Email: lliu@mmri.org, Web: www.mmri.org

Wow! You Did That Map with SAS/GRAPH®?

2011-08-16T17:38:00.000-04:00

原文載點：http://www.nesug.info/Proceedings/nesug07/np/np01.pdf

身為一個專門從事空間時間序列分析的學者我來講，每天跟不同的地圖奮鬥是很稀鬆平常的是。目前市面上最熱門的地圖分析軟體莫過於ArcGIS，但這套軟體僅提供有限的分析方法，而其主要的功能僅僅就是繪製地圖而已。如果你用了別的軟體進行分析，就必須要將結果另存新檔，然後把結果跟地圖檔一起輸入到ArcGIS裡面才能開始繪圖。想當然爾，這裡面也沒有太多強大的資料整理功能，所以一切的一切都還是要經過另外的軟體將資料結果整理好才行。如果分析和繪製地圖的動作可以一次在SAS裡面完成的話將會省下許多時間和精力。NESUG 2007年有一篇教學文件，提供一些在SAS繪製地圖的資訊。

PROC MAPIMPORT

這是SAS專門用來輸入地圖資料的語法，不過必須先知道去哪裡可以找到地圖資料。以下是作者列出的幾個不錯的網頁資源：

(1) ESRI (a GIS software company)
http://arcdata.esri.com/data/tiger2000/tiger_download.cfm

(2) US Census Bureau
http://www.census.gov/geo/www/cob/index.html

(3) CDC (Centers for Disease Control)
http://www.cdc.gov/epiinfo/shape.htm

(4) CIP (International Potato Center, part of the Consultative Group on International Agricultural Research)
http://research.cip.cgiar.org/gis/index.php

其中，ESRI 就是 ArcGIS 的母公司，而美國地圖大都可以在 US Census Bureau 的網站上免費下載。CDC 的網站則是可以下載一個名叫Epi Info的軟體，裡面除了有地圖外還有其他的流病資料。CIP 的網站則是除了地圖以外還提供氣候資料。個人的經驗是 US Census Bureau 網站就足夠提供所有需要使用美國資料來進行研究的地圖。如果想進行台灣的研究，則可以到交通部的運研所下載：

http://www.iot.gov.tw/ct.asp?xItem=154948&ctNode=1091

如果要別的國家的地圖資料，有兩個網站可以下載：

http://www.diva-gis.org/gData
http://www.vdstech.com/map_data.htm

接下來的範例是用北卡（North Carolina）的資料畫地圖。先去 US Census Bureau 找北卡的 shape 檔：

North Carolina - zt37_d00_shp.zip (1,672,705 bytes)

解壓縮後存在自定目錄，然後執行下列程式：

/*Version 9.1 */
proc mapimport 
       datafile="path/zt37_d00.shp" 
       out=mymap91;
run;
/* Version 9.2 */
proc mapimport;
       datafile="path/zt37_d00.shp" 
       out=mymap92;
       id zcta;
run;

記得把上面的 path 改為自己的路徑。如此一來 SAS 就可以把 shapefile 裡面的資料都輸入進去。

MAPSONLINE

這是 SAS 官網提供的線上工具，裡面也有地圖資料和一些其他資訊。網址如下：

http://support.sas.com/rnd/datavisualization/mapsonline/html/home.html

重點是如果你不小心把 sashelp 資料庫裡面一個叫做 ZIPCODE 的檔案刪掉的話，可以去上面這個網址重新下載，因為這個資料裡面包含美國所有的地理資訊，要畫美國地圖一定得連結到這個檔案。裡面還有一個很威的資料名為 WORLDCTS，包含全世界各大城市的基本資料，可以到這個連結直接下載：

http://support.sas.com/rnd/datavisualization/mapsonline/html/sampledata.html

CHOOSING MAP COLORS

在繪製地圖時，定義顏色是相當重要的一件事情。如何讓地圖上的顏色不會眼花撩亂並且能夠使讀者一目了然，這已經無關統計分析能力，而是一門藝術了。SAS 裡面對於顏色的指定可以使用顏色本來的英文名稱、RGB值、HLS值以及HEX值。所有的定義都可以去 SAS TS-688 的網站去看：

http://support.sas.com/techsup/technote/ts688/ts688.html

下面用一個簡單範例介紹：

%let color1=lightgreen ; 
%let color2=cornsilk ;
proc template;
   define style styles.grade;
   parent=styles.listing;
   style twocolorramp / startcolor=&color1 endcolor=&color2;
   end; 
run;

DISTANCE CALCULATIONS

計算點和點之間的距離是空間分析裡面的基本功，當知道兩個地點的經緯度（不知道的可以去用google map查詢）後，利用下面這個巨集程式便可輕鬆計算出來：

%macro geodist(lat1,long1,lat2,long2);
%let pi180=0.0174532925199433;
    7921.6623*arsin(sqrt((sin((&pi180*&lat2-&pi180*&lat1)/2))**2+
    cos(&pi180*&lat1)*cos(&pi180*&lat2)*(sin((&pi180*&long2-    
    &pi180*&long1)/2))**2));
%mend

如果是做美國的地理研究，則只需要知道每個 county 的 ZIP code 即可。在SAS V9.2版裡面有一個 ZIPCITYDIST() 函式就是專門用來計算兩個 ZIP code 之間的距離。範例如下：

%let zip1=27513;
%let zip2=21202;
/* latitude and longitude coordinates */
proc sql;
create table zip1 as select * from sashelp.zipcode where zip=&zip1;
create table zip2 as select * from sashelp.zipcode where zip=&zip2;
select x into :long1 from zip1;
select y into :lat1 from zip1;
select x into :long2 from zip2;
select y into :lat2 from zip2;
quit; run;
%let city1=%sysfunc(zipcity(&zip1));
%let city2=%sysfunc(zipcity(&zip2));
/* two new 9.2 functions */
%let zipdist=%sysfunc(zipcitydistance(&zip1,&zip2));
%let geodist=%sysfunc(geodist(&lat1, &long1, &lat2, &long2, 'M'),comma10.5);

GRAPHIC OVERLAYS

上面這張圖是利用 PROC GCHART 畫好柱狀圖，然後用 PROC GMAP 畫好地圖，之後把柱狀圖疊在地圖上面。首先把地圖畫好：

pattern v=me c=gray98;
proc gmap data=maps.us (obs=1) map=maps.us all anno=region_outline;
      id state;
      choro state / discrete nolegend coutline=gray98;
run;
quit;

再來把柱狀圖下面的文字設定好：

axis1 label=('NEW ENGLAND'        j=c '%CHANGE 48.7') ;
axis2 label=('MIDDLE ATLANTIC'    j=c '%CHANGE 21.1') ;
axis3 label=('SOUTH ATLANTIC'     j=c '%CHANGE 18.7') ;
axis4 label=('WEST NORTH CENTRAL' j=c '%CHANGE 63.3') ;
axis5 label=('EAST NORTH CENTRAL' j=c '%CHANGE 33.3') ;
axis6 label=('EAST SOUTH CENTRAL' j=c '%CHANGE 48.4') ;
axis7 label=('WEST SOUTH CENTRAL' j=c '%CHANGE 22.5') ;
axis8 label=('MOUNTAIN'           j=c '%CHANGE 50.9') ;
axis9 label=('PACIFIC'            j=c '%CHANGE 23.9') ;
axis10 order=0 to 20 by 5 label=none value=none style=0 major=none minor=none;

其中，axis1~axis9 自然就是九個柱狀圖 X 軸的文字和統計數據（這自己另外算好先），而 axis10 則表示用來控制 axis1~axis9 的 Y 軸共同設定。然後依序畫出每張柱狀圖，以右上角的 New England 為例：

proc gchart data=tpay;
goptions horigin=8.75 in vorigin=6.75 in;
where region eq 1;
vbar time / type=mean discrete sumvar=x mean raxis=axis10 maxis=axis1;
run;
quit;

資料集 tpay 裡面因為有所有 region 的資料，所以要用一個 where statement 限制只使用 New England (region = 1) 的數據。而 raxis 和 maxis 就是設定 Y 軸和 X 軸的選項，所以把 axis10 和 axis1 分別填入。當要畫別區的柱狀圖時，只要把 maxis 後面的設定換掉即可。至於 goptions 則是用來控制柱狀圖出現的位置。這就必須自己一直更換不同的數據來看位置對不對。你也可以用一個巨集程式來簡化步驟：

%macro bars(reg,h,v);
goptions horigin=&h in vorigin=&v in;
where region eq &reg;
vbar time / type=mean discrete sumvar=x raxis=axis10 maxis=axis&reg mean;
run;
%mend;

因此九張柱狀圖就可以用下面短短的程式碼一口氣完成：

proc gchart data=tpay;
%bars(1,8.75,6.45) %bars(2,7.75,4.90) %bars(3,7.50,3.00) 
%bars(4,4.00,5.00) %bars(5,6.00,4.25) %bars(6,6.10,2.50) 
%bars(7,4.00,1.85) %bars(8,1.95,4.20) %bars(9,0.00,2.00)
quit;

最後一個步驟就是要把這些柱狀圖疊到地圖上。我們先用 ODS 開一個空白的 pdf 檔：

ods pdf file='z:\map_bars.pdf' notoc startpage=never;

後面的 notoc 和 startpage 主要是抑制 pdf 檔的書籤出現以及讓所有圖形都出現在同一頁。然後進行 pdf 檔裡面的環境設定：

goptions reset=all ftext='Helvetica/bo' htext=3.5 gunit=pct ctext=blue lfactor=3;

這些都是來設定字形、字色還有字體大小。這樣就大功告成了。整個流程可以用下列四步驟來簡單描述：

步驟一：用 GOPTIONS 進行圖形整體環境設定。
步驟二：用 ODS 開啟一個空白的 pdf 檔。
步驟三：用 PROC GMAP 畫地圖。
步驟四：用 PROC GCHART 畫柱狀圖

最後用 ODS CLOSE 把整個圖形打包起來。

REFERENCES & RECOMMENDED READING

以下是一些不錯的線上參考資料：

http://support.sas.com/documentation/onlinedoc/index.html
http://support.sas.com/rnd/datavisualization/mapsonline/html/home.html
http://support.sas.com/rnd/papers
http://support.sas.com/samples
http://ftp.sas.com/techsup/download/sample/graph/other-color.html
http://arcdata.esri.com/data/tiger2000/tiger_download.cfm
http://www.census.gov/geo/www/cob/index.html
http://www.cdc.gov/epiinfo/shape.htm
http://research.cip.cgiar.org/gis/index.php
http://support.sas.com/techsup/technote/ts688/ts688.html
http://www.colorbrewer.org
http://robslink.com/SAS/Home.htm

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the authors at:
Robert Allison
robert.allison@sas.com
Louise Hadden
louise_hadden@abtassoc.com
Mike Zdeb
msz03@albany.edu

MCPOWER: a Flexible Macro Suite for Generating Monte Carlo Power Estimates for Linear, Logistic, and Poisson Regression Models

2011-07-09T10:17:00.001-04:00

原文載點：http://support.sas.com/resources/papers/proceedings11/430-2011.pdf

哈佛的一位統計分析師分享了一個名為 %mcpower 的巨集程式，能夠輕鬆的完成在線性迴歸、羅吉斯迴歸以及卜瓦松迴歸的蒙地卡羅檢定力模擬與估計。

%mcpower 主要包含 %glmpower 和 %manypower 兩個巨集程式。程式可從此處下載：

http://www.kenkleinman.net/home/index.php/sas-and-r-code/sas-macros.html

其中，%glmpower 是無法直接被執行的，所有的參數都是定義在 %manypower 裡面，格式如下：

%macro manypower(alt1=., alt2=., alt3=.,
                 n1=.,n2=.,n3=.,n4=.,n5=.,n6=.,n7=.,n8=.,n9=.,n10=.,
                 outcomedist = "NORMAL", null=0, scaley=1,
                 xdist="NORMAL", xloc=., xscale =.,
                 nreps=,plotrows=2, descending = desc, intercept = 0,
                 plotstyle=graphic, outputpower=YES, outputnsubs=YES, test=FALSE);

不意外地，巨集變數很多，但因為模擬本來就是一件即為複雜的工作，所以參數多也是很正常的。這些參數主要可以分類為三大群：

假設巨集變數：

null：所要檢定之虛無假設的值。預設值為 0。
alt1, alt2, alt3：所要檢定之對立假設(Ha)的數值。若是線性迴歸，則檢定的項目是斜率，若是羅吉斯迴歸，則檢定的項目是odds ratio。若是卜瓦松迴歸，則檢定的項目是risk ratio。無預設值。
n1, n2, ..., n10：每個模擬資料集裡面的樣本個數。最多可設定十個。無預設值。

分配巨集變數：

outcomedist：依變數的分配。可填入 "NORMAL"、"BERNOULLI" 和 "POISSON"。預設值為 "NORMAL"。(注意，雙引號要寫上）
scaley：殘差變異數。若使用線性迴歸(即outcomedist = "NORMAL")，則可以省略這個參數的設定。預設值為 1。
intercept：模擬的模型裡面的截距項。預設值為 0。
xdist：應變數的分配。預設值是 "NORMAL"。
xloc：應變數分配的第一個參數值。若 xdist = "NORMAL" 時，此值表示平均數。若 xdist = "POISSON" 時，則此值表示 lambda。若 xdist = "BERNOULLI" 時，則此值表示 p。
xscale：應變數分配的第二個參數值。若 xdist 為 "POISSON" 或 "BERNOULLI" 時，則此值要省略，因為這兩個分配只有一個參數值。若 xdist = "NORMAL" 時，則可設定變異數。

輸出巨集變數：

outputpower：輸出檢定力表格。預設值為 YES。
outputnsubs：輸出樣本數表格。預設值為 YES。
plotstyle：輸出表格的格式。如果設定為 sg，則圖形會用 PROC SGPLOT 程序來畫。如果設定為 graphic，則圖形會用 PROC GPLOT 來畫。如果沒有設定，則圖形就不會輸出。預設值為 graphic。
plotrows：若 plotstyle = sg 時，則可設定圖形要有幾行。預設值為 2。

範例：

%manypower(null = 0, scaley=1,xdist = "bernoulli", xloc=.5, xscale=.,
           nreps = 100, n1=1000,n2=1100,n3=1200,n4=1400,n5=1700,n6=2000,
           alt1 = 1.15, alt2 = 1.20, alt3 = 1.25,
           plotrows=2,plotstyle=graphic, outcomedist="poisson", intercept = 1);

結果：

CONTACT INFORMATION
Your comments and questions would be welcome and valued. Contact the author at:
Ken Kleinman
Harvard Medical School and Harvard Pilgrim Health Care
133 Brookline Ave.
Boston, MA 02215
ken.kleinman@gmail.com

Fundamental Diagnostics for Two-Level Mixed Models: The SAS ® Macro MIXED_DX

2011-07-02T10:57:00.001-04:00

原文載點：http://support.sas.com/resources/papers/proceedings10/201-2010.pdf

在進行 mixed model 分析時，模型的檢定是相當重要但也是相當煩人的一環。幾位南卡大和南佛大的學者寫了一個巨集程式來簡化相關程式的寫作，雖然目前這個程式只能使用在 2-level model，但也是具有相當貢獻。

這個名叫 %mixed_dx 的巨集程式如下所示：

%mixed_dx (ModelI = ,
                    Dims = ,
                    solnF = ,
                    solnR =,
                    Level1 = ,
                    Influence = ,
                    min_nj = ,
                    pr = ,
                    GroupVar = );

每個參數的設定如下：

ModelI：讀取由ODS OUTPUT的ModelInfo選項所生成的包含所有模型資訊的資料
Dims：讀取由ODS OUTPUT的Dimensions選項所生成的資料維度數據
solnF：讀取由ODS OUTPUT的SolutionF選項所生成的固定效應參數估計值
solnR：讀取由ODS OUTPUT的SolutionR選項所生成的隨機效應參數估計值
Level1：讀取由model statement的outp選項所生成的Level-1殘差值資料
Influence：讀取由model statement的influence選項所生成的模型檢定數據
min_nj：設定Level-2 unit的最小組內樣本數
pr：定義檢定內所需要用到的threshold percentile rank(預設值=90)
GroupVar：若有使用名義(nominal)類別的Level-2 unit的因子，則必須要把每個名義類別的定義另存一個新檔並用此巨集變數來呼叫。

首先，利用PROC MIXED程序去估計一個模型，並且利用ODS OUTPUT的功能把%mixed_dx要用到的資料給另存出來：

ods exclude influence SolutionR;
title;
proc mixed data = temp covtest noclprint NAMELEN=32;
      class schoolid;
      model mathach = SIZE FEMALE SES MEANSES SES*MEANSES SES*SIZE FEMALE*SIZE/ solution outp=L1Resid influence(effect=schoolid iter=5);
      random  intercept FEMALE SES / sub=schoolid solution type=un;
      ods output SolutionR=L2Resid SolutionF=Fixed ModelInfo=ModelStuff Dimensions=DatStuff Influence=influence;
run;

其中，outp選項所輸出的資料L1Resid是level 1的殘差值，這是要用在Level1參數裡面的。橘色高亮部分即為原始PROC MIXED內定的模型檢定指令，若配合ODS HTML的話可以輸出高解析度的圖表，但就是因為裡面包含太多圖了，所以才需要使用%mixed_dx來把幾個重要的圖整理出來。而紅色部分即為ODS OUTPUT輸出的指令。此處是缺一不可，即為重要。

接下來就來執行巨集程式：

%mixed_dx (ModelI = modelstuff, Dims = datstuff, solnF = Fixed, solnR = L2Resid, Level1 = L1Resid, Influence = influence, min_nj = 5, pr = 90, GroupVar = None);

在生成的圖形當中，主要以殘差的箱型圖、柱狀分布圖以及散布圖為主。另外也有殘差變異數的散布圖（用來檢定同質性）還有離群值的柱狀分布圖。如下所示：

表格方面，會整理出level-2 unit的常態統計量表、Levene’s同質性檢定表、離群值檢定和各種檢定統計量：

目前這段程式碼可以在他們的校方網站上面下載：http://www.ed.sc.edu/bell/

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author at:
Bethany A. Bell
University of South Carolina
College of Education, Wardlaw 133
Columbia, SC 29208
Work Phone: 803-777-2387
E-mail:babell@sc.edu
Web: http://www.ed.sc.edu/bell/

Expanding PROC EXPAND

2011-06-25T10:35:00.000-04:00

原文載點：http://support.sas.com/resources/papers/proceedings11/417-2011.pdf

PROC EXPAND 程序是 SAS/ETS 裡面其中一個語法，提供許多跟時間相關資料的簡便處理。比方說，將日資料累加成月資料，月資料累加成季資料或年資料，也可以倒過來把季/年資料變成月資料，月資料變成日資料。若遇到缺漏值，該程序可用內插法將資料補齊。也可進行諸如moving average的資料轉換。Bruce Gilsen 在 2011 年的 SAS Global Forum 發表了一篇技術文件，介紹了兩個用 PROC EXPAND 運作的例子。

Calculate an un-weighted quarterly average from monthly data

第一個例子介紹如何將月資料轉成季資料。使用的範例資料如下(部分)：

  Obs     date         id   income
    1    20070101    1      1
    2    20070201    1      2
    3    20070301    1     15
    4    20070401    1      4
    5    20070501    1      5
    6    20070601    1     12
    7    20070701    1      7
    8    20070801    1      8
    9    20070901    1      9

程式如下：

proc expand data=one out=two from=month to=quarter ;
      id date;
      by id;
      convert income / method=aggregate observed=average;
run;

紅色部分的 from=month to=quarter 顧名思義就是把資料從月(month)轉成季(quarter)。藍色部分表示轉換的變數為 income，方法為累加(aggregate)，最後求出"非加權"平均值(average)。結果如下：

   Obs    id      date     income    income
                                               (expected result)
    1      1    2007:1     6.1333       6
    2      1    2007:2     6.9780       7
    3      1    2007:3     7.9891       8

最後一欄income(expected result)便是所求，不過程式也會把加權後的結果一併列出，即為倒數第二欄的結果。加權的方式是用天數為基準。以第三季(七～九月)結果來看，前三個月的天數分別是31, 31, 30。總計為92天，因此加權平均後的結果是

((31 * 7) + (31 * 8) + (30 * 9)) / 92 = 7.9891

Convert weekly data to a monthly average when the weekly values represent values for seven days

先看一下範例：

   Obs   date       id   income
    1    20070530    1     11
    2    20070606    1     30
    3    20070613    1     50
    4    20070620    1     70
    5    20070627    1     90
    6    20070704    1     50
    7    20070711    1     40
    8    20070718    1    120
    9    20070725    1     80
   10    20070801    1    100
   11    20070808    1      5
   12    20070815    1      7
   13    20070822    1      9
   14    20070829    1     11
   15    20070905    1     24

將週資料轉成月資料的方法同第一個範例：

proc expand data=one out=two from=week to=month ;
      id date;
      by id;
      convert income / observed=average;
run;

結果如下：

   Obs    id      date            income    income
                                                       (expected result)
    1      1    JUN2007      57.093     60
    2      1    JUL2007       75.905     80
    3      1    AUG2007     18.852     12

只要把 from=quarter 改成 from=week 即可。但是加權平均income的數據，是用插補法先將週四至下週二的數據進行插補，然後再做平均。作者提供另一種作法：

data two;
     set one;
     do dailydate = date-6 to date;
       output;  /* transform each weekly observation to 7 daily observations */
     end;
   run;
proc expand data=two out=three from=day to=month ;
      id dailydate;
      by id;
      convert income / observed=average;
run;

這個作法是先將原始的週資料用複製的方法將沒有數字的天數補齊，如下所示：

   Obs   dailydate    id     income
    1    20070524      1       11
    2    20070525      1       11
    3    20070526      1       11
    4    20070527      1       11
    5    20070528      1       11
    6    20070529      1       11
    7    20070530      1       11
    8    20070531      1       30
    9    20070601      1       30
   10    20070602      1       30
   11    20070603      1       30
   12    20070604      1       30
   13    20070605      1       30
   14    20070606      1       30
   15    20070607      1       50
   16    20070608      1       50
   17    20070609      1       50
   18    20070610      1       50
   19    20070611      1       50
   20    20070612      1       50
   21    20070613      1       50
   22    20070614      1       70
   23    20070615      1       70
   24    20070616      1       70
   25    20070617      1       70

然後再把這個日資料換算成月資料即可。

Contact information
Bruce Gilsen
Federal Reserve Board
Mail Stop 157
Washington, DC 20551
e-mail: bruce.gilsen@frb.gov
phone: 202-452-2494.

MORE CUSTOMIZATION?: CREATING SYMBOLS IN RTF FILES USING ODS

2011-06-18T10:49:00.001-04:00

原文載點：http://www.lexjansen.com/pharmasug/2004/technicaltechniques/tt12.pdf

有時候想要在 SAS 輸出報表裡面寫一些標題或註解，需要加上一些科學符號、希臘字母、甚至是邏輯運算子等特殊符號時，得使用特別的方法顯示出來。這一篇 2004 年的 PHARMASUG 技術文件可供參考。

首先必須先去查這些特殊符號的 Unicode decimal number，因為 SAS 必須藉由這些代碼才能對應到特殊符號。查詢的網址如下：

http://www.alanwood.net/unicode/index.html

http://www.unicode.org/

以"alpha"這個希臘字母當做範例，程式如下：

ods escapechar="^";
%let symbol_alpha=%bquote(^R/RTF"\u945\") ;

第一行用ods escapechar來定義一個轉義字母，在此設定為"^"。SAS 看到轉義字母才知道後面的 Unicode 代碼是要轉成什麼特殊符號。接著用 %bquote 函式來把特殊符號做出來。把上面的程式複製到 SAS，再把紅色的 Unicode 換成你想要的 Unicode 即可。

一些常用的特殊符號的 Unicode 代碼：
希臘字母：http://www.alanwood.net/unicode/greek.html
數學符號運算子：http://www.alanwood.net/unicode/mathematical_operators.html
箭頭：http://www.alanwood.net/unicode/arrows.html
羅馬數字和分數：http://www.alanwood.net/unicode/number_forms.html
上標與下標：http://www.alanwood.net/unicode/superscripts_and_subscripts.html
其他：http://www.alanwood.net/unicode/letterlike_symbols.html

原文提供一個範例。假設要製作下列表格的註解：

程式如下：

footer f1 ;
define f1 ;
text
   " +Above Normal Range" "^1n"
   "-Below Normal Range" "^1n"
   "A reading &symbol_ge 40 for SGPT/ALT is considered Above Normal"  "^1n"
   "A reading &symbol_le 30 for SGPT/ALT is considered Below Normal"   "^1n" 
   "An example of &symbol_alpha." "^2n "
   "&font9.Output:" " &tab.&outdir.\&pgmonly..rtf" "^1n"
    "Source:&tab.&source."  ;
   style={font_face=Arial font_size=2.5 font_weight=medium background=white} ;
   just=left ;
end ;

如果覺得這樣太麻煩，就送去Word編輯吧！

CONTACT INFORMATION

Richard Rowell
Connetics Corporation
3290 West Bayshore Road
Work: (650) 843-2876
rrowell@connetics.com
Jim Lenihan
15 Moonlight Ct
South San Francisco, CA 94080
Work: (650) 742-0131
jameslenihan@sbcglobal.net

Creating Forest Plots from Pre-computed Data using PROC SGPLOT and Graph Template Language

2011-06-04T00:53:00.000-04:00

原文載點：http://support.sas.com/resources/papers/proceedings10/195-2010.pdf

Forest plot 是研究環境流行病學的人經常會用到的一種圖表，主要是可以將許多變數或樣本的參數估計值以及信賴區間畫成一條一條的橫線，以供其他人易於判讀相對大小以及有無顯著。我在自己的博士論文裡面也放了相當多這種圖型，但之前我都是自己用 R 寫程式畫的，過程挺複雜的，沒想到 SAS 裡面有現成的語法可以畫這種圖。
假設用 PROC LOGISTIC 來估計一個模型，如下所示：

ods output "Odds Ratios"=orci; 
proc logistic data=uis descending; 
      model dfree=age beck ivhx ndrugtx race treat site ; 
run;

估計好後把所有參數的 Odds Ratio 另存成 orci 這個新資料集，內容如下所示：

然後用 PROC SGPLOT 直接來畫 Forest plot：

title "SGPLOT: Forest Plot"; 
proc sgplot data=orci; 
  scatter x=oddsratioest y=effect / xerrorlower=lowercl                
                                    xerrorupper=uppercl  
                                    markerattrs=or(symbol=DiamondFilled size=8); 
  refline 1 / axis=x;   
  xaxis label="OR and 95% CI " min=0; 
  yaxis label="Covariates"; 
run;

陽春一點的圖型只需要四個語法：

scatter：後面的 x= 放估計參數變數名稱，y= 放參數本身變數的名稱。後面三個選項分別是：

xerrorlower = 信賴區間下限
xerrorupper = 信賴區間上限
markerattrs = 設定估計量的符號形式及大小

refine：設定在哪一個軸畫基準線
xaxis：設定 x 軸
yaxis：設定 y 軸

做出來的圖型如下所示：

原文內還有其他的教學，主要著重於如何調整圖型品質，但根據自己的經驗是上述這種圖就差不多可以拿去放在 paper 裡面了，所以就不多做介紹。有興趣的人可以再去參閱原文後段的進階內容。

CONTACT INFORMATION
Zoran Bursac, PhD, MPH
Associate Professor
Biostatistics
Fay W. Boozman College of Public Health
University of Arkansas for Medical Sciences
4301 W. Markham, Slot 781
Little Rock, AR 72205
Work Phone: (501) 526-6723
Fax: (501) 526-6729
E-mail: zbursac@uams.edu
Web: www.uams.edu/biostat/bursac/

Automagically Copying and Pasting Variable Names

2011-05-28T06:00:00.001-04:00

原文載點：http://support.sas.com/resources/papers/proceedings10/046-2010.pdf

如果使用者從第三處拿到一筆資料，想要建立coding book或目錄，或者需要寫新的程式必且在keep或drop指令輸入大量變數名稱，因而需要大量複製和貼上變數名稱的話，有什麼比較好的方法呢？Arthur S. Tabachneck 等人提供了一個方法讓使用者能夠一次複製且貼上全部的變數名稱。

首先，執行下面這段程式碼：

FILENAME _cb CLIPBRD;
DATA _NULL_;
   WINDOW DSN rows=8 columns=80
   irow=1  icolumn=2 color=black
   #2  @3  'Enter 1 or 2 level data set name: ' 
            color=gray dsn $41.  required=yes
            attr=underline       color=yellow;
   DISPLAY DSN blank;
   FILE _cb;
   length name $32;
   do dsid = open(dsn,'I') while(dsid ne 0);
     do i = 1 to attrn(dsid,'NVARS');
       name = varname(dsid,i);
       put name @;
     end;
     dsid = close(dsid);
   end;
RUN;
FILENAME _cb CLEAR;

切記裡面一個指令都不用改，執行完後 SAS 裡面會跳出一個黑色的視窗，如下所示：

裡面有一個空行，是要你輸入資料所在的 library 名稱還有資料本身名稱，如：「Orz.wahaha」，亦即表示你要複製 Orz library 裡面的 wahaha 資料內的所有變數名稱。如果沒有指定 library，則程式會自行去 Work library 裡面搜尋。輸入好後按 Enter，要進行最後"關掉DSN視窗"的動作，這時程式才會將你指定的資料集裡面的所有變數名稱複製好，然後你就可以用 Ctrl + v 的熱鍵把複製好的變數名稱任意地貼到任何地方或任何軟體裡面。

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the authors at:
Arthur Tabachneck, Ph.D. Randy Herbison,
Director, Data Management Senior Systems Analyst
Insurance Bureau of Canada Westat
2235 Sheppard Ave. East 1650 Research Boulevard
Toronto, ON L3T 5K9 Canada Rockville, MD 20850
E-mail: atabachneck@ibc.ca E-mail: RandyHerbison@westat.com
John King Andrew Clapson
Ouachita Clinical Data Services, Inc. Ottawa, ON Canada
Mount Ida, AR E-mail: andy_clapson@hotmail.com
ouachitaclinicaldataservices@gmail.com
Roger DeAngelis Tom Abernathy
CompuCraft Inc Pfizer, Inc.
1770 Via Petirrojo Apt A 235 E. 42nd Street
Newbury Park CA 91320 New York, NY 1001
E-mail: xlr82sas@aol.com E-mail: tom.abernathy@pfizer.com

Top Ten SAS® Sites for Programmers: A Review

2011-05-21T20:57:00.003-04:00

原文載點: http://support.sas.com/resources/papers/proceedings11/055-2011.pdf

一位相當資深的 SAS 程式設計師在 SAS Global Forum 2011 發表了一篇文章推薦了十個值得瀏覽的 SAS 技術網站。這十個網站是根據三百位SAS使用者的抽樣調查而來，並非原作者主觀意見，所以算是滿客觀的一個結論。雖然不是說每個我都很喜歡，但在此分享給大家。

Support.SAS.com
SAS 客服網站。SAS 的線上手冊也在裡面。我個人也推薦這個網站是因為除了程式上的技術問題外，使用者有時候也會遇到軟體上的問題，通常這些問題都可以透過 SAS 客服來解決。我以前經常寫信去請教 SAS 客服解決一些怪異的 error message 或記憶體問題，而對方總是相當快速地給我回覆。當然，前提是你用的 SAS 必須是正版的
LexJansen.com
相當老牌的一個 SAS 教學網站，裡面集合了所有 SAS 官方或非官方的論壇文章，本站側欄的 Good Sites 就是參考這個網站而設計的。
sasCommunity.org
類似 SAS 版的維基百科，不過我沒用過所以不熟。
SAS-L archive
也是相當老牌的 SAS 討論群組，可以在上面發問或回答別人的問題。缺點是介面老舊，彷彿像 90 年代的風格，也不是很人性化。
SAS.com
SAS 官方網站，裡面有相當多 SAS 官方提供的巨集程式可以使用。
NESUG.org
北美東 SAS 使用者論壇，我投過一次 poster 被接受但因為生病所以就放棄沒參加了。
blogs.SAS.com
一堆 SAS 大頭在裡面寫的部落格，個人是覺得沒有甚麼用處。
SESUG.org
南美東 SAS 使用者論壇。我曾經發表過一篇技術文件在那邊。
SAS Global Forum
本站參考文獻主要來源，強力推薦!!!
Google.com
估狗大神，無須多言。

Comments and suggestions can be sent to:

Charles Edwin Shipp

JMP 2 Consulting, Inc.

E-mail: CharlieShipp@aol.com

Kirk Paul Lafler

Software Intelligence Corporation

E-mail: KirkLafler@cs.com

Practicalities of Using ESTIMATE and CONTRAST Statements

2011-05-21T01:26:00.001-04:00

原文載點：http://support.sas.com/resources/papers/proceedings10/269-2010.pdf

在 PROC GLM 和 PROC MIXED 裡面，最令一般 SAS 使用者感到頭痛的不是模型的檢定而是事後分析(post hoc test)。這牽涉到 contrast 和 estimate 兩個語法。在我接觸過所有非統計科班出身但需要用到線性模型以及事後分析的人，沒有人不對這兩個語法感到聞風喪膽。多年人我曾經試圖用各種方法來教人如何使用這兩種語法，但效果都不怎麼顯著。簡單來講，Contrast 和 estimate 無論是在 PROC GLM 還是 PROC MIXED 裡面，功能完全一樣，語法內的選項也一樣，而兩者之間唯一的差別是，contrast 可以只能算出 p-value，而 estimate 可以把指定的點估計量及其 p-value 一次輸出。David J. Pasta 在 SAS Global Forum 2010 發表了一篇技術文件將這兩種語法做一個詳盡的介紹，堪稱一大貢獻。在本文中，我會介紹幾個比較實用的範例。
原文作者套用一個醫學資料並詳述其背景，但為了簡化初學者的時間，我將會以程式來主。第一個範例資料只有三個變數，分別是作為應變數的離散型變數 sex (兩層) 和 race (五層)，而依變數則為一個連續型變數 y1。有了對資料的基本概念，首先來看一個最簡單的例子：

proc glm data=anal;
      class sex;
      model y1 = sex / solution;
      lsmeans sex / stderr tdiff e;
run;

這個例子只單獨分析 sex 對 y1 的影響，而 lsmeans 是用來算出每個性別對 y1 的最小平方估計量。後面的 tdiff 則是去求兩個性別的差異。這段程式會跟用 PROC TTEST 做出來的結果一模一樣。如果不用 lsmeans 而改用 estimate 語法，則可以這樣寫：

proc glm data=anal;
      class sex;
      model y1 = sex / solution;
      estimate 'male' intercept 1 sex 0 1;
      estimate 'female' intercept 1 sex 1 0;
      estimate 'female-male' sex 1 -1;
run;

頭兩個 estimate 語法可算出不同性別在 y1 上的平均值，後面接著 intercept 和 sex 是為了要告訴 SAS 估計值是從哪些變數來的。一般使用者最常忘記的就是把 intercept 1 加上去，因為除非在 model 的選項上面加上 noint 抑制掉截距項的計算，否則一切用 estimate 求算的點估計量都要加上 intercept 1。不過你也不用擔心不加會怎麼樣，因為 log 視窗裡面自然會出現紅色的錯誤訊號提醒你加上。而 sex 後面跟著「0 1」和「1 0」則是要準確地將要計算的性別定位出來，當然你也要對自己如何將 sex 編碼要非常了解，而不是亂安插一通。在本例中，女性的編碼順位是排在男生前面（若真的不清楚編碼順位，看一下輸出報表一開始的 Class Level Information 裡面就會清清楚楚地把所有類別變數的編碼順位列出來），所以當 sex 後面寫「1 0」時，自然是去算女生的值，反之則是算男生的值。

而當要比較女生和男生的差異時，則有兩點要注意。一是 intercept 可以不用寫了，因為兩兩相減的情況下，截距項自然就被減掉了。二是 sex 後面要改成「1 -1」，很明顯地這就是指用女生的估計值去減男生的估計值。如果要算男生的估計值減女生的估計值，把「1 -1」改成「-1 1」即可。偷懶的方法是不改，但把 output 輸出的估計值變號即可，至於 p-value 是不會改變的。

當考慮兩個應變數時，程式如下：

proc glm data=anal;
      class race sex;
      model y1 = race sex / solution e;
      lsmeans race sex / stderr tdiff e;
      lsmeans race sex / stderr tdiff e om;
run;

第一個 lsmeans 語法會一次把兩個變數的最小平方估計量通通算出來，並且會把成對比較的結果也列印出來。第二個 lsmeans 語法和第一個的不同是後面多了 om 這個選項。它的功用是去計算不同類別的權重，並在計算平均數的時候考量到權重大小。簡單來講就是去計算加權的估計量。以 sex 為例，在沒有考量加權的情況下，每個 race 的權重都會是 0.2。而考量權重的話，每個 race 的權重則等於每個種族在全部樣本裡面所佔的比例。

若用 estimate 語法則程式如下：

proc glm data=anal;
      class race sex;
      model y1 = race sex / solution e;
      estimate 'male' intercept 1 sex 0 1;
      estimate 'female' intercept 1 sex 1 0;
      estimate 'female-male' sex 1 -1;
run;

寫法跟第一個程式範例一樣，但我會習慣加上「race 0.2 0.2 0.2 0.2 0.2」，只是說不加的話也沒關係，因為程式會自動用等同的權重去平均。寫到這邊只有兩點要強調：一是類別變數有幾層，estimate 語法內變數後面要跟的數字就有幾個。有人習慣將 1 後面可能還會出現的 0 都省略掉，但我強烈建議補上，因為多寫幾個 0 不會花太多時間，而且有可好處是可以把所有的 estimate 語法排列的工工整整，要 debug 也比較方便。二是腦筋一定要清楚知道變數的編碼順位，在寫之前最好再去對一下 Class Level Information 表格內的訊息。

當考慮交互作用項時，estimate 語法就沒有辦法簡化了：

proc glm data=anal;
      class sex race;
      model y1 = race sex race*sex / solution;
      estimate 'male'   intercept 1 sex 0 1 race .2 .2 .2 .2 .2  race*sex 0 0 0 0 0  .2 .2 .2 .2 .2;
      estimate 'female' intercept 1 sex 1 0 race .2 .2 .2 .2 .2  race*sex .2 .2 .2 .2 .2  0 0 0 0 0;
      estimate 'female-male' sex 1 -1 race*sex .2 .2 .2 .2 .2  -.2 -.2 -.2 -.2 -.2;
run;

在 sex 的點估計中，除了「race 0.2 0.2 0.2 0.2 0.2」一定要加上外，交互作用項 race*sex 後面要寫的數字向量變將整個程式複雜程度提昇了許多。由於 race 有五層，sex 有兩層，所以交互作用項就會是 5X2 = 10 層。因此開始寫這段數字向量時，先要知道到底該寫幾個數字進去，這是第一個要注意的事情。其次，將 sex 的向量和 race 的向量用 Kronecker product 的方式去算。如果不清楚 Kronecker product，我在這邊提供一個小程式碼：

proc iml;
      a = {0 1};
      b = {0.2 0.2 0.2 0.2 0.2};
      ab = a@b;
      print ab;
quit;

自行替換 a 和 b 裡面的數字就可以算出兩個變數的交互作用項的向量。然後就直接從 output 視窗裡面複製貼上即可。在成對比較上面，交互作用的數字向量寫法也是一樣可以用這個小程式去算，只要把{0 1}改成{1 -1}即可。但在寫程式時，「race 0.2 0.2 0.2 0.2 0.2」因為在相減的過程中被扣掉了，所以不用寫上去。一個比較白痴的方法是把第一個 estimate 和第二個 estimate 排好，然後下面減上面，自然就得到最後一個的數字向量。只是當你有大量變數要比較時，這種作法就有點慢了。

可是當考慮權重時，程式如下：

proc glm data=anal;
      class sex race;
      model y1 = race sex race*sex / solution e;
      estimate 'male' intercept 1  sex 0 1  race  .104 .142 .128 .047 .579 race*sex 0 0 0 0 0  .104 .142 .128 .047 .579;
      estimate 'female' intercept 1  sex 1 0 race  .104 .142 .128 .047 .579 race*sex .104 .142 .128 .047 .579  0 0 0 0 0;
      estimate 'female-male' sex 1 -1 race*sex .104 .142 .128 .047 .579 -.104 -.142 -.128 -.047 -.579;
      estimate 'female-male (off)' sex 1 -1 race*sex .105 .142 .128 .047 .579 -.105 -.142 -.128 -.047 -.579;
run;

紅色部分就是 race 在五個不同總族下的權重，這可以用 PROC FREQ 額外算出來。只要這個部分解決，後面交互作用項就簡單了，因為可以直接套用上面我提供的小程式，之後便是一連串的剪貼。兩兩比較的技巧也如同之前所述。最後一個標記為 'female-male (off)' 的程式，只是要提醒，在計算權重時，我們可能會使用到進位的權重，但進位後可能會造成無法估計的結果。本例中，由於進位的關係，race*sex 後面的數字總和雖然還是零，但是 race 權重的總和變成 1.001 而不是 1.000 了。因此你就會看到一段紅色錯誤訊息出現在 log 視窗告訴你這個事後比較無法估計。解決的方法沒有特別，就是手動去微調裡面的數據，東加一點西減一點，調到總合是 1 即可。
有時候事後比較的目的是要看某離散型的變數是否有線性趨勢，則程式如下：

proc glm data=anal;
      class sex educat;
      model y1 = sex educat / solution;
      lsmeans sex educat / stderr e;
      lsmeans sex educat / stderr e om;
      estimate 'educat linear (no divisor)' educat -2 -1 0 1 2;
      estimate 'educat linear (divisor=10)' educat -2 -1 0 1 2 / divisor=10;
      estimate 'educat linear (divisor=40)' educat -4 -2 0 2 4 / divisor=40;
      estimate 'educat by year (not centered)' educat 8 12 14 16 20;
      estimate 'educat by year (centered)' educat -6 -2 0 2 6 / divisor=80;
run;

此程式把 race 換成 educat，但一樣是五層。線性趨勢的事後檢定所使用的數字向量為「-2 -1 0 1 2」，這個比例大致上是不會改變的，當然不同數目的層次有不同的寫法，請自行孤狗。有時會有比較特殊的寫法，就比方'educat by year (not centered)'的寫法，不過情況比較少。大部分還是會寫成像'educat by year (centered)'所示，也就是讓數字向量內的數字是對秤的，而不是單調遞升或遞減的。至於 divisor 的用法，在此例沒有特殊意義，只是當做是分母把要用來估計的數字向量一次全部除掉，換句話說，以第二個'educat linear (divisor=10)'為例，你也可以寫成

estimate 'educat linear (divisor=10)' educat -0.2 -0.1 0 0.1 0.2;

但是 divisor 的真正威力，是在下面這個範例：

proc glm data=anal;
      class sex race;
      model y1 = race sex race*sex / solution e;
      estimate 'male (333)' intercept 1  sex 0 1  race  .333 .333 .333 race*sex 0 0 0  .333 .333 .333; 
      estimate 'female (333)' intercept 1  sex 1 0 race  .333 .333 .333 race*sex  .333 .333 .333  0 0 0; 
      estimate 'female-male (333)' sex 1 -1 race*sex  .333 .333 .333  -.333 -.333 -.333; 
      estimate 'male (334)' intercept 1  sex 0 1  race  .333 .333 .334 race*sex 0 0 0  .333 .333 .334; 
      estimate 'female (334)' intercept 1  sex 1 0 race  .333 .333 .334 race*sex  .333 .333 .334  0 0 0; 
      estimate 'female-male (334)' sex 1 -1 race*sex  .333 .333 .334  -.333 -.333 -.334; 
      estimate 'male (d)' intercept 3  sex 0 3  race  1 1 1  race*sex 0 0 0  1 1 1 / divisor=3; 
      estimate 'female (d)' intercept 3  sex 3 0 race  1 1 1  race*sex 1 1 1  0 0 0 / divisor=3; 
      estimate 'female-male (d)' sex 3 -3 race*sex 1 1 1   -1 -1 -1 / divisor=3; 
      estimate 'male (o)' intercept 849  sex 0 849  race  142 128 579 race*sex 0 0 0  142 128 579 / divisor=849; 
      estimate 'female (o)' intercept 849  sex 849 0 race  142 128 579 race*sex 142 128 579  0 0 0 / divisor=849; 
      estimate 'female-male (o)' sex 849 -849 race*sex 142 128 579 -142 -128 -579 / divisor=849; 
run;

此例中，race 不再是五層，而是三層！所以當要計算 sex 內男女分別的點估計量，SAS 裡面是不允許寫分數而是得直接寫小數上去，造成 race 後面的數字向量不能是{1/3 1/3 1/3}而是{0.333 0.333 0.333}。這種情況使得三個數字的總和變成 0.999 而不是 1.000，然後就會看到 SAS 的 log 視窗給你錯誤的訊息了。有一種解決方法是改成{0.333 0.333 0.334}，但一旦層數變多了，如七層、九層、十一層、二十一層....，則不但整條 estimate 又臭又長，每次都要去微調裡面的數字也很麻煩，所以此時就要靠 divisor 這個選項幫我們一次除盡所有數。首先，先把原本每個變數後面的數字都乘上 3，接著在最後面加上 / divisor = 3 即可。如果考量到 race 的權重，請回憶國小國中學過的最小公倍數。你可能會問，一對小數點怎麼去求最小公倍數呢？在此額外提供一個 SAS 小程式方便你的計算：

data _null_;
    x = lcm(a1,a2,....,aN);
    put x;
run;

把你算出來的權重全部丟到 lcm function 裡面，這個函數顧名思義就是最小公倍數 least common multiple 的縮寫，馬上就可以算出 divisor 該填什麼數字進去了。

以上就是很基本但卻很實用的 estimate 寫法。contrast 的用法一模一樣，所以沒有多做介紹。原文後面還有一些更進階的語法，但我個人的經驗是極少數的情況會使用到，因此大家熟練以上的程式就夠了。不過有時候離散變數內的層數很多，在建立 estimate 語法的時候，後面會跟著一長串的數字向量，非常驚人（恐怖）。原文作者建議將一些會一直重複用到的數字向量用 %let 先定義好，等到要用的時候再呼叫即可，比方說：

%let S_1    = 1 0 0 0 0 0 0 0 0 0;
%let S_2    = 0 1 0 0 0 0 0 0 0 0;
%let S_3    = 0 0 1 0 0 0 0 0 0 0;
%let S_4    = 0 0 0 1 0 0 0 0 0 0;
%let S_5    = 0 0 0 0 1 0 0 0 0 0;
%let S_6    = 0 0 0 0 0 1 0 0 0 0;
%let S_7    = 0 0 0 0 0 0 1 0 0 0;
%let S_8    = 0 0 0 0 0 0 0 1 0 0;
%let S_9    = 0 0 0 0 0 0 0 0 1 0;
%let S_10   = 0 0 0 0 0 0 0 0 0 1;
%let S_none = 0 0 0 0 0 0 0 0 0 0;
%let P_U    = .1 .1 .1 .1 .1 .1 .1 .1 .1 .1;
%let P_O    = 0.0666048524 0.0776077886 0.0888270746 
              0.0955339206 0.1017462525 0.1059496214 
              0.1124710246 0.1164580436 0.1168907433 
              0.1179106784 ;
  *** Intercept at the end of the BEFORE and at the beginning of the AFTER ***;
  estimate 'LS:1_Before'   intercept             1
                           decile                &S_1.
                           decile*tafter         &S_none.;      
  estimate 'LS:1_After'    intercept             1
                           decile                &S_1.
                           decile*tafter         &S_1.;   
  estimate 'LS:2_Before'   intercept             1
                           decile                &S_2.
                           decile*tafter         &S_none.;      
  estimate 'LS:2_After'    intercept             1
                           decile                &S_2.
                           decile*tafter         &S_2.;      
 . . . 
  estimate 'LS:10_Before'  intercept             1
                           decile                &S_10. 
                           decile*tafter         &S_none.;      
  estimate 'LS:10_After'   intercept             1
                           decile                &S_10.
                           decile*tafter         &S_10.;      
  estimate 'LS:Before_U'   intercept             1
                           decile                &P_U.
                           decile*tafter         &S_none.; 
  estimate 'LS:After_U'    intercept             1
                           decile                &P_U.
                           decile*tafter         &P_U.;  
  estimate 'LS:Before_O'   intercept             1
                           decile                &P_O.
                           decile*tafter         &S_none.;   
  estimate 'LS:After_O'    intercept             1
                           decile                &P_O.
                           decile*tafter         &P_O.;

試想，如果不用 %let 去額外定義巨集變數，而把那些 0 和 1 或 -1 等數字全部填回到 estimate 後面，則整串程式碼絕對會讓人暈眩。因此，我強烈建議使用 %let 來定義好 estimate 後面該用到的數字向量。此外，當要考量權重時，雖然可以自行算好寫上，但當層數多的時候（如本例），一個個複製貼上也是很麻煩，所以原文作者提供了一段程式碼讓 SAS 自動幫你把把這些權重貼到巨集變數裡面。

proc sort data=temp01(keep=patid_age decile) out=temp02 nodupkey;
  by patid_age decile;
run;
data _null_;
  set temp02;
  by patid_age decile;
  if not(first.patid_age and last.patid_age) 
then error 'ERROR: decile varies within patid_age';
run;
ods listing close;
ods output onewayfreqs = temp03(keep=decile frequency cumfrequency);
proc freq data = temp02;
  tables decile;
run;
ods listing;
data _null_;
  set temp03;
  call symput('N_'||left(trim(put(decile,2.))),frequency);
run;

data _null_;
  set temp03;
  if (decile eq 10) then call symput('Ntot',cumfrequency);
run;
data _null_;
  call symput('P_1', &N_1. / &Ntot.));
  call symput('P_2', &N_2. / &Ntot.);
  call symput('P_3', &N_3. / &Ntot.);
  call symput('P_4', &N_4. / &Ntot.);
  call symput('P_5', &N_5. / &Ntot.);
  call symput('P_6', &N_6. / &Ntot.);
  call symput('P_7', &N_7. / &Ntot.);
  call symput('P_8', &N_8. / &Ntot.);
  call symput('P_9', &N_9. / &Ntot.);
  call symput('P_10', &N_10. / &Ntot.);
run;
*** Macro variable for model specifications ***;
%let P_O = &P_1. &P_2. &P_3. &P_4. &P_5. &P_6. &P_7. &P_8. &P_9. &P_10.;

這段程式碼簡而言之，就是先用 PROC FREQ 把次數分配表的結果算出來，然後把這必要的變數（也就是 frequency 和 cumfrequency）存出來計算每一層的比例，再把這比例結果用 symput 函數存成不同的巨集變數，然後一次丟進 %let 裡面整合成一個巨集變數供 estimate 語法使用。不過我覺得原文作者寫的這段程式稍嫌複雜些，可能因為他用的資料結構比較複雜的因素。在一般的情況下，用 PROC FREQ 可以立刻求出類別變數內各層的比例：

proc freq data=xxx;
      tables varcat;
      ods output onewayfreqs = temp(keep=percent);
run;

用 ods output 另存的資料裡面有個 percent 變數就是每一層的比例。然後用 proc transpose 轉置：

proc transpose data=temp out=temp;
      var percent;
run;

假設有十層的話，轉置後的 percent 會被存在 COL1~COL10 裡面。然後用 symput 函數分別將他們存成不同的巨集變數：

data _null_;
      set temp;
      call symput('P_1', COL1);
      call symput('P_2', COL2);
      call symput('P_3', COL3);
      call symput('P_4', COL4);
      call symput('P_5', COL5);
      call symput('P_6', COL6);
      call symput('P_7', COL7);
      call symput('P_8', COL8);
      call symput('P_9', COL9);
      call symput('P_10', COL10); 
run;

最後再全部丟到 %let 後面用一個巨集變數整合他們：

%let P = &P_1 &P_2 &P_3;

如此一來，在寫 estimate 時需要用到這串數字向量時，輕鬆填上&P 即可！

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author at:
David J. Pasta
Vice President, Statistics & Data Operations
ICON Clinical Research
188 Embarcadero, Suite 200
San Francisco, CA 94105
+1.415.371.2111
david.pasta@iconplc.com
www.iconplc.com

Updates to SAS® Power and Sample Size Software in SAS/STAT® 9.2

2011-05-14T00:43:00.000-04:00

Link: http://www2.sas.com/proceedings/forum2008/368-2008.pdf

PROC POWER程序從SAS 9.1版發表出來後，讓許多power和樣本計算的過程不需要靠使用者自行寫程式才能完成。但一開始的功能並不完備，不過在SAS 9.2版釋出後，PROC POWER 內建的功能就逐漸被補齊了。Wayne Watson 於 SAS Global Forum 2008 發表了一份技術文件，整理了 SAS 9.2版的 PROC POWER 一些新增的功能。

(1) LOGISTIC REGRESSION
羅吉斯迴歸的 power analysis 在此已經可以進行了。使用前要先寫 logistic 告知 SAS 要跑這種模型的 power analysis，並把要檢定的預測變數、勝算比以及反應變數 Y=1 的機率都寫上去，重點是要把檢定的預測變數的分配設定好（藍色部分）。

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;

(2) CONFIDENCE INTERVAL FOR ONE PROPORTION
二項變數的信賴區間也可以算了，但這結果其實跟用 PROC FREQ 算出來的結果是一樣的，所以如果你已經會用 PROC FREQ 來算二項變數的信賴區間的話，就不用多學這一段了。

proc power;
    onesamplefreq ci = Wilson
        alpha = 0.05
        proportion = 0.3
        halfwidth = 0.1
        ntotal = 70
        probwidth = .;
run;

(3) EQUIVALENCE TESTS FOR ONE PROPORTION
假設檢定中的比例檢定的 power analysis 也可以做了。重點就是要把紅色那段程式碼寫上去，其餘的參數設定跟以往都差不多。

proc power;
    onesamplefreq test = equiv_exact
        alpha = 0.05
        proportion = 0.35
        lower = 0.2
        upper = 0.4
        ntotal = 500
        power = .;
run;

(4) WILCOXON MANN-WHITNEY TEST FOR TWO INDEPENDENT GROUPS
無母數 Wilcoxon Mann-Whitney 檢定的 power analysis 可以用 twosamplewilcoxon 呼叫出。比較麻煩的是類別變數的分配要自己將每一層的 coding 和機率值打上去（藍色部分），注意後面機率值的總和要等於一。

proc power;
    twosamplewilcoxon
        alpha = 0.05
        vardist(’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;

最後原文有介紹一個SAS專門拿來做 Power analysis 的副產品Power and Sample Size 3.1，簡稱 PSS。不過這個軟體並沒有附在 SAS/STAT 裡面，需要另外安裝。我自己的 SAS 9.2 光碟片裡面好像也沒有這個附加軟體，研判應該是需要另外付費才有。我想應該也不會有很多人有 PSS，所以在此略過不提。把 PROC POWER 熟用應該就可以應付大部分的 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

Creating High-Quality Scatter Plots: An Old Story Told by the New SGSCATTER PROCEDURE

2011-05-07T06:03:00.003-04:00

原文載點：http://support.sas.com/resources/papers/proceedings10/057-2010.pdf

散佈圖算是統計圖表裡面一個很基本的呈現資料分布型態的表示方法，但是用 PROC GPLOT 畫出來的散佈圖不甚美觀。一位辛辛那提大學數學系的中國學生寫了這一篇技術文件，採用了 PROC SGSCATTER 程序，將單純的散佈圖的品質提昇到另一個境地。

資料來源是採用 SAS 裡面現成的資料檔 sashelp.cars，但作者只採用部分的內容：

ods html style=harvest;
data cars;
  set sashelp.cars;
  where make in ('Jeep' 'Chevrolet' 'Ford' 'Chrysler');
run;

變數的定義如下：

接下來所有的範例都會被 ODS 輸出成 html 的格式，採用了 harvest 的風格。如果你不喜歡 harvest 的風格，還有其他三種可以選：

範例一
以往用 PROC GPLOT 要畫 Y*(X1 X2) 的圖時，會產生兩張獨立的圖給你。現在用 PROC SGSCATTER 會自動給你合併成一張的圖。

proc sgscatter data=cars;
  plot invoice*(weight length); 
run;

範例二
如果想共用同一條軸，只要把 y= 和 x= 加上去就可以了。這樣可以讓散佈圖可視面積稍大一些。

proc sgscatter data=cars;
  compare y=invoice x=(weight length);
run;

範例三
以前介紹過的散佈矩陣圖。在 PROC SGSCATTER 裡面用 matrix 語法來執行。

proc sgscatter data=cars;
  matrix invoice weight length;
run;

範例四
用 rows= 和 columns= 來指定合併圖時的行列數目。

proc sgscatter data=cars;
  plot invoice*(weight length) / rows=2 columns=1;
run;

範例五
想要知道不同汽車製造商的散佈位置，首先在 plot 語法後面加上一個 group = 的選項，把汽車製造商的變數名稱 make 放入，這樣 PROC SGSCATTER 就會知道要針對那個變數做分組，並且自動替不同的製造商加上不同的顏色和點型。

proc sgscatter data=cars;
  plot MPG_city*weight / group=make;
  where make in ('Ford' 'Chrysler' 'Chevrolet');
  title 'Scatter Plot by Make';
run;

範例六
想要在散佈圖上加上廠商的名稱，則是用 datalabel= 這個選項來處理。

proc sql;
  create table cars2 as
select origin, make, mean(MSRP) as MSRP,
       mean(MPG_city) as MPG_city,
       mean(MPG_highway) as MPG_highway
from sashelp.cars
group by origin, make
order by origin, make;
quit;
proc sgscatter data=cars2;
  plot MSRP*MPG_highway / datalabel=make group=origin grid;
  title 'Averaged MSRP vs. Highway MPG for Car Makers by Origin';
  format MSRP dollar6.0;
  label MSRP='Manufacturer Suggested Retail Price' MPG_highway='Highway MPG';  
run;

範例七
若要加上一條迴歸線以及信賴區間，則用 reg= 的選項來畫。

proc sgscatter data=cars2;
  plot MSRP*MPG_highway / datalabel=make group=origin grid reg=(degree=2 clm nogroup);
  title 'Averaged MSRP vs. Highway MPG for Car Makers by Origin';
  title2 '-- with quadratic regression fitting and conf. intervals --';
  format MSRP dollar6.0;
  label MSRP='Manufacturer Suggested Retail Price' MPG_highway='Highway MPG';  
run;

範例八
要話95%的預測橢圓，則用 ellipse = 的選項來畫。

proc sgscatter data=cars2;
  compare y=MSRP x=(MPG_highway MPG_city) / group=origin ellipse=(alpha=0.05 type=predicted);
  title 'Averaged MSRP vs. Highway/City MPG for car makers by Origin';
  title2 '-- with 95% prediction ellipse --';
  format MSRP dollar6.0;
  label MSRP='Manufacturer Suggested Retail Price'
          MPG_highway='Highway MPG' MPG_city='CITY MPG';  
run;

範例九
回到散佈矩陣圖，如果想要在圖的對角線畫上每個變數的次數分配柱狀圖以及常態曲線，則可用 diagnoal = 的選項處理。

title 'Scatter Plot Matrix with Histograms and Normal Fitting Curves';
proc sgscatter data=cars;
  matrix invoice weight length / diagonal=(histogram normal);
run; quit;

範例十
最後，用 ODS GRAPHICS 設定圖型內部的參數，然後在 ODS HTML 後面加上路徑來讓生出來的圖存到指定位置去。

ods html gpath='C:\' style=harvest; 
ods graphics / reset=all width=12in height=6in border=off imagename='example' imagefmt=png;
proc sgscatter data=cars2;
  plot MSRP*(MPG_highway MPG_city) 
       / datalabel=make group=origin 
         grid reg=(degree=2 clm nogroup);
  title 'Averaged MSRP vs. Highway/City MPG for Car Makers by Origin';
  title2 '-- with quadratic regression fitting and conf. intervals --';
  format MSRP dollar6.0;
  label MSRP='Manufacturer Suggested Retail Price'
          MPG_highway='Highway MPG'
          MPG_city='City MPG';  
run; 
ods html close;

幾個重要的 ODS GRAPHICS 參數在此介紹一下：

WIDTH=, HEIGHT= ：設定圖型的長度和寬度
IMAGENAME=, IMAGEFMT= ：設定圖型的名稱和格式
BORDER=ON|OFF ：設定要不要畫圖型的邊界
RESET=ALL ：繪圖結束後重設所有參數

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author at:
Xiangxiang Meng
Department of Mathematical Science University of Cincinnati
mengxa@mail.uc.edu

A Grad Student ‘How-To’ Paper on Grant Prep: Preliminary Data, Power Analysis, and Presentation

2011-04-30T01:28:00.001-04:00

原文載點：http://support.sas.com/resources/papers/proceedings10/274-2010.pdf

SAS 語法百百種，身為研究生的你/妳，有沒有想過有什麼是一定要學起來並且用在報告或論文上面？這篇技術文件整合了許多你一定得用在報告和論文的SAS程式！

1. PRELIMINARY ANALYSIS I

Code1.1：將變數製作成聚集變數
這一步驟是用於一次要算幾百個變數的平均數的研究，讓使用者免去一個個輸入變數名稱的痛苦。但如果你只要算兩三個變數，那這段可以跳過。此程式所有的變數會被 PROC SQL 程序定義在一個名為 LIST 的巨集變數裡面。

PROC SQL noprint;
 select distinct NAME
 into :LIST
 separated by " "
 from work.cells_1;
quit;

Code1.2：計算平均數並用ODS OUTPUT輸出
此例的變數名稱就直接用 &LIST 把上面整合好的變數名稱一次叫進來。同樣地，如果你只有少量變數，那直接輸入即可。此程式會將算好的基本統計量另存在 basic_measures 的新資料集裡面。

ODS OUTPUT BASICMEASURES=work.basic_measures;
proc univariate data=work.&input.;
class day;
var &LIST.;
run;
ODS OUTPUT CLOSE;

Code1.3：重整結果以生成繪圖所需的資料
重新編輯資料集 basic_measures 以算出額外兩個新變數 High 和 Low。這兩個變數被定義為平均數正負一個變異數的大小。但我建議正規的作法應該要去算 95% 信賴區間比較好。這部份各位可以自行修改。

data work.basic_measures_all;
set work.basic_measures (rename=( LocValue=Mean));
if LOCMEASURE in ("Mean");
High=mean+varvalue;
Low=mean-varvalue;
Keep varname day mean high low;
run;

做出來的繪圖資料如下圖所示：

Code1.4：將變數和標籤製作成聚集變數
這一段是要將剛剛做出來的繪圖資料裡面的變數以及標籤製作成巨集變數以用於繪圖程式當中。如果你打算手動輸入，此段也可以省略。

proc sql noprint;
  select NAME as VAR into:VAR1- :VAR&SYSMAXLONG.
  from work.cells_1;
  quit;
proc sql noprint;
  select LABEL as LABEL into:LABEL1- :LABEL&SYSMAXLONG.
  from work.cells_1;
  quit;

Code1.5：繪圖並輸出
用 PROC SGPLOT 進行曲線圖的製作，並將結果用 ODS RTF 輸出到外部文件檔案。如果你想要製作別種圖，那就需要對 PROC SGPLOT 部分進行調整。

ods output;
ods RTF file="C:\documents and settings\desktop\GRAPHS_MEANS.RTF"  ;
%do j=1 %to &SQLOBS.;
data work.values;
set work.basic_measures_all;
if varname="&&VAR&J.";
run;
proc sgplot data=work.values noautolegend;
 xaxis type=discrete;
 series x=day y=mean;
 scatter x=day y=mean /
 markerattrs=(size=0)
 yerrorlower=low
 yerrorupper=high;
run;
%end;

製作出來的圖型如下：

Code1.6：用 PROC MIXED 程序進行分析並估計所需的組間差異
此例由於是在做重複觀測值的模型估計，所以用到 PROC MIXED 程序。這部份可視使用者的情況來調整。

ods proclabel "&&VAR&K. &&LABEL&K.";
ods output  Estimates=work.est  SolutionF=work.sol tests3=work.tests3 ;
proc mixed data=work.mixed1 COVTEST cl NOCLPRINT NOITPRINT ;
  class patient_id day_numeric;
  model &&var&k. =day_numeric /solution noint cl residual ;
  repeated day_numeric / subject=patient_id type=ar(1) rcorr r;
  estimate '-4 vs 0'      day_numeric   -1 1 0 0 0 0 0 0;
  estimate '-4 & 0 vs 4'  day_numeric   -0.5 -0.5 1 0 0 0 0 0;
  estimate '-4 & 0 vs 8'  day_numeric   -0.5 -0.5 0 1 0 0 0 0;
  estimate '-4 & 0 vs 12' day_numeric   -0.5 -0.5 0 0 1 0 0 0;
  estimate '-4 & 0 vs 15' day_numeric   -0.5 -0.5 0 0 0 1 0 0;
  estimate '-4 & 0 vs 16' day_numeric   -0.5 -0.5 0 0 0 0 1 0;
  estimate '-4 & 0 vs 22' day_numeric   -0.5 -0.5 0 0 0 0 0 1;
run;
quit;
ods output close;

Code1.7：格式化結果
這段只是單純把上面用 PROC MIXED 跑出來的結果整理一下。

data work.est;
length variable $32 variable_label $32 variable_subset $32;
set work.est;
variable="&&var&k.";
run;
PROC APPEND BASE=work.est_all_&source. DATA=work.est force;
RUN;

最後的表格就可以直接輸出到報告和論文上，如下所示：

2. PRELIMINARY ANALYSIS II

Code2.1：繪製個人曲線圖
注意此處的巨集變數都是從 Code1.4 生成出來的。如果跳過那步驟的人，就自己手動輸入變數名稱吧！

axis1 label=(angle=90);
symbol i=j repeat=100 w=5;
proc gplot data=work.mixed2;
   plot &&var&i.*week_numeric=patient_id/vaxis=axis1 nolegend;
run;
quit;

圖型如下所示：

Code2.2：PROC SGPANAL
這段是在畫另一種圖，但我覺得不是很重要，可以省略。

proc sgpanel data=work.data_transpose_2;
  where _NAME_ in ("&&VAR&I.");
  panelby _LABEL_/novarname ;
  colaxis type=discrete values=(-7,0,1,3,7,10,14,21,28);
  rowaxis label="log (Count)";
  series x=day y=_1 /markers LINEATTRS = (THICKNESS = 3);
  series x=day y=_2 /markers LINEATTRS = (THICKNESS = 3);
  series x=day y=_3 /markers LINEATTRS = (THICKNESS = 3);
  series x=day y=_4 /markers LINEATTRS = (THICKNESS = 3);
  series x=day y=_5 /markers LINEATTRS = (THICKNESS = 3);
  series x=day y=_6 /markers LINEATTRS = (THICKNESS = 3);
run;
quit;

圖型如下：

Code2.3：用 PROC MIXED 估計參數並用 ODS OUTPUT 輸出共變異數矩陣裡的參數

ODS OUTPUT "Covariance Parameter Estimates"=work.parameters;

proc mixed data=mixed2 COVTEST cl NOCLPRINT NOITPRINT ;
  class patient_id ;
  model &&var&i. =day_numeric /solution cl ;
  random intercept / subject=patient_id ;
  run;
  quit;
ODS OUTPUT CLOSE;

Code2.4：整理共變異數並計算 ICC
這段程式看起來非常冗長，但其實簡單來講就是把 Code2.3 做出來的資料整理後算出 ICC（藍色部分），不過原作者還算了很多其他的估計值。如果你只是想要求出 ICC，那綠色的程式碼可以不用寫。

data work.parameters_between;
 set work.parameters;
if covparm="Intercept" then do;
   Between_Subject_Variance=Estimate;
   Between_p_value=ProbZ; 
   Between_Lower=Lower; 
   Between_Upper=Upper; 
 end;
 if covparm="Intercept" ;
 keep Between_Subject_Variance Between_P_Value Between_Lower Between_Upper;
 run;
 data work.parameters_within;
 set work.parameters;
 if covparm="Residual" then do;
   Within_Subject_Variance=Estimate;
   Within_p_value=ProbZ; 
   Within_Lower=Lower; 
   Within_Upper=Upper; 
   end;
 if covparm="Residual" ;
 keep Within_Subject_Variance WITHIN_P_VALUE Within_Lower Within_Upper;
 run;
 data work.parameters_all;
 merge work.parameters_between work.parameters_within;
 ICC=between_subject_variance/(within_subject_variance+between_subject_variance); 
 variable="&&var&i."; variable_label="&&label&i."; data_source="&source.";
 run;
 quit;
 proc print data=work.parameters_all noobs;
 run;

最後用 PROC PRINT 列印出來的結果如下：

3. POWER SIMULATION
最後重頭戲是 power analysis。不過此處用到模擬的原因是因為，作者想要做 power analysis 的對象是 ICC 和組間變異數，而非我們一般常看到的假設檢定。而由於 ICC 和組間變異數的 power 計算並沒有封閉形式(close form)的解，所以必須用模擬的方式來算出。
Code3.1：生成種子(seed)
用亂數去生成種子讓每一次的模擬都能隨機。種子的數量取得於母體大小乘上需要模擬的次數。本例的種子數目高達一千五百萬。

%MACRO INITIALIZE ();
%Global Initial;  
%Let Initial=0;
data Initialseed (keep=id InitialSeed);
      do I = 1 to 15000000;   
            InitialSeed =int(10000*(ranuni(0)));
            id=I;
            output;
      end;
run;
data SIM.InitialSeed;
 set work.InitialSeed;
run;
%MEND INITIALIZE;

Code3.2：製作依序取得種子的巨集程式

%LET Initial=35;
%MACRO SEED ();
%let Initial=%EVAL(&Initial. +1);
 data work.seed ;
  set sim.initialseed (Firstobs=&Initial. obs=&Initial.);;
 run;
 data _null_;
  set work.seed;
  call symputx ('InitialSeed',InitialSeed);
 run;
%put Seed is # &Initial. and is valued &InitialSeed. ;
%MEND SEED;

Code3.3：模擬資料的巨集程式
原文中 Code3.3~Code3.9 都是在講資料的生成。我將他們放在一起。簡單來講，先生成 baseline 資料(work.test)，再生成其他時間點的資料(work.test_2)，然後用 PROC TRANSPOSE 把生成好的結果轉置(直的變成橫的)，接著就是用 PROC MIXED 不斷重複地去估計共變異數矩陣的參數，把這重複估計出來的參數存成一個新檔案(work.mixed)，然後去算有多少顯著的比例。最後把結果整理成報表。

%MACRO TEST(mean=,between_var=,within_var=,simulation=,population=,timecount=);

data work.test;
*first loop defines the simulation number*;
do i= 1 to &simulation.;
*second loop defines the number of people in the study*;
 do j=1 to &population.;
   TIME_0=&mean.+sqrt(&between_var.)*rannor(&initialseed.);
  Simulation_ID=i;
  Person_ID=J; 
  output;
 end;
end;
drop I j;
run;

data work.test_2;
set work.test;
array time (*) TIME_1-TIME_&timecount.;
  do k=1 to &timecount.;
   Time[k]=TIME_0 + sqrt(&within_var.)*rannor(&initialseed.);
 end;
  output;
drop k;
run;

proc transpose data=work.test_2 out=work.test_3 prefix=value;
 by simulation_ID person_id;
run;
data work.test_3;
set work.test_3;
 TIME= input(scan(_NAME_,2,'_'),2.);
 TIME2=TIME;
run;

ODS LISTING CLOSE;
ODS OUTPUT COVPARMS=work.mixed;
proc mixed data=work.test_3 covtest;
  by simulation_ID;
  model value1=time;
  random intercept/ sub=person_id;
run;
ODS OUTPUT CLOSE;

data work.mixed_1;
set work.mixed;
if Covparm="Intercept";
if probz<0.1 then P_point1=1;
if probz<0.05 then P_point05=1;
if probz<0.01 then P_point01=1;
if probz<0.001 then P_point001=1;
if probz=. then do;
 p_point1=.;
 p_point05=.;
 p_point01=.;
 p_point001=.;
  end;
run;

ODS OUTPUT onewayfreqs=work.mixed_freqs_1;
ODS LISTING;
proc freq data=work.mixed_1 ;
table p_point1 p_point05 p_point01  p_point001/missing;
run;
ODS LISTING CLOSE;
ODS OUTPUT CLOSE;

proc append Base=work._mixed_freqs_all data=work._mixed_freqs_all_&initial. FORCE;
run;
%MEND;

Code3.10：開始模擬
最後就是呼叫前面寫好的巨集程式，輸入不同的組間變異數和組內變異數的數值，定義總共要模擬的次數、樣本大小、時間點數量以及平均數。

%test(mean=-0.2740, between_var=0.01,within_var=0.67,simulation=1000,population=20,timecount=12);
%test(mean=-0.2740, between_var=0.02,within_var=0.66,simulation=1000,population=20,timecount=12);
%test(mean=-0.2740, between_var=0.03,within_var=0.65,simulation=1000,population=20,timecount=12);
%test(mean=-0.2740, between_var=0.04,within_var=0.64,simulation=1000,population=20,timecount=12);
%test(mean=-0.2740, between_var=0.05,within_var=0.63,simulation=1000,population=20,timecount=12);
%test(mean=-0.2740, between_var=0.06,within_var=0.62,simulation=1000,population=20,timecount=12);
%test(mean=-0.2740, between_var=0.50,within_var=0.16,simulation=1000,population=20,timecount=12);

Code3.11：整理繪圖資料
模擬完後的資料存在 work._mixed_freqs_all 裡面，把不同組間變異數相對應的 ICC 放進去，並且把變數的標籤寫好。

data work.plots ;
set work._mixed_freqs_all;
 if between_var=  0.01 then ICC =  0.01 ;
 if between_var=  0.06 then ICC =  0.09 ;
 if between_var=  0.11 then ICC =  0.16 ;
 .……
label percent="Power" population="Number of Subjects" timecount="Number of Time  
            Points:";
run;

Code3.12：繪圖

proc sort data=work.plots;
 by population;
run;
symbol i=j width=5;
ODS LISTING;
proc gplot data=work.plots;
 by population ;
 where table ="P_point01" and mixed=2;
 plot percent*ICC=timecount;
run;
quit;

最後圖型如下：

CONTACT INFORMATION

Your comments and questions are valued and encouraged. Contact the author at:
Elisa L. Priest, MPH
elisapriest@hotmail.com

Improving Your Statistical Consulting Prowess by Adding R to Your SAS® Repertoire

2011-04-09T15:56:00.002-04:00

原文載點：http://support.sas.com/resources/papers/proceedings10/290-2010.pdf

我不清楚有多少人和我一樣，用SAS分析資料，但是用R來做圖。其實我一開始堅持資料分析和做圖都在SAS內一次完成，但後來我終於明白SAS/GRAPH要做出等同用R畫出來的圖形品質，需要花費更多時間，於是我終於將做圖的工作完全交付給R，也展開了我要一直切換兩個軟體視窗的日子。但我一直相信可以在SAS裡面呼叫R來完成畫圖的工作，終於讓我在SAS Global Forum 2010上面發現到一篇技術文件，在此跟大家分享。

首先要讓SAS執行R的程式所需的第一件事情，是要去SAS官網下載SAS/IML Studio 3.2，這是一套龐大但免費的軟體，總容量是整整1GB，下載前需要去SAS官網註冊帳號，註冊好就可以直接下載了。網址是 http://support.sas.com/rnd/app/da/workshop.html

原文內雖然是說用SAS/IML Workshop 2.2，但我實際去官網查詢的結果，發現只有 SAS/IML Workshop 2.1的連結，而且沒有辦法用在最新的SAS 9.2版上面。所以請直接去下載SAS/IML Studio 3.2版即可。經過我的測試是可以執行以下的範例的。

下載完畢並安裝好後，就可以在SAS裡面寫R的程式並呼叫R來執行，而結果會顯示在SAS的視窗裡面。不過根據本人的經驗，事情沒有那麼簡單，因為當我裝好並執行一段官網上的範例時，視窗一直彈跳出軟體找不到R的錯誤訊息，經過與SAS官網客服人員的聯繫，才知道SAS/IML Studio 3.2目前只能與R的2.9.1和2.11.1這兩個版本相容。如果你已經安裝了最新的R 2.12.2版，目前唯一的解決方法就是移除最新版然後降級安裝成2.11.1版。如果你有某些特殊的函式比需要在R 2.12.2版才能執行，則必須再次安裝最新版的R才行。

所有的R程式，必須要用"submit /R;"開頭，用"endsubmit;"結尾。而以往在SAS裡面常用的注解符號"*"也必須換成井字符號"#"。在進行資料分析前，得先將SAS資料轉成R可以用的資料格式，可利用下面這段程式碼進行轉檔；

run ExportDataSetToR("libname.dataname", "R_dataname");

其中libname是寫SAS library的名稱，而dataname自然就是SAS資料的名稱。R_dataname則是填入資料集在R裡面所用的名稱。作者建議儘量少用底線當做用在R裡面的資料名稱，因為有可能會出問題。此外，如果是用自定的library，則必須照舊先在一開始用libname語法宣告自定library的路徑。底下的範例因為都是用SAS內部已經建好的Sashelp，所以沒有必要在一開始用libname來宣告。在此需特別注意這一點。以下用幾個範例來說明：

範例一；表格
這個範例是教如何在SAS裡面使用一些R的指令來整理表格，雖然我覺得SAS在這方面就比R強太多了，但這就是一個範例，程式如下；

run ExportDataSetToR("Sashelp.cars", "cars");
submit /R;
library(Hmisc, T); 
sink("c/sugi2010/car_summaries.txt");
car.sum <- summary.formula(DriveTrain ~ Type + Origin + MPG_City + MPG_Highway + Weight + Wheelbase + Length, data = cars, method = "reverse", overall = FALSE, test = TRUE,digits = 3, prmsd=TRUE, exclude1 = FALSE, long = TRUE);
print(car.sum, test = TRUE,digits = 3, prmsd=TRUE, exclude1 = FALSE, long = TRUE);
sink();
endsubmit;

為了顧及一些人不懂R的語法，所以在此稍微解釋一下，之後的範例會省略對R的介紹，有需要的人請自行參照R的線上手冊。

第一行就是作資料轉檔的語法，剛剛已經介紹過了。本例是拿Sashelp裡面的一個cars的檔案來做表格。"submit /R;"之後開始就全部都是R的語法了，裡面不會有任何SAS程式，直到最後用"endsubmit;"才是回到SAS語法。library() 是R裡面呼叫某個package的語法，Hmisc就是專門拿來做summary table的package。sink("....txt") 是先定義接下來所有的報表是要輸出到某txt文件裡面，記得結束要補上"sink();"在最後面才會完成存檔動作。summary.formula()就是製作的函數，裡面要放很多參數，在此就不一一介紹了，總之做完的結果是放在car.sum這個物件裡面。而print()是把car.sum這個物件印進剛剛用sink函數建立的純文字檔。而那個純文字檔打開後如下所示；

R繪製出來的圖形相當漂亮，但表格就很陽春，非常像古早PE2打出來的效果。因此這個時候就是SAS展現威能的時候了。用一個data step把那個純文字檔叫進來；

data cars2;
  infile "("c/sugi2010/car_summaries.txt"" dsd dlm="|"   
         lrecl=1024 firstobs=5 truncover;
  format junk $20. label $34. all $34. front $35. rear $35. test $35. ;
  input junk $ label $ all $ front $ rear $ test $ ; 
  if substr(junk,1,1)="+" then delete;
  drop junk;
run;

其中：

用infile把純文字檔呼叫進SAS，用dlm="|"將每一行分段，用lrecl來定義整行最大長度，並用firstobs指定SAS從第五行開始讀。
用format把剛剛用"|"將每一行切成數段，所以每一段就變成一個變數，因此需要定義變數的格式。在此將所有變數定義成文字格式，因為接下來不用再做分析了，所以只要能列印出及可。
用input設定各變數的名稱，這些名稱要和format裡面所寫的變數名稱一樣。
用junk變數把不必要的隔線刪除掉。

把不要的內容清乾淨後，用ods和proc report把結果印出來；

ods rtf file="c:\sugi2010\sas_r\car_summaries2_09oct20b.rtf" ;
title "Summaries by treatment group";
proc report data=cars2 nowindows;
  columns label all front rear test;
run; ods rtf close;

結果如下；

範例二；繪圖
在SAS裡面用R畫圖是我比較關心的，但和上面這個範例比起來並沒有太多不同。程式如下；

run ExportDataSetToR("Sashelp.CARS", "cars");
submit /R;
library(Hmisc, T);  
## Create data frame – cars2 – dropping cars with 3, missing, and 5 cylinders; 
cars2 <- cars[cars$Cylinders>3 & !is.na(cars$Cylinders) & cars$Cylinders != 5,] 
## place cars2 in search path 
attach(cars2); 
## Create summary data sets for Highway and city gas mileage MPG by number of 
cylinders and origin.  ; 
## Bootstrap resampling used to calculate 95% confidence intervals of mean; 
mean.hw <- summarize(MPG_Highway, llist(Cylinders, Origin), smean.cl.boot); 
mean.city <- summarize(MPG_City, llist(Cylinders, Origin), smean.cl.boot); 
## rename MPG_Highway and MPG_City to common name:  mean.mpg.; 
mean.hw2 <- upData(mean.hw, rename=list(MPG_Highway="mean.mpg")); 
mean.city2 <- upData(mean.city, rename=list(MPG_City="mean.mpg")); 
##create new variable type to label type of mean.mpg; 
mean.hw2$type[mean.hw2$mean.mpg>0]<- "Highway" 
mean.city2$type[mean.city2$mean.mpg>0]<- "City" 
## Vertical combine of the 2 data frames – new dataframe is named all; 
all <- rbind(mean.hw2, mean.city2);
##  Use Hmisc function xYplot to plot mean gas mileage and 95% confidence intervals; 
## by cylinders with separate panels for each of 3 origins.  ; 
## groups – separate curves by type within each panel.; 
## col = c(1,2) sets colors for 2 groups (1=black, 2=red). Lty sets line types; 
## layout=c(3,1)  determines number of rows and columns – c(3,1) is 3 columns and 1 row; 
xYplot(Cbind(mean.mpg, Lower, Upper) ~ Cylinders | Origin, data=all, groups=type, 
col=c(1,2), lty=c(1,2), layout=c(3,1), type='o'); 
endsubmit;

基本上你只要把前兩行和最後一行先在SAS裡面寫好，然後在第三行前面貼上R的程式即可。不過上面這段用R繪圖的語法，其實在SAS裡面是可以直接做的，下面是兩個軟體畫出來的圖形的比較；

左圖是用R畫出來的圖，右邊是用SAS畫的。乍看之下好像右圖比較好，但那是因為新版的SAS剛好有PROC SGPANEL這個程序把所有繪圖參數都最佳化了，如果遇到PROC SGPANEL沒有支援的圖形，那就不要太指望SAS能給你甚麼好圖了。而根據本人實際測試的結果，其實還是R做的圖比較好看，所以就大膽地在SAS/IML Studio底下執行R的程式吧！

CONTACT INFORMATION
Name: Jeff Gossett
University of Arkansas for Medical Sciences. Department of Pediatrics.
One Children's Way
Little Rock, AR 72202
Work Phone: 501-364-6631
Fax: 501-364-1431
E-mail: GossettJeffreyM@uams.edu
Web: http://www.arpediatrics.org/research/biostatistics

So Many Variables; Too Many Labels, Moving Labels from One Variable to Another

2011-04-03T01:23:00.001-04:00

原文載點：http://support.sas.com/resources/papers/proceedings10/047-2010.pdf

在 SAS 資料處理的過程中，使用者會發現，當複製一個變數並給予新的變數名稱時，舊變數的格式雖然會一併複製到新變數裡面，但舊變數的 label 卻沒有辦法複製過去。若只有一兩個變數時還可以用手動的方式在 data step 裡面重打，但若變數過多時，這就會成為一個相當棘手的問題。John Ladds 提供了一個簡易的方式來解決這個 label 複製的問題。

以下面這個資料集當作範例：

proc format;
   value OriginalOrderAgree
       1 = 'Strongly agree'
2 = 'Somewhat agree'
3 = 'Neither agree nor disagree'
4 = 'Somewhat disagree'
5 = 'Strongly disagree'
6 = 'Not applicable'
      97 = "Don't know"
99 = 'Not stated';
   value ReverseOrderAgree
       5 = 'Strongly agree'
4 = 'Somewhat agree'
3 = 'Neither agree nor disagree'
2 = 'Somewhat disagree'
1 = 'Strongly disagree'
       other='Missing';
run;
** Create some test data;
data work.theData;
   input @01 id    $10.
         @11 Q001    2.
         @13 Q002    2.
         @15 Q003    2.
         @17 Q004    2.
         @19 Q005    2.
         @21 Q006    2.
         @23 Q007    2.
         @25 Q008    2.
         @27 Q009    2.
         @29 Q010    2.;
   label Q001 = "Have material's and equipment I need";
   label Q002 = "Material and tools avail in lang choice";
   label Q003 = "I use the language of my choice";
   label Q004 = "My jobs fits my interests.";
   label Q005 = "I have support at work.";
   label Q006 = "I am satisfied, current work arrangement";
   label Q007 = "I can claim overtime";
   label Q008 = "Overall, I like my job.";
   label Q009 = "I get satisfaction from my work.";

   label Q010 = "I know how my work contributes";
format Q001-Q010 OriginalOrderAgree.;
datalines4;
1234567890 1 2 3 4 5 6 1 2 3 4;
1234567892 2 3 4 5 1 1 2 3 4 5;
1234567893 3 4 5 1 2 2 2 3 5 1;
1234567894 4 5 1 2 3 3 3 3 1 2;
1234567895 5 1 2 3 4 4 4 4 2 3;
1234567896 1 2 3 4 5 5 5 5 3 4;
1234567897 2 3 4 5 1 1 1 1 4 5;
1234567898 3 4 5 1 2 2 2 2 5 1;
1234567899 4 598 2 3 3 3 3 1 2;
123456781099 199 3 4 4 4 5 2 3;
;;;;
run;

根據這個資料集，裡面有十個變數，從Q001~Q010，每個變數為數據「1,2,3,4,5」。現在假設要製造新變數R_Q001~R_Q010，而這些新變數是將就舊變數的數據倒轉過來，也就是 1 變成 5，2 變成 4，...，5 變成 1。於是作者利用下面這個 macro 來處理這段：

%macro reverseValueOrder(Q=,Qfmt=);
   select(&Q.);
     when(1)   R_&Q.=5;
     when(2)   R_&Q.=4;
     when(3)   R_&Q.=3;
     when(4)   R_&Q.=2;
     when(5)   R_&Q.=1;
     otherwise R_&Q.=.;
   end;
   ** Transfer the questions label to the new variable;
   if eof then call execute("label R_&Q.="||'"'||strip(vlabel(&Q.))||'";');
   ** Add the new value labels to the new variable;
   format R_&Q. &Qfmt..;
%mend reverseValueOrder;

其中，「select(&Q);.....end;」主要就是在做變數倒轉的工作，而下面有一行「if...then...」搭配「call execute」便是在處理 label 複製的動作。裡面的 vlabel 函數的作用是將舊變數裡面的 label 給列出來，放在 strip 函式後面主要是移除不必要的空格，並確保 vlabel 函式叫出來的 label 是字串形式。

接著就可以在新的 data step 裡面一邊進行數據倒轉的工作，一邊進行搬移 label 的工作，如下所示：

data work.thedata;
   set work.theData end=eof;
** Open the post DATA STEP Processing;
if eof then call execute('data &syslast.; set &syslast.;'); 
**A. My Job World                              ;
   %reverseValueOrder(Q=Q001, Qfmt=ReverseOrderAgree);
   %reverseValueOrder(Q=Q002, Qfmt=ReverseOrderAgree);
   %reverseValueOrder(Q=Q003, Qfmt=ReverseOrderAgree);
   %reverseValueOrder(Q=Q004, Qfmt=ReverseOrderAgree);
   %reverseValueOrder(Q=Q005, Qfmt=ReverseOrderAgree);
   %reverseValueOrder(Q=Q006, Qfmt=ReverseOrderAgree);
   %reverseValueOrder(Q=Q007, Qfmt=ReverseOrderAgree);
   %reverseValueOrder(Q=Q008, Qfmt=ReverseOrderAgree);
   %reverseValueOrder(Q=Q009, Qfmt=ReverseOrderAgree);
   %reverseValueOrder(Q=Q010, Qfmt=ReverseOrderAgree);
;... ... ... ... ... ... ... ... ... ... ... **  
** Close the post DATA STEP Processing;
if eof then call execute ('run;');
run;

特別提醒的地方有三個。一是 set 後面一定要加上「end=eof」字樣。二是執行 macro 前加上「if eof then call execute('data &syslast.; set &syslast.;'); 」，目的是確保程式執行的資料是當前資料。所以記得在運作此程式之前，別去跑其他不相干的資料。三是結束前加上「if eof then call execute ('run;');」，這樣就大功告成了。

CONTACT INFORMATION
Name: John Ladds
Enterprise: Statistics Canada – SDD SAS Technology Centre
Address: 14-R-90 R.H. Coats Building,
100 Tunney’s Pasture Driveway
Ottawa, Ontario, Canada
K1A 0T6
Phone: 613 951-1767
E-mail: John.Ladds@statcan.gc.ca
Web: www.statcan.ca

ROC Hard? No, ROC Made Easy!

2011-03-15T23:12:00.005-04:00

原文載點: http://support.sas.com/resources/papers/proceedings10/222-2010.pdf

SAS 裡面並沒有內建繪製 ROC 曲線的程序和語法，這和目前其他統計軟體都已經內建繪製 ROC 曲線的功能的情況比較起來，顯得相當奇怪。使用者在 SAS 的環境底下，必須先用 PROC LOGISTIC 程序將模型估計出來，再將一些報表另存新的檔案，然後才能根據這些新的資料來繪製 ROC 曲線並且計算 AUC 面積。SAS 官網有一段教學，詳情可參照 http://support.sas.com/kb/25/018.html，但繪製出來的圖形不甚美觀。Kriss Harris 在 SAS Global Forum 2010 發表了一篇技術文件，提供了一個相當實用的巨集程式 %ROC_CUTOFF 讓使用者可以在 SAS 裡面迅速且美觀地畫出 ROC 曲線。

首先，利用一段模擬的資料來當作本篇的範例：

data blood_for_roc;
CALL STREAMINIT(400);  
do i = 1 to 40;
Sputum_Eosinophils = RAND('NORMAL',3, 1); 
cht = 1;
output;
end;
do i = 41 to 80;
Sputum_Eosinophils = RAND('NORMAL',10, 5);  
cht = 2;
output;
end;
run;

%ROC_CUTOFF 的格式如下所示：

%macro roc_cutoff(DATASET = blood_for_roc ,
                  OUTCOME = cht , 
                  OUTCOME_LEV = 2 /* Number */ , 
                  XVAR = Sputum_Eosinophils,
                  XVAR_LABEL = Sputum Eosinophils,
                  AUC_TABLE = Y,
                  CONFUSION = N,
                  CUTOFF = );

每個參數的定義如下：

DATASET = SAS 資料集名稱
OUTCOME = 類別型的依變數名稱
OUTCOME_LEV = 依變數內的類別數
XVAR = 連續型的應變數名稱
XVAR_LABEL = 應變數的標籤
AUC_TABLE = 要不要顯示 AUC 表 (Y/N)
CONFUSION = 要不要顯示 confusion 表 (Y/N)
CUTOFF = 指定 cutoff 的數字(可留空白不填)

以模擬的資料來跑 ROC：

%macro roc_cutoff(DATASET = blood_for_roc ,
                  OUTCOME = cht , 
                  OUTCOME_LEV = 2 /* Number */ , 
                  XVAR = Sputum_Eosinophils,
                  XVAR_LABEL = Sputum Eosinophils,
                  AUC_TABLE = Y,
                  CONFUSION = N,
                  CUTOFF = );

結果如下所示：

若將CUTOFF 設定成 5，結果如下所示：

若將 AUC_TABLE 和 CONFUSION 都設定為 N，則一開頭的兩個表格就不會被秀出來：

原始碼如下：

%macro roc_cutoff(DATASET = blood_for_roc ,
                  OUTCOME = cht , 
                  OUTCOME_LEV = 2 /* Number */ , 
                  XVAR = Sputum_Eosinophils,
                  XVAR_LABEL = Sputum Eosinophils,
                  AUC_TABLE = Y,
                  CONFUSION = N,
                  CUTOFF = );
;*****************************************************************%
%* GETTING ROC OUTPUT FROM MODEL *;
;*****************************************************************%
%LET OUTCOME_LEV = %SYSFUNC(TRANWRD(&OUTCOME_LEV,%STR(%"),));
%LET OUTCOME_LEV = %SYSFUNC(TRANWRD(&OUTCOME_LEV,%STR(%'),)); 
ods output Association = roc_association (WHERE = (UPCASE(LABEL2) = "C"));
proc logistic data=&dataset outest = outest;
model &OUTCOME. (event = "&OUTCOME_LEV.") = &XVAR./EXPB outroc=ROCData;
output out = roctest p = p ;
run;
DATA ROCData;
SET ROCData; 
INDEX =1;
RUN;
DATA Outest (KEEP=INTERCEPT &XVAR. INDEX);
SET Outest; 
INDEX =1;
RUN;
/* Calculating X value */
DATA _merge;
MERGE ROCData Outest;
BY INDEX;
SPEC = 1-_1MSPEC_;
YI = (_SENSIT_ + SPEC) - 1;
N = _POS_+ _NEG_+ _FALPOS_+ _FALNEG_;
TA = (_POS_ + _NEG_)/N;
X_VALUE = (LOG(_PROB_/(1-_PROB_))-INTERCEPT)/&XVAR.;
IF (_NEG_+_FALNEG_) NE 0 THEN DO;
NPV = _NEG_/(_NEG_ +_FALNEG_);
END;
IF (_POS_+_FALPOS_) NE 0 THEN DO;
PPV = _POS_/(_POS_ +_FALPOS_);
END;
IF
((_POS_+_FALPOS_)*(_POS_+_FALNEG_)*(_NEG_+_FALPOS_)*(_NEG_+_FALNEG_)) NE 0
THEN DO;
MCC = ((_POS_*_NEG_)-(_FALPOS_*_FALNEG_))/
SQRT(((_POS_+_FALPOS_)*(_POS_+_FALNEG_)*(_NEG_+_FALPOS_)*(_NEG_+_FALNEG_)));
END;
RUN;
/* Merging in the original data, so that this can be plotted to */
proc sql;
create table merge2 as
select a.*, b.&XVAR. as XVAR1, b.&OUTCOME., b.p,
case when b.&OUTCOME. = &OUTCOME_LEV. then 3
else 4 end as Index1
from _merge as a right join Roctest as b
on a._PROB_ = b.p ;
quit;
data merge3;
set merge2;
string_check = anyalpha(&OUTCOME.);
if string_check = 0 then &OUTCOME._char = strip(put(&OUTCOME., best.))
;
else &OUTCOME. = &OUTCOME._char;
run;
;*****************************************************************%
%* CALCULATING X Value CUT-OFFS for Particular Specificities and 
Sensitivities and Vice Versa *;
;*****************************************************************%
proc sql;
create table One_Hundred_PC_Sens as
select max(X_VALUE) as max_x format 10.2
from _merge
where _SENSIT_ = 1;
quit;
proc sql;
create table sens_spec as
select *, abs(_SENSIT_ - SPEC) as abs_sens_spec, min(calculated abs_sens_spec) as min
from _merge;
quit;

proc sql;
create table sens_spec_intersection as
select X_VALUE format 10.2
from sens_spec 
where abs_sens_spec = min;
quit;
proc sql;
create table One_Hundred_PC_Spec as
select min(X_VALUE) as min_x format 10.2
from _merge
where spec = 1;
quit;
/* Fetching the values */
%let dsid=%sysfunc(open(Roc_association,i));
%let num_AUC=%sysfunc(varnum(&dsid,nValue2));
%let rc=%sysfunc(fetch(&dsid,1));
%let AUC=%sysfunc(getvarn(&dsid,&num_AUC));
%let AUC_format=%sysfunc(putn(&AUC,bestd10.2));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(One_Hundred_PC_Sens,i));
%let num_Sens=%sysfunc(varnum(&dsid,max_x));
%let rc=%sysfunc(fetch(&dsid,1));
%let Sens=%sysfunc(getvarn(&dsid,&num_Sens));
%let Sens_format=%sysfunc(putn(&Sens,bestd10.2));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(One_Hundred_PC_Spec,i));
%let num_Spec=%sysfunc(varnum(&dsid,min_x));
%let rc=%sysfunc(fetch(&dsid,1));
%let Spec=%sysfunc(getvarn(&dsid,&num_Spec));
%let Spec_format=%sysfunc(putn(&Spec,bestd10.2));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(sens_spec_intersection,i));
%let num_section=%sysfunc(varnum(&dsid,X_VALUE));
%let rc=%sysfunc(fetch(&dsid,1));
%let section=%sysfunc(getvarn(&dsid,&num_section));
%let rc=%sysfunc(close(&dsid));
%if &CUTOFF = %NRSTR() %then %do;
%let section2 = §ion;
%let section_format=%sysfunc(putn(§ion2,bestd10.2));
%let rc=%sysfunc(close(&dsid));
%end;
%else %do;
%let section2 = &CUTOFF;
%let section_format = &CUTOFF;
%end;
;*****************************************************************%
%* 2*2 Tables *;
;*****************************************************************%
data freq_input;
set merge3;
if XVAR1 ne .;
if XVAR1 <=  §ion2 then new_result = 1;
else new_result = 0;
run;
ods output  CrossTabFreqs =  CrossTabFreqs;
proc freq data = freq_input;
tables new_result * &OUTCOME._char / NOROW NOCOL;
run;
data crossTabFreqs_ref;
set CrossTabFreqs;
if &OUTCOME._char ne "";
if new_result ne .;
run;

data TP;
set CrossTabFreqs_ref;
if new_result = 0 and &OUTCOME._char = &OUTCOME_LEV.;
run;
data FP;
set CrossTabFreqs_ref;
if new_result = 0 and &OUTCOME._char ne &OUTCOME_LEV.;
run;
data FN;
set CrossTabFreqs_ref;
if new_result = 1 and &OUTCOME._char = &OUTCOME_LEV.;
run;
data TN;
set CrossTabFreqs_ref;
if new_result = 1 and &OUTCOME._char ne &OUTCOME_LEV.;
run;
/* Fetching the values */
%let dsid=%sysfunc(open(TP,i));
%let num_TP=%sysfunc(varnum(&dsid,Frequency));
%let rc=%sysfunc(fetch(&dsid,1));
%let TP=%sysfunc(getvarn(&dsid,&num_TP));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(FP,i));
%let num_FP=%sysfunc(varnum(&dsid,Frequency));
%let rc=%sysfunc(fetch(&dsid,1));
%let FP=%sysfunc(getvarn(&dsid,&num_FP));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(FN,i));
%let num_FN=%sysfunc(varnum(&dsid,Frequency));
%let rc=%sysfunc(fetch(&dsid,1));
%let FN=%sysfunc(getvarn(&dsid,&num_FN));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(TN,i));
%let num_TN=%sysfunc(varnum(&dsid,Frequency));
%let rc=%sysfunc(fetch(&dsid,1));
%let TN=%sysfunc(getvarn(&dsid,&num_TN));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(TP,i));
%let num_TP_per=%sysfunc(varnum(&dsid,Percent));
%let rc=%sysfunc(fetch(&dsid,1));
%let TP_per=%sysfunc(getvarn(&dsid,&num_TP_per));
%let TP_format=%sysfunc(putn(&TP_per,best3.0));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(FP,i));
%let num_FP_per=%sysfunc(varnum(&dsid,Percent));
%let rc=%sysfunc(fetch(&dsid,1));
%let FP_per=%sysfunc(getvarn(&dsid,&num_FP_per));
%let FP_format=%sysfunc(putn(&FP_per,best3.0));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(FN,i));
%let num_FN_per=%sysfunc(varnum(&dsid,Percent));
%let rc=%sysfunc(fetch(&dsid,1));
%let FN_per=%sysfunc(getvarn(&dsid,&num_FN_per));
%let FN_format=%sysfunc(putn(&FN_per,best3.0));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(TN,i));
%let num_TN_per=%sysfunc(varnum(&dsid,Percent));
%let rc=%sysfunc(fetch(&dsid,1));
%let TN_per=%sysfunc(getvarn(&dsid,&num_TN_per));
%let TN_format=%sysfunc(putn(&TN_per,best3.0));
%let rc=%sysfunc(close(&dsid));
/* Selecting distinct OUTCOME values */

proc sql;
create table distinct_values as
select distinct &OUTCOME.
from &dataset;
quit; 
data LEVEL NOTLEVEL;
  set distinct_values;
   if  &OUTCOME. =  UPCASE(&OUTCOME_LEV.) then output LEVEL;
else OUTPUT NOTLEVEL;   
run;
/* Turning it to character if it's not */
data level_char;
set level;
string_check = anyalpha(&OUTCOME.);
if string_check = 0 then &OUTCOME._char2 = strip(put(&OUTCOME., best.))
;
else &OUTCOME. = &OUTCOME._char2;
run;
data notlevel_char;
set notlevel;
string_check = anyalpha(&OUTCOME.);
if string_check = 0 then &OUTCOME._char2 = strip(put(&OUTCOME., best.))
;
else &OUTCOME. = &OUTCOME._char2;
run;
/* Fetching the distinct levels */
%let dsid=%sysfunc(open(level_char,i));
%let num_level=%sysfunc(varnum(&dsid,&OUTCOME._char2));
%let rc=%sysfunc(fetch(&dsid,1));
%let level=%sysfunc(getvarc(&dsid,&num_level));
%let rc=%sysfunc(close(&dsid));
%let dsid=%sysfunc(open(notlevel_char,i));
%let num_notlevel=%sysfunc(varnum(&dsid,&OUTCOME._char2));
%let rc=%sysfunc(fetch(&dsid,1));
%let notlevel=%sysfunc(getvarc(&dsid,&num_notlevel));
%let rc=%sysfunc(close(&dsid));
/* Plot */
ods html style = statistical;
proc template;
define statgraph sgplot;
mvar OUTCOME_TITLE OUTCOME_LEV OUTCOME AUC Sens Section2 Specif level notlevel confusion 
auc_table
XVAR_LABEL  TP FP FN TN TP_per FP_per FN_per TN_per Sensitivity Specificity; 
dynamic XVAR1 YVAR1;
begingraph;
EntryTitle "Operating Characterisitcs for " OUTCOME_TITLE " = " OUTCOME_LEV /;
if (confusion = "Y") 
EntryFootnote textattrs=(color = red )"*" textattrs=()"Confusion Matrix is for the "
XVAR_LABEL " at the cutoff = " textattrs=(color = blue ) Section2 ; 
*EntryFootnote "At cutoff, Sensitivity = " textattrs=(color = blue ) Sensitivity " and 
Specificity = " textattrs=(color = blue ) Specificity  ;
endif;
layout lattice / columns=1 rows=2  rowgutter=2px columngutter=2px
                        rowweights=(.7 .3 );
       columnaxes;
         columnaxis / label=XVAR_LABEL;
endcolumnaxes;
    sidebar / align=right; 
      layout gridded / border = false; 
layout overlay;

DiscreteLegend "SERIES" "SERIES1" / location = outside ACROSS = 1 down 
= 2; 
endlayout; 
layout overlay;
*DiscreteLegend "SERIES3" / location = outside;
endlayout; 
endlayout; 
     
    endsidebar; 
layout overlay / cycleattrs=true xaxisopts = (Label = XVAR_LABEL) yaxisopts=( 
Label="Proportion" type=linear linearopts=( tickvaluelist=( 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0.9 1 ) viewmin=0 viewmax=1 ) );
    StepPlot X=X_VALUE Y=SPEC / primary=true LegendLabel="Specificity" NAME="SERIES";
StepPlot X=X_VALUE Y=_SENSIT_ / LegendLabel="Sensitivity" NAME="SERIES1";
Referenceline x =  Section2 / LINEATTRS = (color = black pattern = 2);
if (confusion = "Y") 
/* 2 by 2 */
layout gridded /  columns=3 rows = 3
         halign = CENTER valign = top border=true 
opaque=true backgroundcolor=GraphWalls:color location = outside; 
entry halign=center " " / textattrs=GraphValueText; 
       entry halign=center level  / textattrs=(WEIGHT=bold) /* Setting the titles to bold */  
 ;
entry halign=center notlevel  / textattrs=(WEIGHT=bold); 
   entry halign=center "Event" / textattrs=(WEIGHT=bold) ; 
entry halign=center TP textattrs=(STYLE=ITALIC color = blue )"(" TP_per "%)" /* 
Formatting the percentages */ ; 
entry halign=center FP textattrs=(STYLE=ITALIC color = blue )"(" FP_per "%)" ; 
entry halign=center "Non Event" / textattrs=(WEIGHT=bold); 
      entry halign=center FN  textattrs=(STYLE=ITALIC color = blue )"(" FN_per "%)"; 
       entry halign=center TN textattrs=(STYLE=ITALIC color = blue )"(" TN_per "%)"; 
    
endlayout;
endif;
if (auc_table = "Y") 
layout gridded / columns=4 rows = 1
         halign = center valign = top border=true 
opaque=true backgroundcolor=GraphWalls:color location = outside; 
entry halign=center "AUC " / textattrs=(WEIGHT=bold); 
       entry halign=center "100% Sensitivity: " / textattrs=(WEIGHT=bold); 
entry halign=center "Intersection: "  / textattrs=(WEIGHT=bold); 
   entry halign=center "100% Specificity: " / textattrs=(WEIGHT=bold) ; 
entry halign=center AUC; 
      entry halign=center Sens; 
       entry halign=center Section2; 
entry halign=center Specif ; 
endlayout;

endif;
*DiscreteLegend "SERIES" "SERIES1" / location = outside;
endlayout;
layout lattice / columns=1 rows=2  rowgutter=2px columngutter=2px
                        rowweights=(.5 .5 ) COLUMNDATARANGE = Union;
columnaxes;
         columnaxis / label=XVAR_LABEL;
       endcolumnaxes;
    layout gridded / AUTOALIGN=NONE  columns=3 rows = 3
         halign = right valign = top border=true 
opaque=true backgroundcolor=GraphWalls:color location = outside; 
endlayout;
layout overlay /   yaxisopts=( Label=" " );
Scatterplot x=XVAR1 y = OUTCOME / group = OUTCOME INDEX = Index1 NAME = 
"SERIES3";
Referenceline x =  Section2 / LINEATTRS = (color = black pattern = 2);
endlayout;
layout overlay /   yaxisopts=( Label=" " );
Boxplot y = XVAR1 x = OUTCOME / ORIENT = horizontal;
Referenceline x =  Section2 / LINEATTRS = (color = black pattern = 2);
endlayout;
endlayout;
endlayout;
endgraph;
end;
run;
%let OUTCOME_TITLE = &OUTCOME;
%let OUTCOME_LEV = &OUTCOME_LEV;
%let OUTCOME =  &OUTCOME._char;
%let AUC =  &AUC_format;
%let Sens =  &Sens_format;
%let Section2 = &Section_format;
%let Specif =  &Spec_format;
%let level = &level;
%let notlevel = ¬level;
%let confusion = &confusion;
%let auc_table = &auc_table;
%let XVAR_LABEL = &XVAR_LABEL;
%let TP = &TP;
%let FP = &FP;
%let FN = &FN;
%let TN = &TN;
%let TP_per = &TP_format;
%let FP_per = &FP_format;
%let FN_per = &FN_format;
%let TN_per = &TN_format;
ods listing sge = on;
proc sgrender data=merge3 template=sgplot;
dynamic XVAR1 = "XVAR1" YVAR1 = "YVAR1";
run;
ods listing close;
proc sql;
drop table roc_association, Roctest, ROCData, Outest, _merge, merge2, merge3, 
One_Hundred_PC_Sens, 
sens_spec, sens_spec_intersection, One_Hundred_PC_Spec, freq_input,
CrossTabFreqs, crossTabFreqs_ref, TP, FP, FN, TN, distinct_values, 
LEVEL, NOTLEVEL, level_char, notlevel_char;
quit;
%mend roc_cutoff;

CONTACT INFORMATION

Your comments and questions are valued and encouraged. Please contact the author at:

Kriss Harris

GlaxoSmithKline

Mail Stop : 1B106003

Gunnels Wood Road

Stevenage SG1 2NY

Phone: 01438 76 4904

Email: Kriss.5.Harris@gsk.com

Web: http://www.krissharris.co.uk

A Macro to Create a Batch Submit SAS Program

2011-03-07T12:09:00.006-04:00

原文載點: http://support.sas.com/resources/papers/proceedings10/092-2010.pdf

在一些很特殊的情況下，我們需要大量執行不同的 SAS 程式。若一個一個從 SAS 裡面開啟再執行，不但相當消耗時間，也不沒有太大的效率。Helen Sun 和 Cindy Wong 在 SAS Global Forum 2010 發表了一個簡易的巨集程式，讓 SAS 使用者在 PC 的視窗系統底下可以輕鬆批次地將同一個目錄底下的 SAS 程式碼給一次執行完畢。

假設有一堆 SAS 程式被儲存在 C:\temp 的目錄底下，如下圖所示：

先執行以下的巨集程式一次：

    %MACRO BATCH(dir=); 
      /*read in file names*/ 
      FILENAME dirfile PIPE "dir ""&dir"" /b" LRECL=200; 
      DATA in; 
        INFILE dirfile LENGTH=len; 
        INPUT char $varying200. len; 
      RUN; 
  
     /*only choose sas program names*/ 
     DATA sasprg; 
       SET in; 
       IF INDEX(char,'.sas '); 
     RUN; 
     /*create statement*/ 
     DATA line; 
       LENGTH name $200; 
       SET sasprg; 
       name="%include '"||trim(char)||"';"; 
       DROP char; 
     RUN; 
     /*write out to a file*/ 
     FILENAME fileref "&dir.batch.sas"; 
     DATA _null_; 
       SET line; 
       FILE fileref; 
       PUT name; 
     RUN; 
     FILENAME fileref clear; 
   %MEND;

然後再把路徑輸進 %BATCH 裡面：

%BATCH(dir=C:\temp\);

接著就會在這個路徑看到一個名為 batch.sas 的程式，打開來一看如下所示：

也就是說 %BATCH 幫使用者把該路徑底下的程式名稱一一放進 %include 裡面，然後我們就可以直接執行 batch.sas 這個程式來批次執行該目錄底下所有 SAS 程式。

CONTACT INFORMATION
Helen Sun
SAS Programmer
hsun@robarts.ca
Robarts Clinical Trials
Robarts Research Institute
P.O. Box 5015, 100 Perth Drive
London, Ontario, Canada
N6A 5K8
Tel: 519-931-5700x24254
Fax: 519-931-5284
www.robartsclinicaltrials.ca

Disagreement on Agreement: Two Alternative Agreement Coefficients

2011-02-25T17:02:00.003-04:00

原文載點：http://www2.sas.com/proceedings/forum2007/186-2007.pdf

Kappa一致性係數是用來評量許多評分者對某一類別變數評分是否具有一致性的指標。根據公式，它是實際評定一致的次數百分比與評分者理論上評定的最大可能次數百分比的比率。但是近年來有越來越多的學者發現這個指標的不完美處。於是有許多替代性的指標被研究出來，以克服Kappa一致性係數的缺點。然而SAS並沒有對這些新的一致性指標的程式語法更新至新版的SAS中。也許我們可以期待新的一致性係數會在不久的將來被納入到SAS內建的語法，但畢竟SAS改版的頻率並沒有很高，所以還是用自製的巨集程式比較有效率。本篇教學文件主要是介紹由Kilem Gwet所發表的一階和二階一致性係數(first-order and second order agreement coefficient，簡稱AC1和AC2)在SAS裡面如何操作。其巨集程式由 Emily Blood 和 Kevin Spratt 編寫，並發表於 2007 SAS Global Forum 上。

AC1和AC2的分別在於，AC1是使用任意的類別評分工具下的一致性係數，而AC2是使用已經排序好的類別平分工具下的一致性係數。其理論背景可參照原文。本文主要介紹巨集程式 %AC1AC2 的使用方式。

巨集語法如下:

%AC1AC2(dataset=,numcategories=,betadata=)

內含三個參數，分別為:

dataset: 使用SAS資料檔的名稱
numcategories: 評量工具內的等級數量，假設評分標準是1,2,3,4,5，則此參數並須設定為"5"
betadata: 錯誤分類機率矩陣的SAS資料檔名稱，主要是拿來計算AC2之用

用一個範例來看執行結果:

資料:

DATA work.one;
 INPUT r1 r2 r3 r4 r5 r6;
CARDS;
 4 4 4 4 4 4
 2 2 5 2 5 5
 3 3 5 2 3 3
 5 5 5 5 5 5
 2 4 2 4 4 2
 1 3 3 3 1 3
 3 5 3 3 5 3
 1 1 3 3 3 4
 4 4 4 4 1 1
 5 5 5 5 5 5
 1 4 4 4 4 4
 1 4 2 4 4 4
 2 3 2 2 3 3
 4 1 4 4 4 4
 2 2 4 4 4 5
 3 3 5 3 3 3
 5 5 1 1 1 4
 1 1 1 1 2 1
 2 2 4 4 4 4
 1 3 3 5 5 5
 5 5 5 5 5 5
 4 4 2 4 4 4
 5 2 5 5 4 2
 1 4 4 4 1 4
 5 4 4 4 4 1
 2 4 2 2 2 2
 1 5 1 1 1 5
 4 2 4 4 4 2
 1 3 3 3 3 3
 5 5 5 5 5 5
;
RUN;

data work.beta;
input cat1 cat2 cat3 cat4 cat5 ;
cards;
.9  .9 .2 .1 0
.05 .1 .8 .7 0
.03  0  0 .1 0
.01  0  0 .1 0
.01  0  0  0 1
;
RUN;

第一筆資料『one』為原始資料，包含六名心理學家對患有某種心理疾病的病人所做的評量，分數從1~5。第二筆資料『beta』裝有錯誤分類的機率資料，代表了不一致性發生時的嚴重程度。其矩陣表示法如下所示：

其中『a|b』表示當原先評分為 b 時，之後評分變成 a 。因此，很明顯的，對角線的數據即表示具有一致評分的機率，而非對角線的數據即表示不一致評分的機率。特別要注意的是，設定這個矩陣的SAS資料集時，每一欄的總和一定要等於1。這個矩陣只有拿來運算AC2，在AC1的公式裡面並不包含這個矩陣。

設定好參數並執行巨集，則會在SAS output視窗上面看到 AC1 和 AC2 各有五個數據出現，從左到右分別是一致性係數值(Agree Stat)、條件變異數(CondVar)、非條件變異數(UncondVar)、評分者評定一致的百分比(PA)以及評分者理論上評定一致的百分比(PE)。

原始巨集碼可至這個網址下載：http://mcrc.hitchcock.org/SASMacros/Agreement/AC1AC2.txt

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the authors at:
Emily A. Blood or Kevin F. Spratt
Department of Orthopedic Surgery
Dartmouth Medical School
One Medical Center Drive
Hanover, NH 03756
Work Phone: 603-653-6019
E-mail: Emily.a.blood@dartmouth.edu or Kevin.F.Spratt@dartmouth.edu

Automatically Run Fisher's Exact Test When the Chi-square Test Might Not Be Valid

2010-08-02T10:06:00.001-04:00

Link: http://support.sas.com/resources/papers/proceedings10/096-2010.pdf

在使用PROC FREQ執行卡方檢定時，如果看到列聯表下面出現warning的訊息，就表示這個列聯表有50%以上的期望值小於五，因此程式會建議使用Fisher's exact test。然後使用者必須回到程式裡面，在tables statement後面加上一個exact的選項，再重跑一次PROC FREQ才能得到Fisher's exact test的結果。Wei Xu於SAS GLOBAL FORUM 2010發表了一個macro程式，讓電腦自動幫使用者判斷要不要進行Fisher's exact test，若需要的話會順便把Fisher's exact test的結果跑出來，使用者不需再去重新執行PROC FREQ一次。

這個macro名叫%RUN_FISHERS，裡面只有四個參數須要設定:

Data= 要進行卡方檢定的資料名稱
Row= 放在行的變數名稱
Col= 放在列的變數名稱
Count= 如果資料是已經整理好的列聯資料，則會多一個計數變數來表示某一欄位的個數，此時需要把這個變數定義在這個參數裡面。如果不是列聯資料的話，可在此留空白。

範例:

 data test; 
   input r c ct; 
   datalines; 
   1 1 12 
   1 2 5 
   1 3 6 
   2 1 13 
   2 2 2 
   2 3 4 
   ; 
 run;   
%run_fishers(data=test, row=r, col=c, count=ct);

由於此資料會出現warning訊息，所以Fisher's exact test的結果就會直接顯示在報表最下方。如果今天資料變成:

data test; 
input r c ct; 
datalines; 
1 1 12 
1 2 50 
1 3 60 
2 1 13 
2 2 20 
2 3 40 
; 
run;   
%run_fishers(data=test, row=r, col=c, count=ct);

報表最下方就不會出現Fisher's exact test的結果。

最後是這個macro的原始碼:

%macro run_fishers (version, data=, row=, col=, count=); 
%let _version=1.0; 
%if &version ne %then %put RUN_FISHERS macro Version &_version; 
%let opts = %sysfunc(getoption(notes))  
            _last_=%sysfunc(getoption(_last_)); 
%if &version ne debug %then %str(options nonotes;); 
 
/* Check for newer version */ 
%if %sysevalf(&sysver >= 8.2) %then %do; 
%let _notfound=0; 
 
filename ver url 'http://ftp.sas.com/techsup/download/stat/versions.dat' termstr=crlf; 
 data _null_; 
    infile ver end=_eof; 
    input name:$15. ver; 
    if upcase(name)="&sysmacroname" then do; 
       call symput("_newver",ver); stop; 
    end; 
 
 
    if _eof then call symput("_notfound",1); 
run; 
  %if &syserr ne 0 or &_notfound=1 %then 
    %put &sysmacroname: Unable to check for newer version; 
  %else %if %sysevalf(&_newver > &_version) %then %do; 
    %put &sysmacroname: A newer version of the &sysmacroname macro is available.; 
    %put %str(         ) You can get the newer version at this location:; 
    %put %str(         ) http://support.sas.com/ctx/samples/index.jsp; 
  %end; 
 %end; 
  
proc freq data=&data noprint; 
  %if &count ne %then %str(weight &count / zeros;); 
  tables &row*&col / sparse outexpect out=_out1; 
run; 
 
proc means data=_out1 noprint; 
  var count; 
  output out=_out2; 
run; 
 
data _null_; 
  set _out1; 
  if expected<=5 then warn+1; 
  if _n_=1 then set _out2; 
  pct_lt5=warn/_freq_; 
  if _freq_=_n_; 
  warning=(pct_lt5>=.2); 
  call symput('warning',warning); 
run; 
 
options &opts; 
 
proc freq data=_out1; 
  weight count / zeros; 
  tables &row*&col / chisq; 
  %if &warning=1 %then %do; 
  exact fisher; 
  %end; 
run; 
 
%mend;

CONTACT INFORMATION
Wei Xu
Boston Scientific
100 Boston Scientific Way
Marlborough, MA 01752-1234
(508) 683-4264
wei.xu@bsci.com

Stupid Human Tricks with PROC EXPAND®

2010-07-23T17:51:00.000-04:00

Link: http://support.sas.com/resources/papers/proceedings10/093-2010.pdf

鮮少有人知道 SAS/ETS 裡面有個叫做 PROC EXPAND 的程序。這個程序主要是用來整理時間序列資料，比方說可以自己定義時間區間後再來做一些次數量表的製作。David L. Cassell 利用這個程序進行了一些平常使用者可以用來做資料處理的功能，並發表了一篇教學文件在 SAS GLOBAL FORUM 2010 上。讓我們來看看如何用 PROC EXPAND 程序來取代一些複雜的資料操作。

TRICK 1: LAGS AND LEADS
在 data step 裡面，我們可用 LAG() 或 LAGn() 函數來求出某個變數的前1~n天的數據，但如果反過來要找某變數的後1~n天數據，雖然有許多方法可以做，但並沒有一個簡便的 call function 可以處理。但是，用 PROC EXPAND 可以輕鬆地製作出這兩種數據出來。範例程式如下:

proc expand data=YourData method=none;
by pt;
convert dosedt = lead1_dt / transformout = (lead 1);
convert dosedt = lead2_dt / transformout = (lead 2);
convert dosedt = lag1_dt / transformout = (lag 1);
convert dosedt = lag2_dt / transformout = (lag 2);
run;

其中，我們打蘇案要轉換的變數是dosedt，而在 PROC EXPAND 裡面，只需要用 convert 指令，把新變數名稱定義好先(lead1_dt, lead2_dt, lag1_dt, lag2_dt)，然後在後面斜線後加上一個 transformout= 的選項，然後在等號括弧裡面寫上 lead n 或 lag n 便可以造出後 n 天或指是前 n 天的變數。
TRICK 2: MAX OR MIN OF THE LAST K RECORDS
如果想要知道某一個變數在前n筆或後n筆資料的最大值或最小值，一般的做法是會在PROC MEANS或PROC SUMMARY後加上一個where statement來限制程序要讀取的資料範圍，但在PROC EXPAND裡面，可以直接指定前n筆或後n筆資料，如下所示:

proc expand data=YourData out=YourMax method=none;
by factory;
convert x = max_x / transformin=(movmax 50);
convert x = min_x / transformin=(movmin 50);
run;

同樣利用convert statement，在後面加上transformin選項，並用"movmax n"或"movmin n"來讓SAS去讀前n筆或後n筆資料的最大值最小值。實際上SAS還是會從第一筆資料開始讀，每讀一筆資料就會算一次最大最小值，直到讀進來的資料數量大於movmax和movmin所指定的n值後，SAS才會開始去計算前n筆或後n筆資料的最大值最小值，而所有的結果都會完整呈現在新的資料裡面，我們只需要取最後一筆資料的x_max和x_min數據即為所求。
TRICK 3: AGGREGATING CHUNKS OF RECORDS
我曾經在網路上看過很多人詢問過如何讓資料每跳幾筆資料算一個平均數，或者是每跳幾筆資料取一個數據。許多人用數個PROC程序和data step來完成這兩個動作，但在PROC EXPAND裡面，只需要兩個指令便可完成。

假設資料如下:

data temp1;
input obs vol price;
datalines;
1 2 11
2 2 11
3 2 12
4 3 13
5 3 11
6 3 12
7 4 14
8 4 12
9 4 12
10 5 11
11 5 16
12 5 14
13 6 10
;
run;

總計十三筆資料，兩個變數(vol, price)。如果我們想要每三筆資料算一次vol的總值，以及每三筆資料抓出price的數據，則程式如下:

proc expand data=temp1 out=exp1 factor=(1:3);
convert vol=aggrvol / observed=total;
convert price=new_price / observed=end;
run;

結果如下:

TIME AGGRVOL NEW_PRICE
0 6 12
3 9 12
6 12 12
9 15 14

如何讓成是每三筆資料做一次運算，關鍵完全在PROC EXPAND後面的factor(n:m)選項。以此例來看，factor(1:3)便是要求SAS每三筆資料做一次運算。至於運算的方法，則是定義在convert statement後面的observed=選項。當observed=total時，則會求出每三筆資料的總合，當observed=end時，則會求出每三筆資料的最後一筆數據。實際上，observed=選項有許多參數可以設定，比方要算平均數就只要用observed=average即可。所有可使用的選項可以參考SAS線上手冊。

TRICK 4: THE MOVING AVERAGE
最後一個是算moving average，而這也是PROC EXPAND在處理時間序列資料時最常使用到的計算。假設我們要算每筆資料含前四筆資料的平均，程式如下:

proc expand data=temp1 out=YourOut method=none;
convert number1=mean1 / transformout=(movave 5);
run;

如同前面所述，我們只要在convert statement後面加上一個transformout=的選項，並於括弧內寫上"moveave 5"便可輕鬆算出moving average。

CONTACT INFORMATION
David L. Cassell
Design Pathways
3115 NW Norwood Pl.
Corvallis, OR 97330
DavidLCassell@msn.com
541-754-1304

Dropping Automatically Variables with Only Missing Values

2010-07-23T15:14:00.005-04:00

Link: http://support.sas.com/resources/papers/proceedings10/048-2010.pdf

在進行資料分析前，有些人會習慣把一些含有missing data的樣本給清除掉，雖然這不會影響到分析結果，因為大部分的SAS程序都是用CCA(Complete Case Analysis)來處理含有missing data的數據，不過若要用手動的方法來清掉missing data的話，在遇到龐大數量的變數時，需要消耗很多時間在輸入遍數名稱上。美國人口普查局的Selvaratnam Sridharma發表了一個macro程序於SAS Global Forum 2010，讓這個程式撰寫的過程只需要幾秒鐘的時間就可以完成。

這個macro如下所示:

%DROPMISS(DSIN, DSOUT, NODROP);

裡面只需要定義兩個macro參數:

DSIN: 原始資料的名稱
DSOUT: 新資料的名稱
NODROP: 不要處理missing data的變數名稱，此為optional選項。

其中，DSIN和DSOUT應該不用多做解釋。NODROP如果沒有指定特定的變數名稱的話，這個macro就會針對所有在DSIN所指定的資料裡面的變數進行missing data的處理。但如果想要讓這個macro不處理所有文字型變數或數值型變數，則可以用 __NUMERIC_ 或_CHARACTER_ 來限制。如下所示:

%DROPMISS (DSIN=olddata,DSOUT=newdata, nodrop= _NUMERIC_ );
%DROPMISS (DSIN=oldedata,DSOUT=newdata, nodrop= _CHARACTER_ );

此macro的原始碼如下:

/******************/
options nomprint noSYMBOLGEN MLOGIC;
/****************************/
%macro DROPMISS( DSNIN /* name of input SAS dataset
*/
, DSNOUT /* name of output SAS dataset
*/
, NODROP= /* [optional] variables to be omitted from dropping even if
they have only missing values */
) ;
/* PURPOSE: To find both Character and Numeric the variables that have only
missing values and drop them if
* they are not in &NONDROP
*
* NOTE: if there are no variables in the dataset, produce no variables
processing code
*
*
* EXAMPLE OF USE:
* %DROPMISS( DSNIN, DSNOUT )
* %DROPMISS( DSNIN, DSNOUT, NODROP=A B C D--H X1-X100 )
* %DROPMISS( DSNIN, DSNOUT, NODROP=_numeric_ )
* %DROPMISS( DSNIN, DSNOUT, NOdrop=_character_ )
*/
%local I ;
%if "&DSNIN" = "&DSNOUT"
%then %do ;
%put /------------------------------------------------\ ;
%put | ERROR from DROPMISS: | ;
%put | Input Dataset has same name as Output Dataset. | ;
%put | Execution terminating forthwith. | ;
%put \------------------------------------------------/ ;
%goto L9999 ;
%end ;
/*###################################################################*/
/* begin executable code
/*####################################################################/
/*===================================================================*/
/* Create dataset of variable names that have only missing values
/* exclude from the computation all names in &NODROP
/*===================================================================*/
proc contents data=&DSNIN( drop=&NODROP ) memtype=data noprint out=_cntnts_( keep=
name type ) ; run ;
%let N_CHAR = 0 ;
%let N_NUM = 0 ;
data _null_ ;
set _cntnts_ end=lastobs nobs=nobs ;

if nobs = 0 then stop ;
n_char + ( type = 2 ) ;
n_num + ( type = 1 ) ;
/* create macro vars containing final # of char, numeric variables */
if lastobs
then do ;
call symput( 'N_CHAR', left( put( n_char, 5. ))) ;
call symput( 'N_NUM' , left( put( n_num , 5. ))) ;
end ;
run ;
/*===================================================================*/
/* if there are no variables in dataset, stop further processing
/*===================================================================*/
%if %eval( &N_NUM + &N_CHAR ) = 0
%then %do ;
%put /----------------------------------\ ;
%put | ERROR from DROPMISS: | ;
%put | No variables in dataset. | ;
%put | Execution terminating forthwith. | ;
%put \----------------------------------/ ;
%goto L9999 ;
%end ;
/*===================================================================*/
/* put global macro names into global symbol table for later retrieval
/*===================================================================*/
%LET NUM0 =0;
%LET CHAR0 = 0;
%IF &N_NUM >0 %THEN %DO;
%do I = 1 %to &N_NUM ;
%global NUM&I ;
%end ;
%END;
%if &N_CHAR > 0 %THEN %DO;
%do I = 1 %to &N_CHAR ;
%global CHAR&I ;
%end ;
%END;
/*===================================================================*/
/* create macro vars containing variable names
/* efficiency note: could compute n_char, n_num here, but must declare macro names
to be
global b4 stuffing them
/*
/*===================================================================*/
proc sql noprint ;
%if &N_CHAR > 0 %then %str( select name into :CHAR1 - :CHAR&N_CHAR from
_cntnts_ where type = 2 ; ) ;
%if &N_NUM > 0 %then %str( select name into :NUM1 - :NUM&N_NUM from
_cntnts_ where type = 1 ; ) ;
quit ;
/*===================================================================*/
/* Determine the variables that are missing

/*
/*===================================================================*/
%IF &N_CHAR > 1 %THEN %DO;
%let N_CHAR_1 = %EVAL(&N_CHAR - 1);
%END;
Proc sql ;
select %do I= 1 %to &N_NUM; max (&&NUM&I) , %end; %IF &N_CHAR > 1 %THEN %DO;
%do I= 1 %to &N_CHAR_1; max(&&CHAR&I), %END; %end; MAX(&&CHAR&N_CHAR)
into
%do I= 1 %to &N_NUM; :NUMMAX&I , %END; %IF &N_CHAR > 1 %THEN %DO;
%do I= 1 %to &N_CHAR_1; :CHARMAX&I,%END; %END; :CHARMAX&N_CHAR
from &DSNIN;
quit;
/*===================================================================*/
/* initialize DROP_NUM, DROP_CHAR global macro vars
/*===================================================================*/
%let DROP_NUM = ;
%let DROP_CHAR = ;
%if &N_NUM > 0 %THEN %DO;
DATA _NULL_;
%do I = 1 %to &N_NUM ;
%IF &&NUMMAX&I =. %THEN %DO;
%let DROP_NUM = &DROP_NUM %qtrim( &&NUM&I ) ;
%END;
%end ;
RUN;
%END;
%IF &N_CHAR > 0 %THEN %DO;
DATA _NULL_;
%do I = 1 %to &N_CHAR ;
%IF "%qtrim(&&CHARMAX&I)" eq "" %THEN %DO;
%let DROP_CHAR = &DROP_CHAR %qtrim( &&CHAR&I ) ;
%END;
%end ;
RUN;
%END;
/*===================================================================*/
/* Create output dataset
/*===================================================================*/
data &DSNOUT ;
%if &DROP_CHAR ^= %then %str(DROP &DROP_CHAR ; ) ; /* drop char variables
that
have only missing values */
%if &DROP_NUM ^= %then %str(DROP &DROP_NUM ; ) ; /* drop num variables
that
have only missing values */
set &DSNIN ;
%if &DROP_CHAR ^= or &DROP_NUM ^= %then %do;

%put /----------------------------------\ ;
%put | Variables dropped are &DROP_CHAR &DROP_NUM | ;
%put \----------------------------------/ ;
%end;
%if &DROP_CHAR = and &DROP_NUM = %then %do;
%put /----------------------------------\ ;
%put | No variables are dropped |;
%put \----------------------------------/ ;
%end;
run ;
%L9999:
%mend DROPMISS ;

CONTACT INFORMATION

Selvaratnam Sridharma
Economic Planning and Coordination Division
U.S. Bureau of the Census
Address
Washington, DC 20233-6100
301-763-6774
Email: selvaratnam.sridharma@census.gov

PROC MIXED: Macro to assess fixed and random effects for significance using the Likelihood Ratio test and the approximate Mixture Method

2010-03-26T00:21:00.006-04:00

原文載點：http://analytics.ncsu.edu/sesug/2009/SC013.Bardenheier.pdf

若要從兩個nested mixed model挑選較佳的模式，比較專業的人會推薦使用Likelihood Ratio Test（以下簡稱LRT）來處理這類的問題。PROC MIXED程序並沒有語法來進行這個檢定，但這個檢定需要使用到的參數，卻全部都可以從PROC MIXED程序的輸出報表內找到。懂得運作的使用者可以從報表拿擷取相關參數，然後自己算一算再查個表大概就知道結果了。我自己曾經在準備博士班資格考時寫了一個簡易的macro來做LRT，但沒想到真的有人把這個觀念發表出來。我稍微看了一下這份技術文件內的macro寫法，與我當年寫的大同小異，但作者有將輸出報表格式設計的比較好看一點，所以我就來介紹一下這個方便進行LRT的macro。

假設兩個nested mixed model分別稱為reduced model和full model，通常reduced model會被設定在虛無假設H0下，full model會被設定在對立假設H1下。而整個LRT會從這兩個模式各取出兩個數據出來，一個是他們的自由度，另一個是他們-2log likelihood function的值。這兩個數據全部都可以從PROC MIXED的報表內找到。我以PROC MIXED線上手冊所附的一個範例為例，請參照這個網頁：

http://www.tau.ac.il/cc/pages/docs/sas8/stat/chap41/sect30.htm

這兩個數據可以在這個地方找到，圖示如下（紅色箭頭處）：

至於公式，引用 wiki 的資料：

白話一點，就是把 reduced model 的 -2 log likelihood function 的值減掉 full model 的 -2 log likelihood function 的值，就可以算出上面那個 D。然後兩個模式的自由度數量的差異套進卡方分配裡面，就可以算出 p-value 了。如果 p-value < 0.05，就表示代表 reduced model 的 H0 被拒絕了，也就是說 full model 比 reduced model 好。反之亦然。

這種計算其實根本不需要使用到電腦的，但懶人永遠可以找藉口，所以原文作者還是花了將近兩頁的篇幅寫了這個 macro 來服務大眾。

%macro LRMixtureTest(fullmodel=,redmodel=,DFfull=,DFred=); 
/* 
This macro calculates the Likelihood Ratio Test for fixed effects in a 
mixed model 
See Singer Applied Longitudinal Data Analysis pp 116-120. 
and also calculates the mixture method pvalue for random effects in a 
mixed model -see KKNM chapter 26. 
The macro argument &fullmodel refers to the output dataset produced by 
running the full model in proc mixed and using the statement: 
    ods output FitStatistics=fm; 
The macro argument &DFfull refers to the output dataset produced by 
running the full model in proc mixed and adding to the statement: 
    ods output FitStatistics=fm SolutionF=SFfm; The number of 
parameters being tested (ie, fixed effects) determine the degrees of 
freedom in the LRT and come from this dataset. 
The argument &redmodel refers to the output dataset produced by running 
the reduced model in proc mixed using the statement: 
   ods output FitStatistics=rm  SolutionF=SFrm ; 
Maximum Likelihood should be used when comparing fixed effects because 
it maximizes the likelihood of the sample data. Using the Restricted 
Maximum likelihood should only be used when assessing a random term 
(ie, the fixed effects should be the same in both models)  
because in RML we maximize the likelihood of the residuals. The degrees 
of freedom would not be correct if the number of fixed effects were 
different in each model when using RML -which is the default in PROC 
MIXED. 
*/ 
data &fullmodel; set &fullmodel; 
 if descr='-2 Res Log Likelihood' then chisqfullRML=value;  *RML; 
 else if descr='-2 Log Likelihood' then chisqfullFML=value; *FML; 
run; 
data &redmodel; set &redmodel; 
  if descr='-2 Res Log Likelihood' then chisqredRML=value; *RML; 
  else  if descr='-2 Log Likelihood' then chisqredFML=value; *FML; 
run; 
proc sql; 
CREATE TABLE flDF AS SELECT  effect, count(*) AS fullDF from 
work.&DFfull; quit; *Count number of effects to get degrees of freedom; 

proc sql; 
create table rdDF as select  effect, count(*) as redDF from 
work.&DFred; quit; 

data degfree (drop=effect); 
  merge  work.flDF work.rdDF; 
if _n_=1 then  LRTDF= fullDF  - redDF; run; 

data likelihood; 
 merge &fullmodel &redmodel degfree; 
 testintRML=abs((chisqredRML)-(chisqfullRML)); **Models can yield 
negative LLs, those that are smaller in absolute value -ie, closer to 0 
fit better (pg 116-117, Singer);  

7 
testintFML=abs((chisqredFML)-(chisqfullFML)); **Models can yield 
negative LLs, those that are smaller in absolute value -ie, closer to 0 
fit better (pg 116-117, Singer); 
 pvaluemixture=(( .5*(1-probchi(testintRML,2)) + .5*(1-
probchi(testintRML,1)) ));    *for random terms; 
  pvalueLRT=(1-probchi(testintFML,LRTDF));      *For fixed terms; 
  
proc print data=likelihood split='*' noobs; 
   var  testintFML  pvalueLRT; 
 format  pvalueLRT 6.4; 
 where testintFML ne .; 
   label testintFML='likelihood ratio*test statistic* 
         comparing*reduced model to full model' 
                   pvalueLRT='p-value LRT'; 
      title "Likelihood Ratio Test for fixed effects"; run; 
  
proc print data=likelihood split='*' noobs; 
   var  testintRML  pvaluemixture; 
 format pvaluemixture  6.4; 
 where  testintRML ne . ; 
   label testintRML='likelihood ratio*test statistic* 
         comparing*reduced model to full model' 
                   pvaluemixture='mixture method p-value'; 
      title "Mixture method Test for random effects"; 
 run;  
%mend LRMixtureTest;

這個 macro 裡面綠色的部分完全都是作者的註釋和說明，只有藍色的部分才是真正 SAS 會去執行的程式碼。其中只有四個參數要說明：

fullmodel：full model 的 -2 log likelihood function 值
redmodel：reduced model 的 -2 log likelihood function 值
DFfull：full model 的自由度
DFred：reduced model 的自由度

四個參數設定就和上一段的敘述相同。全部都是必要的參數，缺一不可。因此，只要呼叫 %macro LRMixtureTest 再套上四個參數值，電腦就會幫你做 LRT。我們來看原文內的一個範例。

範例模式：

這模式只使用一個因變數叫做 PCTBLK，但他同時是固定效應也是隨機效應。假設這是個 full model，而 reduced model 就設定為拿掉 PCTBLK 固定效應的模式，也就是拿掉上面劃紅線的那一段。

因此，先跑兩個程式把那四個要用到的參數算出來：

%include '\\cdc\private\mixture method pvalue macro1.sas';
proc mixed data=bhb7.allUS2levge10blk method=ml covtest empirical 
noclprint ;
class  state_id;
model diffcov=pctblk/ ddfm=contain s;
random intercept pctblk /subject=state_id type=ar(1) s ;
ods output FitStatistics=fm  SolutionF=SFfm ;
run;

proc mixed data=bhb7.allUS2levge10blk method=ml covtest empirical
noclprint;
class  state_id;
model diffcov=/ddfm=contain   s;
random intercept  pctblk /subject=state_id type=ar(1) s;
ods output FitStatistics=rm SolutionF=SFrm ;
run;

跑完後，兩個模式的 -2 log likelihood function 值和自由度會被 ods output 輸出到不同的資料裡面，接著使用那個 macro 去把從這四個資料裡面抓數值出來算：

%LRMixtureTest(fullmodel=fm,redmodel=rm,DFfull=SFfm,DFred=SFrm);

最後這個 macro 就會輸出下面這個報表給你：

結果顯示 p-value < 0.05，所以 full model 比 reduced model 好。

從這個範例中，特別要強調一點是，如果兩個 nested mixed model 差異只有在固定效應時，一開始那兩個 PROC MIXED 執行時，後面的 method option 必須強制設定為 ML，因為如果沒有特別去設定的話，PROC MIXED 預設的參數估計方法是 REML，但這只有在兩個模式的差異是在隨機效應或者是 covariance structure 時才能使用。原文中的第一個範例就是用 REML 再估計參數，因為兩個模式的差別只有在一個 random intercept 而已。

Contact Information:
Barbara Bardenheier, MPH, MA
Centers for Disease Control and Prevention
1600 Clifton Rd, MS E-52
Atlanta, GA 30333
Tel : (404) 639-8789
Fax : (404) 639-8614
Email : bfb7@cdc.gov

A Macro for Computing the Mantel-Fleisś Criterion

2010-03-21T21:32:00.001-04:00

原文載點：http://analytics.ncsu.edu/sesug/2009/PO008.Welch.pdf

在作2乘2的卡方分析時，我們通常會看每一格裡面的樣本期望值是不是大於五，如果不大於五的比例超過一半的話，便需要使用 Fisher's exact test 的結果來取代原本的卡方分析結果。但這種判斷方法，在2乘2乘H的表格下並不適用。Mantel 和 Fleisś 在 1980 年代發明了一種判斷方法來專門處理這種情況，不過這個方法並沒有內建在 SAS 現存的程序內，因此 Brandon Welch, Jane Eslinger 和 Rob Woolson 在 2009 年的 SESUG 發表了一篇技術文件，提供一個簡便的 macro 供使用者使用。

有關 Mantel-Fleisś criterion 的理論，可詳見原文。以下就直接切入這個 macro 的使用方法：
%MantelFleiss(In,Strat,Var1,Var2,Count)

這個名為 MantelFleiss 的 macro 只有五個參數，定義如下：

In：資料集名稱
Strat：層數，亦即2乘2表格的數量(H)
Var1：行變數
Var2：列變數
Count：此變數是加權變數（weight），是非必要的變數。如果省略的話，程式會自動定義為「1」。

原始碼如下：

%macro MantelFleiss(In,Strat,Var1,Var2,Count);

%*-----------------------------------------------------------;
%* &In  - input data set                 ;
%* &Strat - stratification variable (e.g., site, race, etc.);
%* &Var1 - row variable                  ;
%* &Var2 - column variable                 ;
%* &Count - optional weight variable            ;
%*-----------------------------------------------------------;

%*reset output data sets;
PROC DATASETS nodetails nolist;
delete _counts _expect _mantel:;
RUN;
QUIT;

%*determine number of strata levels;
PROC SQL noprint;
select count(distinct &strat) into: stratflag from &In;
QUIT;
%put **Number of strata levels &stratflag;

%*get expected counts and totals in each partial table;
PROC FREQ data = &In;
tables &Strat*&Var1*&Var2 / expected outexpect out = _expect;
%if %nrbquote(&Count.) ne %then weight &Count.;;
ODS output CrossTabFreqs = _counts;
RUN;

%*get sum of expected cell counts in cell n11 in each partial table
(assign to macro var SUMEXP);
DATA _null_; 
set _expect end = eof;
by &Strat;
retain sumexp 0;
if first.&Strat then do;
sumexp = sumexp + expected;
end;
if eof then call symput('sumexp',put(sumexp,8.4));
RUN;
%put Sum of expected values in cells nh11: sum(mh11) = %cmpres(&sumexp);

PROC SORT data = _counts; 
by &Strat &Var1 &Var2;
RUN;


%*Subset to the needed totals output by PROC FREQ
number the &Strat, &Var1, and &Var2 for use below
0 - totals
1 - first &Strat/&Var1/&Var2
2 - second &Strat/&Var1/&Var2
...;

DATA _mantelf1;
%*only include marginal totals from ODS;
set _counts(where = (_type_ not in ('100' '111'))); 
by &Strat &Var1 &Var2;
retain &Strat._ 0 &Var1._ &Var2._;

%*reset for change in strata;
if first.&Strat then do;
&Var1._ = 0; &Var2._ = 0;
end;

%*get numeric version of strata variable;
if first.&Strat then &Strat._ = &Strat._ + 1;

%*define zero as total for any &Var1/&Var2;
if missing(&Var1) then &Var1._ = 0; 
else &Var1._ + 1;

if missing(&Var2) then &Var2._ = 0; 
else &Var2._ + 1;

keep &Strat._ &Var1._ &Var2._ &Strat. &Var1. &Var2. frequency;
RUN;

DATA _mantelf2;
set _mantelf1;
by &Strat._;

retain n11_ n1_1 n1_2;

if first.&Strat._ then do;
n11_ = 0;
n1_1 = 0;
n1_2 = 0; 
end;

%*get each total;
if &Var1._ = 1 then n11_ = frequency; %*for row  one in the hth strata;
if &Var2._ = 1 then n1_1 = frequency; %*for column one in the hth strata;
if &Var2._ = 2 then n1_2 = frequency; %*for column two in the hth strata;


%*define (nh11)L and (nh11)U for final computation;
max_ = max(0,n11_-n1_2);
min_ = min(n1_1,n11_);

if last.&Strat._;

RUN;

%*sum acoss min and max to get lower/upper bounds of criterion;
DATA _mantelf3(keep = mf_suml mf_sumu mflendpt mfrendpt mf_r);
set _mantelf2 end = eof;

retain mf_suml mf_sumu;

if _n_ = 1 then do;
mf_suml = 0; mf_sumu = 0;

end;

mf_suml = mf_suml + max_;
mf_sumu = mf_sumu + min_;

if eof then do;
mflendpt = &sumexp-mf_suml;
mfrendpt = mf_sumu-&sumexp;
mf_r   = min(mflendpt,mfrendpt);

put "Sum of lower bound (sum[(nh11)L]) = " mf_suml;
put "Sum of upper bound (sum[(nh11)U]) = " mf_sumu;
put "Left end-point: sum[mh11] - sum[(nh11)L] = " mflendpt;
put "Right endpoint: sum[(nh11)U] - sum[mh11] = " mfrendpt;
put "Mantel-Fleiss R = min[" mflendpt "," mfrendpt "] = " mf_r;

if mf_r>5 then put "Mantel-Fleiss criterion satisfied (since R > 5)";
else      put "Mantel-Fleiss criterion NOT satisfied (since R <= 5)";
output; 
end;

label
mf_suml = "Summation of Maximum (sum[(nh11)L])"
mf_sumu = "Summation of Minumum (sum[(nh11)U])"
mflendpt = "Left End-Point: sum[mh11] - sum[(nh11)L]"
mfrendpt = "Right End-Point: sum[(nh11)U] - sum[mh11]"
mf_r   = "R Value: min[sum[mh11] - sum[(nh11)L],sum[(nh11)U] -  sum[mh11]]";
RUN;
%mend;

範例一：
資料檔：

 PROC FORMAT;
value trtf
1 = 'TreatmentA'
2 = 'TreatmentB'
;

value sitef
1 = 'Arkansas'
2 = 'Indiana'
3 = 'Illinois'
;

value respf
1 = 'Present'
2 = 'Absent'
;
RUN;

DATA test;
do site = 1 to 3;
do treat = 1 to 2;
do resp = 1 to 2;
input counts @@;
output;
end;
end;
end;
label site = "Site" treat = "Treatment" resp = "Response";
format treat trtf. site sitef. resp respf.;
cards;
12 5 7 3
1 2 5 2 
0 9 5 6
;
RUN;

上述程式在表格內呈現的結果如下：

接著執行 macro：
%MantelFleiss(test,site,treat,resp,counts) ;

要特別注意的一點是，這個 macro 的結果並不是呈現在 output 視窗裡面，而是直接呈現在 log 視窗裡面，所以送出這行程式碼後，便可以直接在 log 視窗上看到結果，如下所示：

Sum of expected values in cells nh11: sum(mh11) = 20.5130
Sum of lower bound (sum[(nh11)L]) = 9
Sum of upper bound (sum[(nh11)U]) = 25
Left end-point: sum[mh11] - sum[(nh11)L] = 11.513
Right endpoint: sum[(nh11)U] - sum[mh11] = 4.487
Mantel-Fleiss R = min[11.513 ,4.487 ] = 4.487
Mantel-Fleiss criterion NOT satisfied (since R <= 5)

原文作者很好心地幫使用者做出該資料是否有滿足 Mantel-Fleiss criterion 的判斷句，所以就不用管上面的數據所代表的意義了，但若在寫 paper 時，可能還是要把一些數據呈現出來，但其實這些數據並沒有解釋上的明顯意義。

這個例子並沒有滿足 Mantel-Fleiss criterion，也就是說我們得使用 Fisher's exact test 來取代原本的卡方檢定。

我們再來看另一個有滿足 Mantel-Fleiss criterion 的例子。

範例二：

我們直接來看資料的表格結構。

所有參數設定都和範例一同樣，只是數據改變而已。執行 macro 後的結果如下：

Sum of expected values in cells nh11: sum(mh11) = 19.3630
Sum of lower bound (sum[(nh11)L]) = 9
Sum of upper bound (sum[(nh11)U]) = 30
Left end-point: sum[mh11] - sum[(nh11)L] = 10.363
Right endpoint: sum[(nh11)U] - sum[mh11] = 10.637
Mantel-Fleiss R = min[10.363 ,10.637 ] = 10.363
Mantel-Fleiss criterion satisfied (since R > 5)

從最後一句話得知這份資料是有滿足 Mantel-Fleiss criterion 的，所以我們便可以直接使用卡方檢定的結果。

CONTACT INFORMATION
Brandon Welch
Rho®, Inc.
Statistical Programming
6330 Quadrangle Dr., Ste. 500
Chapel Hill, NC 27517
Email: Brandon_Welch@rhoworld.com
Phone: 919-595-6339

User-Written DATA Step Functions

2009-11-24T23:23:00.000-04:00

原文載點：User-Written DATA Step Functions

這是一個 SAS V9.2 版新釋出的功能，類似 macro，可自行定義 SAS 沒有的函式（function），並可在不用重複呼叫的情況之下，自由使用於每個 data step 裡面。這個自行定義函式的功能，是由 PROC FCMP 來完成。

直接來看一個範例：

proc fcmp outlib=sasuser.funcs.trial;
function study_day(intervention_date, event_date);
if event_date < intervention_date then
return(event_date – intervention_date);
else
return(event_date – intervention_date + 1);
endsub;
run;

options cmplib=sasuser.funcs;
data _null_;
start = '15Feb2006'd;
today = '27Mar2006'd;
sd = study_day(start, today);
put sd=;
run;

這個範例可分為兩部分。第一部份使用 PROC FCMP 去製作一個名叫 study_day 的函式，而裡頭包含了兩個參數，一個叫 intervention _date，另一個叫 event_date。裡面的內容包含了一個 if-else- 的條件判斷式，並依照不同的條件傳回 event_date 減去 intervention_date 的天數或 event_date 減去 intervention_date 加一的天數。結束後使用 endsub 把 function 裡面的內容包起來，並將結果存到 sasuser 底下的 funcs 這個資料裡面，並命名為 trials。

等到要使用時，必須用 options cmplib 把 sasuser.funcs 叫進來，之後便可以無限次的在任何 data step 裡面使用 study_day 這個函式。即便重新開啟 SAS，也不用再跑一次 PROC FCMP 程序，只要執行 options cmplib=sasuser.funcs 這一行後，便可以直接使用 study_day 這個函式了。

如果想要製作 call routine 的函式，則需要在 PROC FCMP 裡面使用 subroutine statement 來製作。範例如下：

proc fcmp outlib=sasuser.funcs.math;
subroutine subA();
x = 5;
call subB();
put 'In subA:' x=;
endsub;
subroutine subB();
x = 'subB';
put 'In subB:' x=;
endsub;
run;

options cmplib=sasuser.funcs;
data _null_;
x = 99;
call subA();
put 'In DATA step: ' x=;
run;

在這個範例中，PROC FCMP 程序先用 subroutine 製作一個叫做 subA() 的 call routine，其功能是會定義 x=5，並且呼叫另一個 call routine 名為 subB()，最後會在 log 視窗列印出「In subA: x=」的字樣。但 subB() 這個 call routine 也是自創的，所以需要用另一個 subroutine statement 把他的功能寫好。完成後兩個 subroutine statement 都要用 endsub 把 subA() 和 subB() 打包好，如此一來，這個自創的 subA() 就可以在第二部份的 data step 裡面使用。

原文裡面還有一些更進階的寫法，但我覺得上面兩個範例已經足夠一般使用者來使用，所以就不花篇幅寫進階的教學了，有興趣的人可以自行下載原文來研究。

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author:
Jason Secosky
SAS Institute Inc.
SAS Campus Drive
Cary, NC 27513
919-677-8000
Jason.Secosky@sas.com

Stopping stepwise: Why stepwise and similar selection methods are bad, and what you should use

2009-11-24T20:22:00.000-04:00

原文連結：Stopping Stepwise: Why Stepwise Selection Methods are Bad and What you Should use Instead

從事統計諮詢工作多年，經常遇到這樣的問題：為什麼經過stepwise過後的線性迴歸結果與預期的相差很大，無法有效解釋。這類問題也反應出理論統計和實際應用上的出入。SAS 在選擇變數的過程中，只能藉由讀進去的數字來套用公式，但每個變數的背後意義，在 SAS 還沒演化成「天網」Skynet 之前，應該是猜不出來。以往的解決方式，通常是希望使用者根據經驗法則，從電腦挑出的變數裡面刪除研究內不會去討論到的變數，或者是從刪除掉的變數裡面強迫加入研究內會討論到的變數。但這些比較主觀的方法，總是有點失去科學性。這篇教學文件主要是回答諸多常用線性迴歸做研究的人內心的疑惑，並使用 PROC GLMSELECT 來解決這個問題。

早在2001年，Dr. Frank Harrell 就提出傳統的線性迴歸模式選擇法（含backward, forward and stepwise）的盲點，諸如：
。R平方有偏誤（過高）
。參數估計值的標準差太小
。F值或卡方統計量並沒有真正服從F分配或卡方分配
。參數估計值的信賴區間太小
。p-value有偏誤（太小）
。參數估計值（絕對值）太大
。共線性

有人可能懷疑會不會是因為樣本數不多導致 power 不夠，但這跟模式選擇是沒有太大的關係。原文作者甚至利用模擬的資料來比較 100 subjects 和 1000 subjects 在同一個模式下進行模式選擇，發現 1000 subjects 的結果也沒有來得特別好。

為了要克服模式選擇的問題，原文作者提出幾個方案：

(1) A FULL(ER) MODEL
顧名思義，就是完全不要去做模式選擇，直接用 full model 去分析資料。當然使用這種有點逃避現實的方法很明顯的會出現嚴重的共線性問題。以我個人的經驗來看，為了要規避模式選擇卻讓共線性影響分析的方式並沒有高明到哪裡去。因此我得建議是，共線性分析一定要做，並把他視為模式選擇的替代品。若某個變數不是研究的重點，也沒有統計上的顯著性，又跟其他獨立變數產生高度的相關性時，便可以直接砍除了。

(2) EXPERT KNOWLEDGE
這個方法就是我在第一段提到的，使用個人對研究主題和資料背景的認知，主觀地決定什麼變數該留下，什麼變數開刪除。雖然這並不是非常科學的方法，但我一直認為統計分析必須適度跳脫教科書的範疇，在整個統計分析的過程中，完全交由電腦去決定並不是一個百分之百有保障的作法，真正的決策還是得依靠人腦才行。

(3) MODEL AVERAGING
這個方法的觀念有點像是在稀釋掉使用 full model 時會造成的問題，其概念是建立數個可能的模式，再利用 AIC 去計算每個變數在每個模式的權重，最後用加權的方法把每個變數在不同模式下的估計值給平均起來。詳細的作法可以參考這本書：

Burnham, K. P. & Anderson, D. R. (2002), Model selection and multimodel inference, Springer, New York.

(4) PARTIAL LEAST SQUARES
這個方法主要是用來解決用主成分分析法（Principal Component Analysis, PCA）建立 PCA regression 以達成變數減少目的所產生不易解釋 PCA 變數的問題。

(5) LASSO
這個方法是借用脊迴歸（ridge regression）的參數估計方法來重新估計參數，詳細過程可見：

Trevor Hastie, R. T. & Friedman, J. (2001), The elements of statistical learning, Springer-Verlag, New York.

(6) LAR
Least Angle Regression 的縮寫，其概念是將所有參數中央化，讓所有的參數從零開始，然後再慢慢加上從殘差相關係數所產生的數據。詳細方法可見：

Bradley Efron, Trevor Hastie, I. J. & Tibshirani, R. (2004), ‘Least angle regression’, Annals of Statistics 32, 407–499.

(7) CROSS VALIDATION
類似 bootstrap 和 jackknife 的 resampling 方法，先決定要重複抽出 K 組資料，每一組資料裡面用 (K-1)/K 筆資料去估計模式，剩下的 1/K 筆資料拿來做驗證。

上述的方法中，LASSO 和 LAR 已經可以透過 PROC GLMSELECT 來完成，此程序的基本語法如下：

PROC GLMSELECT (options);
CLASS variable;
MODEL variable = (effects)/(options);
SCORE (data=dataset) (out=dataset)

此程式在執行模式選擇的主要關鍵語法是 MODEL statement 後面的 selection=。透過 selection= 這個 option 來指定選模方法和設定變數選取的準則（AIC, AICC, BIC,...,etc）。如果要計算 LASSO 和 LAR，只需要簡單地寫上「selection=LAR」和「selection=LASSO」即可。傳統的 backward、foreward 以及 stepwise 法仍舊可以使用在 PROC GLMSELECT 程序中，但提供額外參數供使用者自行設定，例如：

selection=forward(select=SL stop=validation)

selection=backward(select=SL)

selection=stepwise(select=SL SLE=0.1 SLS=0.08 choose=AIC)

這些參數其實就跟 PROC REG 裡面會用到的模式選擇參數一樣，在此不額外補充（其實原文裡面也沒有，是我自己從 PROC GLMSELECT手冊裡面添加的）。

使用 LASSO 和 LAR 法可提供傳統的模式選擇結果出現問題的替代方案，但這兩個方法也不是一本萬利，因為他們也各自有一些假設需要去檢驗。由於原文裡面沒有提到要檢定哪些假設，所以請自行參考原著。但原文作者也提到，這兩個方法仍舊比傳統的 stepwise 要好很多。若今後我有看到有專門在探討這兩種方法的 SAS 技術文件，會再以專文提出。

CONTACT INFORMATION
Peter L. Flom
515 West End Ave
Apt 8C
New York, NY 10024
peterﬂomconsulting@mindspring.com
(917) 488 7176

David L. Cassell
mathematical statistician
Design Pathways
3115 NW Norwood Pl.
Corvallis OR 97330

SAS® Abbreviations Are Your Friends; Use a Template Method to Code!

2009-11-18T19:01:00.004-04:00

原文載點：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」。

假設我們要在一個 data step 裡面使用 tranwrd function，結果打到一半.....

發現忘記 tranwrd function 所要用到的第二個和第三個語法，結果只好趕快去翻書或是求救於孤狗大神。但畢竟這個函式並不是每天在使用，所以隔了一個星期後突然又要用到，又得再一次去翻書或上網搜尋。為了避免這種情況，便可以開啟 SAS Abbreviation的功能，來把這些參數的定義註明起來，若下次再遇到時，SAS 便可以馬上提醒該函式的參數用法。方法如下：

首先，點選 SAS 工具列的「Tools」，再點選裡面的「Add Abbreviation」：

接著會開啟一個新的對話框：

其中，「Abbreviation」欄位可填上函式名稱 tranwrd，「Text to insert for abbreviation」欄位可填上 tranwrd function 的參數註解。如下所示：

輸入好點選 OK 就大功告成了。

下次當你要再次使用 tranwrd function 時，打出 tranwrd 則 SAS Abbreviation 會自動啟動，彈出之前新增過的 tranwrd function 註解，如下圖所示：

如果你是照著 tranwrd function 的格式去加註解，則看到註解跳出來時，可以直接按 Enter 或 Tab 鍵，SAS 會把註解直接輸入到 Editor 裡面。如果你只是要看看，或已經記得了正確完整的寫法，那就可以繼續輸入，不用管拿個字動彈跳出來的註解。

若你已經完全記住了這個語法，並且保證永遠都不會忘記，那就可以刪除這個註解。刪除的方法是，一樣到 SAS 工具列點選「Tools」，然後在下拉選單裡面點選「Keyboard Macros」，再點選「Macros...」：

此時會有一個新的對話視窗彈出，可以看到 tranwrd在一個框框裡面，點選 tranwrd 後再按 Delete 按鈕即可刪除。如果想要編輯的話，點 Edit 按鈕即可進行編輯，以下是一個例子。

假設我們今天有一個模版程式會經常使用，但還沒有需要用到 macro 來寫，如下所示：

我們也可以使用上述的方法把這個模版寫進 SAS Abbreviation 裡面，並且命名為 DataReport。但用了一陣子之後，發現 PROC REPORT 程序裡面的 Define statement 還滿常用到的，因此想要加入這一行，則一樣操作「Tools -> Keyboard Macros -> Macros -> 選擇 DataReport -> 點 Edit」，之後會彈出下面這個視窗：

視窗下方的 Modify 按鈕會亮起來，點擊下去會彈出第二個視窗：

使用者便可以在裡面將 Define的語法輸入到 PROC REPORT 裡面，完畢後按 OK 鍵。此後，當再一次於 Editor 視窗輸入 DataReport 並按下 Enter 鍵時，便會出現更新後的模版程式：

如果你經常需要使用到某個 SAS Abbreviation，而且你又是個超級懶人的話，SAS 提供熱鍵功能，讓你免於重複輸入，而只需要按個熱鍵便可把模版整個貼上數次。

定義熱鍵的方法，一樣是「Tools -> Keyboard Macros -> Macros」，然後在視窗內選擇想要設定熱鍵的 SAS Abbreviation，接著點擊右邊的「Assign Keys」按鈕，會彈出一個新的視窗：

接著在「Press new shortcut key」底下的空白欄位填上你愛用的熱鍵，如「Alt+E」等等。輸入完畢後，旁邊的「Assign」按鈕會亮起來，點擊下去後，便會看到這個熱鍵出現在上方的「Current keys」的空白欄位裡面。點擊下方的 OK 按鈕，便可以在 Keyboard Macros 視窗裡面看到設定好的熱鍵出現在 Keys 的欄位裡面：

關閉這個對話視窗後，該熱鍵就算是正式可以啟用了。

原文內還有教導如何在 SAS Enterprise 4.2 裡面使用 SAS Abbreviation 的功能，請有需要的人自行參考原文。

CONTACT INFORMATION
Elizabeth Ceranowski
Student Programs Manager
SAS Institute Inc.
SAS Campus Drive
Cary, NC 27513
Work Phone: (919) 531-9347
E-mail: Elizabeth.Ceranowski@sas.com

Examining Mediator and Moderator effect using Rural Women HIV Study

2009-06-18T16:57:00.002-04:00

原文載點：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;

第一個 PROC REG 得到 β=-0.98 (p-value=0.02)，因此 C1 達成。第二個 PROC REG 程序得到 β=0.143 (p-value=0.003)，因此 C2 也達成。至於 Eq3 的配適結果，得到 (β1, β2)=(-0.79, -0.44)，其 p-value 分別為 0.055 和 0.02，表示 M 在 Eq3 依舊顯著，但 IV 在 Eq3 變成不顯著了。所以基本上 C3 是達成了，但由於 β1 沒有降到接近 0，所以 C4 不算達成，因此 Spiritual Activity 只能稱是個 partial mediator，而不是 full mediator。

不過本篇技術文件並沒有繼續去做 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;

其中 sscopesa 是 tssqav 和 tcopesa 的交互作用項，這是由於 PROC REG 程序裡面不能使用 tssqav*tcopesa 這種語法來代表交互作用。因此在跑這個程式之前，一定要先用一個 data step 把交互作用項用另一個變數名稱給建立起來。

最後我們得到 IV*M 的估計參數 β=0.00175，其 p-value=0.5172 並不顯著，因此可以測得 Spiritual Activity 並不是 Moderator。

Using PROC SGPLOT for Quick High-Quality Graphs

2009-06-11T15:23:00.010-04:00

原文載點：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 DATA=sg.countries;
PANELBY Region;
REG X=NumParticipants Y=TotalMedals;
TITLE 'Number of Participants by Total Medals Won for Each Country';
RUN;

首先呼叫 PROC SGPANEL，再用 PANELBY statement 定義 Region 是類別變數。這功能就很像 CLASS statement。其餘的程式都和之前的 PROC SGPLOT 一樣，結果如下：

這份技術文件裡面有製作幾個表格，分別說明每種圖形要在 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

Updates to SAS® Power and Sample Size Software in SAS/STAT® 9.2

2009-04-27T16:44:00.003-04:00

原文載點：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;

語法的使用方法如下：

logistic：呼叫 logistic regression 下的 power analysis。
alpha：設定顯著水準。
vardist：設定 testpredictor 所指定的因變數 X 的分配，記得要再 vardist 後面用 ('....') 寫上在 testpredictor 一樣的變數名稱。此例是設定變數 Duration 的分配是 Normal(4, 1.5)。
testpredictor：因變數名稱
testoddsratio：因變數的 OR
responseprob：反應變數 Y 的機率
ntotal：樣本總數
power：power 設定值。

用法和一般的 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;

語法解釋如下：

onesamplefreq ci：計算係賴區間的方式。此例用 Wilcox score。另外五種分別為：Wald, continuity-corrected Wald, exact, Agresti-Coull, and Jeffreys
alpha：顯著水準
proportion：該二項變數的比例。
halfwidth：設定信賴區間一半的寬度。
ntotal：樣本總數。
probwidth：信賴區間

此例是計算一個比例為 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;

語法解釋如下：

twosamplewilcoxon：呼叫兩獨立樣本檢定的 power analysis。
alpha：顯著水準。
vardist：指定 group 變數的屬性。以此為例，兩個變數都是 ordinal，並且每個類別（-3 -2 -1 0 1 2 3）對應到的機率（.01 .04 .20 .50 .20 .04 .0）都必須要寫上。由於有兩個 group 要比較，所以必須宣告 vardist 兩次。
variables：指定兩個 group variable 的變數名，中間記得用「|」隔開來。
sides：標明該檢定單尾還是雙尾，此例用「u」表單尾。
ntotal：樣本總數。
power：power設定值。

此例兩個要比較的群組分別是 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

Combining, Combining, Combining, Splitting, Splitting, Splitting

2008-11-10T03:47:00.000-04:00

原文載點：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

Funny ^Stuff~ in My Code – Using ODS ESCAPECHAR

2008-11-01T01:17:00.001-04:00

原文載點：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

Performing Iterative Processes with the Macro Facility

2008-10-23T15:18:00.003-04:00

原文載點：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

A macro for nearest neighbor imputation

2008-10-21T00:37:00.007-04:00

原文載點：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;

語法使用方式如下：

INDATA：原始資料的名稱，包含其所用的 library
MISSVAR：含有 missing data 的變數（只能放一個）
RESPVAR：沒有 missing data 的變數（也只能放一個）
IDVAR：放 id number 或 case number，如果沒有的話必須要自己造一個
OUTDATA：插補過後輸出的完整資料

這邊使用一個實際範例來說明如何使用 %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

A Legend is Not Just a Legend

2008-10-14T17:30:00.003-04:00

原文載點：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

Names, Names, Names - Make Me a List

2008-09-11T18:14:00.006-04:00

原文載點：http://www2.sas.com/proceedings/forum2007/052-2007.pdf

如果你曾經讀過去年的這一篇文章「Text Utility Macros for Manipulating Lists of Variable Names」，而且也親自使用過該原文作者所提供的 text utility macro 的話，想必一定會對其強大的批次 rename 功能印象深刻。不過這個 macro 有個缺點，那就是在使用前要先定義這個東西：

%let orig_vars = a01a a01b d01 d02 t1_r t2_r;

這個 %let 的作用是把所有原始變數集中定義成 orig_vars 這個變數，然後再讓 text utility 來使用。不過，如果當變數很多時，我們需要手動去輸入，仍舊是相當耗時。而之前我曾經發過這篇文章「Automatically Renaming Common Variables Before Merging」，裡面是利用 PROC SQL 來抓出舊變數名稱。不過，如果舊變數名稱是有規則可尋，比方說 y1~y1000 時（ %let 不允許使用 y1-y1000 這種寫法），則這篇在 SAS Global Forum 2007 由 Ian Whitlock 發表的技術文件可以有很大的幫助。

這篇技術文件一共發佈了五個有用的 macro 讓使用者可以輕易操作大量變數名稱。以下就一一介紹：

。%ZIP
這個 macro 顧名思義就和壓縮檔 zip 一樣，可以組合數個不同的變數或 index 來產生有規則的變數。語法如下：

%zip (l1= , lv1= , sep1=, l2= , l2v= , sep2= , osep= )

變數解釋如下：
。l1= 第一列字串
。lv1= 取代第一列字串的字串（可省略不用）
。sep1= 用來區隔第一列字串的字元，預設值是空白（%str()）
。l2= 第二列字串
。l2v= 取代第二列字串的字串（可省略不用）
。sep2= 用來區隔第二列字串的字元，預設值是空白（%str()）
。osep= 用來區隔輸出字串的字元，預設值是空白（%str()）

比方說想要生成「Ax By Cz」這三個字串，可使用下列程式：

%zip(l1=A B C, l2=x y z);

l1 所表示的「A B C」會去抓在 l2 相對應位置的「x y z」，這便是 %zip 所呈現的功能。

。XPROD
這是可以拿來產生連續字尾或字首的 macro。語法如下：

%xprod ( l1= , lv1= , l2= , lv2=, sep1=, sep2=, osep=) )

使用方式和 %zip 一模一樣，但產生出來的效果完全不同。比方說想要產生「lib.w1 lib.w2 lib.w3」這三個字串，可發現這三個字串只有尾數不同，所以可以用 %xprod 來產生：

%xprod(l1=lib.w, l2=1 2 3);

換句話說，在 %xprod 裡面，l1 是拿來放共同字串，l2 是拿來放會改變的字串。他甚至可以做出更複雜的效果。假設想要產生「a_x a_y b_x b_y c_x c_y」這六個字串，我們可以想像這是一個 3X2 的矩陣：

      x   y
a_
b_
c_

而實際上 %xprod 就是用這種原理去做字串合併，因此我們只需要輸入：

%xprod(l1=a_ b_ c_, l2=x y);

就可自動產生那六個字串出來。

。RANGE
這是拿來生成連續數字的程式。語法如下：

%range ( to=, from= , step= , sep= )

變數解釋：
。to= 連續數字最後一位數，預設值是 1
。from= 連續數字的第一位數，預設值是 1
。step= 跳號的間格數，預設值是 1
。sep= 分隔連續數字的字元，預設值是空格（%str( )）

顧名思義，如果要生出「1 2 3 4 5」這串連續數字，則可使用：

%range(to=5);

。REPLACE
這是拿來取代字元的 macro，特別是使用在重新命名的動作裡面。語法如下：

%replace ( l= , lv= , code= , key= , lsep=, osep=)

變數解釋：
。l= 輸入字串
。lv= 取代上列字串的字串（可省略不用）
。code= 包含 key 所代表的符號變數的字串，通常是在程式中不會被變動的字串
。key= 用來當作取代「l」代表的字串的符號變數，預設值是 #
。lsep= 用來區隔輸入字串的字元，預設值是空白（%str()）
。osep= 用來區隔輸出字串的字元，預設值是空白（%str()）

舉例來說，如果在套入一組資料時要順便改變某些變數名稱，則程式碼為：

data ww;
 set w(rename=(x=__x y=__y z=__z);
run;

則可使用 %replace 來簡化：

data ww;
 set w(rename=(%replace(l=x y z, code=#=__#));
run;

這個程式會先把 l 所定義的 x y z 用依序一個 # 來取代，而 code 的作用就是把 x y z 一一套入「#=__#」的框架裡面，進而產生「x=__x y=__y z=__z」這三個字串出來。

另外還有更進階的寫法，就是把 %replace 放進另一個自創的 %rename 裡面，程式如下：

%macro rename ( list, pref=__) ;
%replace ( l=&list, code = # = &pref# )
%mend rename ;

%rename 的第一個參數就等於 %replace 的第一個參數，而 pref 就等於輸出字串裡面共有的字串（此例預設值是「__」）。因此原來的程式可以更簡化為：

data ww;
 set w(rename=(%rename(x y z)));
run;

這時「x y z」就等於「list」，也等於「l」，而 pref 所表示的「__」被丟進 %replace 的 code 參數裡面，則 code 就變成「code= #=__#」，這就便回原來只用 %replace 的那個程式碼了。

%replace 也可以使用在變數轉換上面。比方說上面這個例子，在 rename 後要把新變數轉換成 best 32. 的格式。如果不用額外的 macro 則語法是：

data ww;
  set w(rename=(%rename(x y z)));
 x=input(__x, best32.);
 y=input(__y, best32.);
 z=input(__z, best32.);
run;

同樣地，我們可以另外自製一個命為 %CHAR2NUM 的 macro，然後套用 %replace 去寫。方法如下：

%macro char2num ( list , pref = __ ) ;
%replace ( l=&list
, code= %str(# = input(&pref#,best32.);)
)
%mend char2num ;

這個新 macro 的參數用法和 %rename 一樣，但我們來看看 %replace 在裡面發會什麼功用。l 去呼叫 list 所定義的字串，而 code 生出「#=input(&pref#,best32.);」，其中 # 就是暫時替代 l/list 參數所定義的字串，pref 是利用預設值「__」，則 code 其實就是產生下面這三行：

x=input(__x, best32.);
y=input(__y, best32.);
z=input(__z, best32.);

因此原本的程式可以改寫成：

data ww;
   set w(rename=(%rename(x y z)));
 %char2num(x y z);
run;

。QT
這個 macro 只是很簡單地幫每個字元加上雙引號。語法如下：

%qt ( l=, lv= , lsep= %str( ), qt = %str(%"), osep=%str( ) )

變數解釋：
。l= 輸入字串
。lv= 取代上列字串的字串（可省略不用）
。lsep= 用來區隔輸入字串的字元，預設值是空白（%str()）
。qt= 幫每個字元左右加上的符號，預設值是雙引號（"）
。osep= 用來區隔輸出字串的字元，預設值是空白（%str()）

因此，如果要生成出「"a" "b" "c"」只需要簡單地使用：

%qt(l=a b c);

以下是各 macro 的原始碼：
。ZIP

%macro zip
( l1= /* first list */
, lv1= /* external variable override for first list */
, sep1=%str( ) /* separator between the joined elements */
, l2= /* second list */
, lv2= /* external variable override for second list */
, sep2=%str( ) /* separator between the joined elements */
, osep=%str( ) /* separator between new elements */
) ;
/* %zip ( l1= a b , l2= c d ) produces ac bd
so does
%let list1 = a b ;
%let list2 = c d ;
%zip (lv1=list1, lv2=list2)
If lists do not have same length shorter length used and
warning to the log. Empty lists result in empty list and
no message.
If the LV options are used then L1, L2, and ZIP_: should be
avoided for variable names.
*/
%local zip_i zip_1 zip_2 zip_list ;
%if %length(&lv1) = 0 %then
%let lv1 = l1 ;
%if %length(&lv2) = 0 %then
%let lv2 = l2 ;
%do zip_i = 1 %to &sysmaxlong ;
%let zip_1 = %qscan(%superq(&lv1) , &zip_i, &sep1 ) ;
%let zip_2 = %qscan(%superq(&lv2) , &zip_i, &sep2 ) ;
%if %length(&zip_1) = 0 or %length(&zip_2) = 0 %then
%goto check ;
%if &zip_i = 1 %then
%let zip_list = &zip_1&zip_2 ;
%else
%let zip_list = &zip_list&osep&zip_1&zip_2 ;
%end ;
%check:
%if %length(&zip_1) > 0 or %length(&zip_2) > 0 %then
%put WARNING: Macro ZIP - list lengths do not match - shorter used. ;
%unquote(&zip_list)
%mend zip ;

。XPROD

%macro xprod
( l1= /* first list */
, lv1= /* external variable override for first list */
, sep1=%str( ) /* separator between elements of first list */
, l2= /* second list */
, lv2= /* external variable override for second list */
, sep2=%str( ) /* separator between elements of second list */
, osep=%str( ) /* separator between elements of new list */
) ;
%local xp_i xp_j xp_1 xp_2 xp_list ;
%if %length(&lv1) = 0 %then
%let lv1 = l1 ;
%if %length(&lv2) = 0 %then
%let lv2 = l2 ;
%do xp_i = 1 %to &sysmaxlong ;
%let xp_1 = %qscan(%superq(&lv1), &xp_i, &sep1) ;
%if %length(&xp_1) = 0 %then %goto endloop1 ;
%do xp_j = 1 %to &sysmaxlong ;
%let xp_2 = %qscan(%superq(&lv2), &xp_j, &sep2) ;
%if %length(&xp_2) = 0 %then %goto endloop2 ;
%if &xp_i = 1 and &xp_j = 1 %then
%let xp_list = &xp_1&xp_2 ;
%else
%let xp_list = &xp_list&osep&xp_1&xp_2 ;
%end ;
%endloop2:
%end ;
%endloop1:
%unquote(&xp_list)
%mend xprod ;

%RANGE

%macro range /* second more efficient version due to Chang Chung */
( to=1 /* end integer value */
, from=1 /* starting integer value */
, step=1 /* increment integer */
, osep=%str( ) /* sparator between integers */
) ;
/*
return sequence of integers starting at &FROM going to &TO
in steps of &step
%range(to=5) produces 1 2 3 4 5
*/
%local rg_i ;
%do rg_i = &from %to &to %by &step ;
%if &rg_i = &from %then
%do;&rg_i%end ;
%else
%do;&osep&rg_i%end ;
%end ;
%mend range ;

。REPLACE

%macro replace
( l= /* value list */
, lv= /* external variable override for value list */
, lsep=%str( ) /* separator between values */
, code= /* block of code containing symbolic variable */
, key=# /* symbolic variable to replace (#abc# etc.) */
, osep=%str( ) /* separator between new elements */
/* may be %str(;) when code is statement */
/* if so remember to add closing semicolon */
) ;
%local rg_i rg_w rg_list ;
%if %length(&lv) = 0 %then
%let lv = l ;
%if %length(%superq(&lv)) = 0 /*or %index(%superq(code),&key) = 0*/ %then
%do ;
%let rg_list = %superq(code) ;
%goto mexit ;
%end ;
%do rg_i = 1 %to &sysmaxlong ;
%let rg_w = %qscan(%superq(&lv),&rg_i,&lsep) ;
%if %length(&rg_w) = 0 %then %goto mexit ;
%if &rg_i = 1 %then
%let rg_list = %sysfunc(tranwrd(%superq(code),&key,&rg_w)) ;
%else
%let rg_list =
&rg_list&osep%sysfunc(tranwrd(%superq(code),&key,&rg_w)) ;
%end ;
%mexit:
%unquote(&rg_list)
%mend replace ;

Contact information
Ian Whitlock
29 Lonsdale Lane
Kennett Square, PA 19348
Ian.Whitlock@comcast.net

Using 22 Easy Tricks with ODS to Generate Colorful Reports

2008-08-23T16:13:00.002-04:00

原文載點：http://www.lexjansen.com/mwsug/2007/Tutorials/T10-2007.pdf

ODS 系統在 SAS V9.1 版中已經逐漸發揮他的威力，到了 V9.2 版後更是往上提升了一個階層。不過很多人以為用用 ODS 把圖或報表直接輸出就已經是最方便的，殊不知他可以結合最簡單的 PROC PRINT 程序將輸出結果強化。這一篇在 NESUG 2007 發表的技術文件詳細記載了如何讓 ODS 搭配 PROC PRINT 來產生更美觀的報表，也同時揭露了一些在 SAS 手冊裡都沒有介紹的指令。

簡單來講，如果要用 PROC PRINT 程序加上 ODS 來列印報表，在完全沒有額外設定，只使用預設值的情況下，語法大概就諸如此款：

然後得到的結果是：

如果你覺得上面這種呈現方式有點像是墓碑上刻出來的文字，而想要換成別的樣式或顏色，那麼就得靠 PROC PRINT 程序手冊裡面都很少提到（至少 SAS online documentation 沒有寫）的 style option！

style option 基本的撰寫格式如下：

以這個圖例說明：

利用四個 style 指令將 Obs 標題顏色、表格內資料背景顏色、Obs 列顏色以及底部 total 底色一次進行修改。

接下來又一步一步來介紹不同 style 指令的用法。

1. 更改標題列顏色

如果想要將標題列著成黃色，如下圖所示：

則可以用下面這道程式來實現：

其中，黃色標記處 style(header)={background=yellow}; 即是讓標題列變成黃色的關鍵。

2. 更改總和列顏色和字型和大小

PROC PRINT 裡面有個 SUM statement 可以計算每個欄位的總和，如想要更改此處的背景顏色、字型樣式和字型大小，如下圖所示：

則可以使用這道程式：

黃色標明處 style (grandtotal)={background=grayee font=(Arial) font_size=6}; 即是更改由橘色標明處 SUM statement 所算出來的兩個總和欄位的背景顏色(grayee)、字型(Arial)和字型大小(font_size)。如果沒有 SUM statement，則黃色標明處的語法將不會有任何作用。

3. 更改表格資料欄背景顏色

如想要改變表格資料的背景顏色，如下圖所示：

則語法如下：

黃色標明處 style(data)={background=grayee} 即是用來設定資料欄位背景顏色。附帶一提，標題列會自動斷成兩行是利用紅色標明處的 split='*' 進行設定，然後在橘色標明處 label statement 裡面將標題要斷行的地方打上 * 即可。

4. 更改分類總和欄位背景顏色

如果有使用 by statement 來計算不同分類欄位的總和，則其背景顏色也可以修改，如下圖所示：

語法如下：

此程式由於有使用 by origcity; 和 sum statement 來計算每個程式 sales2005 和 sales2006 的總和，所以要修改分類總和欄位的背景顏色，就得用黃色標明處的 style(total)={background=light green}; 來設定。

5. 更改特定欄位字體顏色

若想指定某個欄位的字型顏色，如下圖所示：

可使用下列程式：

黃色標明處的 var month 後面加上 /style={foreground=red}; 就可以更改 month 這個欄位字型顏色。此外，一開頭的 SAS logo 是用 title3 後面的 "^S=preimage='c:\ben\qp\SAS_logo.bmp'}" 所貼上。

總結：作者有坦承這篇文件沒有把 22 個 trick 都用上，不過會在 conference 裡面說明。不過由於沒有她在會議裡面的任何影音檔，所以無法得知其他沒有在文件裡面的小技巧是什麼。但上述這些語法已經足夠讓表格變的更漂亮，因此有興趣的人可以玩玩不同組合最出現什麼情況。

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author at:
Ben Cochran
The Bedford Group
Raleigh, NC 27607
Office: 919.831.1191
Fax: 919.831.1191
Email: bedfordgroup@nc.rr.com
Web: www.bedford-group.com

Adding One Value to All Observations

2008-08-14T16:33:00.005-04:00

原文載點：http://www.nesug.info/Proceedings/nesug07/cc/cc45.pdf

由於最近在工作上太多人問到這個問題，所以我特別找出這一篇文章來說明。在 SAS 裡面，要在每筆資料後面加上一個特定值，依據資料排列格式的不同，有許多種作法。有得很簡單，只需要在 data step 加上一行即可，有得很複雜，得動用到 merge statement 或 proc sql 才能解決。這一篇在 NESUG 2007 年發表的技術文件列出了一些解法可供參考。

本例使用這個樣本資料：

data temp;
input Obs ID EARNHR @@;
cards;
1 1 .
2 2 5.75
3 3 5.75
4 4 6.30
5 5 6.75
6 6 8.00
7 7 8.25
8 8 9.60
9 9 10.05
10 10 12.80
;

如果打算算出 EARNHR 的中間值，然後讓他變成一個獨立的變數放在原資料 temp 的最後面。最常使用的方法是先使用 PROC MEANS 把中間值的結果另存到新的資料檔裡面：

proc means data = temp noprint ;
  var EARNHR ;
  output out = add (drop = _TYPE_ _FREQ_) median(EARNHR) = MEDIANEARNHR;
run ;

然後再用 merge 的方法放回原資料檔：

data temp3 ;
  set temp ;
  if _N_ = 1 then set add ;
run ;

接下來要教兩個進階用法。由於第一個方法會生出一個無用的資料檔 MEDIANEARNHR，雖然他只是放在暫存檔，關掉 SAS 後就會消失。但使用 PROC SQL 就根本連這個檔都不會產生。方法如下：

proc sql ;
  select mean(EARNHR) as AVGEARNHR format = 8.2 into :MEANEARNHR
  from temp ;
quit;

data temp ;
  set temp ;
  mean= &MEANEARNHR ;
run ;

本例改成放 EARNHR 的平均值 MEANEARNHR。先用 PROC SQL 把 AVEARNHR 算出來，放進一個叫做 MEANEARNHR 的巨集變數裡面，然後在把他命名成mean放進原來的資料內。這還有很多其他的用途。比方說如果想要挑出 EARNHR 小於平均值的樣本，則可以用這個 data step 來挑出：

data temp2 ;
  set temp ;
  where EARNHR lt &MEANEARNHR ;
run ;

挑出來的結果是：

Obs ID EARNHR
1 2 5.75
2 3 5.75
3 4 6.30
4 5 6.75
5 6 8.00

第三種方法則是用 SYMPUT 函式來定義巨集變數。此例，先用 PROC MEANS 去計算 EARNHR 的中間值，並存在 add2 裡面。

proc means data = temp noprint ;
  var EARNHR ;
  output out = add2 (drop=_TYPE_ _FREQ_) median(EARNHR) = MEDIANEARNHR;
run ;

然後開一個虛無的資料集「_null_」（這個資料連在work library都不會出現），然後把 add2 裡面的中間值變數名稱用 SYMPUT 定義成巨集變數。

data _null_ ;
  set add2 ;
  call symputx ('MDNEARNHR', MEDIANEARNHR) ;
run ;

最後就可以利用這個巨集變數在新的 data step 裡面做很多動作。無論是新增變數還是去挑資料都相當容易。

data temp;
  set temp ;
  median=&MDNEARNHR;
  where EARNHR ge 0 and EARNHR lt &MDNEARNHR ;
run ;

個人傾向使用 PROC SQL 是因為他能夠減少無謂的資料檔輸出，只是需要相當熟悉他的語法才能夠用上手。

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author at:
Anne W. Warren
MDRC
16 East 34th Street, 19th floor
New York, NY 10016
Work Phone: 212.340.8661
Fax: 212.684.0832
Email: anne.warren@mdrc.org
Web: www.mdrc.org

Reliability analysis: Calculate and Compare Intra-class Correlation Coefficients (ICC) in SAS

2008-06-23T17:16:00.001-04:00

原文載點：http://www.nesug.info/Proceedings/nesug07/sa/sa13.pdf

在 reliability analysis 裡面，intra-class correlation coefficients (ICC) 是一個相當重要的統計指標。不過長久以來，SAS 內部並沒有可以直接計算 ICC 的程序。雖然有許多線上的 ICC 計算器，不過仍有許多人研發出在 SAS 上運作的 ICC macro 供 SAS 使用者使用。Li Lu 和 Nawar Shara 在 2007 年的 NESUG 發表了一篇自行撰寫的 ICC macro 技術文件，不但可用來計算 ICC 本身，還可以求出其信賴區間。

原始碼如下：

%macro Icc_sas(ds, response, subject);
 ods output OverallANOVA =all;
 proc glm data=&ds;
     class &subject;
     model &response=&subject;
 run;
 data Icc(keep=sb sw n R R_low R_up);
     retain sb sw n;
     set all end=last;
     if source='Model' then sb=ms;
     if source='Error' then do;sw=ms; n=df; end;
     if last then do;
     R=round((sb-sw)/(sb+sw), 0.01);
     vR1=((1-R)**2)/2;
     vR2=(((1+R)**2)/n +((1-R)*(1+3*R)+4*(R**2))/(n-1));
     VR=VR1*VR2;
     L=(0.5*log((1+R)/(1-R)))-(1.96*sqrt(VR))/((1+R)*(1-R));
     U=(0.5*log((1+R)/(1-R)))+(1.96*sqrt(VR))/((1+R)*(1-R));
     R_Low=(exp(2*L)-1)/(exp(2*L)+1);
     R_Up=(exp(2*U)-1)/(exp(2*U)+1);
     output;
     end;
 run;
 proc print data=icc noobs split='*';
     var r r_low r_up;
     label r='ICC*' r_low='Lower bound*' r_up='Upper bound*';
     title 'Reliability test: ICC and its confidence limits';
 run;
%mend;

包含三個參數：

ds：資料名稱（含 library）
response：反應變數名稱
subject：ID 名稱

使用方式相當簡單。如：

%Icc_sas(ds=one, response=cimt, subject=subject_id);

由於原文沒有提供 one 這個資料集的原始資料，所以不知道跑出來的結果如何。有興趣的人可以拿自己的資料測試，並與其他線上 ICC 計算器的結果進行比對。如有發現異常現象請回報給原作者。

CONTACT INFORMATION
Your code requests, comments and questions are valued and encouraged. Contact the author at
Li Lu
Dept of Epidemiology and Statistics
MedStar Research Institute
6495 New Hampshire Avenue
Hyattsville, MD 20783
(301) 560-7313
li.lu@medstar.net

ALPHA_CI: A SAS® Macro for Computing Confidence Intervals for Coefficient Alpha

2008-06-05T22:09:00.001-04:00

原文載點：http://www2.sas.com/proceedings/forum2008/230-2008.pdf

在進行效度分析的時候，我們通常都會需要去求一個叫做 Cronbach's alpha的值，而這個值越大代表整個問卷設計的可信度越高。在 SAS 中可利用 PROC CORR 程序輕鬆地計算出這個數值出來，但卻沒有辦法去估計他的信賴區間。事實上，已有不少的學者發表相關文獻來估計 Cronbach's alpha 的信賴區間，而 Jeffrey D. Kromrey, Jeanine Romano 和 Susan T. Hibbard 收集了八個不同的估計方法，並將其整合在一個 macro 程式裡面。

相關的理論請參考原文以及其附錄的參考文獻。

整個 macro 的原始程式碼如下：

%macro ALPHA_CI (dsn = _last_,n_items = 100,confidence = .95);
proc iml symsize = 500;
start Bonett(N1,items,rxx,confidence,lowerCI,upperCI);
* +--------------------------------------------------------------------------------------+
Bonett 2002 method of CI computation
Inputs to subroutine are number of examinees (N1), number of items (items),
sample value of coefficient alpha (rxx), and desired level of confidence
(Confidence).
Outputs are the upper and lower limits of the confidence interval (upperCI, lowerCI)
+--------------------------------------------------------------------------------------+;
critZ = -1#probit((1-Confidence)/2);
transBon_z = log(1-abs(rxx));
if rxx <0 transbon_z =" transBon_z" var_z =" (2#items)/((items" se_bon_z =" sqrt(var_z);" lowerci =" 1" upperci =" 1" lowerf =" (1-Confidence)/2;" upperf =" Confidence" fb =" FINV(upperF,N1-1,(N1-1)#(items-1),0);" fa =" FINV(lowerF,N1-1,(N1-1)#(items-1),0);" lowerci =" 1" upperci =" 1" critz =" -1#probit((1-Confidence)/2);" transfish_z =" 0.5#log((1+abs(rxx))/(1-abs(rxx)));" transfish_z =" transfish_z" se_fish_z =" sqrt((n1-3)##-1);" upperci =" (exp(2#(transfish_z" lowerci =" (exp(2#(transfish_z" critz =" -1#probit((1-Confidence)/2);" transhawh_z =" (1-rxx)**(1/3);" var_z =" ((18#items)#(n1" se_hawh_z =" sqrt(var_z);" c_star =" ((9#n1" lowerci =" 1" upperci =" 1" critz =" -1#probit((1-Confidence)/2);" bigone =" j(N1," means =" ((bigone`)*out3pl)/N1;" score3pld =" out3pl-(bigone*means);" itemcov =" (1/(N1-1))*((SCORE3PLd`)*SCORE3PLd);" one =" j(items,1);" jtphij =" (one`)*itemcov*one;" trphisq =" trace(itemcov*itemcov);" trsqphi =" (trace(itemcov))**2;" jtphisqj =" (one`)*(itemcov*itemcov)*one;" omega =" jtphij*(trphisq+trsqphi);" omega =" omega-(2*(trace(itemcov))*jtphisqj);" omega =" (2/(jtphij**3))*omega;" s2 =" (items**2)/((items-1)**2);" s2 =" s2*omega;" se =" sqrt(s2/N1);" lowerci =" rxx-(CritZ*se);" upperci =" rxx+(CritZ*se);" lowerf =" (1-Confidence)/2;" upperf =" Confidence" fb =" FINV(upperF,n1#(items-1),n1,0);" fa =" FINV(lowerF,n1#(items-1),n1,0);" upperci =" 1" lowerci =" 1" critz =" -1#probit((1-Confidence)/2);" upperci =" 1" lowerci =" 1" parm =" J(1,7,.);" nsub =" NROW(data);" nvar =" NCOL(data);" nv2 =" nvar*(nvar+1)/2;" rsub =" 1." rs1 =" 1." vrat =" nvar" mean =" data[+,]" quant =" -1#probit((1-prob)/2);" cov =" rs1" vvar =" vecdiag(cov);" summat =" cov[+,+];" sumvars =" vvar[+];" sumcovs =" .5*(summat" wcv =" cov" wvv =" vecdiag(wcv);" sumvar2 =" wvv[+];" summat2 =" wcv[+,+];" alpha =" vrat" tmp =" 2." t1 =" summat" t2 =" summat" t3 =" 2.*sumvars" nase =" tmp*(t1" nase =" sqrt(rsub*nase);" dwrtvar =" -2." dwrtcov =" vrat" jac =" J(nvar,nvar,dwrtcov);" j=" 1" trac =" 0.;" isub="1" v =" data[isub,]" wcv =" jac" tmp =" wcv[+,+];" trac =" trac" nnase =" sqrt(rsub*rs1*trac);" i =" 1" i =" 1" item_scores =" temp_x;"> 1 then do;
Item_Scores = Item_Scores || temp_x;
End;
%end;
confidence = &confidence;
* +---------------------------------+
Computation of Cronbach Alpha
+---------------------------------+;
N_items = ncol(Item_Scores);
N_obs = nrow(Item_Scores);
mu1 = J(N_items,1,0);* ;
var = J(1,N_items,0);* ;
do k = 1 to N_items;
do i=1 to N_obs;
mu1[k,1] = mu1[k,1] + Item_Scores[i,k];* ;
end;
var[1,k]=(mu1[k,1]/N_obs)*(1 - mu1[k,1]/N_obs);
end;
sumvar=0;
do k = 1 to N_items;
sumvar = sumvar + var[1,k];* ;
end;
rowsum = J(N_obs,1,0);
do p = 1 to N_obs;
do k = 1 to N_items;
rowsum[p,1]=rowsum[p,1] + Item_Scores[p,k];
end;
end;
sumscore = 0;
sumscore2 = 0;
do p = 1 to N_obs;
sumscore = sumscore + rowsum[p,1];
sumscore2= sumscore2 + rowsum[p,1]##2;
end;
vartotal= (sumscore2-(sumscore##2/N_obs))/(N_obs);
* +------------------------------------------------------------+
Be sure we have some score variance before going any further
* +------------------------------------------------------------+;
if vartotal <=0 then do; print 'Total Score Variance =' vartotal; print 'Check the Data!'; end; if vartotal > 0 then do;
rxx = (N_items/(N_items -1))*((vartotal- sumvar)/vartotal);
if rxx < rxx =" .00001;"> .99 then rxx = .99;
* +-------------------------------------------+
Call subroutines for confidence intervals
+------------------------------------------+;
run Bonett(N_obs,N_items,rxx,Confidence,lowerBonett,upperBonett);
run Feldt(N_obs,N_items,rxx,Confidence,lowerFeldt,upperFeldt);
run Fisher(N_obs,N_items,rxx,Confidence,lowerFisher,upperFisher);
run HW(N_obs,N_items,rxx,Confidence,lowerHw,upperHW);
run ID(N_obs,N_items,rxx,Item_Scores,confidence,lowerID,upperID);
run KF1(N_obs,N_items,rxx,Confidence,lowerKF1,upperKF1);
run KF2(N_obs,N_items,rxx,Confidence,lowerKF2,upperKF2);
parm = scalpha(Item_Scores,confidence);
lowerADF = parm[1,6];
upperADF = parm[1,7];
end; * end the 'if vartotal > 0' conditional;
rxx = round(rxx,.001);
* +---------------------------------+
Printed macro output
+---------------------------------+;
file print;
put @1 'Confidence Intervals for Coefficient Alpha' /
@1 '-------------------------------------------------------' /
@1 'Level of Confidence:' @30 Confidence /
@1 'Number of Observations:' @30 N_obs /
@1 'Number of Items:' @30 N_items /
@1 'Sample Value of Alpha:' @30 rxx //
@1 'Method' @40 'Lower' @50 'Upper' /
@1 '------------------------------' @40 '-----' @50 '-----' /
@1 'Bonett' @40 lowerBonett 5.3 @50 upperBonett 5.3 /
@1 'Feldt' @40 lowerFeldt 5.3 @50 upperFeldt 5.3 /
@1 'Fisher' @40 lowerFisher 5.3 @50 upperFisher 5.3 /
@1 'Hakstain & Whalen' @40 lowerHW 5.3 @50 upperHW 5.3 /
@1 'Iacobucci & Duchachek' @40 lowerID 5.3 @50 upperID 5.3 /
@1 'Koning & Frances Exact' @40 lowerKF1 5.3 @50 upperKF1 5.3 /
@1 'Koning & Frances Asymptotic' @40 lowerKF2 5.3 @50 upperKF2 5.3 /
@1 'Asymptotic Distribution Free' @40 lowerADF 5.3 @50 upperADF 5.3 /
@1 '-------------------------------------------------------';
quit;
%mend;

以上長達數百行的原始碼，我們不需要特別瞭解內部構造，只要知道該怎樣用即可。整個 macro 包含三個參數：

dsn：使用的 data set 名稱。如果沒有特別指定，則這個 macro 會直接使用最後一次使用過的 data set。
n_items：變數數量。預設值是 1000。
confidence：信賴區間大小。預設值是 0.95。

以下為一個小小使用範例。假設有個問卷裡面有五道題目，而每題的選項只有兩種：1 和 0。可利用下面這個 data step 完成資料輸入：

data one;
input @1 (X1 - X5)(1.);
cards;
10000
11000
11100
11110
11111
00000
10000
11000
11100
11110
11111
00000
;

欲算出 95%、99% 和 90% 的信賴區間，則只要輸入下面這三行即可：

%ALPHA_CI (dsn = one,n_items = 5);
%ALPHA_CI (dsn = one,n_items = 5,confidence=.99);
%ALPHA_CI (dsn = one,n_items = 5,confidence=.90);

以第一個 95% 信賴區間為例，報表如下：

由於沒有方法去檢定哪個信賴區間是最適合的，所以作者列出下面這個表格，讓使用者根據自己的樣本大小去評斷哪種方法所算出來的信賴區間有較高的覆蓋機率（coverage probability）。

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Please contact Jeff Kromrey at:
University of South Florida
4202 East Fowler Ave. EDU 162
Tampa, FL 33620
Work Phone: 813-974-5739
Fax: 813-974-4495
Email: kromrey@tempest.coedu.usf.edu