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