4

漏斗プロットを使用して実際のデータの信頼帯を作成するために、ブートストラップを介してデータをシミュレートしようとしています。以前の質問に対する受け入れられた回答の戦略に基づいて構築して います。データをシミュレートするために単一の確率分布を使用する代わりに、シミュレートされるデータの部分に応じて異なる確率分布を使用するように変更したいと考えています。

質問に答えるのを手伝ってくれる人、または質問をより明確に表現するのを手伝ってくれる人に感謝します。

私の問題は、適切な R コードを記述して、より複雑な形式のデータ シミュレーションを行うことです。

現在のコードは次のとおりです。

n <- 1e4
set.seed(42)
sims <- sapply(1:80, 
               function(k) 
                 rowSums(
                   replicate(k, sample((1:7)/10, n, TRUE, ps))) / k)

このコードは、各データ ポイントが1:80観測間の平均である値を持つデータをシミュレートします。たとえば、データ ポイントの値が 10 回の観測値の平均 ( =10) である場合、確率分布にk基づいて 10 個の値 (0.1、0.2、0.3、0.4、0.5、0.6 または 0.7 のいずれか) をランダムにサンプリングします。psこれにより、各値の確率が得られます (経験的分布全体に基づく)。

ps は次のようになります。

ps <- prop.table(table((DF$mean_score)[DF$total_number_snps == 1]))
#        0.1         0.2         0.3         0.4         0.5         0.6         0.7 
#0.582089552 0.194029851 0.124378109 0.059701493 0.029850746 0.004975124 0.004975124

たとえば、観測値が である確率は0.1です0.582089552。

ここで、すべてのシミュレーションに 1 つの度数分布を使用する代わりに、各データポイントの基礎となる観測の数に応じて、条件付きで異なる度数分布を使用したいと考えています。

cond_probs実際のデータ ポイントごとに 1 つの行を含むテーブル を作成しました。オブザベーションの数を示す列と、total各オブザベーションの各値の頻度を示す列があります。

cond_probs テーブルの例:

gene_name   0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 total
A1  0.664   0.319   0.018   0.000   0.000   0.000   0.000   0.000   0.000   113.000
A2  0.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   1.000

したがって、データ ポイントについては、値が である観測A2のみがあります。したがって、観測頻度は です。については、観察があり、それらの大部分 ( ) は値を持っています。考え方は に似ていますが、すべてのデータに対して 1 つの確率分布ではなく、各データ ポイントに対して確率分布があります。10.10.11A11130.6640.1cond_probspscond_probs

上記のコードを変更して、頻度分布のcond_probs代わりに使用するようにサンプリングを変更したいと思います。また、どの行からサンプリングするかを選択する際の基準として、ps観測数 を使用します。したがって、次のように機能します。kcond_probs

k観測数のあるデータ ポイントの場合:

テーブルを見て、観測値の数が k: と同じサイズcond_probsの行をランダムに選択します。そのような行が存在しない場合は、続行します。total0.9k-1.1k

データポイントが選択されると、元のコードで使用されているのとcond_probs同じように、その行の確率分布を使用して、観測数をランダムにサンプリングし、これらの観測値の平均を出力します。psk

nの反復ごとに、の値が( )の現在の値と類似しているすべての行replicateから、 から新しいデータ ポイントを置換してランダムにサンプリングします。cond_probstotalk0.9k-1.1k

このデータセットでは、データ ポイントの基礎となる観測の数に基づいて、どの確率分布を使用するかを調整する必要があります。これは、このデータセットでは、観察の確率が観察数の影響を受けるためです (SNP が多い遺伝子は、遺伝的連鎖と背景選択により、観察ごとのスコアが低くなる傾向があります)。

以下の回答を使用して更新:

以下の回答を使用してみましたが、例のシミュレートされた cond_probs データでは機能しますが、実際の cond_probs ファイルでは機能しません。cond_probs ファイルをインポートしてマトリックスに変換しました

cond_probs <- read.table("cond_probs.txt", header = TRUE, check.names = FALSE)
cond_probs <- as.matrix(cond_probs)

最初の例の 10 行 (~20,000 行のうち) は次のようになります。

>cond_probs
       total   0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1.0
[1,]     109 0.404 0.174 0.064 0.183 0.165 0.009 0.000 0.000 0.000 0.000
[2,]     181 0.564 0.221 0.144 0.066 0.006 0.000 0.000 0.000 0.000 0.000
[3,]     289 0.388 0.166 0.118 0.114 0.090 0.093 0.028 0.003 0.000 0.000
[4,]     388 0.601 0.214 0.139 0.039 0.008 0.000 0.000 0.000 0.000 0.000
[5,]     133 0.541 0.331 0.113 0.000 0.008 0.008 0.000 0.000 0.000 0.000
[6,]     221 0.525 0.376 0.068 0.032 0.000 0.000 0.000 0.000 0.000 0.000
[7,]     147 0.517 0.190 0.150 0.054 0.034 0.048 0.007 0.000 0.000 0.000
[8,]     107 0.458 0.196 0.252 0.084 0.009 0.000 0.000 0.000 0.000 0.000
[9,]      13 0.846 0.154 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

