Analysis of Method Comparison Studies Using SAS®

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

這一篇是 Mohamed Shoukri 在 2003 年在 SUGI 28 發表的一篇關於用 SAS 檢測兩種不同測量方法差異的文章。在某些跟 bio 相關的研究領域中，通常會有不同的測量方式來針對某一種特殊的計量。每一種測量方式當然都各有其發展的時空背景和優缺點，但使用的人總希望這些測量方式對結果來說不會造成太大的差異。最簡單最直覺的方法就是用 T-test 和 Pearson correlation，但陸續有學者發現這並不適用所有的情況，Altman 和 Bland 分別在 1983 年和 1986 年發表兩篇論文來指出一些問題，並提出新的統計方法來檢定。

Bland 和 Altman 的基本概念是去算兩種不同方法所測出的差異（Vi=Xi-Yi），然後去繪製其和兩種方法的平均值（Ui=(Xi+Yi)/2）的散佈圖。如果描繪點都在零附近打轉的話，那就表示兩種方法應該沒太大的偏差。至於比較有系統的統計檢定方法來測兩種是否有顯著不同，則根據 Shukla 在 1973 年發表的方法，利用 t 分配來檢定下面這個值：



上述的值服從自由度為 k-2 的 t 分配，而 r 是 X 和 Y 的相關係數。

此外，V 的上下限可用下列公式來表示：

Upper Limit(V) = mean(V) + 2*SD(V)
Lower Limit(V) = mean(V) - 2*SD(V)

Shoukri 利用 Ontario Veterinary College 裡的一份資料來跑這個分析。共有下列四個變數：

slide = 鑑定樣本（15 個）
rater = 臨床醫師（2 人）
reading = 醫師判讀鑑定樣本的次數（2 次）
count = 量測結果

由此可知，每位醫生判讀十五份鑑定樣本，每份鑑定樣本都被判讀兩次。

THE TRANSPOSE PROCEDURE

第一步驟需要將資料從單變量架構轉置成多變量架構。也就是從：



轉置成



可利用 PROC TRANSPOSE 程序來完成，如下所示：

proc transpose data=rater out=ratertrans (rename=(rater1=rater11 rater2=rater12 rater3=rater21 rater4=rater22) drop=_NAME_) prefix=rater;
var count;
by slide;
run;

接著從轉置後的資料來計算不同醫生對兩次判讀結果的平均值以便計算 U 和 V。程式如下：

data newrater;
set ratertrans;
* 1. Compute the averages of the readings made by each rater;
x = (rater11 + rater12)/2;
y = (rater21 + rater22)/2;
* 2. Compute the sum and the difference of the averages of the readings;
u = x + y;
v = x - y;
run;

THE CORR PROCEDURE

接著，就可以用 PROC CORR 來計算 U 和 V 的相關係數。結果如下：



報表顯示相關係數是 0.362，p-value 是 0.184 大於 0.05 顯著水準。因此可以解釋為兩次的判讀有相當程度的準確性。

THE MEANS PROCEDURE

為了要計算 V 的上下限，可以用 PROC MEANS 加上一個 Data procedure 來完成：

proc means data=newrater;
var v;
output mean=meanv std=sdv out=reflines;
run;

data lines;
set reflines;
* 1. Compute the upper and lower limits of v;
up = meanv + (2*sdv);
lo = meanv - (2*sdv);
* 2. Convert the variables, up and lo, into macro variables;
call symput('upper',up);
call symput('lower',lo);
run;

特別一提的是，Data procedure 內的 symput 指令是將 upper 和 lower 變數轉成 macro 變數並重新命名為 up 和 lo。

THE GPLOT PROCEDURE

最後，使用 PROC GPLOT 指令把點和上下限描繪出來就大功告成了。

goptions reset=all;
symbol1 v=plus c=blue;
proc gplot data=newrater;
plot v*u / vref=&upper &lower lvref=2 cvref=red;
run;
quit;



