Ann Jessup,一個護理系即將畢業的博士班學生,其 committee member 之一正是我老闆 Dr. Bangdiwala。由於這一陣子我老闆跑去智利不知道去幹啥,因此 Ann 的博士論文中的統計部分就落到我的頭上了。



她的論文題目裡面用到所謂的 GEE model,這對我來說並不是太困難的事情。但之後卻發生了連我 RA 老闆 Mark 都不知道該怎樣解決的麻煩。第一個麻煩是,Ann 要從下面這三個 GEE model 中挑選最好的一個





Y=X1+X2+X3+W

Y=X1+X2+X3+V

Y=X1+X2+X3+Z



這三個模式除了最後一個獨立變數不一樣之外,其餘都一樣。簡單的作法就是分別做三個模式,如果 W,V,Z 其中一個顯著而另兩個不顯著,自然就可以認定含有那一個顯著變數的模式是比較好的。



但是,如果有兩個以上顯著呢?



當然,我沒有一開始就直接去跑這三個模式,所以不知道結果會是如何,但總是要假設一下各種可能性。如果是在 mixed model,這種 non-nested 的情況就直接用 AIC 或 BIC 這兩個數據來判定,越小的越好。那 GEE model 有沒有類似的準則可以判定呢?答案是有的!GEE model 也可以用 AIC,但不幸的是,目前 SAS 9.1.3 版的 PROC GENMOD 沒有內建這個功能。(聽說9.2版會內建了,但不知啥時才會 release 出來)。於是我就上網去看看有沒有其他人發明相關的方法,果真發現一個叫做 QIC 的數據,而且 SAS 公司已經把他的 macro 程式放在網上任意下載。所以我就拿 SAS 提供的程式搞定了這一部份。



第二個麻煩是,當選好了 GEE model,要如何做模式檢測來挑出可能的離群值或影響值呢?其他的統計模式都有些方法可以來做模式檢測,就唯獨 GEE model 沒看過什麼模式檢測。我翻了以前上課的講義,老師也沒講這一段。跑去問 Mark,他居然也不知道。此時,又得去拜一下 Google 大神了。沒想到不只有人把這段理論給弄了出來,連程式也寫好了。不過更幸運的是,發表這些 paper 的學者中,有一個共同名字,叫做 John Pressier。而這位仁兄,正巧就是我的 committee member 之一。



於是乎當然就很高興的把 paper 稍微讀一下,接著就馬上下載他寫好的 macro 來用。中間遇到一個問題,那就是他的程式裡面並沒有辦法設定某些變數是離散的。後來直接寫信去問他,他叫我把那些變數改成 dummy variable 就搞定了。



所以,就在這麼許多巧合和運氣之下,Ann 的兩大問題就這樣被我解決了。同時間自己也學到了很多,其實我應該要更感謝 Ann 才對,否則我永遠也不會知道 QIC 和 GEE model diagnostic 的存在。來護理系工作真是惠我良多啊。












創作者介紹

ToTo 奇妙の冒險

cchien 發表在 痞客邦 PIXNET 留言(11) 人氣()


留言列表 (11)

