关于原始值和引导值差异的说明答案

【问题标题】：Explanation about difference in original and bootstrap value关于原始值和引导值差异的说明
【发布时间】：2016-06-07 01:41:43
【问题描述】：

我有以下功能；（1）计算我拥有的每个变量的偏差差异，（2）对我在第一步中计算的偏差差异进行引导

set.seed(1001)

xfunction <- function(d,i)
{
glm.snp1 <- glm(PHENOTYPE~d[i], family="binomial", data=training1)
null <- glm.snp1$null.deviance
residual <- glm.snp1$deviance
dDeviance <- null-residual
return(dDeviance)
}

myboot <- function(d)
{
boot(d,xfunction, R=1000)
}

result <- lapply(training1,function(x)myboot(x))

所以基本上从结果中我得到了原始 dDeviance 的值（没有引导程序），我可以从引导程序计算平均值（dDeviance）。我需要帮助来解释为什么原始引导值和平均引导值差异太大？例如，对于其中一个变量，dDeviance 的原始值为 0.024，而 dDeviance 的引导平均值为 0.000412。

【问题讨论】：

这似乎不太可能，因为自举统计的分布主要围绕原始数据中的统计。查看您的代码，我认为您的问题在于索引。您确实使用了 d[i]，但仍使用原始训练集和原始表型。
感谢@Heroka 的评论。我也同意你的观点，该值应该在原始值附近。我试图摆脱 xfunction 中的数据，但它无法识别我的 PHENOTYPE 来自哪里。如果您能对我的职能提出修改建议，我将不胜感激。供您参考，我的结果是 PHENOTYPE，预测变量由 1500 个 SNP 和 1000 个样本组成。
根据我的引导经验，我有一个函数可以获取一个包含我需要的所有内容和索引的数据集，然后使用它。因此，表型应该在训练中1，然后您可以使用 data=d[i,] 拟合您的模型。
我同意@Heroka，您应该在引导期间而不是在公式中子集 data.frame。我也不明白你用lapply 循环在那里做什么。
感谢两位cmets。我会努力解决这个问题。好吧，我提出了 lapply 的想法，以便它将计算 training1[,2:1501] 中每个列的 dDeviance 并使用引导函数引导该值。

标签： r resampling statistics-bootstrap

【解决方案1】：

正如 cmets 中所指出的，最好对 data.frame 索引进行子集化，以使其与 boot 兼容。如果你想遍历不同的变量，并应用这个函数，最好在迭代中进行。

对于每个引导程序，您都可以拟合不同的模型并提取值。

所以我们重写了函数，它接受了一个额外的公式参数indep_var，它指定了要针对 PHENOTYPE 回归的列：

xfunction <- function(d,ind,indep_var){

    dDeviance = sapply(indep_var,function(V){
         f = reformulate(V,response="PHENOTYPE")  
         glm.snp1 <- glm(f, family="binomial", data=d[ind,])
         glm.snp1$null.deviance-glm.snp1$deviance
    })
return(dDeviance)
}

我们可以建立一个示例数据集：

set.seed(111)
training1 = data.frame(PHENOTYPE=rbinom(100,1,0.5),matrix(rnorm(500),ncol=5))

head(training1)
  PHENOTYPE         X1          X2         X3         X4         X5
1         1  0.1916634 -0.09152026 -0.9685094 -1.0503824 -0.9137690
2         1  1.5525443 -1.87430581  0.9119331  0.3251424  0.1126909
3         0  0.9142423 -0.66416202  0.0928326 -2.1048716 -2.3597249
4         1  0.3586254  0.20341282 -0.5329309 -0.9551027 -1.5693983
5         0  0.1750956 -2.59444339 -1.6669055 -0.5306399  1.2274569
6         0 -0.8472678 -0.09373393 -0.5743455  0.8274405  0.7620480

然后在数据上试试这个：

xfunction(training1,1:nrow(training1),indep_var=c("X1","X2","X3"))
       X1        X2        X3 
0.1189847 0.2512539 0.2684927

然后使用引导程序：

library(boot)
boot(training1,xfunction,R=1000,indep_var=c("X1","X2","X3"))

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = training1, statistic = xfunction, R = 1000, indep_var = c("X1", 
    "X2", "X3"))


Bootstrap Statistics :
     original   bias    std. error
t1* 0.1189847 1.033286    1.564971
t2* 0.2512539 1.047503    1.863012
t3* 0.2684927 1.005062    1.747959

【讨论】：