[程式] R 的 factor 設定

看板Statistics作者 (統計的巴比倫塔)時間14年前 (2010/03/31 15:47), 編輯推噓1(1推 0噓 24→)
留言25則, 5人參與, 5年前最新討論串1/1
我看了一篇 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)
→
03/31 21:17, , 1F
我通常儘量避免用as.factor,因為不是很確定dummy variable
03/31 21:17, 1F
→
03/31 21:18, , 2F
會被怎麼設定, 自己設比較安全
03/31 21:18, 2F
→
03/31 21:38, , 3F
那我把as.factor()都改成factor()會有效嗎?還是跟本一樣?
03/31 21:38, 3F
→
04/01 14:39, , 4F
他的結論要你用type III的test....
04/01 14:39, 4F
→
04/01 14:40, , 5F
try Anova() in library("car") with type="III"
04/01 14:40, 5F
→
04/01 14:42, , 6F
block & v 設反了? 可以用 tmp ~ trt*block + v
04/01 14:42, 6F
→
04/01 14:43, , 7F
block,v,tr t都是 factor.
04/01 14:43, 7F
我同時把 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") -----------------------------------------------
推
04/01 15:02, , 8F
那R的兩個anova結果不同的原因是?? 我用lm跑估計係數是一樣
04/01 15:02, 8F
→
04/01 15:02, , 9F
應該不是factor不同的問題@@?
04/01 15:02, 9F
我前面好像沒講清楚 這個 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)
→
04/01 15:57, , 10F
我是自己隨機產生資料去跑 lm(a~factor(b)+c) 跟lm(a~c+fact
04/01 15:57, 10F
→
04/01 15:58, , 11F
or(b) 結果是一樣的 但用anova()還真的不一樣...
04/01 15:58, 11F
※ 編輯: gsuper 來自: 140.113.239.247 (04/01 16:07)
→
04/01 17:31, , 12F
因為SSE有type1,2,3種. type1是最直覺, type3不會因解
04/01 17:31, 12F
→
04/01 17:32, , 13F
釋變數的次序不同而影響.
04/01 17:32, 13F
→
04/01 17:33, , 14F
但原po的問題應該不是這個, 直覺上.
04/01 17:33, 14F
→
04/01 17:35, , 15F
不同原因在於計算先後及不對稱樣本數.
04/01 17:35, 15F
→
04/01 17:36, , 16F
問題可能也不在facor, 所有的解釋變數都是factor...
04/01 17:36, 16F
→
04/01 17:39, , 17F
v 跟 blcok 設反了?
04/01 17:39, 17F
設反應該沒差 因為只是 factor 命名而已 trt = trt factor v = block factor block = pair factor
→
04/01 17:40, , 18F
可以請問是哪篇paper嗎?
04/01 17:40, 18F
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)
→
04/02 09:09, , 19F
line78: res$Pr[1:4] for 2-way ANOVA
04/02 09:09, 19F
→
04/02 09:10, , 20F
sorry, line36.
04/02 09:10, 20F
→
04/02 09:11, , 21F
line39以下全都得改. block和v不可換.否則給改comP[,3]
04/02 09:11, 21F
→
04/02 09:11, , 22F
line24也得改.
04/02 09:11, 22F
→
04/02 09:13, , 23F
BMC的IF看似不錯, 但paper看看即可.
04/02 09:13, 23F
→
04/02 09:25, , 24F
要確定拿到的p-value是trt和probset的交互項.
04/02 09:25, 24F
細節方面沒問題 我先前已經把 function 拆開 step by step 的執行確認 作者原本用 v 就是當 block factor 只是新加入的 pair factor 用 "block" 命名 可以參考上面的資料矩陣 (這篇文章的6th頁) 現在是 linear model 不曉得怎麼設才正確 ※ 編輯: gsuper 來自: 140.113.239.247 (04/02 14:26)
→
01/02 15:05, 5年前 , 25F
我通常儘量避免用as. https://noxiv.com
01/02 15:05, 25F
‣ 更多 gsuper 的文章
文章代碼(AID): #1BiluL2v (Statistics)