如何在R中执行嵌套重采样 - 多样性指数不确定性

Question

我目前正在尝试应用 Herrmann 等人中描述的重采样方法。 (2022) 具有外部和内部重采样循环： 我的数据集由 100 行 x 60 列组成。行对应于尝试 ID（每行一次尝试），列对应于物种（每列一个物种）。每个单元格代表每次尝试的每个物种 i 的个体数量 K = NiK

NiK 可能会因重采样而变化，但给定尝试中所有混杂物种的个体总数（称为 NK，对应于行总和）需要保持不变。

让我们生成一个类似的数据框：

df <- matrix(sample(0:15, 6000, replace = T), ncol = 60)
df

这是我用于外循环的代码：

#Outer loop
#Resampling among the 100 tries 1000 times
seed(1234)
outer = matrix(replicate(1000, sample(1:nrow(df), nrow(df), replace = T)), ncol = 100)

现在添加内部循环，我被困住了。 内部循环在每次尝试中对给定物种的个体数量进行重新采样，而外部循环则在内部循环生成的新数据集中对尝试进行重新采样：

i = 物种数

r = 外部重采样 ID

inner <- NULL
complete <- NULL
for (r in 1:1000) {
  for (i in 1:ncol(df)) {
    inner = cbind(inner,sample(df[,i], nrow(df), replace = T))
  }
  complete = rbind(complete, inner[outer[r,],])
  inner <- NULL
}
complete

该循环可以很好地通过内部和外部循环生成 100 次尝试的 1000 次重采样，但我目前不知道如何限制每一行的总和为 NK，这可以针对给定的尝试进行计算K :

rowSums(df[K,])

有谁知道一个函数或方法可以允许在使用内部循环在每列中重新采样时保持 NK 恒定吗？

预先感谢您的帮助，希望我的英语和我的解释是可以理解的。

Answer 1

使用

rowSums(df)

和

rmultinom

，我们可以对一行重新采样并最终得到相同的总数。示范：

rs <- rowSums(df)
all(replicate(100, sum(rmultinom(1, rs[1], df[1,])) == rs[1]))
#> [1] TRUE

由于要对内部样本求和，因此对于每个外部样本，我们可以对多个行重采样进行行求和。

set.seed(1234)
K <- 100L # "tries"
N <- 60L # species
R <- 1e3L # resamples
df <- matrix(sample(0:15, N*K, replace = 1), K, N)
rs <- rowSums(df)

outer <- mapply( # outer sampling
  \(x) rowSums(
    # inner sampling
    mapply(\(i) rmultinom(1, rs[i], df[i,]), sample(nrow(df), K, 1))
  ),
  1:R
)

dim(outer)
#> [1]   60 1000

outer

的每一列现在对应于外部重采样。

Answer 2

谢谢jblood94的快速答复！

因为我已经找到了解决方案，所以我也可以为回答我自己的问题做出贡献。 我尝试用错误的数据框解决问题。

确实，我的原始数据集看起来更像：

ID <- sort(rep(1:100,times=c(sample(1:10, 100, replace = T))))
df <- data.frame(ID, Species.name= sample(letters[1:10], length(ID), replace = T))
df

with ID = 尝试 ID 和 Species.命名物种的标识。

使用

acast()

我获得了与我原来的帖子中描述的类似的数据框：

 library(reshape2)
 acast(df, ID~Species.name)

适合计算多样性指数。

尽管如此，使用 df，在保持每次尝试的个体数量的同时更容易执行嵌套采样，因为一个个体等于一行！

工作代码：

library(reshape2)
library(dplyr)
             # First extract number of individuals per try to keep it constant
            groupsize = df %>%
              group_by(ID)%>%
              summarise(n = n())
            groupsize
            
            #Test inner resampling on new dataframe
            df2 = df %>%
              left_join(groupsize, by = "ID") %>%
              group_by(ID) %>%
              mutate(samp = sample(Species.name, replace = T)) %>%
              filter(length(samp) == n) %>%
              ungroup()
            df2
            #Seems to work. For a given try, only the initially observed species are resampled 
            #to constitute the new resampled try, while keeping the number of individual.
            #Now we only need to cast df2 and select row (tries) from a resampled list of tries (outer resampling)
            
            #Lets creating the loop with inner + outer resampling
            
            inner <- NULL
            outer <- NULL
            complete <- NULL
              for (i in 1:1000) {
                inner = df %>%  left_join(groupsize, by = "ID") %>%
                                group_by(ID) %>%
                                mutate(Sp = sample(Species.name, replace = T)) %>%
                                filter(length(Sp) == n) %>%
                                ungroup() #inner resampling
                inner = data.frame(ID = rownames(acast(inner, ID ~ Sp, fun.aggregate = length)), 
                acast(inner, ID ~ Sp, fun.aggregate = length), bloc = i) #rotate df from inner resampling
                outer = sample(1:100,100, replace = T)
                complete = rbind(complete, inner[outer,]) #outer resampling
                }
            complete
            nrow(complete)

然后我们获得一个 1000 x 100 行的数据框，其中列数对应于物种，非常适合计算多样性指数。 烦人的是来自

acast()

的错误消息，但这并不影响结果。

遵循 Herrmann 等人。 （2022），我们现在需要：1）对每个块的列进行求和以获得 1000 行，重采样尝试中的每个物种池计数，2）然后计算每个多样性指数的 1000 个值，并获得 Efron 95% CI。

希望这个帖子能帮助人们进行类似的分析！