我有大约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()
的示例不适合我的情况。我摆弄了一下,但无济于事。我不确定如何编写函数,这可能也是我无法弄清楚的原因。
这里有一个刺
...# 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)
的数据帧。
这次我明白了吗?