[問題] bootstrap regression test + parallel

看板R_Language作者skylikewater (choc.)時間11年前 (2013/09/26 01:32)推噓1(1推 0噓 4→)

留言5則, 3人參與討論串1/2 (看更多)

[問題類型]: 程式諮詢(我想用R 做某件事情，但是我不知道要怎麼用R 寫出來) [軟體熟悉度]: 使用者(已經有用R 做過不少作品) [問題敘述]: 這個月開始做多變項(Xi)線性迴歸的拔靴，但幾乎是看文件自己摸的，對結果沒什麼信心。由於涉及到package，經考慮還是來 R 版發，希望能徵詢到有經驗的前輩。我的資料是「篩選過」的 MRI 參數，例如 n 個區域，每個區域內含的體素(voxel)不等，全部經內插轉為100個voxel。也就是 Y = 100*n，個案有 200 人以上，Xi 都是常態分配的資料。對 Y=X1+X2+...+e ，我希望做拔靴 1. 增加power 2. 變相增加SNR 3. 老闆"建議" 來做個別 X 變項在這個迴歸模型中顯著與否的考驗，作為「該 X 變項是否在這個 voxel 有顯著解釋力」的說明參考本篇教學：　　短網址:http://ppt.cc/ilty 　　http://socserv.mcmaster.ca/jfox/Books/Companion/appendix/ 　　Appendix-Bootstrapping.pdf 1. 計算個別變項 t-test 的 p-value？　我對他方法的理解是　「把抽出來的參數分配，對原資料對迴歸求得的參數做 z-test」　但對它第十五頁 bootp 計算很疑惑：　看起來像是數一數「有幾個抽出來的參數絕對值，比實際的參數絕對值大」？　而且也看不懂分母為何設 2000 (twice boot time?) 　 John Fox 有一篇同標題2002年的文件，精神類似，　但他則建議 1. 不用絕對值 2. 因為是雙尾所以要乘 2 (?) 說到底這兩個統計量，　跟其他網路上可以找到的 bootstrapping regression test 教學，　都寫得不是很清楚。　像這兩篇 p 值的計算都不涉及標準誤。　但我實際用 pnorm((tboot-t)/se) 算出來的結果會過於顯著，　網路上可以找到的資源都看了，請各位指點怎麼實作、或推薦我補強什麼統計概念。 2. 多變項的 sample() 實作？　有可能直接的用 sample() 直接取出放回「多個變項」嗎？　經過實測結果很怪（p 都等於 1）。說實話我也不確定在迴歸的架構下　取出放回是不是「需要以同行/同個案/同樣本，來取出放回」　還是可以自由地每個變項分開抽 3. multicore？ boot 有多核選項，分為 parallel 跟 snow 　但以我的資料，在 Linux (CentOS 64bit) 上實測沒什麼加速的感覺。　如果可以用 sample 純粹自己手作，　我會想自己動手用 Shell script + R 分配不同次 boot 的結果　但我可能得用 I/O .rdata 的方法實現，不是很實際　想請問 R 目前有類似 matlab parallel toolbox 之類的東西嗎？ [程式範例]: Betas=function(formula,data,indices,maxit=20) { D=data[indices,] LMFit=rlm(formula=formula,data=D,maxit=maxit) X1D=deltaMethod(LMFit,"X1") X2D=deltaMethod(LMFit,"X2") X3D=deltaMethod(LMFit,"X3") X1ZDiff=X1D[1,1]/X1D[1,2] X2ZDiff=X2D[1,1]/X2D[1,2] X3ZDiff=X3D[1,1]/X3D[1,2] ZDiff=c(X1ZDiff,X2ZDiff,X3ZDiff) return(ZDiff) } BootTime=1000 set.seed(666) # for reproducibility for (依次取不同Voxel) { Voxel=該點Voxel (1*SubjNum vector) Data=data.frame(Voxel=Voxel,X1=X1,X2=X2,X3=X3) Results=boot(data=Data,statistic=Betas,R=BootTime,maxit=100, parallel="multicore",formula=Voxel~X1+X2+X3) X1P[該點]=(1+sum(abs(Results$t[,1])>abs(Results$t0[1])))/(2*BootTime) X2P[該點]=(1+sum(abs(Results$t[,2])>abs(Results$t0[2])))/(2*BootTime) X3P[該點]=(1+sum(abs(Results$t[,3])>abs(Results$t0[1])))/(2*BootTime) } -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 140.112.121.113