作者skylikewater (choc.)
看板R_Language
標題[問題] bootstrapping regression.
時間Thu Sep 26 01:32:21 2013
[問題類型]:
程式諮詢(我想用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
1F：推 Wush978:(無關)看到這個選模問題，讓我想到這篇文章 09/26 23:03
3F：→ Wush978:簡單來說，就是不推薦用檢定方式選擇 X 09/26 23:04
4F：→ celestialgod:可以用snowfall 可以參考本版#1Hocf0cY 09/28 21:29
5F：→ skylikewater:感謝回答的幾位前輩 我如果比較成功了再跟大家分享 09/29 19:55