私が実行した場合：

sampleSize <- 20
set.seed(42)
#replace 1:80 with 1: max number of SNPs in gene in dataset
sims_test <- sapply( 1:50, simulateData, sampleSize )

x個の観測値でサンプリングの平均を見てください。20個あるはずなのに、1つの結果しか得られません。

例えば：

> sims_test[[31]]
[1] 0.1

Andsims_testは と同じ順序ではありませんsims:

>sims_test
   [,1] [,2]      [,3]  [,4] [,5]      [,6]      [,7]   [,8]      [,9]
 [1,]  0.1  0.1 0.1666667 0.200 0.14 0.2666667 0.2000000 0.2375 0.1888889
 [2,]  0.1  0.1 0.1333333 0.200 0.14 0.2333333 0.1571429 0.2625 0.1222222
 [3,]  0.1  0.1 0.3333333 0.225 0.14 0.1833333 0.2285714 0.2125 0.1555556
 [4,]  0.1  0.1 0.2666667 0.250 0.10 0.1500000 0.2000000 0.2625 0.2777778
 [5,]  0.1  0.1 0.3000000 0.200 0.16 0.2000000 0.2428571 0.1750 0.1000000
 [6,]  0.1  0.1 0.3666667 0.250 0.16 0.1666667 0.2142857 0.2500 0.2000000
 [7,]  0.1  0.1 0.4000000 0.300 0.12 0.2166667 0.1857143 0.2375 0.1666667
 [8,]  0.1  0.1 0.4000000 0.250 0.10 0.2500000 0.2714286 0.2375 0.2888889
 [9,]  0.1  0.1 0.1333333 0.300 0.14 0.1666667 0.1714286 0.2750 0.2888889

更新 2

を使用cond_probs <- head(cond_probs,n)して、コードが n = 517 まで機能すると判断した後、これよりも大きいすべてのサイズに対して、上記と同じ出力が生成されます。これがファイル自体の問題なのか、メモリの問題なのかはわかりません。行 518 を削除し、行を数回前に複製してより大きなファイルを作成すると、行自体が問題を引き起こしていることがわかりました。518 行目は次のようになります。

9.000   0.889   0.000   0.000   0.000   0.111   0.000   0.000   0.000   0.000   0.000

私は別の4つの問題のある行を見つけました:

9.000   0.444   0.333   0.111   0.111   0.000   0.000   0.000   0.000   0.000   0.000

9.000   0.444   0.333   0.111   0.111   0.000   0.000   0.000   0.000   0.000   0.000

9.000   0.111   0.222   0.222   0.111   0.111   0.222   0.000   0.000   0.000   0.000

9.000   0.667   0.111   0.000   0.000   0.000   0.222   0.000   0.000   0.000   0.000

私はそれらについて何も異常に気づきません。それらはすべて「合計」で 9 つのサイトを持っています。これらの行を削除して、これらの前の行のみを含む「cond_probs」ファイルを実行すると、コードが機能します。しかし、「cond_probs」全体がまだ機能しないため、他にも問題のある行があるはずです。

これらの問題のある行を小さな「cond_probs」ファイルに戻そうとしたところ、このファイルは機能するため、行に本質的に問題があるようには見えないため、非常に混乱しています。一方、全部で 9 つのサイトがあるという事実は、何らかの原因パターンを示唆しています。

トラブルシューティングのために次に何をすべきかわからないので、それが役立つ場合は、ファイル全体を非公開で共有できれば幸いです.

発生するもう 1 つの問題は、コードが期待どおりに機能しているかどうかわからないことです。「合計」が「1」の観測値を持つ 2 つのデータ ポイントがあるダミーの cond_probs ファイルを作成しました。

total   0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.000   0.000   0.000   0.000   0.000   0.000   1.000   0.000   0.000   0.000   0.000
1.000   0.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000

したがって、それらは両方とも「1」の観測値でデータポイントに対してサンプリングされるため、平均値が「0.2」で観測値の約 50%、平均が「0.6」で 50% になると予想されます。ただし、平均は常に 0.2 です。

sims_test[[1]]
 [1] 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2

10000 回サンプリングしても、すべての観測値は 0.2 であり、決して 0.6 ではありません。私のコードの理解では、観察ごとに同様のサイズの cond_probs から新しい行をランダムに選択する必要がありますが、この場合はそうではないようです。コードを誤解していますか、それとも入力が正しくないという問題がありますか?

cond_probs ファイル全体は、次のアドレスにあります。

cond_probs

更新 3

シミュレーションの実行時に変更sapplyすると、この問題が修正されました。lapply

もう 1 つの理由は、そのままcond_probsにして分布sampleSizeを何度も選択することが最善の解決策である可能性があることです。分布を選択する確率は、cond_probs. 分布を組み合わせると、分布を選択するオッズは、これらの合計の観測total 9数10に依存しなくなります。例: と を含むディストリビューションがある場合、を90含むディストリビューションを選択する機会があるはずです。分布を組み合わせると、'total'= 9 または 10 の分布を選択する確率は 50/50 になりませんか (これは理想的ではありません)。total=1010total=990%total=10

