作者coccolegacy (我要告解)
看板R_Language
標題Re: [問題] bootstrap regression test + parallel
時間Sat Sep 28 02:26:36 2013
很久以前碰的
有說錯的話還起各位大大指正
※ 引述《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)