使用pRoc软件包进行ROC分析后如何获得p值?

问题描述 投票:1回答:1

ROC分析一组数据后,如何计算p值?使用相同的统计信息,我看到可以在SPSS中输出p值。示例代码如下:

library(pROC)
data(aSAH)
head(aSAH)
roc(aSAH$outcome, aSAH$s100b, plot=T)

> head(aSAH)
   gos6 outcome gender age wfns s100b  ndka
29    5    Good Female  42    1  0.13  3.01
30    5    Good Female  37    1  0.14  8.54
31    5    Good Female  42    1  0.10  8.09
32    5    Good Female  27    1  0.04 10.42
33    1    Poor Female  42    3  0.13 17.40
34    1    Poor   Male  48    2  0.10 12.75
> roc(aSAH$outcome, aSAH$s100b,plot = T)
Setting levels: control = Good, case = Poor
Setting direction: controls < cases

Call:
roc.default(response = aSAH$outcome, predictor = aSAH$s100b,     plot = T)

Data: aSAH$s100b in 72 controls (aSAH$outcome Good) < 41 cases (aSAH$outcome Poor).
Area under the curve: 0.7314
r spss roc
1个回答
0
投票

没有在roc中获得p值的选项,您可以设置选项ci=TRUE以获得置信区间。抓取roc的不可见输出,然后使用str进行检查。

library(pROC)
rr <- roc(aSAH$outcome, aSAH$s100b, plot=T)
str(rr, 1)
# List of 16
#  $ percent           : logi FALSE
#  $ sensitivities     : num [1:51] 1 0.976 0.976 0.976 0.976 ...
#  $ specificities     : num [1:51] 0 0 0.0694 0.1111 0.1389 ...
#  $ thresholds        : num [1:51] -Inf 0.035 0.045 0.055 0.065 ...
#  $ direction         : chr "<"
#  $ cases             : num [1:41] 0.13 0.1 0.16 0.12 0.44 0.71 0.49 0.07 0.33 0.09 ...
#  $ controls          : num [1:72] 0.13 0.14 0.1 0.04 0.47 0.18 0.1 0.1 0.04 0.08 ...
#  $ fun.sesp          :function (thresholds, controls, cases, direction)  
#  $ auc               : 'auc' num 0.731
#   ..- attr(*, "partial.auc")= logi FALSE
#   ..- attr(*, "percent")= logi FALSE
#   ..- attr(*, "roc")=List of 15
#   .. ..- attr(*, "class")= chr "roc"
#  $ call              : language roc.default(response = aSAH$outcome, predictor = aSAH$s100b, ci = TRUE, plot = TRUE)
#  $ original.predictor: num [1:113] 0.13 0.14 0.1 0.04 0.13 0.1 0.47 0.16 0.18 0.1 ...
#  $ original.response : Factor w/ 2 levels "Good","Poor": 1 1 1 1 2 2 1 2 1 1 ...
#  $ predictor         : num [1:113] 0.13 0.14 0.1 0.04 0.13 0.1 0.47 0.16 0.18 0.1 ...
#  $ response          : Factor w/ 2 levels "Good","Poor": 1 1 1 1 2 2 1 2 1 1 ...
#  $ levels            : chr [1:2] "Good" "Poor"
#  $ ci                : 'ci.auc' num [1:3] 0.63 0.731 0.833
#   ..- attr(*, "conf.level")= num 0.95
#   ..- attr(*, "method")= chr "delong"
#   ..- attr(*, "boot.n")= num 2000
#   ..- attr(*, "boot.stratified")= logi TRUE
#   ..- attr(*, "auc")= 'auc' num 0.731
#   .. ..- attr(*, "partial.auc")= logi FALSE
#   .. ..- attr(*, "percent")= logi FALSE
#   .. ..- attr(*, "roc")=List of 15
#   .. .. ..- attr(*, "class")= chr "roc"
#  - attr(*, "class")= chr "roc"

str演示了如何访问ci

rr$ci
# 95% CI: 0.6301-0.8326 (DeLong)

所以至少您有一个置信区间。

但是,要获得roc分析的p值,您可能需要阅读Cross Validated: https://stats.stackexchange.com/a/386506/163114

上的有趣答案。

使用的数据:

