使用R引导多列数据

Question

我在R方面比较新，我想建立一个函数，它可以循环遍历导入的表格中的列，并生成由均值和95％置信区间组成的输出。理想情况下，应该能够用不同样本大小对列进行自举，但是首先我想让迭代工作起来。我有一些好像有点用的东西，但我无法完全实现它。代码如下，包括示例数据和输出：

#cdata<-read.csv(file.choose(),header=T)#read data from selected file, works, commented out because data is provided below
#cdata #check imported data

#Sample Data
#   WALL NRPK CISC WHSC LKWH YLPR
#1    21    8    1    2    2    5
#2    57    9    3    1    0    1
#3    45    6    9    1    2    0
#4    17   10    2    0    3    0
#5    33    2    4    0    0    0
#6    41    4   13    1    0    0
#7    21    4    7    1    0    0
#8    32    7    1    7    6    0
#9     9    7    0    5    1    0
#10    9    4    1    0    0    0

x<-cdata[,c("WALL","NRPK","LKWH","YLPR")] #only select relevant species

i<-nrow(x) #count number of rows for bootstrapping 
g<-ncol(x) #count number of columns for iteration

#build bootstrapping function, this works for the first column but doesn't iterate

bootfun <- function(bootdata, reps) {

  boot <- function(bootdata){

    s1=sample(bootdata, size=i, replace=TRUE)
    ms1=mean(s1)
    return(ms1)

  } # a single bootstrap

  bootrep <- replicate(n=reps, boot(bootdata))

  return(bootrep)

} #replicates bootstrap of "bootdata" "reps" number of times and outputs vector of results

cvr1 <- bootfun(x$YLPR,50000) #have unsuccessfully tried iterating the location various ways (i.e. x[i])
cvrquantile<-quantile(cvr1,c(0.025,0.975))
cvrmean<-mean(cvr1)
vec<-c(cvrmean,cvrquantile) #puts results into a suitable form for output
vecr<-sapply(vec,round,1) #rounds results
vecr

      2.5% 97.5% 
 28.5  19.4  38.1 

#apply(x[1:g],2,bootfun) ##doesn't work in this case

#desired output:

#Species    Mean LowerCI UpperCI
#WALL       28.5    19.4      38.1
#NRPK       6.1 4.6    7.6
#YLPR       0.6 0.0    1.6

我也尝试过使用boot包，很好地通过均值进行迭代，但是我无法在置信区间上做到同样的效果。上面的“普通”代码还有一个优点，就是您可以轻松地检索引导结果，这些结果可能用于其他计算。为了完整起见，这是boot代码：

#Bootstrapping using boot package
library(boot)
#data<-read.csv(file.choose(),header=TRUE) #read data from selected file
#x<-data[,c("WALL","NRPK","LKWH","YLPR")] #only select relevant columns
#x #check data

#Sample Data

#  WALL NRPK LKWH YLPR
#1    21    8    2    5
#2    57    9    0    1
#3    45    6    2    0
#4    17   10    3    0
#5    33    2    0    0
#6    41    4    0    0
#7    21    4    0    0
#8    32    7    6    0
#9     9    7    1    0
#10    9    4    0    0

i<-nrow(x) #count number of rows for resampling 
g<-ncol(x) #count number of columns to step through with bootstrapping
boot.mean<-function(x,i){boot.mean<-mean(x[i])} #bootstrapping function to get the mean

z<-boot(x, boot.mean,R=50000) #bootstrapping function, uses mean and number of reps
boot.ci(z,type="perc") #derive 95% confidence intervals
apply(x[1:g],2, boot.mean) #bootstrap all columns

#output:
#WALL NRPK LKWH YLPR 
#28.5  6.1  1.4  0.6 

我已经查看了所有可以找到的资源，但似乎无法让事情顺利进行。我的输出需求是每列的加权平均值和相应的置信区间。谢谢！

- JK101

---

当你说一个输出由均值和置信区间组成时，你的意思是什么？你是指想要计算出的统计量的平均值、0.025和0.975百分位数吗？ - alexwhitworth

此外，对于50k个副本*n列，您可能希望使用boot::boot，因为它提供了与parallel包的集成。 - alexwhitworth

当我说输出时，我的意思是我希望结果显示在一个表格中，该表格由每列的列名、自举均值和相关置信区间组成。感谢您对boot::boot的建议。 - JK101

你的回复根本没有回答我的问题，但很高兴你找到了需要的答案。 - alexwhitworth

很抱歉，我可能误解了你的问题。在我的代码示例中，我提供了所需输出的示例，其中显示了带有95％置信区间的自助法平均值。请原谅我，我感冒了（并且已经服用了感冒药），所以显然出了点差错。感谢您的尝试帮助，非常感激。 - JK101

Answer 1

注意：apply(x[1:g],2, boot.mean) #bootstrap all columns并没有进行任何自助法。你只是计算每列的平均值。

要进行自助平均和置信区间，请尝试以下内容：

apply(x,2,function(y){ 
   b<-boot(y,boot.mean,R=50000); 
   c(mean(b$t),boot.ci(b,type="perc", conf=0.95)$percent[4:5])
})

- fishtank

---

这个完美地运作了。我的代码还是有一些奇怪的小错误（不是你的问题），但是你提供的内容产生了我所期望的完全正确的结果。非常感谢你！ :) - JK101