禁止留言
  • 路人
  • COMMENT:
    很喜歡你寫的blog,

    問個小問題,

    一般linear regression, logistic regression, cox regression都有model selection

    (stepwise)的指令,

    但若在GEE中要做model selection,

    你知道要怎麼做嗎?

    謝謝
  • cchien
  • COMMENT:
    SAS沒有在PROC GENMOD中設計model selection的指令。最完整的作法就

    是用Likelihood Ratio Test來檢定兩個nested model是否有顯著差異。



    Ex:

    H0:(Reduced model)Y=b1*X1+b2*X2

    H1:(Full model)Y=b1*X1+b2*X2+b3*X3



    上述就等於:



    H0:b3=0

    H1:b3=!0



    當然一開始就要從saturated model當作full model,然後把其中一個變

    數拿掉,當作reduced model。若經過LRT顯示H0 is not rejected,就表

    示那個變數不顯著。然後下一步驟把沒有這個變數的模式變成full

    model,再拿掉一個變數當成reduced model,反覆操作上面的流程即可。



    過程很麻煩,但這是最精確的作法。
  • 路人
  • COMMENT:
    謝謝你!



    我也是這麼做,

    希望SAS在9.1.3之後的版本能趕快增加這方面的功能,

    就不用人工的方式了~~



    請問你有做過count data (Poisson regression)在time series方面的分析嗎?
  • cchien
  • COMMENT:
    你可以看我相簿內放的一篇 paper,可以給你一點幫助。
    -----
  • 請問....
  • 您好
    不好意思打擾 在搜尋GEE
    恰好看到您的網誌
    像是看到救星
    想請問GEE只能用SAS跑嗎? SPSS可以嗎?例如說第14.0或16.0版有這個功能嗎?
    謝謝
  • SPSS V16.0(英文版)的 GEE 操作流程:

    工具列->Analyze->Generalized Linear Models->Generalized Estimating Equations

    但我沒用過SPSS跑GEE,因為我都用SAS跑
    不過我還是順便進去看了一下SPSS的GEE功能
    發現沒有我在文章中提到的model selection criterion
    也就是QIC
    所以你要用QIC做model selection的話
    還是得用我的committee member, Dr. John Pressier, 所寫的QIC macro來跑

    cchien 於 2009/01/08 03:29 回覆

  • 請益
  • 哈囉您好:
    因為搜尋GEE而來到您的網誌,
    想問個有關GEE的問題:請問目前GEE是否能跑mutiple outcome的迴歸?
    我的意思是,outcome變項是三類以上的類別變項。
    謝謝^^
  • proc genmod的model statement後面的dist option用multi即可。

    cchien 於 2009/04/23 08:22 回覆

  • 正在寫護理論文的阿惠
  • 請問mixed model

    您好,看到您對於linear mixed model(LMM)也相當了解,在此希望您能救救我...
    請問進行LMM前變項需進行的檢定除常態分佈外還有其他的嗎?
    另外,於LMM的模式檢定,不看-2LL(-2 Restricted Log Likelihood)嗎,其與AIC( Akaike's Information Criterion)及BIC( Schwarz's Bayesian Criterion)這三個的差別?....感謝
  • 你要做的檢定不是在進行GLMM前,而是在進行GLMM後
    請去搜尋GLMM的model diagnositcs和influential analysis
    如果你用sas跑程式的話,就把下面這篇SUGI的技術文件讀一讀

    http://www2.sas.com/proceedings/sugi29/189-29.pdf

    至於你說的模式檢定,其實應該是「模式選擇」(Model Selection)
    -2 Restricted Log Likelihood只是拿來計算AIC等criteria的元素
    本身並沒有辦法直接決定模式的選擇
    所以只要看AIC/BIC即可
    AIC/BIC的差別google一下就有
    通常用任一個criterion應該都會得到一致的結果
    比方說A模式的AIC小於B模式的AIC
    則A模式的BIC大都也會小於B模式的BIC
    不過我通常會建議使用LRT(Likelihood Ratio Test)來進行GLMM底下的model selection
    用LRT是最正統的方法,也可以避免遇到AIC和BIC結果不同調的情況

    cchien 於 2009/05/15 00:56 回覆

  • joanna
  • 對不起又要打擾了

    又是我, 不好意思可以請教GEE的問題嗎? 由於這個研究case crossover study 除了用conditional logistic regression外, 我也想用GEE 看drug A,B的OR有沒有與conditional logit一致. 所以我也放了先前跟您提過的,將drug A,B可能使用的情形設成contrast"pairwise" x, 若model y=x, 可以跑出GEE result, 但再加上其他變項,GEE有出現parameter estimate, 但均無95%CI, Z,p,想不透到底錯在哪呢? 可以幫我看一下嗎? 非常謝謝!!
    對象均為有發生event的, 取case period,3 control periods prior to the event. 所以一人有四筆觀察資料. 請看我只放一個x 的program: proc genmod data=insu_fin2 descend;
    class x(ref='4') id;
    model outcome=x / dist=bin;
    repeated subject=id / corr=unstr corrw;
    contrast "pairwise" x 1 0 0 -1,
    x 1 0 -1 0,
    x 1 -1 0 0 ,
    x 0 1 -1 0,
    x 0 1 0 -1,
    x 0 0 1 -1;run;
    接著放多個covariates, 有的是(0,1), 有的是連續var
    proc genmod data=insu_fin2 descend;
    class id x(ref='4');
    model outcome=x drug_3 drug_4 drug_5 age pr sex c_index MPR period*x/ dist=bin expected;
    repeated subject=id / corr=unstr corrw;
    run;
    另外參考書的作法用macro 也是做不出來~_~
    %GEE1 (data=insu_fin2,
    YVAR=outcome,
    XVAR=intercpt x drug_3 drug_4 drug_5 age pr sex c_index MPR,
    LINK=3,
    VARI=3,
    CORR=6);
    run;
  • 要看log裡面的error message比較快。

    cchien 於 2009/06/11 03:15 回覆

  • JOANNA
  • 它寫error in computing variance function. error in parameter estimate. 這樣是只無法收斂嗎? 可是資料無missing data, 變項也覺得沒錯, 不知哪裡錯了呢. 謝謝~~
  • 上面看了一些論壇上面的文章
    發現有不少人也出現同樣的問題
    但卻沒有看到一個比較正確的回覆
    個人是推測你用unstructure covariance structure這邊可能有問題
    因為UN的未知參數很多,所以有時候沒有辦法完全估計出來
    你改成CS看看
    不過嚴謹一點的作法得去做model selection
    也就是得把每一個covariance structure都套進去看看
    然後挑一個最好的
    有時候如果repeated index是time的話
    會自動假設covariance structure是AR(1)
    不過當然會有例外
    就看你做研究的嚴謹度如何了

    cchien 於 2009/06/11 12:02 回覆

  • joanna
  • 我也用過exchangeable, 還是不行. 你說要每個covar structure 套的意思是Corr 1-6 都分別看看嗎? 有可能有autoregressive, 所以我試看看corr 的所有可能看看. (有一點看不太懂你的意思, 硬做看好了:p)
  • 一些基本的cov structure如CS, AR(1), TOEP都可以試,不過用AR(1)的時候必須是要在時間變數的間隔是完全一樣才行。這些都試不出來的話再用其他進階的。

    cchien 於 2009/06/11 13:38 回覆

  • JJ
  • 找了GEE MODEL 相關的資料~還是不知怎麼跑。重複量數到底要怎麼放在SPSS的GEE 變項??
    可以請妳教我GEE在SPSS的流程嗎? 感激不盡
  • 老實講我不會用SPSS跑GEE,但我知道怎樣用google。

    Enjoy it!

    http://www.ats.ucla.edu/stat/spss/library/gee.htm

    cchien 於 2009/06/19 15:27 回覆