公告

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

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


2010年3月26日 星期五

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

原文載點: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 的資料:


\begin{align}
D & = -2(\ln(\text{likelihood for null model}) - \ln(\text{likelihood for alternative model})) \\
& = -2\ln\left( \frac{\text{likelihood for null model}}{\text{likelihood for alternative model}} \right).
\end{align}

白話一點,就是把 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);  


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。我們來看原文內的一個範例。

範例模式:

doodle.png

這模式只使用一個因變數叫做 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
CODE { display: block; /* fixes a strange ie margin bug */ font-family: Courier New; font-size: 8pt; overflow:auto; background: #f0f0f0 url(http://klcintw.images.googlepages.com/Code_BG.gif) left top repeat-y; border: 1px solid #ccc; padding: 10px 10px 10px 21px; max-height:200px; height:200px; // for IE6 line-height: 1.2em; }