R> 25分層的分層引導

Question

我有大約25個不同組的數據。爲了瞭解如果我有不同的樣本規模，每個組的方差會發生怎樣的變化，我正在嘗試進行分層自舉。例如，在樣本大小爲5時，它應該爲每個組生成1000個5個重採樣點的集合。我喜歡根據需要收集最小樣本量，可能範圍爲每組5至30個。

我遇到的問題是我必須對每個組進行子集分類，然後在各個組上運行bootstrapping，然後將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次...

example here使用函數boot()不符合我的情況。我已經擺弄了一下，但無濟於事。我不知道如何編寫函數，這也可能是我無法弄清楚的原因。

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)的數據框。

這次我收到了嗎？