作者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)