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
没有在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")