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 个观测值的平均值 ( k=10) 时，它会根据概率分布随机抽样 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。

现在，我不想对所有模拟使用一个频率分布，而是根据每个数据点的观察数量有条件地使用不同的频率分布。

我制作了一个表格，cond_probs其中包含我的每个真实数据点的一行。有一列包含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观察值，其值为0.1。0.1因此，观察的频率是1。因为A1, 有113观察结果，其中大多数 ( 0.664) 具有价值0.1。这个想法cond_probs类似于ps，但cond_probs每个数据点都有一个概率分布，而不是所有数据的一个概率分布。

我想修改上面的代码，以便修改采样以使用cond_probs而不是ps频率分布。并且在选择从哪一行中采样时，使用观察的数量k作为标准。cond_probs所以它会像这样工作：

对于具有k观察次数的数据点：

查看cond_probs表格并随机选择total观察数量与 k: 大小相似的行0.9k-1.1k。如果不存在这样的行，请继续。

一旦选择了一个数据点，就使用该行的概率分布，cond_probs就像ps在原始代码中使用的一样，随机抽样k观察的数量并输出这些观察的平均值。

对于 的每次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)

第一个示例十行（约 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 个。

例如：

> sims_test[[31]]
[1] 0.1

并且sims_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 个站点，这表明某种原因模式。

如果这有帮助，我很乐意私下分享整个文件，因为我不知道下一步该做什么来进行故障排除。

出现的另一个问题是我不确定代码是否按预期工作。我制作了一个虚拟 cond_probs 文件，其中有两个数据点的“总数”为“1”观察：

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”观察值的数据点，因此得到大约 50% 的观察值的平均值为“0.2”，50% 的观察值的平均值为“0.6”。但是平均值始终为 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修复了此问题。

cond_probs我认为保持原样并选择分布sampleSize次数的另一个原因可能是最好的解决方案：选择分布的概率应该与其在cond_probs. 如果我们将分布结合起来，选择一个分布的几率total 9将10不再取决于这些总数的观察数量。示例：如果存在90带有total=10和的分布10，total=9则应该有90%机会选择带有 的分布total=10。如果我们结合分布，选择 'total'= 9 或 10 的分布的几率不会变成 50/50（这不是理想的）吗？

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

合适的分布是第三种：

> 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 回答