CONCLUSION

一、一般來說，如果有點超出來上下限，就可以判定是離群值。
二、如果（U,V）沒有特殊的趨勢，可以解釋兩個醫生的判讀具有一致性。
三、如果（U.V）有特殊趨勢，則可能要將原始資料做轉換，然後重新跑上述步驟。通常 log-transform 就可以滿足。

CONTACT INFORMATION
Your comments and questions are valued and encouraged.
Contact the author at:

Samia Hashim
King Faisal Specialist Hospital & Research Center
Biostatistics, Epidemiology & Scientific Computing Department
MBC 03-BESC
PO Box 3554
Riyadh 11211
KSA
Work Phone: (+9661) 464-7272 Ext. 39229
Fax: (+9661) 442-4542
Email: samia@kfshrc.edu.sa
Website: www.kfshrc.edu.sa
Mohamed Shoukri
King Faisal Specialist Hospital & Research Center
Biostatistics, Epidemiology & Scientific Computing Department
MBC 03-BESC
PO Box 3554
Riyadh 11211
KSA
Work Phone: (+9661) 464-7272 Ext. 32526
Fax: (+9661) 442-4542
Email: shoukri@kfshrc.edu.sa
Website: www.kfshrc.edu.sa

SAS-L – A VERY POWERFUL RESOURCE FOR SAS USERS WORLDWIDE

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

這一篇由 JoAnn Matthews 等人於 2003 年發表在 SUGI 28 的文章，並不是在講關於 SAS 的用法，而是在介紹一個很棒的 SAS 討論群組，叫做：SAS-L。

這有點像是一個新聞群組，也可以說是網頁論壇，由 University of Georgia 所提供的平台，讓眾多 SAS 的愛好者在這邊互相交換心得，並且提出自己遇到的問題。雖然 SAS-L 的介面並不太具有親和力，不過這邊的確吸引了很多高手前來。如同一般的論壇一樣，如果沒有註冊，則只能觀看留言。註冊後就可以開始貼文章了。 JoAnn 在文內有詳細介紹如何註冊，如何用新聞群組軟體來觀看，或搭配 Google 搜尋。有興趣的人可以瞧一下她的講解。

Splitting a Large SAS® Data Set

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

當處理龐大檔案，若沒有很強的硬體設備，通常我們會建議將檔案稍做一些切割以便分別處理。但 SAS 中沒有任何指令是可以任意切割資料，即便最簡單的「均分」都沒有指令提供。Selvaratnam Sridharma 於 2003 年的 SUGI 28 發表了一篇關於切割資料的文章，並寫了兩個小 macro 分享給大家。


The %split1 Macro

第一個 macro 程式可供使用者將資料切割成數塊，每一塊包含 n 個觀測值，SAS 會自動命名新的小資料集，而切剩的資料會放進最後一個小資料集中（此資料的觀測值總數就會小於 n）。

%macro split1(num);
data _null_;
if 0 then set orig nobs=count;
call symput('numobs',put(count,8.));
run;
%let m=%sysevalf(&numobs/&num,ceil);
data %do J=1 %to &m ; orig_&J %end; ;
set orig;
%do I=1 %to &m;
if %eval(&num*(&i-1)) <_n_ style="color: rgb(51, 51, 255);">%mend split1;

範例：

data orig;
do i =1 to 84;
output;
end;
run;
%split1(18);

資料集 orig 會被切割成五塊（orig_1~orig_5）。前四塊各含 18 個觀測值，最後一個則只有 15 個觀測值。

The %split2 Macro

第二個 macro 會將檔案切成 n 塊，每塊的觀測值會相等。同樣地，切剩的資料仍舊會被丟到最後一個資料集中。

