→ andrew43:你的問題2是不是想求「一種因子」的綜合檢定? 06/22 14:20
→ andrew43:Anova(..., test="...", type=...) 06/22 14:22
→ andrew43:至於問題1, 你自己容易解釋即可, 倒沒什麼要緊的事. 06/22 14:24
我手上有 114 個因子 (factors , features , SNPs)
想做 Stepwise Insert , multiple logistic regression
建 prediction model
因此要先 filter 掉一些完全無用的因子
必須初步做 single logistic regression , 檢定 "單項因子是否有價值"
====================================================================
而在觀察 summary(glm(簡單回歸)) 的結果時
發現有以下兩種現存的檢定 , 但不知到哪種檢定符合我的要求
1. 看 summary(glm(單項回歸式))$coefficients 中 , 是否認一項顯著
2. null devaince - model deviance 的 LRT
主要想問的是這個...
※ 編輯: gsuper 來自: 140.113.239.247 (06/22 22:35)
→ andrew43:我被你的 "因子" 搞糊塗了. 從你最先的例子,只有一個因子 06/23 00:01
→ andrew43:不知你是否可以重新舉一個簡化的實例, 並附上r原碼 06/23 00:02
→ andrew43:方便之後討論? 06/23 00:02
→ andrew43:在你原文中所談到的deviance相減, 是檢驗該迴歸式 06/23 00:36
→ andrew43:是否存在任一個以上的有顯著效果的變數, 似乎不是你要的. 06/23 00:36
Anova(model)$p.value
= 1-pchisq(mod$null.deviance-mod$deviance, 2)
= LR test
至少在簡單邏輯回歸上述應該是成立的
LR test 在 複邏輯回歸
1. 觀察整個 model 是否有效 (compare with null deviance)
2. 新增 factor , 對舊 model 是否有顯著的進步 (compare with old deviance)
我是這樣理解的
也不知道對不對
→ gsuper:我正在努力寫 等我一下 06/23 01:00
先釐清名詞
x變項 = factor(因子) , factor內含有 3 個 levels (AA,AG,GG)
===============================================================
===============================================================
主要目的有兩項
1. 先決定 x 是否有用
(是否具備基礎的預測能力)
2. x 是三元變項 , 但是否轉換成二元變項 , 預測能力會更佳
(e.g. AA vs AG+GG)
==============================================================
==============================================================
先設定資料
y = 1 = case , 共 20 人
x = AA , AG , GG
y <- c(1,1,1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,0,0)
x <- c("AA","AA","AA","AA","AA","AA","AA","AA","AA","AG",
"AG","AG","GG","GG","GG","GG","GG","GG","GG","GG")
初步觀察資料
AA 與 case 相近
AG 稍微偏向 control
GG 與 control 相近
==============================================================
==============================================================
組四種 models , 第一個3元 , 後3個2元
m123 <- glm(y~factor(x,level=c("AA","AG","GG")),family=binomial)
m1_23 <- glm(y~factor(x=="AA") ,family=binomial)
m2_13 <- glm(y~factor(x=="AG") ,family=binomial)
m3_12 <- glm(y~factor(x=="GG") ,family=binomial)
==============================================================
==============================================================
從 Deviance 來看
最小代表最好 (因子解釋力最強)
> summary(m123)$deviance
[1] 3.819085
> summary(m1_23)$deviance
[1] 6.701994
> summary(m2_13)$deviance
[1] 27.32723
> summary(m3_12)$deviance
[1] 10.81347
#結論 1 > 2 > 4 > 3
==============================================================
==============================================================
從 LR test P.value 來看 (至少有一個 level 顯著)
> library(car)
> Anova(m123)[[3]]
[1] 6.437302e-06
> Anova(m1_23)[[3]]
[1] 4.535914e-06
> Anova(m2_13)[[3]]
[1] 0.5277844
> Anova(m3_12)[[3]]
[1] 3.914466e-05
#結論 : 2 > 1 > 4 >>>>>> 3
#結論 : model 3 不可用 , 僅 model 1 2 4 有預測能力
=============================================================
=============================================================
觀察顯著的 effect size (斜率)
以max(顯著的斜率群)互相比較
> summary(m123) $coefficients[,c(1,4)]
Estimate Pr(>|z|)
(Intercept) 21.56607 0.9982341
factor(x, level = c("AA", "AG", "GG"))AG -22.25922 0.9981773
factor(x, level = c("AA", "AG", "GG"))GG -43.13214 0.9975772
> summary(m1_23)$coefficients[,c(1,4)]
Estimate Pr(>|z|)
(Intercept) -2.302585 0.02813286
factor(x == "AA")TRUE 22.868654 0.99691267
> summary(m2_13)$coefficients[,c(1,4)]
Estimate Pr(>|z|)
(Intercept) 0.1177830 0.8084737
factor(x == "AG")TRUE -0.8109302 0.5382557
> summary(m3_12)$coefficients[,c(1,4)]
Estimate Pr(>|z|)
(Intercept) 1.609438 0.03773005
factor(x == "GG")TRUE -21.175506 0.99555629
#結論 2 > 4 >>>>> (3與1都沒有顯著的斜率)
=============================================================
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:09)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:10)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:12)
整理上述結果
Deviance : 1 > 2 > 4 > 3
LR test : 2 > 1 > 4 >>>>>>>> (3不顯著)
斜率觀察 : 2 > 4 >>>>> (3與1都沒有顯著的斜率)
----------------------------------------------------
1. 先決定 x 是否有用
首先我知道了 model 3 不可用 (LR test)
為何從斜率來看 , 3元的模型 (model 1) , 三個因子都不顯著?
照理說 AA 和 GG 的預測能力都應該很好才對
2. 四個 model 裡面 , 如何判定哪種最佳?
照理說該看 Deviance , 所以 model 1 最佳
但從 LR test 來看 , model 2 稍佳
(LR test 和 Deviance 為何會不一致?)
而從斜率來看 , model 1 沒有斜率顯著 , 所以不可用 model 1
總覺得不能用 model 1 非常奇怪...
3. 三元模型的 Deviance 是否衡大於 二元模型?
p.s. 總覺得一定沒人看得懂我在說啥了 0rz...
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:34)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:38)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:38)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:42)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:48)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:52)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:54)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 01:58)
推 andrew43:LR test 你的理解應該沒錯. 06/23 02:01
→ andrew43:一直讓你補充內容, 辛苦了. 06/23 02:03
→ gsuper:大大累了先去睡吧 ~~ 我撐不住了 0rz.. 06/23 02:03
→ andrew43:認為deviance越小越好不總是個好方法. 06/23 02:12
→ andrew43:這概念可以類比成: 一般線性中 R^2 最大就是最好的模型? 06/23 02:13
→ andrew43:答案往往是否定的. 06/23 02:13
所以天秤的兩端
一邊是解釋力 ( 降低殘差 )
另一端是什麼? ( 提高 LR test 顯著性?
提高 AUC of ROC ? )
→ andrew43:在合併不同level時, 也請多加考慮是否有實際意義. 06/23 02:17
→ andrew43:就常見選擇模型的方法, 後者是正確的方式. 06/23 02:19
所謂"後者"是什麼?
是指一般是用 LR test 選有效因子, 而非 devaiance 最小者?
→ andrew43:不過你不是在 "挑因子" 而是在 "併水準", 這我也不敢多說 06/23 02:19
→ andrew43:如何為正解了. 06/23 02:20
→ andrew43:註: 你的第二解和第三解其實有相同意義, 只是檢定量不同. 06/23 02:22
第二解和第三解是在說下面這個 table 嗎?
Deviance : 1 > 2 > 4 > 3
LR test : 2 > 1 > 4 >>>>>>>> (3不顯著)
斜率觀察 : 2 > 4 >>>>> (3與1都沒有顯著的斜率)
我後來發現三元的 model , 全斜率不顯著的原因
是因為細格有0所導致 (odds ratio 分子分母都不可為0)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 16:09)
※ 編輯: gsuper 來自: 140.113.239.247 (06/23 16:36)
→ andrew43:我似乎有誤解你「比較斜率的方法」。先不討論之。 06/25 02:57
→ andrew43:光比解釋力是因為有最好的解釋力可能是過份多自變數造成. 06/25 03:00
→ andrew43:為了模型穩定等原因, 應避色使用這個方法. 06/25 03:01
→ andrew43:但因為你的例子不是挑因子而且修改因子, 所以我不敢確定. 06/25 03:02
→ andrew43:目前我能說你的方法二是常見的合理做法. 06/25 03:03
→ andrew43:抱歉. 方法二不太對. 應拿二個模型相比. 06/25 03:04
→ andrew43:也就是例如m1m2和m2比來決定m2是否足夠. 06/25 03:05
→ andrew43:但你目前做的是只和僅有常數項的模型相比. 06/25 03:06
→ andrew43:以上我敢說的都說完了. 還有待高手相助了. 06/25 03:08
補個程式筆記 (自用,與本篇內容無關)
str <- paste(paste("FACTOR[[",1:5,"]]",sep=""),collapse="*")
code <- paste("glm(PHENO~",str,",family=binomial)")
mod <- eval(parse(text=code))
※ 編輯: gsuper 來自: 140.113.239.247 (07/08 17:54)