R中的分层引导，具有> 25层

Question

我有大约 25 个不同组的数据。为了了解如果我有不同的样本量，每组的方差会如何变化，我正在尝试进行分层引导。例如，在样本大小为 5 时，它应该为每组生成 1000 个包含 5 个重采样点的集合。我喜欢在每组 5 到 30 个可能的范围内尽可能收集最小的样本量。

我遇到的问题是我必须对每个组进行子集化并在各个组上运行引导，然后将 R 输出复制并传递到 excel 中。 （我在 R 和如何编码方面相当熟悉）。时间太长了。我需要自动化引导以识别组，并以某种方式将 1000 个组的集合的统计信息保存到数据框中。这有意义吗？

这是我到目前为止的代码：....

#sample data
set.seed(1234)
df <- data.frame(g.name = as.factor(sample(c(LETTERS),100, replace = T)),
            C.H = as.numeric(sample(c(1:9),100, replace=T)))

#subset data by group... here only a three examples
Agroup=subset(df,C.H=='A')
Bgroup=subset(df,C.H=='B')
Cgroup=subset(df,C.H=='C')

#Bootstrap selecting a sample size of "i", "B" number of times. i.e. I am 
selecting sample sizes from 5 to 30, 1000 times each. I then apply var() to 
the sample, and take the multiple variances(or the variance of the 
variances). C.H is the measurement ranging from 1 to 9.  

B=1000
cult.var=(NULL)
for (i in 5:30){
 boot.samples=matrix(sample(Agroup$C.H,size=B*i, 
replace=TRUE),B,i)
  cult.var[i]=var(apply(boot.samples,1,var))
}
print(cult.var)

这可行，但需要大量复制和粘贴。我想我需要使用 for 循环来按组进行引导或找出其他方法。我确实找到了一种无需引导即可自行进行分层抽样的方法。所以也许我可以弄清楚如何以某种方式重复 1000 次......

使用函数boot() 的example here 不适合我的情况。我已经摆弄了一点，但无济于事。我不确定如何编写函数，这也可能是我无法弄清楚的原因。

Answer 1

这是一个尝试......

# generating data
set.seed(1234)
df <- data.frame(g.name = as.factor(sample(c(LETTERS),100, replace = T)),
                 C.H = as.numeric(sample(c(1:9),100, replace=T)))

boot.samples <- with(df, tapply(C.H, g.name, function(x) lapply(5:30, function(i) replicate(1000, sample(x,size=i,replace=T)))))

str(boot.samples$A)
## List of 26
##  $ : num [1:5, 1:1000] 7 7 3 7 7 7 3 3 2 7 ...
##  $ : num [1:6, 1:1000] 7 2 2 2 3 7 7 2 2 7 ...
##  $ : num [1:7, 1:1000] 2 3 2 7 2 3 7 2 3 3 ...
##  $ : num [1:8, 1:1000] 7 7 3 3 3 7 2 7 7 3 ...
##  $ : num [1:9, 1:1000] 2 2 2 7 2 7 3 3 3 7 ...
## ...and so on

variances <- lapply(boot.samples, function(y) sapply(y, function(x) apply(x, 2, var)))
    str(variances)
## List of 26
##  $ A: num [1:1000, 1:26] 3.2 5.8 6.2 3.2 0.3 4.8 5 5.8 6.7 3.2 ...
##  $ B: num [1:1000, 1:26] 3.2 0.8 4.7 5.3 5.3 5.3 1.2 4.7 4.2 3.8 ...
##  $ C: num [1:1000, 1:26] 9 4.8 2.7 9.8 8.3 9.8 10.2 10.2 9 12.3 ...
##  $ D: num [1:1000, 1:26] 8.3 7.5 9.8 3.8 3.5 3.5 5.7 3.7 6.7 3.2 ...
## ...and so on

variancesvariances <- lapply(variances, function(x) apply(x, 2, var))
str(variancesvariances)
## List of 26
##  $ A: num [1:26] 3.15 2.27 1.53 1.3 1.03 ...
##  $ B: num [1:26] 4.32 3.54 2.83 2.46 2.09 ...
##  $ C: num [1:26] 13.06 10.08 8.46 6.98 5.59 ...
##  $ D: num [1:26] 4.9 3.7 3.02 2.39 2.07 ...
## ...and so on

似乎随着广告中样本量的增加而下降......让我们拍一张漂亮的照片

cols <- rainbow(26)
plot(NA, xlim=c(1,26), ylim=c(0,max(unlist(variancesvariances))))
for(i in 1:26) {
  lines(variancesvariances[[i]], col=cols[i])
  text(1,variancesvariances[[i]][1],names(variancesvariances)[i],col=cols[i])
}

请注意，这可以转换为带有as.data.frame(variancesvariances) 的data.frame。

这次我收到了吗？