aSAH <- structure(list(gos6 = structure(c(5L, 5L, 5L, 5L, 1L, 1L, 4L, 
1L, 5L, 4L, 1L, 5L, 3L, 1L, 5L, 5L, 1L, 5L, 3L, 1L, 3L, 5L, 3L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 5L, 5L, 1L, 3L, 5L, 3L, 
3L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 3L, 3L, 5L, 5L, 4L, 1L, 
5L, 5L, 5L, 1L, 5L, 1L, 5L, 4L, 1L, 1L, 1L, 5L, 3L, 1L, 5L, 1L, 
3L, 1L, 1L, 5L, 5L, 5L, 1L, 3L, 5L, 1L, 5L, 5L, 5L, 5L, 5L, 5L, 
1L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 4L, 5L, 1L, 1L, 1L, 
3L, 1L, 1L, 5L, 5L, 4L, 1L, 5L, 5L, 5L), .Label = c("1", "2", 
"3", "4", "5"), class = c("ordered", "factor")), outcome = structure(c(1L, 
1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 
1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
2L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 
2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L
), .Label = c("Good", "Poor"), class = "factor"), gender = structure(c(2L, 
2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 
2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 
2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 
2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L
), .Label = c("Male", "Female"), class = "factor"), age = c(42L, 
37L, 42L, 27L, 42L, 48L, 57L, 41L, 49L, 75L, 37L, 49L, 49L, 71L, 
64L, 49L, 63L, 33L, 51L, 54L, 78L, 47L, 56L, 76L, 58L, 50L, 36L, 
43L, 61L, 31L, 58L, 61L, 76L, 54L, 55L, 57L, 64L, 67L, 71L, 74L, 
72L, 62L, 65L, 50L, 51L, 50L, 35L, 44L, 45L, 48L, 61L, 36L, 53L, 
66L, 73L, 58L, 31L, 28L, 32L, 36L, 47L, 29L, 46L, 73L, 45L, 40L, 
50L, 57L, 37L, 60L, 53L, 57L, 56L, 31L, 59L, 66L, 50L, 54L, 55L, 
51L, 44L, 45L, 47L, 57L, 29L, 52L, 18L, 68L, 47L, 45L, 34L, 60L, 
33L, 38L, 48L, 67L, 64L, 70L, 24L, 74L, 42L, 63L, 81L, 62L, 33L, 
51L, 29L, 68L, 53L, 58L, 32L, 39L, 34L), wfns = structure(c(1L, 
1L, 1L, 1L, 3L, 2L, 5L, 4L, 1L, 2L, 5L, 2L, 5L, 5L, 1L, 2L, 5L, 
2L, 2L, 5L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 3L, 4L, 2L, 
1L, 4L, 5L, 4L, 1L, 4L, 2L, 5L, 5L, 1L, 1L, 2L, 4L, 1L, 1L, 2L, 
5L, 5L, 2L, 1L, 2L, 5L, 1L, 2L, 2L, 2L, 1L, 5L, 1L, 1L, 2L, 5L, 
2L, 1L, 4L, 5L, 2L, 4L, 4L, 1L, 4L, 1L, 2L, 1L, 4L, 2L, 1L, 2L, 
3L, 2L, 4L, 2L, 1L, 4L, 5L, 1L, 1L, 5L, 4L, 2L, 2L, 1L, 1L, 1L, 
3L, 1L, 2L, 5L, 2L, 2L, 5L, 5L, 5L, 2L, 4L, 4L, 5L, 1L, 1L, 1L
), .Label = c("1", "2", "3", "4", "5"), class = c("ordered", 
"factor")), s100b = c(0.13, 0.14, 0.1, 0.04, 0.13, 0.1, 0.47, 
0.16, 0.18, 0.1, 0.12, 0.1, 0.44, 0.71, 0.04, 0.08, 0.49, 0.04, 
0.07, 0.33, 0.09, 0.09, 0.07, 0.11, 0.07, 0.17, 0.07, 0.11, 0.13, 
0.19, 0.05, 0.16, 0.41, 0.14, 0.34, 0.35, 0.48, 0.09, 0.96, 0.25, 
0.5, 0.46, 0.16, 0.07, 0.43, 0.45, 0.11, 0.08, 0.09, 0.86, 0.52, 
0.08, 0.06, 0.13, 2.07, 0.1, 0.14, 0.15, 0.07, 0.06, 0.77, 0.05, 
0.09, 0.3, 0.03, 0.09, 0.04, 0.23, 0.7, 0.09, 0.27, 0.71, 0.08, 
0.26, 0.08, 0.16, 0.09, 0.13, 0.1, 0.08, 0.11, 0.33, 0.11, 0.28, 
0.07, 0.1, 0.32, 0.22, 0.07, 0.05, 0.24, 0.38, 0.1, 0.15, 0.08, 
0.14, 0.1, 0.07, 0.04, 0.19, 0.56, 0.14, 0.58, 0.32, 0.82, 0.74, 
0.15, 0.47, 0.17, 0.44, 0.15, 0.5, 0.48), ndka = c(3.01, 8.54, 
8.09, 10.42, 17.4, 12.75, 6, 13.2, 15.54, 6.01, 15.96, 17.86, 
5.18, 8.9, 13.41, 20.75, 11.6, 16.11, 32.37, 54.82, 32.41, 49.94, 
40.34, 9.47, 6.29, 12.53, 6.54, 6.3, 80.3, 12.8, 9.8, 9.81, 9.85, 
18.21, 5.03, 14.04, 21.93, 8.02, 7.42, 8.38, 6.59, 9.63, 13.12, 
7.96, 14.34, 41.43, 7.63, 7.06, 12.59, 13.56, 3.87, 9.44, 7.66, 
12.98, 419.19, 27.19, 22.27, 9.95, 21.22, 11.73, 10.4, 58.83, 
6.39, 11.09, 12.22, 13.67, 17.21, 22.63, 12.9, 9.7, 21.57, 8.23, 
72.57, 15.54, 10.51, 10.6, 14.57, 5.19, 9.63, 4.61, 21.48, 17.3, 
12.71, 9.44, 11.07, 19.46, 10.83, 5.37, 11.97, 7.75, 9.83, 10.55, 
28.49, 8.53, 12.57, 12.9, 46.83, 11.68, 12.67, 9.01, 9.57, 34.06, 
11.72, 14.26, 47.61, 11.67, 24.58, 10.33, 13.87, 15.89, 22.43, 
6.79, 13.45)), row.names = 29:141, class = "data.frame")
© www.soinside.com 2019 - 2024. All rights reserved.