這一篇其實一個多月前就該發的,但時間一拖自己也熊熊忘了。剛剛才想到所以稍微寫一下發表。

我過去這一年來在北卡大護理系接到最主要的統計諮詢案子是一個在護理系念了八年的超級大老級博士生,我們都叫她 Sue,老家在田納西。她本身就已經很有年紀了,可能比我媽還老,女兒好像也長大嫁人了。我剛進到護理系工作時,就曾被一個當年也在護理系念博班的學姐,同時也是她的「同梯」,帶去她家吃飯。結果現在那個學姐都回北醫當 faculty 好幾年了,但 Sue 卻還一直陷在自己論文的無底深淵裡面。

Sue 的研究是有關於人體疼痛後的反應,而資料是她自己一個一個慢慢找自願者接受測試收集得來的,樣本大概才六十幾人,但變數超多,所以好像收了好幾年才收完。收完後我光是幫她清資料就花了個把月,接著幫她架模式,中間也是折騰很久才弄完。結果在去年十月還是十一月的時候,她就偷偷地完成了 final defense 了。她口試這件事護理系聽說沒幾個人知道,連我跟她那麼熟的人都是之後才曉得。

但事情就這樣結束了嗎?當然不是,否則就不會有這篇文章了。在她號稱非常機密的 final defense 結束後,居然有口試委員不肯簽名,而有些已經簽名的口試委員則要求她要改善統計模式。OK,於是我就教她怎樣做 mixed model 的模式鑑定。結果做出來後,果然發現原先使用的模式違反了一些模式假設。而要解決的最簡單的方法就是去做 log transformation。

經過 log transformation 之後,模式果然好看了一點,雖然還不是相當完美。不過實證研究的統計分析本來就很難做出跟教科書裡面的範例一樣漂亮的結果,所以我並不意外。不過那段時間那些口試委員意見還是很多,所以 Sue 花了很多時間改模式。但改到後來,我發現她的模式已經很難拿去解釋她論文裡面的 research questions。我一直提醒她,過度配飾統計模式的結果,就已經不是在「用」統計,而是在「玩」統計了。而且,她現在要拿的是「護理博士」,不是「統計博士」,把統計分析的部分弄到那麼細那麼複雜,對她真正應該著重於「護理研究分析」上沒太大幫助。不過她還是沒有辦法抵抗口試委員的要求,畢竟我沒有辦法代替她的口試委員簽名,所以她只能聽口試委員的。不過後來我看有些口試委員要她補做的統計分析居然是錯誤地使用統計工具在不適用的資料上,當然程式還是會幫你跑出一堆報表出來,但在我的眼中都是無用的啊!

Whatever,她持續努力了一個多月,本來打算在去年十二月畢業的她,結果得拖到今年五月,這意味著她還要多付一個學期的學費,就只為了改論文。總之,當她終於把所有的模式作 log transformation 後,新的問題出現了。

在她博士論文裡面許多命題都是建立在實驗組和對照組的比較上,在一開始沒有經過變數變換的模式裡面,只需要輕鬆地在 SAS 的 PROC MIXED 程序裡面加上 contrast statement,並把 group 1 -1 這個比較係數帶入即可。這樣 SAS 會幫忙算出 E(Y_group1)-E(Y_group0) 的統計檢定。若用 estimate statement 則連組間比較的結果都算出來。但經過變數變換後,同一個 contrast statement 或 estimate statement 算出來的結果不再是 E(Y_group1)-E(Y_group2) 的假設檢定,而變成:

E(log(Y-group1))-E(log(Y_group0))

如果還記得這種簡單的運算的人應該知道,上式並不等於

E(log(Y_group1 - Y_group0))

而是

E(log(Y_group1/Y_group0))

白話之,就是 log(A)-log(B) 不等於 log(A-B) 而是等於 log(A/B)。因此,Sue 想要算出「組間差異值」(即 Y_group1 減 Y_group0)的命題是完全沒有辦法從這個經過 log transformation 後的模式解出。

針對這個問題,我和小老闆 Todd 討論之後得出來一個稍微可行的替代方案來解決這類問題。

在使用任何變數變換過後的統計模式做組間差異檢定,仍舊用 group 1 -1 這個比較係數帶入 SAS 運算,之後得到 E(log(Y_group1))-log(Y_group0))=E(log(Y_group1/Y_group0))=M 的統計檢定值。如果 p-value 小於顯著水準時,則可以下一個結論:Y_group1 的期望值顯著地高於 Y_group0 期望值的 M 倍。之後再用兩個 estimate statement 算出 E(log(Y_group1)) 和 E(log(Y_group0)) 兩個點估計量。假設分別是 A 和 B,則 Y_group1 減去 Y_group0 便是 Exp(A)-Exp(B)。但這個值並沒有辦法直接求得其統計量和 p-value,所以無法下結論說 Y_group1 是顯著地大於 Y_group0 多少。但由於之前已經下了「Y_group1 的期望值顯著地高於 Y_group0 期望值的 M 倍」的結論,因此也可以引伸出說 Y 在兩組間是有差距的(但沒有統計量和 p-value 輔助證實)。在製作表格方面,也要拆成兩個表格各自表述。

