如何通过向量计算？

Question

Updated：现在可以了，但还是不知道其他方法是如何工作的。

 cuts <- seq(from=3, to=36, by=0.01)

    for (i in cuts) {
      cut_off<- i
      set.seed(666)
      samp_h <-rnorm(1000,mean=12,sd=3)
      samp_d <-rnorm(1000,mean=18,sd=6)
      a <- sum(samp_h <= cut_off)
      c <- sum(samp_h > cut_off)
      b <- sum(samp_d <= cut_off)
      d <- sum(samp_d > cut_off)
      sens <- a / (a+c)
      spci <- d / (d+b)
      assign(paste("ss",as.character(cut_off),sep = ""), sens)
      assign(paste("sp",as.character(cut_off),sep = ""), spci)}

    ss_v<- unlist(
      lapply(               
        paste0("ss",cuts), 
        get)              
    )

    sp_v<- unlist(
      lapply(               
        paste0("sp",cuts), 
        get)              
    )

    plot(1-sp_v, ss_v)

大家好： 我试图使用不同的“cut_off”来获得不同的“sens”（敏感）和“spci”（特异性）。上面代码的问题是，对于 34 个“剪切”，我可以获得结果。但如果我将剪辑更改为：

cuts <- seq(from=3, to=36, by=0.01)

此方法无法返回结果。问题是我计算每个向量中的数字，所以我问如何使用向量直接计算“ss_v”和“ss_p”。非常感谢。

背景资料： 假设“健康”患者的抗体水平分布正常（12,32），而“患病”患者的抗体水平分布正常（18,62）。请注意，这些是“编造”的数字，并不现实。 模拟大量患病和健康患者的抗体计数（例如，每个患者 1000 人）——使用 R 中的“rnorm”函数。如果选择 15 的截止值，灵敏度和特异性是多少？ 记录 3 到 36 之间的临界值范围（例如 3、3.01、3.02、……、35.98、35.99、36）的灵敏度和特异性。提示：使用 R 中的“seq”函数生成截止值，然后使用“for”循环或矢量化计算计算灵敏度和特异性。 绘制 x 轴为“1-Specificity”，y 轴为“Sensitivity”的图。

Answer 1

您的代码是尝试将 R 用作宏语言的示例。更好的是学习如何正确使用 R 向量。由于您使用了for 循环，因此您应该预先分配sens 和spci，而不是分配给sens 和spci 作为索引向量。 （因此，我赞同您对结果向量的请求，这是明智的做法。）然后给向量命名，而不是在您的工作区中乱扔大量单独的、不连贯的命名对象。试试这个：

cuts <- seq(from=3, to=36, by=1)
sens <- numeric(length(cuts)); spci=numeric(length(cuts))
for (i in cuts) {
  cut_off<- i
  set.seed(666)
  samp_h <-rnorm(1000,mean=12,sd=3)
  samp_d <-rnorm(1000,mean=18,sd=6)
  hth <- table(samp_h)
  dis <-table(samp_d)
  a<-length(hth[names(hth) <= cut_off])
  c<-length(hth[names(hth) > cut_off])
  b <-length(dis[names(dis) <= cut_off])
  d <-length(dis[names(dis) > cut_off])
  sens[i] <- a / (a+c)
  spci[1] <- d / (d+b)
} 
 names(sens) <- paste0("ss",cuts)  
 names(spci) <- paste0("sp",cuts)

我不认为在每次循环迭代中处理一个新的模拟数据集的概念真的给我留下了深刻的印象，但如果你用 diff 模拟某些东西可能会这样。我也不确定您是否正确构建了sens 和spci 作为敏感性和特异性，但至少您现在可以看到结果是什么样的。有几个包可以构建 ROC 曲线。

这就是我怀疑循环内的算法是否正确的原因：

> sens
  ss3   ss4   ss5   ss6   ss7   ss8   ss9  ss10  ss11  ss12  ss13 
0.000 0.000 0.745 0.747 0.752 0.764 0.792 0.836 0.895 0.000 0.123 
 ss14  ss15  ss16  ss17  ss18  ss19  ss20  ss21  ss22  ss23  ss24 
0.239 0.374 0.485 0.593 0.661 0.700 0.721 0.736 0.744 0.745 0.745 
 ss25  ss26  ss27  ss28  ss29  ss30  ss31  ss32  ss33  ss34  ss35 
0.745 0.745 0.745 0.745 0.745 0.745 0.745 0.747 0.747 0.747 0.747 
 ss36  <NA>  <NA> 
0.747 0.747 0.747 

这看起来不像我所期望的敏感结果。我可能使用过abcd <-table( samp_h >= cut_off, samp_d >= cutoff) 之类的代码来生成您对a、b、c、d 的值。然后，您可以对该表结果使用矩阵索引。另一种选择可能是跳过你的表工作并使用这个代码块：

  a <- sum(samp_h <= cut_off)
  c <- sum(samp_h > cut_off)
  b <- sum(samp_d <= cut_off)
  d <- sum(samp_d > cut_off)

sens-itivity 结果现在看起来更合理，但spci 结果却不是这样。`（因为我的索引错误，现在在下面的代码中修复。）

cuts <- seq(from=3, to=36, by=1)
sens <- numeric(length(cuts)); spci=numeric(length(cuts))
  set.seed(666)
  samp_h <-rnorm(1000,mean=12,sd=3)
  samp_d <-rnorm(1000,mean=18,sd=6)
#Only need to make the test data.frame once
 dfrm <- data.frame( vals = c(samp_h, samp_d), 
                     grp = c( rep("H", 1000), rep("D",1000) ) )

for (i in seq_along(cuts) ) {
  cut_off<- i

  abcd <- with(dfrm, 
    table(Test_res = vals > cut_off, 
          status=grp ) )
  sens[i] <- abcd["TRUE","D"] / sum( abcd[, "D"])
  spci[i] <- abcd["FALSE", "H"] / sum( abcd[, "H"])
} 
 names(sens) <- paste0("ss",cuts)  
 names(spci) <- paste0("sp",cuts)

plot(  1-spci, sens, type="b")
text( 1-spci[c(TRUE,FALSE,FALSE,FALSE,FALSE)]+.05, 
      # hack to print every 5th cutoff value
      sens[c(TRUE,FALSE,FALSE,FALSE,FALSE)], 
      label=(3:36)[ c(TRUE,FALSE,FALSE,FALSE,FALSE)] )