3

我在 R 中运行一个简单的遗传漂移模拟。

# Population size
N<-5000
# Number with focal allele
X1<-(N/2)
# Number of generations
ngens<-(2000)
# Number of replicates
nreps<-10

# Drift function
drift <- function(N, X1, ngens, nreps) {
# Makes a matrix of NA's of nreps columns, and ngen rows
       p <- matrix(NA, nrow=ngens, ncol=nreps)
  # Set base population
       p[1,] <- X1/N
  # Repetitive sampling function, each generation sample 10 times from the generation before (gen-1)  
       for(gen in 2:ngens)
         p[gen,] <- rbinom(n=nreps, size=N, prob=p[gen-1,]) / N
       p
}
# Run function "drift" & output as data frame
p <- data.frame(drift(N, X1, ngens, nreps))
# Plot
matplot(p, type="l", ylim=c(0,1), lty=1, xlab="Generations", ylab="Proportion Focal",col="grey")
# Mean value
p$mean<-apply(p[,c(1:10)],1,mean)
matplot(p$mean, type="l", ylim=c(0,1), lty=1, xlab="Generations", ylab="Proportion Focal",col="black",add=T)

我想：

  1. 围绕每一代的平均值计算置信区间（通过自举）
  2. 在我的数据框中添加一对具有上下置信区间估计值的列，然后可以将其绘制在 matplot 上，就像我对平均值所做的那样

谁能建议一种方法来做到这一点？我知道我需要 Boot 包并且大致知道如何使用它，但指导会很好。

问题（对我来说）是得到一个循环，该循环为每一代模拟生成一个 CI 并将其粘贴到“p”数据帧

编辑：

我已经按照@bakyaw 的部分建议尝试了这个，并且“for”循环改编自我曾经使用过的旧脚本。

myBootFunction<-function(x){
  b <- boot(x, function(u,i) mean(u[i]), R = 999)
  boot.ci(b, type = c("norm"))
}

meanList<-apply(p[c(2:ngens),c(1:nreps)],1,function(x)myBootFunction(x))
for(i in 1:49) {
  low<-meanList[[i]][[4]][[2]]
  high<-meanList[[i]][[4]][[3]]
  CIMatrix<-cbind(high,low)}

注意添加了 c(2:ngens)，没有它就会出现这个错误。

[1] “t 的所有值都等于 0.5 \n 无法计算置信区间”

然而，这仍然只是将 CIMatrix 创建为一个 1x2 双矩阵，而不是每一代都在其中。

4

1 回答 1

0

如果您已经解决了问题，这可以作为答案。开机功能：

    library(boot)

    myBootFunction<-function(x){
        b <- boot(x, function(u,i) mean(u[i]), R = 999)
        boot.ci(b, type = c("norm", "basic", "perc"))
    }

然后代替代码中的这一行：p$mean<-apply(p[,c(1:10)],1,mean)

您可以使用 ：

meanList<-apply(p[,c(1:10)],1,function(x)myBootFunction(x))

获得具有置信区间的列表后，您可以将其转换为数据框，然后对其进行处理。

于 2012-11-13T16:05:45.900 回答