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;此例的變數名稱就直接用 &LIST 把上面整合好的變數名稱一次叫進來。同樣地，如果你只有少量變數，那直接輸入即可。此程式會將算好的基本統計量另存在 basic_measures 的新資料集裡面。
ODS OUTPUT BASICMEASURES=work.basic_measures;
proc univariate data=work.&input.;
class day;
var &LIST.;
run;
ODS OUTPUT CLOSE;重新編輯資料集 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;
這一段是要將剛剛做出來的繪圖資料裡面的變數以及標籤製作成巨集變數以用於繪圖程式當中。如果你打算手動輸入，此段也可以省略。
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;用 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;
此例由於是在做重複觀測值的模型估計，所以用到 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;圖型如下所示：
這段是在畫另一種圖，但我覺得不是很重要，可以省略。
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;圖型如下：
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.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;
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;最後就是呼叫前面寫好的巨集程式，輸入不同的組間變異數和組內變異數的數值，定義總共要模擬的次數、樣本大小、時間點數量以及平均數。
%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
沒有留言:
張貼留言
要問問題的人請在文章下方的intensedebate欄位留言，請勿使用blogger預設的意見表單。今後用blogger意見表單留言的人我就不回應了。