這是個相當特殊的案例,在此分享我們後來所使用的折衷方法給大家參考。

創作者介紹

ToTo 奇妙の冒險

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


留言列表 (10)

禁止留言
  • 阿民
  • 真是普渡眾生
    功德無量啊
    阿彌陀佛...
    :P
  • 你從中國廚藝學院逃出來了嗎?XD

    cchien 於 2009/01/14 00:35 回覆

  • joanna
  • 請問我要比較A,B藥兩者有無差異,用logistic reg y=A+B+C, 得到ORa, ORb, 這樣無法說ORa>ORb, 表示A比B藥更高風險, 那可以設計成: y=(A-B)/2+(A+B)/2+C
    這樣就得到A vs B的差異嗎? 還是要設dummy variable才可以? 請問您這篇有沒有參考什麼文獻呢?
  • 你的y是什麼?C又是什麼?

    cchien 於 2009/06/03 03:38 回覆

  • joanna
  • 急~

    Y是發生event(1,0),C是另一個藥(1,0).請問用standardized logit coefficients將A,B的beta除以各別報表裡的SE,數字絕對值大表示y越重要(是這樣做嬤), 但仍不知這兩者有無差異呢? 而A,B藥非互斥, 個人可能同時用,這樣設dummy 就不適合了嗎? 請順便回信箱好嗎?非常謝謝
  • 抱歉我的原則是絕對不回信箱,因為我曾經受到騷擾過。有問題就留言討論。

    我大概瞭解為什麼你要把A,B分開來,而不是用一個二項變數來表示。但先無論模式,有人同時服用A,B兩藥,這樣就沒有辦法單純地測出哪種藥效會比較好,因為一定會互相影響。

    你可以用一個3-level變數,然後去測A,B和A+B這三種藥程哪一種比較好,否則單比較A,B絕對會被老闆打槍。

    另外standardized parameter不是那種用法,別誤用。你如果一開始沒有把整個study design弄正確的話,後面的模式是對你一點幫助都沒有。

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

  • joanna
  • 沒回信箱沒關係, 謝謝你的指正. 請問有人也可能兩藥都未用, 那是否設成4-level 像這樣 x=0 (A0B0),x=1(A1B0),x=2(A0B1), x=3(A1B1)?
    然後以x=0為reference再跑OR, 但若有意義的話, 怎麼知道是x裡哪一個才會增加風險呢? 可以做事後比較嗎? 回到上一次問題, standardized estimate不是在比較不同Xi對y的影響嗎? 所以我想到這種做法看A or B對Y的影響, 但困擾的是drug A,B的單位是暴露有無, 這樣表示單位相同能用standardized estimate嗎? 謝謝!!
  • 如果你有(A,B)=(0,0)這種情況,SAS自然會判讀為4-level變數,這點你不用過於擔心,只要把coding做好即可。

    至於你想知道哪個Odds最高,得需要由SAS報表所給的估計值去算,SAS只會給你任兩種組合比較的Odds Ratio,但這個表格也是很有意義的,所以別忽略了。至於怎麼計算,建議你去圖書館借這一本書來看:

    Categorical Data Analysis Using The SAS System. 2nd Edition.
    Author: Maura Stokes, Charles Davis & Gary Koch.

    你沒時間看的話就直接翻到p.191看Table 8.4,那邊有個簡單範例公式。

    關於standardization的問題,在logistic regression裡面,這個動作的目的和在linear regression是一樣的,主要是消除不同變數的單位。但是這種目的在類別變數是沒什麼意義的。建議你可以看下面這個連結的第二段,裡面有詳盡解釋:

    http://goliath.ecnext.com/coms2/gi_0199-762729/Six-approaches-to-calculating-standardized.html

    以你的研究目的來講,標準化也許不是必要的。而這篇文章的第一段最後面也有說:there is as yet no widely accepted definition for a standardized coefficient in logistic regression.

    Hope it helpful.

    cchien 於 2009/06/04 03:47 回覆

  • 阿民
  • 你被騷擾啊...
    你什麼時候變正咩了...
    :P
    很羨慕啊
    我只收過一堆色情簡訊而已
  • 可以說你帶賽嗎?
    我今天一早起床看到我的另一個講棒球的blog遭到大量色情廣告留言攻擊。。。。。="=
    害我花了好久的時間才慢慢清完。

    cchien 於 2009/06/05 04:41 回覆

  • JOANNA
  • 快要碩士論文口試了!

    正要找書中. 請問照我這樣設計, 看x=1(A1B0)及x=2(A0B1)的OR誰大, 就可以回答這兩者藥有無差異嗎? 還是他們都以x=0為reference, 不能直接這樣比? 本來我假說是寫 AB藥對event risk 有無差異, 後來才想到寫錯, 原本式子是各以A=0, B=0 為reference, 所以只能設假說: 有用A (B)比沒用A(B)的風險無差異.就不能做comparative effectiveness study, 謝謝!
  • 請自行參考這個範例:
    http://www.tau.ac.il/cc/pages/docs/sas8/stat/chap39/sect46.htm

    用contrast statement。

    cchien 於 2009/06/04 11:59 回覆

  • JOANNA
  • 對不起我對統計很笨, 請問按我想x=0,1,2,3這樣代表A,B使用情況,照您給的網站設contrast"pairwise後要怎麼寫x呢? 不好意思現在改成藥物名稱為LAIA, ILI,
    是否把它轉換為
    A B C P
    LAIA 1 0 1 0
    ILI 0 1 1 0
    而ABCP就是我另外設一個x, 所以proc phreg; class x; model time*censor(0)=x; 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 / estimate=exp; run;
    因為我是用conditional logistic regression, not proc logistic;, 所以class, expb 沒有這個指令該怎麼辦呢?
  • PROC LOGISTIC程序仍舊可以處理conditional logistic regression, 請參照我上次提到那本SAS書的p.278。由於你的研究目的是要做事後比較,就一定要用PROC LOGISTIC程序,因為PROC PHREG沒有contrast事後比較的語法。你不能隨意把PROC LOGISTIC的語法拿到PROC PHREG裡面來用,SAS沒那麼威到可以讓你這樣套用。雖然PROC PHREG程序裡面有個TEST語法,但是那只是用來做迴歸係數的線性檢定。至於contrast statement後面的hypothesis matrix要怎麼架,真的很難用三言兩語來形容。最基本的架構是:

    contrast "comparison statement" var 1 0 0 -1;

    這是假設一個4-level變數裡面,比較1st level和4th level的差異。但其實這裡面的學問真的很大,我自己是花了很長一段時間才慢慢學會,學校的老師則都當作我們都會了,所以也沒有很仔細教。不過身為一個研究生,自我快速學習的能力是很重要的,所以我能教妳的大概就是這樣,其餘的都是case by case,你得自己看書或上網找相關範例來學。我可以再推薦一本書給妳:

    Regression and ANOVA-an integrated approach using SAS software.
    Author: Keith E. Muller & Bethel A. Fetterman

    裡面有很完整詳盡的事後比較分析教學。我知道妳可能沒時間仔細看,但有時候妳沒什麼選擇。

    cchien 於 2009/06/05 04:40 回覆

  • joanna
  • 我看了那本書的10.5.1(p297),跟我的案件很像,它是用proc phreg, 是因為同一人不同periods要用proc phreg not proc logistic? 還是我得把periods再變成不同id, 變更資料結構,就不必在data裡寫time=2-event, 這樣就可以跑logistic,也解決我要的東西呢? 還有曾看網路有人提到proc tphreg能解決proc phreg的categorical var,請問這樣我可以不必變更資料結構下, 把它當logistic在options那邊寫proc logistic的指令, 也能事後比較呢? 非常謝謝!
  • 稍微找了一下資料,但沒有得到太明確的訊息。
    有幫妳問一下我的小老闆,因為我也頗好奇這個問題,但不知道他啥時會回信。
    但一句事後話就是,如果我是妳的話,在一開始就知道要做事後比較的時候,我就會去朝generalized linear model或generalized linear mixed model來做。
    此外,由於不清楚你的資料結構,沒有辦法單純從你的敘述中得知是不是要改變資料結構。不過由於這是做研究生的基本功,而且配適模式時,是要讓資料來配合模式,而不是讓模式來配合資料。所以當選定一個模式後,就要試圖去讓資料去符合模式所能讀進的格式。
    我從妳的IP判斷你是台大公衛學院的學生,所以妳可以去問看看流病所生統組的學生或老師有沒有人能提供進一步的幫助,雖然我曾經聽說流病所的老師都很那個....but you need to try....

    cchien 於 2009/06/06 03:36 回覆

  • joanna
  • 對不起我忘了說,我是1:3 matched, 是否這樣一定得用proc phreg了? 謝謝~
  • joanna
  • 我剛好有生統組的朋友,但未能解答.我一開始不曉得要比較兩物差異,又不能設dummy時, 以為單純用logistic regression可以回答(我們以前課沒教什麼GLM,所以壓根沒想到要弄那個).最早先問(A-B)/2 + (A+B)/2是一流病所畢業的助理教授提供意見,他說一般是在線性model較常這樣做,叫我試試看(可是沒有reference有同此做法,其他老師覺得很怪的方法),昨天我參考你介紹的網站設contrast, 然後用proc tphreg跑proc phreg的寫法加上logistic 寫法,可以出來了! 應該是可以說A比B有無差異,可以心安些. 真的非常謝謝你的指導.
  • ok, 祝口試順利。

    cchien 於 2009/06/07 06:40 回覆