r - r 中的条件概率

Question

问题:

针对影响 0.05% 男性人口的疾病的筛查测试能够在 90% 的个体实际患有该疾病的病例中识别出该疾病。然而，该测试会产生 1% 的假阳性(当个体没有患病时给出阳性读数)。计算检测结果呈阳性的男性患有该疾病的概率。然后，计算一个人在检测结果呈阴性的情况下患有该疾病的概率。

我的错误尝试:

我首先让:• T 指一名男子检测结果呈阳性的事件• 男子检测结果呈阴性的事件• D 是一个人实际上患有该疾病的事件• Dc 是指男性没有患病

因此我们需要找到 P(D|T) 和 P(D|Tc)

然后我写了这段代码:

set.seed(110)
sims = 1000

D = rep(0, sims)
Dc = rep(0, sims)
T = rep(0, sims)
Tc = rep(0, sims)

# run the loop
for(i in 1:sims){
  
  # flip to see if we have the disease
  flip = runif(1)
  
  # if we got the disease, mark it
  if(flip <= .0005){
    D[i] = 1
  }
  
  # if we have the disease, we need to flip for T and Tc, 
  if(D[i] == 1){
    
    # flip for S1
    flip1 = runif(1)
    
    # see if we got S1
    if(flip1 < 1/9){
      T[i] = 1
    }
    
    # flip for S2
    flip2 = runif(1)
    
    # see if we got S1
    if(flip2 < 1/10){
      Tc[i] = 1
    }
  }
}


# P(D|T)
mean(D[T == 1])

# P(D|Tc)
mean(D[Tc == 1])

我真的很挣扎，所以任何帮助将不胜感激!

Answer 1

也许思考此类条件概率问题的最佳方法是通过具体示例。

假设我们对人群中的 100 万人进行了测试。那么预计有 500 人(百万人中的 0.05%)患有该疾病，其中 450 人预计检测呈阳性，50 人检测呈阴性(因为假阴性率为 10%)。

相反，预计 999,500 人不会患有这种疾病(100 万人减去确实患有这种疾病的 500 人)，但由于其中 1% 的检测呈阳性，因此我们预计有 9,995 人(999,500 中的 1%)出现假阳性结果。

因此，如果随机抽取一个阳性检测结果，它要么属于 450 名检测呈阳性的患病者之一，要么属于 9,995 名检测呈阳性的未患病者之一 - 我们不知道是哪一个。

这是第一个问题的情况，因为我们有一个阳性检测结果，但不知道它是真阳性还是假阳性。我们的受试者患有该疾病的概率给定他们的阳性检测结果是他们是 10,445 名阳性结果中的 450 名真阳性之一的概率(9995 假阳性 + 450 真阳性)。简单计算即可得出 450/10,445 或 0.043，即 4.3%。

同样，随机进行的阴性检测要么属于 989505 (999500 - 9995) 名未患病且检测呈阴性的人中的一名，要么属于 50 名患有的人 em> 检测结果呈阴性的疾病，因此患有该疾病的概率为 50/989505，即 0.005%。

我认为这个问题证明了在解释测试结果时需要考虑疾病患病率的重要性，而与编程或 R 关系不大。它只需要一个计算器(最多)。

如果您确实想在 R 中运行模拟，您可以这样做:

set.seed(1) # This makes the sample reproducible

sample_size <- 1000000 # This can be changed to get a larger or smaller sample

# Create a large sample of 1 million "people", using a 1 to denote disease and
# a 0 to denote no disease, with probabilities of 0.0005 (which is 0.05%) and
# 0.9995 (which is 99.95%) respectively.
disease <- sample(x = c(0, 1), 
                  size = sample_size, 
                  replace = TRUE, 
                  prob = c(0.9995, 0.0005))

# Create an empty vector to hold the test results for each person
test <- numeric(sample_size)

# Simulate the test results of people with the disease, using a 1 to denote
# a positive test and 0 to denote a negative test. This uses a probability of
# 0.9 (which is 90%) of having a positive test and 0.1 (which is 10%) of having
# a negative test. We draw as many samples as we have people with the disease
# and put them into the "test" vector at the locations corresponding to the
# people with the disease.
test[disease == 1] <- sample(x = c(0, 1), 
                             size = sum(disease), 
                             replace = TRUE, 
                             prob = c(0.1, 0.9))

# Now we do the same for people without the disease, simulating their test
# results, with a 1% probability of a positive test.
test[disease == 0] <- sample(x = c(0, 1), 
                             size = 1e6 - sum(disease), 
                             replace = TRUE, 
                             prob = c(0.99, 0.01))

现在我们已经运行了模拟，我们可以通过创建列联表来计算真阳性、假阳性、真阴性和假阴性

contingency_table <- table(disease, test)

contingency_table
#>        test
#> disease      0      1
#>       0 989566   9976
#>       1     38    420

并得到这样的阳性测试患上这种疾病的近似概率:

contingency_table[2, 2] / sum(contingency_table[,2])
#> [1] 0.04040015

以及在阴性测试中患病的概率如下:

contingency_table[2, 1] / sum(contingency_table[,1])
#> [1] 3.83992e-05

您会注意到，由于某些采样概率非常小，采样的概率估计值并不那么准确。您可以模拟更大的示例，但您的计算机可能需要一段时间才能运行它。

^{由 reprex package 于 2021 年 8 月 19 日创建(v2.0.0)}