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. 平均値で行ったのと同じように、マットプロットにプロットできる上下の信頼区間の推定値を使用して、データフレームに列のペアを追加します

誰でもこれを行う方法を提案できますか? 私は 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 信頼区間を計算できません」

ただし、これは、すべての世代を含むものではなく、1x2 の double マトリックスとして CIMatrix を作成するだけです。

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 に答える