作者gsuper (统计的巴比伦塔)
看板Statistics
标题[程式] R 的 factor 设定
时间Wed Mar 31 15:47:29 2010
我看了一篇 paper 後
想要重覆原作者的数据分析
但他只有给 unpair data 的 function
所以我就写信问作者要怎麽改成 pair
来来回回通了10几封信後
总算是把 paired data function 写好了
(最後还是他来写 , 我执行後给 report)
问题来了
他没办法重覆他的分析
paper上算出来是 400 , 我却算 2000 多 (多300就算差满多了)
以下是他的理由
不过我看不太懂
想请深入了解 R 的高手解释一下上色的那行是什麽意思?
补充 :
这是在跑 two-way ANOVA 之前建 linear model , 用 anova(lm())来做 ,
unpair function 的 linear model 有2个factor 和一个交互做用
而 pair function , 再多加一个 pair factor )
----------------------------------------------------------------
I found out that using R for anova with more than 2 factors
generate p-values depending on how factors enter into model.
See examples below.
That's an example for my other dataset.
But it applies to the paired sample data analysis.
Basically, we need to add a block effect in the two-way anova model
to account for that effect.
So, the results of the second case study in my paper might not be accurate.
### block effect 就是新的 pair factor ###
Although there are some minor issues,
but the interaction effect reported by R is still correct,
this leads to the corrected pooling of probesets whenever applicable.
This again proves that power of consolidation.
I suggest you not to use the per gene model for the paired samples in R.
If you could implement it in SAS or other softwares
which will give you the right TYPE III test, that should be fine.
Sorry for all the confusion.
--------------------example-------------------------------------
anova(lm(y~as.factor(fcid)+genotype+time, data=y))
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
as.factor(fcid) 4 4.0412 1.01029 5.5987 0.0034300 **
genotype 1 0.1275 0.12748 0.7065 0.4105584
time 4 6.2609 1.56522 8.6740 0.0003138 ***
Residuals 20 3.6090 0.18045
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> anova(lm(y~time+as.factor(fcid)+genotype, data=y))
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
time 4 9.1197 2.27993 12.6347 2.741e-05 ***
as.factor(fcid) 4 1.1806 0.29515 1.6356 0.2044
genotype 1 0.1292 0.12919 0.7159 0.4075
Residuals 20 3.6090 0.18045
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--
他的结论好像是在说
Factors 数量大於 3 的时候
R 会算的不准
而原因是在於 lm() 里面 factors 的置放顺序
会导致 unpredictable 的影响
请问我有误解吗?
---------------------------------------------------
列一下我的资料
Normal_1 与 Tumor_1 为同一病人组织 , 有 pair 关系
Normal_1 Normal_2 Normal_3 | Tumor_1 Tumor_2 Tumor_3
------------------------------------------------------------
block1| 10 20 30 40 50 60
block2| 70 80 90 100 110 120
block3| 130 140 150 160 170 180
转成以下格式跑 2-way ANOVA
anova(lm(tmp~
trt+
v+
trt*v+
block , data=data))
tmp trt v block
--------------------------
10 N 1 1
20 N 1 2
30 N 1 3
40 T 1 1
50 T 1 2
60 T 1 3
70 N 2 1
80 N 2 2
90 N 2 3
100 T 2 1
110 T 2 2
120 T 2 3
130 N 3 1
140 N 3 2
150 N 3 3
160 T 3 1
170 T 3 2
180 T 3 3
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 140.113.239.247
※ 编辑: gsuper 来自: 140.113.239.247 (03/31 15:47)
※ 编辑: gsuper 来自: 140.113.239.247 (03/31 15:49)
※ 编辑: gsuper 来自: 140.113.239.247 (03/31 15:49)
※ 编辑: gsuper 来自: 140.113.239.247 (03/31 15:51)
※ 编辑: gsuper 来自: 140.113.239.247 (03/31 15:55)
※ 编辑: gsuper 来自: 140.113.239.247 (03/31 15:57)
※ 编辑: gsuper 来自: 140.113.239.247 (03/31 16:12)
※ 编辑: gsuper 来自: 140.113.239.247 (03/31 16:13)
※ 编辑: gsuper 来自: 140.113.239.247 (03/31 16:16)
1F：→ bmka:我通常尽量避免用as.factor,因为不是很确定dummy variable 03/31 21:17
2F：→ bmka:会被怎麽设定, 自己设比较安全 03/31 21:18
3F：→ gsuper:那我把as.factor()都改成factor()会有效吗?还是跟本一样? 03/31 21:38
4F：→ clickhere:他的结论要你用type III的test.... 04/01 14:39
5F：→ clickhere:try Anova() in library("car") with type="III" 04/01 14:40
6F：→ clickhere:block & v 设反了? 可以用 tmp ~ trt*block + v 04/01 14:42
7F：→ clickhere:block,v,tr t都是 factor. 04/01 14:43
我同时把 1-way ANOVA 和 2-way 都改了成以下
(这是多重检定 , 要组合两种 Anova)
(factor name 先不改避免乱掉)
-----------------------------------------------
library(car)
Anova(lm(tmp~trt+block),type="III")
Anova(lm(tmp~trt+v+trt*v+block),type="III")
-----------------------------------------------
8F：推 lin15:那R的两个anova结果不同的原因是?? 我用lm跑估计系数是一样 04/01 15:02
9F：→ lin15:应该不是factor不同的问题@@? 04/01 15:02
我前面好像没讲清楚
这个 data set 是多重检定
在 20000 组 ANOVA table 中
paper 上找出 400 多组显着
我却找到 2000 多组
FDR 的调整和细部资料的检查都做过了
( BH method = 0.01 )
※ 编辑: gsuper 来自: 140.113.239.247 (04/01 15:36)
※ 编辑: gsuper 来自: 140.113.239.247 (04/01 15:37)
※ 编辑: gsuper 来自: 140.113.239.247 (04/01 15:39)
※ 编辑: gsuper 来自: 140.113.239.247 (04/01 15:39)
10F：→ lin15:我是自己随机产生资料去跑 lm(a~factor(b)+c) 跟lm(a~c+fact 04/01 15:57
11F：→ lin15:or(b) 结果是一样的 但用anova()还真的不一样... 04/01 15:58
※ 编辑: gsuper 来自: 140.113.239.247 (04/01 16:07)
12F：→ clickhere:因为SSE有type1,2,3种. type1是最直觉, type3不会因解 04/01 17:31
13F：→ clickhere:释变数的次序不同而影响. 04/01 17:32
14F：→ clickhere:但原po的问题应该不是这个, 直觉上. 04/01 17:33
15F：→ clickhere:不同原因在於计算先後及不对称样本数. 04/01 17:35
16F：→ clickhere:问题可能也不在facor, 所有的解释变数都是factor... 04/01 17:36
17F：→ clickhere:v 跟 blcok 设反了? 04/01 17:39
设反应该没差
因为只是 factor 命名而已
trt = trt factor
v = block factor
block = pair factor
18F：→ clickhere:可以请问是哪篇paper吗? 04/01 17:40
PAPER
http://www.biomedcentral.com/1471-2164/9/188
在最後面 method 的部份
有 4 条公式
分别是
1. one-way ANOVA (unpair)
2. two-way ANOVA (unpair)
3. one-way ANOVA (pair)
4. two-way ANOVA (pair)
http://research.stowers-institute.org/hul/affy/perGene.r
然後这是 unpair 的 function
不过别去读
会看很久
-------------------------------------------------------------
改成 type III 後
Anova(lm(tmp~trt+v+block+trt*v, type="III")) #2-way ANOVA
Anova(lm(tmp~trt+block,type="III")) #1-way ANOVA
变成
> x <- mt.rawp2adjp(waoTEST,proc="BH")
> sum(x[[1]]<0.01)
[1]
3091
※ 编辑: gsuper 来自: 140.113.239.247 (04/01 18:19)
※ 编辑: gsuper 来自: 140.113.239.247 (04/01 19:10)
19F：→ clickhere:line78: res$Pr[1:4] for 2-way ANOVA 04/02 09:09
20F：→ clickhere:sorry, line36. 04/02 09:10
21F：→ clickhere:line39以下全都得改. block和v不可换.否则给改comP[,3] 04/02 09:11
22F：→ clickhere:line24也得改. 04/02 09:11
23F：→ clickhere:BMC的IF看似不错, 但paper看看即可. 04/02 09:13
24F：→ clickhere:要确定拿到的p-value是trt和probset的交互项. 04/02 09:25
细节方面没问题
我先前已经把 function 拆开
step by step 的执行确认
作者原本用 v 就是当 block factor
只是新加入的 pair factor 用 "block" 命名
可以参考上面的资料矩阵 (这篇文章的6th页)
现在是 linear model 不晓得怎麽设才正确
※ 编辑: gsuper 来自: 140.113.239.247 (04/02 14:26)