作者gsuper (Logit(odds))
看板Statistics
標題[問題] Logistic regression 的變數設定(2)
時間Fri Jun 22 01:50:33 2012
假設我要在[R]中計算 單變量logistic regression
DATA =
Y X=基因型
------------------
1 AA
1 AA
1 AG
1 AA
1 -- (missing value)
0 GG
0 GG
0 GG
0 AG
0 --
--------------------
attach(DATA)
以下為可能的各種組合 截距 X1 X2 X3
==============================================================================
glm(Y~X) | -- , AA , AG , GG
glm(Y~factor(X=="AA")) | not AA , AA
glm(Y~factor(X=="AG")) | not AG , AG
glm(Y~factor(X=="GG")) | not GG , GG
glm(Y~factor(X=="AA")+factor(X=="AG")) | GG與-- , AA , AG
glm(Y~factor(X=="AA")+factor(X=="GG")) | AG與-- , AA , GG
glm(Y~factor(X=="AG")+factor(X=="GG")) | AA與-- , AG , GG
glm(Y~factor(X=="AA")+factor(X=="AG")+factor(X=="GG"))| -- , AA , AG , GG
==============================================================================
所以說
1. 一般是否習慣將 Low risk genotype 設定為截距項 (Baseline risk)?
2. 若我要定義單變量顯著 , 是要看下表的某項
斜率檢定(見下表)
又或是要做
The likelihood ratio test, G ?
> summary(glm(PHENO~as.factor(C[,5]=="AA")+as.factor(C[,5]=="AG")))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.36390 0.02594 14.026 <2e-16 ***
as.factor(C[, 5] == "AA")TRUE 0.01706 0.07916 0.215 0.829
as.factor(C[, 5] == "AG")TRUE 0.02072 0.03971 0.522 0.602
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.2349203)
Null deviance:
152.29 on 650 degrees of freedom
Residual deviance:
152.23 on 648 degrees of freedom
AIC: 909.47
===========================================================================
G = 152.29 - 152.23 = 0.06 = 卡方 (where
df = #predictors added to the model)
當我計算 glm(Y~X) , 共有四個子變項 , -- , AA , AG , GG
請問 df 等於 3 或 4 或 1?
=======================================================================
排版超亂.....
我盡力了~~
--
祭頌后靈的騎士道與白主教穩守黃金鄉
邊境兵躁動自高自大妄入堡壘
黑主教冷靜人格分裂
籠城王與雙子戰塔一籌莫展
掌握無限的魔女喚醒躺下的靈魂
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.113.239.247
※ 編輯: gsuper 來自: 140.113.239.247 (06/22 02:04)
※ 編輯: gsuper 來自: 140.113.239.247 (06/22 02:05)
※ 編輯: gsuper 來自: 140.113.239.247 (06/22 02:05)
※ 編輯: gsuper 來自: 140.113.239.247 (06/22 02:06)
※ 編輯: gsuper 來自: 140.113.239.247 (06/22 02:06)
1F：→ andrew43:你的問題2是不是想求「一種因子」的綜合檢定? 06/22 14:20
2F：→ andrew43:是的話, 可以使用ptt站內本板 #1DZB3Jf3 文中的 06/22 14:22
3F：→ andrew43:Anova(..., test="...", type=...) 06/22 14:22
4F：→ 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)
5F：→ andrew43:我被你的 "因子" 搞糊塗了. 從你最先的例子,只有一個因子 06/23 00:01
6F：→ andrew43:不知你是否可以重新舉一個簡化的實例, 並附上r原碼 06/23 00:02
7F：→ andrew43:方便之後討論? 06/23 00:02
8F：→ andrew43:在你原文中所談到的deviance相減, 是檢驗該迴歸式 06/23 00:36
9F：→ 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)
我是這樣理解的
也不知道對不對
10F：→ 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)
11F：推 andrew43:LR test 你的理解應該沒錯. 06/23 02:01
12F：→ andrew43:一直讓你補充內容, 辛苦了. 06/23 02:03
13F：→ gsuper:大大累了先去睡吧 ~~ 我撐不住了 0rz.. 06/23 02:03
14F：→ andrew43:認為deviance越小越好不總是個好方法. 06/23 02:12
15F：→ andrew43:這概念可以類比成: 一般線性中 R^2 最大就是最好的模型? 06/23 02:13
16F：→ andrew43:答案往往是否定的. 06/23 02:13
所以天秤的兩端
一邊是解釋力 ( 降低殘差 )
另一端是什麼? ( 提高 LR test 顯著性?
提高 AUC of ROC ? )
17F：→ andrew43:在合併不同level時, 也請多加考慮是否有實際意義. 06/23 02:17
18F：→ andrew43:就常見選擇模型的方法, 後者是正確的方式. 06/23 02:19
所謂"後者"是什麼?
是指一般是用 LR test 選有效因子, 而非 devaiance 最小者?
19F：→ andrew43:不過你不是在 "挑因子" 而是在 "併水準", 這我也不敢多說 06/23 02:19
20F：→ andrew43:如何為正解了. 06/23 02:20
21F：→ 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)
22F：→ andrew43:我似乎有誤解你「比較斜率的方法」。先不討論之。 06/25 02:57
23F：→ andrew43:光比解釋力是因為有最好的解釋力可能是過份多自變數造成. 06/25 03:00
24F：→ andrew43:為了模型穩定等原因, 應避色使用這個方法. 06/25 03:01
25F：→ andrew43:但因為你的例子不是挑因子而且修改因子, 所以我不敢確定. 06/25 03:02
26F：→ andrew43:目前我能說你的方法二是常見的合理做法. 06/25 03:03
27F：→ andrew43:抱歉. 方法二不太對. 應拿二個模型相比. 06/25 03:04
28F：→ andrew43:也就是例如m1m2和m2比來決定m2是否足夠. 06/25 03:05
29F：→ andrew43:但你目前做的是只和僅有常數項的模型相比. 06/25 03:06
30F：→ 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)