如何用向量计算这个值？

Question

​

更新:现在它正在工作，但仍然不知道其他方式是如何工作的。

 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‘(spcificity)。上面代码的问题是，对于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-特异性”的图，在y-轴生成“敏感性”图.

Answer 1

您的代码是尝试使用R作为宏语言的一个例子。最好是学习如何正确地使用R向量。由于您使用的是for-loop，所以应该预先分配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)

我不认为用每次循环迭代来处理一个新的模拟数据集的概念确实会给我留下深刻的印象，但是如果你用不同的方法来模拟某件事情的话，它可能会让我印象深刻。我也不确定您是否正确地构造了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)] )

​

​