作者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