4

1 に答える 1

4

psから適切な分布を選択する関数を単純に書きましたcond_probs:

N <- 10  # The sampled values are 0.1, 0.2, ... , N/10
M <- 8   # number of distributions in "cond_probs"

#-------------------------------------------------------------------
# Example data:

set.seed(1)

cond_probs <- matrix(0,M,N)

is.numeric(cond_probs)

for(i in 1:nrow(cond_probs)){ cond_probs[i,] <- dnorm((1:N)/M,i/M,0.01*N) }

is.numeric(cond_probs)

total <- sort( sample(1:80,nrow(cond_probs)) )
cond_probs <- cbind( total, cond_probs/rowSums(cond_probs) )

colnames(cond_probs) <- c( "total", paste("P",1:N,sep="") )

#---------------------------------------------------------------------
# A function that chooses an appropiate distribution from "cond_prob",
# depending on the number of observations "numObs":

ps <- function( numObs,
                similarityLimit = 0.1 )
{
  similar <- which( abs(cond_probs[,"total"] - numObs) / numObs < similarityLimit )

  if ( length(similar) == 0 )
  { 
    return(NA)
  }
  else
  {
    return( cond_probs[similar[sample(1:length(similar),1)],-1] )
  }
}

#-----------------------------------------------------------------
# A function that simulates data using a distribution that is
# appropriate to the number of observations, if possible:

simulateData <- function( numObs, sampleSize )
{
  if (any(is.na(ps(numObs))))
  {
    return (NA)
  }
  else
  {
    return( rowSums(
              replicate(
                numObs,
                replicate( sampleSize, sample((1:N)/10, 1, prob = ps(numObs))))
            ) / numObs )
  }
}

#-----------------------------------------------------------------
# Test:

sampleSize <- 30
set.seed(42)
sims <- lapply( 1:80, simulateData, sampleSize )

の分布cond_probs:

    total           P1           P2           P3           P4           P5           P6           P7           P8           P9          P10
[1,]    16 6.654875e-01 3.046824e-01 2.923948e-02 5.881753e-04 2.480041e-06 2.191926e-09 4.060763e-13 1.576900e-17 1.283559e-22 2.189990e-28
[2,]    22 2.335299e-01 5.100762e-01 2.335299e-01 2.241119e-02 4.508188e-04 1.900877e-06 1.680045e-09 3.112453e-13 1.208647e-17 9.838095e-23
[3,]    30 2.191993e-02 2.284110e-01 4.988954e-01 2.284110e-01 2.191993e-02 4.409369e-04 1.859210e-06 1.643219e-09 3.044228e-13 1.182153e-17
[4,]    45 4.407425e-04 2.191027e-02 2.283103e-01 4.986755e-01 2.283103e-01 2.191027e-02 4.407425e-04 1.858391e-06 1.642495e-09 3.042886e-13
[5,]    49 1.858387e-06 4.407417e-04 2.191023e-02 2.283099e-01 4.986746e-01 2.283099e-01 2.191023e-02 4.407417e-04 1.858387e-06 1.642492e-09
[6,]    68 1.642492e-09 1.858387e-06 4.407417e-04 2.191023e-02 2.283099e-01 4.986746e-01 2.283099e-01 2.191023e-02 4.407417e-04 1.858387e-06
[7,]    70 3.042886e-13 1.642495e-09 1.858391e-06 4.407425e-04 2.191027e-02 2.283103e-01 4.986755e-01 2.283103e-01 2.191027e-02 4.407425e-04
[8,]    77 1.182153e-17 3.044228e-13 1.643219e-09 1.859210e-06 4.409369e-04 2.191993e-02 2.284110e-01 4.988954e-01 2.284110e-01 2.191993e-02

分布の平均:

> cond_probs[,-1] %*% (1:10)/10
          [,1]
[1,] 0.1364936
[2,] 0.2046182
[3,] 0.3001330
[4,] 0.4000007
[5,] 0.5000000
[6,] 0.6000000
[7,] 0.6999993
[8,] 0.7998670

31 回の観測のシミュレートされたデータの平均:

> sims[[31]]
 [1] 0.2838710 0.3000000 0.2935484 0.3193548 0.3064516 0.2903226 0.3096774 0.2741935 0.3161290 0.3193548 0.3032258 0.2967742 0.2903226 0.3032258 0.2967742
[16] 0.3129032 0.2967742 0.2806452 0.3129032 0.3032258 0.2935484 0.2935484 0.2903226 0.3096774 0.3161290 0.2741935 0.3161290 0.3193548 0.2935484 0.3032258

適切な分布は 3 番目のものです。

> ps(31)
          P1           P2           P3           P4           P5           P6           P7           P8           P9          P10 
2.191993e-02 2.284110e-01 4.988954e-01 2.284110e-01 2.191993e-02 4.409369e-04 1.859210e-06 1.643219e-09 3.044228e-13 1.182153e-17
于 2015-11-04T13:08:40.333 に答える