%macro split2(num);
data _null_;
if 0 then set orig nobs=count;
call symput('numobs',put(count,8.));
run;
%let n=%sysevalf(&numobs/&num,ceil);
data %do J=1 %to &num ; orig_&J %end; ;
set orig;
%do I=1 %to #
if %eval(&n*(&i-1)) <_n_ style="color: rgb(51, 51, 255);">%mend split2;

範例：

data orig;
do i =1 to 85;
output;
end;
run;
%split2(6);

資料集 orig 會被切成六塊（orig_1~orig_6）。前五塊各含 15 個觀測值，最後一個資料集則只有 10 個觀測值。

Contact Information
Selvaratnam Sridharma
U.S. Census Bureau
4700 Silver Hill Road, Stop 8400
Washington, D.C. 20233
301-763-5398

Implementing a Multiple Comparison Test for Proportions in a 2xc Crosstabulation in SAS

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

一般在分析一個 binary 變數在兩個不同群體下的比例是否一致，我們會先繪製一個 2 x 2 表格，然後用卡方分析來檢定兩個比例。如果有 c 個群體，則卡方分析也可以檢定 2 x c 表格。然而，此時的虛無假設變成：

H0: P1=P2=...=Pc

如果拒絕了這個虛無假設，我們只能得到這個結論：至少有兩個比例互不相等。可是有哪些是不相等的，卡方分析並無法提供進一步的報表。於是， Alan C. Elliott 和 Joan S. Reisch 兩人在 2006 年的 SUGI 31 發表了一篇關於在卡方分析底下的多重比較方法，並提供 macro 程式給人使用。

在使用這個 macro 前，需要做一點小小準備。首先，要用 PROC FREQ 製作 2 x c 表格，並且將每個 cell 裡面的 frequency number 另存出來。如圖所示：



如果群組數目不是很多，可以另外新建一個資料，把這些數值（紅色圈選處）寫進去。當這個資料準備好後，立刻接上 Alan C. Elliott, Joan S. Reisch 所提供的 macro，就等於是大功告成了。以他們的程式為例：

data;
do a = 1 to 2;
do b = 1 to 4;
input wt @@;
output;
end;
datalines;
32 43 16 9
55 65 64 16
run;
%compprop(n=4);

結果如下所示（不甚清楚，可詳見原文）：



Macro 程式如下所示，可直接複製使用：

**************************************************************
* SAS® macro to perform multiple comparisons on proportions *
**************************************************************;
%macro compprop(n=);
proc freq;
weight wt;
tables a*b /chisq out=freqout outpct;
run;
*Select variables from data set freqout;
data propmc;set freqout;
if a=1;
* Apply an arcsine transformation to the proportions;
pi180=180/3.141;
n=(count*100)/pct_col;
p=0.5*(arsin(sqrt(count/(n+1)))*pi180 + arsin(sqrt((count+1)/(n+1)))*pi180);
proc sort data=propmc;by p;
proc transpose data=propmc out=propmc1;var n p b;
*Separate out pieces of the data and recombine into the needed data set;
data tmp;set propmc1;
if _N_=1;
%do i=1 %to &ncols;
n&i=col&i;
%end;
data tmp2;set propmc1;
if _N_=2;
data tmp3;set propmc1;
if _N_=3;
%do i=1 %to &ncols;
lab&i=col&i;
%end;
drop col1-col&ncols;
data propmc3;merge tmp tmp2 tmp3;
run;
* The array q05() contains the table lookup values for the q statistic with infinity
degrees of freedom for k groups (1 to 10);
data propmc2;set propmc3;
array p(*) col1-col&ncols;
array na(*) n1-n&ncols;
array lab(*) lab1-lab&ncols;
array q05(10);
q05(1)=.;
q05(2)=2.772;
q05(3)=3.314;
q05(4)=3.633;
q05(5)=3.585;
q05(6)=4.030;
Statistics SUGI 31 and Data Analysis
5
q05(7)=4.170;
q05(8)=4.286;
q05(9)=4.387;
q05(10)=4.474;
icompare=&n;
qcrit=q05(icompare);
* print out the multiple comparison results table;
msg=' ';
file print;
put ' Tukey Style Multiple Comparisons of Proportions';
put ' ';
put ' Compare Diff SE q q(.05) Conclude';
put ' ------------------------------------------------------------';
iskiprest=0;
* as the table is constructed, determine the correct Conclude message;
do i=icompare to 1 by -1;
do j=1 to i;
if i ne j then do;
se=round(sqrt((410.35/(na(i)+.5))+(410.35/(na(j)+.5))),.01);
diff=round(p(i)-p(j),.01);
q=round((p(i)-p(j))/se,.01);
iskip=0;
* Do not test if last one in series was Accept;
if msg='Accept' and j ne 1 then do
msg='Do not test';iskip=1;
put @5 lab(i) 'vs' @10 lab(j) @15 msg;
msg='Accept';
end;
* Do not test last one if previous one was Accept;
if iskip=0 and iskiprest=1 then do
msg='Do not test';iskip=1;
put @5 lab(i) 'vs' @10 lab(j) @15 msg;
end;
if iskip=0 then do;
if q gt qcrit then msg='Reject';else msg='Accept';
put @5 lab(i) 'vs' @10 lab(j) @ 15 diff @25 se @35 q @45 qcrit @55 msg ;
end;
if j=1 and (msg='Accept' or msg='Do no test')then iskiprest=1;
retain msg;
end;
end;
end;
put ' ';
put ' Reference: Biostatistical Analysis, Fourth Edition, Jerrold Zar, 1999, p
564.';
run;
%mend compprop;

