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