Re: [問題] bootstrap regression test + parallel

看板R_Language作者coccolegacy (我要告解)時間11年前 (2013/09/28 02:26)推噓0(0推 0噓 0→)

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

很久以前碰的有說錯的話還起各位大大指正 ※ 引述《skylikewater (choc.)》之銘言： : [問題類型]: : 程式諮詢(我想用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) 算出來的結果會過於顯著， : 　網路上可以找到的資源都看了，請各位指點怎麼實作、或推薦我補強什麼統計概念。 bootstrap最常見的用途，就是用來求一些分配未知或不容易推導的統計量的抽樣分配。可以先參考你這篇文章的第三頁的bootstrap confidence interval，其中的bootstrap percentile interval，就是bootstrap sample所計算出的統計量的empirical quantiles。至於模擬檢定，引用p-value wikipedia的定義如下： In statistical significance testing the p-value is the probability of obtaining a test statistic at least as extreme as the one that was actually observed, assuming that the null hypothesis is true. 他這邊算的，就是從每個bootstrap sample計算出來的檢定統計量大於全部樣本算的檢定統計量的機率 : 2. 多變項的 sample() 實作？ : 　有可能直接的用 sample() 直接取出放回「多個變項」嗎？ : 　經過實測結果很怪（p 都等於 1）。 : 說實話我也不確定在迴歸的架構下 : 　取出放回是不是「需要以同行/同個案/同樣本，來取出放回」 : 　還是可以自由地每個變項分開抽一般的做法會整筆一起抽其中一種做法可以對data.frame的row index抽樣，例如： data(iris) index = sample(1:nrow(iris), 1000, replace = T) bootstrap_sample = iris[index, ] : 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) : } 選項裡的parallel和snow都是R的平行運算套件有賴這裡其他大大們分享了 -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 219.86.165.219 ※ 編輯: coccolegacy 來自: 219.86.165.219 (09/28 02:35) ※ 編輯: coccolegacy 來自: 219.86.165.219 (09/28 03:09)