附帶一提的是，原文內有提到 c 的數目是從 2 到 10，但沒說超過 c 會發生什麼事情。所以在使用時若群組數目超過 10 的話，請特別留意任何錯誤訊息。

CONTACT INFORMATION
Alan C. Elliott
Center for Biostatistics and Clinical Science
UT Southwestern Medical Center Dallas
5323 Harry Hines Blvd
Dallas TX 75390-8822
Phone (214) 648-2712
Fax (214) 648-7673
alan.elliott@utsouthwestern.edu

You Can’t Stop Statistics: SAS/STAT Software Keeps Rolling Along

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

兩位 SAS 內部人員 Maura Stokes 和 Robert Rodriguez在 2006 年的 SUGI 31 發表了一篇相當重要的文章，內容主要在說明最新版本的 SAS 9.2 新功能。其中比較重要的是增加了一些舊版沒有的「內建」的程序，如 GLIMMIX、QUANTREG、GLMSELECT（其實這些可以從 SAS 官網下載外掛程式讓舊版的 SAS 使用）。此外，許多高解析的圖表也一併內建到 ODS 裡面，讓使用者可以輕鬆的繪製一些以往要花很多程式碼才能完成的圖形。由於這篇文章涵蓋的範圍太廣，有新的程式碼，有新的輸出報表和新的高解析圖表，因此就僅拿一些比較「炫」的部分來分享一下。

GENERALIZED LINEAR MIXED MODELS

PROC GLIMMIX 程序在 2005 年就已經發表了，不過僅供外掛，在最新的 SAS 中已經加入了這個相當具有威力的程序。程式的寫法和 PROC MIXED 很像，如下所示：

proc glimmix data=Neuralgia; class Treatment Sex;
model Pain= Treatment Age Treatment*Age /solution oddsratio;
ods select ParameterEstimates OddsRatios;
run;

連輸出結果都很像：




當然報表解讀又是另一回事情了。想更進一步瞭解 PROC GLIMMIX 的使用方法可以參考

Schabenberger, O. (2005). “Introducing the GLIMMIX procedure for Generalized Linear Models,” Proceedings of the Thirtieth Annual SAS Users Group International Conference. Cary, NC: SAS Institute Inc.

MODEL SELECTION

一般在配適線性模式時，如果程序裡面已經有內建一些選擇變數的 option，如 PROC REG 的 selection statement 下的 method 指令，可以用 forward、backward 和 stepwise 來挑。如果程序裡面沒有這類 option，如 PROC GENMOD，就要用手動的方式來嘗試。PROC GLMSELECT 程序則是針對任何標準的廣義線性模式來做最佳變數選擇。特別要注意的是，他只能處理 univariate response 的情況，而無法處理 multivariate responses。此外，其相關的圖形也已經完全和 ODS 整合起來，所以讓我們來看看這個範例：

ods graphics on;
proc glmselect data=baseball plots=all;
class league division; model logSalary = nAtBat nHits nHome nRuns nRBI nBB yrMajor crAtBat crHits crHome crRuns crRbi crBB league division nOuts nAssts nError / details=all selection=stepwise(select=sl) stats=all;
run;
ods graphics off;

此處要特別將 ODS 可以產生的效果提出來一下。此程式啟動了 ODS 功能（紅色標示處），並且在 PROC GLMSELECT 後加上 plots=all 的選項，這樣 SAS 就會自動產生所有在此程序可以生出來的圖形。當然，如果妳只想指定某一張圖形的畫，就要指定特定名稱。這些名稱都有紀錄在 SAS 官網上面，供使用者免費下載。如果指定 plots=criterionpanel ，則 SAS 只會產生 Criteron Panel 這張圖。該圖秀出所有可供判斷最佳模式的圖表。為免仍有人看不懂這張表，SAS 乾脆在圖上標上☆符號，告訴使用者這個就是最佳模式！如下圖所示：



ANNOUNCING NEW SOFTWARE: SAS STAT STUDIO

新版 SAS 還可以呼叫一個新的軟體名為 SAS STAT STUDIO。這整合了包含 SAS 資料、程式、輸出報表和圖形，且可以利用表單點選的方式來完成分析（類似SPSS吧）。但文內並沒有特別敘述如何操作，只秀了一張很 fancy 的圖，不過我們可以拭目以待。



EFFECT PLOT IN PROC LOGISTIC

其實 logistic regression model 並沒有什麼太重要的圖需要特別繪製，因為重點都是在如何解釋那個 Odds Ratio。不過新的 ODS 仍舊聊勝於無地加了幾張圖。只要在 PROC LOGISTIC 後面指定 plots option 就可完成。裡面用了一個範例如下：

ods graphics on;
proc logistic data=one plots=(effect(clband yview=(.5,1)));
class Treatment Diagnosis / param=ref; model Cured/N= Diagnosis Treatment;
ods select effectplot;
run;
ods graphics off;

這個程式可以生出一張 Predicted Probabilities Plot：



BAYESIAN ANALYSIS FOR THE PIECEWISE EXPONENTIAL MODEL IN PROC PHREG

這回連 PROC PHREG 裡面都加入了新的指令可以用 Bayesian analysis 來配適 piecewise exponential model。我在 2005 年修 Bayesian analysis 時有學到這段，可是當時只有 WINBUGS 這個軟體可以比較有效率的估計參數和繪製圖表（但語法還是很難學，而且還要上工作站去跑程式）。現在 SAS 也加入了這個功能，方法相當簡單：

ods graphics on;
proc phreg data=Exposed;
model Days*Status(0)=Treatment Sex;
bayes piecewise=loghazard;
run;
ods graphics off;

令人訝異地，只需要補上紅色那段程式碼，就可以完全搞定。要產生相關圖表甚至不需要在 PROC PHREG 後面加上 plots 選項，只要直接呼叫出 ODS 即可。

在這裡僅列出圖表。



想要進一步地瞭解貝氏理論在倖存分析上的應用，可以參考下面這本教科書：

Ibrahim, J. G., Chen, M., and Sinha, D. (2001). Bayesian Survival Analysis, New York: Springer.

附帶一提的是，本書的作者 Joe Ibrahim 是北卡大生統系的貝氏理論權威，當年教我貝氏理論和高等數理統計學的就是這位大師。Joe 人是不錯，可是他的考試每次沒耗個六七個小時是不可能寫的完的。。。。。

BAYESIAN ANALYSIS FOR THE POISSON REGRESSION MODEL IN PROC GENMOD

同樣地，在 PROC GENMOD 下，貝氏理論也可以輕鬆套用了！^^

ods graphics on;
proc genmod data=liver;
model y = x1-x6 / dist=poisson;
bayes;
run;
ods graphics off;

這回連 bayes statement 後面啥東西都不用加了！Orz

THINGS GO BETTER WITH PROC TTEST

連最簡單的 T-test 都有新的花樣！除了使用 ODS 可以產生長條圖和常態曲線外，一個很炫的「交叉分析」也出現在報表當中。讓我們來先看看程式怎樣寫：

ods graphics on;
proc ttest data=asthma;
var PEF1 PEF2 / crossover= (Drug1 Drug2);
run;
ods graphics off;

紅色那段程式碼就是 SAS 9.2 為 PROC TTEST 程序新增的指令。這個指令新到連目前 SAS 線上手冊的 PROC TTEST 程序指令集都還沒有加入。

他會產生下面這些圖形：



REGRESSION DIAGNOSTICS FOR GEE MODELS IN THE GENMOD PROCEDURE

早期要使用 GEE 在 PROC GENMOD 時，並沒有任何關於模式診斷的指令可供使用。我和一位老師討論過為什麼大家在做 GENMOD 時不重視模式診斷，他自己也說不上來，總之結論就是這一部份的理論很晚才發展出來，因此 SAS 在一開始時沒有加入相關程式。這部分的理論大約是在 90 年代末期才陸續有一些論文發表出來。我們先來看一下程式：

ods graphics on;
proc genmod data = preqaq99 descending plots=(cooksd clustercooksd);
class pract_id ;
model bothered = female age dayacc severe toilet/d=bin itprint ;
repeated sub=pract_id/corr=exch modelse;
run;
ods graphics off;

紅色部分的程式碼可以呼叫大家最經常拿來鑑定 influential data 的 Cooks' D 值。如下所示：



關於這部分的理論，可以參考下面這篇論文：

Preisser JS, Qaqish BF (1999), “Robust Regression for Clustered Data with Application to Binary Responses”, Biometrics, 55, 574–579.

其中 John Preisser 是我的 committee member 之一，而 Qaqish 是他的博士論文指導教授。Qaqish 本人以前是 GEE 的發明者之一梁賡義博士的博士論文指導學生。Preisser 和 Qaqish 目前仍在北卡大生統系執教著。

POWER AND SAMPLE SIZE APPLICATION

SAS 其實已經有可以計算 power analysis 的程序，如 PROC POWER 和 PROC GLMPOWER。其中 PROC GLMPOWER 在 9.2 版前要外掛。文內並沒有特別說 9.2 版會加入這個程序。不過倒是有研發一個視窗介面的程式讓使用者用點選的方式就可以完成 power analysis。這對於一些對程式不熟悉的人來說應該是相當方便。

圖就不另外列了，文章裡面有。

ENHANCEMENTS TO SAS/STAT PROCEDURES

除此之外，還有其他一些林林總總的新增項目，僅列出我覺得比較重要的：

• PROC GENMOD 新增可以做 model selection criteria 的 AIC 和 QIC 值。
• PROC MIXED 新增 method=laplace。
• 在 survery analysis 程序中加入 Jackknife 法。
• PROC PHREG 新增 class statement。
• 新增 PROC QUANTREG 程序。
• PROC NPAR1WAY 新增 conover test 和中位數差的 Hodges-Lehmann 信賴區間。
• PROC GENMOD 和 PROC PHREG 可算出 cumulative residuals。