设为“星标”，精彩可以过

前段时分有群友问有莫得bootstrap ROC的代码，其时莫得，然则兑现起来很毛糙，今天给环球先容4种表率。

这几种王人是通用的表率，包括但不限于单纯二分类数据的bootstrap ROC/AUC及信得过区间，模子里面考据/外部考据获取的多样盘算的bootstrap信得过区间（或ROC/AUC）

在演示前，先说一下这个bootstrap ROC/AUC的念念路。当先你要知说念什么是bootstrap，然后你要知说念在R中奈何绘制ROC弧线。

假如是作念1000次bootstrap，那就会得到1000个自助集，在每一个自助集王人进行1次ROC分析并绘制1条ROC弧线，获取1个AUC值，把这1000条ROC弧线画在一齐，等于bootstrap ROC了，通过这1000个AUC就可以盘算AUC的置信区间了。

念念路了了，底下等于找器用兑现。我选拔R。

演示数据使用aSAH数据集，这是一个动脉瘤性蛛网膜下腔出血的数据集，一共113行，7列。其中：

gos6：格拉斯哥量表评分outcome：死心变量gender：性别age：年岁wfns：寰球神经外科医生纠合会公认的神经学量表评分s100b：生物象征物ndka：生物象征物

data("aSAH"，package = "pROC")str(aSAH)

## 'data.frame': 113 obs. of  7 variables:##  $ gos6   : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 5 5 5 5 1 1 4 1 5 4 ...##  $ outcome: Factor w/ 2 levels "Good"，"Poor": 1 1 1 1 2 2 1 2 1 1 ...##  $ gender : Factor w/ 2 levels "Male"，"Female": 2 2 2 2 2 1 1 1 2 2 ...##  $ age    : int  42 37 42 27 42 48 57 41 49 75 ...##  $ wfns   : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 3 2 5 4 1 2 ...##  $ s100b  : num  0.13 0.14 0.1 0.04 0.13 0.1 0.47 0.16 0.18 0.1 ...##  $ ndka   : num  3.01 8.54 8.09 10.42 17.4 ...

fbroc

先先容一个最毛糙的，用fbroc这个包兑现，因为你在必应省略谷歌搜索bootstrap ROC in R，前几个死心中等于这个包。

library(fbroc)

这个包在使用时需要把死心变量变为逻辑型：

outcome1 <- ifelse(aSAH$outcome == "Good"，FALSE，TRUE)

然后1行代码即可兑现，默许是1000次bootstrap：

set.seed(123)result.boot <- boot.roc(aSAH$s100b， outcome1)result.boot

## ## Bootstraped uncached ROC Curve with 41 positive and 72 negative samples. ##  ## The AUC is 0.73.##  ## 1000 bootstrap samples will be calculated. ## The results use up 0 MB of memory.

获取1000次bootstrap AUC的信得过区间，还同期给出了圭臬误：

set.seed(123)perf(result.boot， "auc"， conf.level = 0.95)

## ## ##                 Bootstrapped ROC performance metric## ## Metric: AUC## Bootstrap replicates: 1000## Observed: 0.731## Std. Error: 0.052## 95% confidence interval:## 0.625 0.824

把这1000条ROC弧线画在一齐，就得到bootstrap ROC了：

plot(result.boot)

图片

这个是我现在找到的最毛糙的表率。

tidyverse

后头的表率等于凭证开端说的念念路，一步一步的兑现了。

先说个tidy的表率，借助tidyverse和tidymodels兑现。

library(yardstick)library(rsample)library(tidyverse)

先说下如安在tidymodels中绘制ROC弧线，细目可参考：tidymodels-yardstick：猜想模子性能

在tidymodels中画一条ROC弧线相配毛糙，当先是盘算绘图需要的数据：

roc_data <- roc_curve(aSAH， outcome， s100b，event_level = "second")

roc_data

## # A tibble: 52 × 3##    .threshold specificity sensitivity##         <dbl>       <dbl>       <dbl>##  1    -Inf         0            1    ##  2       0.03      0            1    ##  3       0.04      0            0.976##  4       0.05      0.0694       0.976##  5       0.06      0.111        0.976##  6       0.07      0.139        0.976##  7       0.08      0.222        0.902##  8       0.09      0.306        0.878##  9       0.1       0.389        0.829## 10       0.11      0.486        0.780## # ℹ 42 more rows

然后是绘图：

autoplot(roc_data)

图片

接下来只好使用bootstrap生成1000个自助集就可以很通俗的绘制1000条ROC弧线了。

生成1000个自助集：

set.seed(123)asb <- bootstraps(aSAH， times = 1000)asb

## # Bootstrap sampling ## # A tibble: 1，000 × 2##    splits           id           ##    <list>           <chr>        ##  1 <split [113/44]> Bootstrap0001##  2 <split [113/43]> Bootstrap0002##  3 <split [113/47]> Bootstrap0003##  4 <split [113/41]> Bootstrap0004##  5 <split [113/37]> Bootstrap0005##  6 <split [113/37]> Bootstrap0006##  7 <split [113/39]> Bootstrap0007##  8 <split [113/38]> Bootstrap0008##  9 <split [113/33]> Bootstrap0009## 10 <split [113/42]> Bootstrap0010## # ℹ 990 more rows

界说一个函数，获取自助集：这是tidymodels中的常见操作，可参考：tidymodels数据差异

ff <- function(split){analysis(split)}

底下等于提真金不怕火1000个自助集的数据，通河县和列复印机有限公司对每个自助集进行1次ROC分析， 上海商展进出口有限公司以获取绘图数据：

plot_data <- asb %>%   mutate(boot_data = map(splits，
首页-九康东仓储有限公司 ff)) %>%   unnest(boot_data) %>%   group_by(id) %>%   roc_curve(outcome， s100b，event_level = "second") dim(plot_data)

## [1] 40007     4

head(plot_data)

## # A tibble: 6 × 4## # Groups:   id [1]##   id            .threshold specificity sensitivity##   <chr>              <dbl>       <dbl>       <dbl>## 1 Bootstrap0001    -Inf         0            1    ## 2 Bootstrap0001       0.04      0            1    ## 3 Bootstrap0001       0.05      0.0779       1    ## 4 Bootstrap0001       0.06      0.143        1    ## 5 Bootstrap0001       0.07      0.195        1    ## 6 Bootstrap0001       0.08      0.312        0.944

临了把1000条ROC弧线画在一齐即可：也等于环球需要的bootstrap ROC：

ggplot()+  # 自助集的ROC弧线，共1000条  geom_path(data = plot_data，            mapping=aes(1-specificity， sensitivity，group=id)，color = "grey")+  # 原始数据的ROC弧线  geom_path(data = roc_data， mapping = aes(1-specificity， sensitivity)，            color="blue"， linewidth=1.5)+  theme_bw()

图片

由于咱们照旧进行了1000次ROC分析，那当然就可以得回1000个AUC，是以凭证这1000个AUC，就可以盘算均值、圭臬差、圭臬误、信得过区间。

先获取1000个AUC：

boot_auc <- asb %>%   mutate(boot_data = map(splits， ff)) %>%   unnest(boot_data) %>%   group_by(id) %>%   roc_auc(outcome， s100b，event_level = "second") #boot_aucdim(boot_auc)

## [1] 1000    4

head(boot_auc)

## # A tibble: 6 × 4##   id            .metric .estimator .estimate##   <chr>         <chr>   <chr>          <dbl>## 1 Bootstrap0001 roc_auc binary         0.799## 2 Bootstrap0002 roc_auc binary         0.721## 3 Bootstrap0003 roc_auc binary         0.774## 4 Bootstrap0004 roc_auc binary         0.707## 5 Bootstrap0005 roc_auc binary         0.743## 6 Bootstrap0006 roc_auc binary         0.701

这1000个AUC基本接近正态散播：

ggplot(boot_auc， aes(x=.estimate))+  geom_density()

图片

盘算置信区间，公式如下（数学常识和统计常识，鸠合搜索省略看讲义王人可以）：

家用电视机 102， 102);">信得过区间下限 = 均值 - z * 圭臬误

信得过区间上限 = 均值 + z * 圭臬误

先盘算圭臬误：

sample_mean <- mean(boot_auc$.estimate)sample_mean

## [1] 0.7315554

sample_size <- nrow(boot_auc)standard_d <- sd(boot_auc$.estimate)se <- standard_d/sqrt(sample_size)se

## [1] 0.001544964

盘算置信区间：

conf_low <- sample_mean - 1.96 * seconf_low

## [1] 0.7285273

conf_high <- sample_mean + 1.96 * seconf_high

## [1] 0.7345836

base R

和tidy的表率莫得实质区别，仅仅兑现步地使用base R语法远程。这让我想起了某个番邦网友对R的辩驳：现在许多东说念主不是纠结于用R如故用Python，而是纠结于用base R如故tidy R。base R和tidy R确凿太割裂了。

先进行1次bootstrap（获取样本编号）望望后果：

set.seed(123)bootset <- sample(nrow(aSAH)， size = nrow(aSAH)， replace = T)bootset

##   [1]  31  79  51  14  67  42  50  43 101  14  25  90  91  69  91  57  92   9##  [19]  93  99  72  26   7  42   9  83  36  78  81  43 103  76  15  32 106 109##  [37]   7   9  41  74  23  27  60  53   7  53  27  96  38  89  34  93  69  72##  [55]  76  63  13  82  97  91  25  38  21  79  41  47  90  60  95  16  94   6##  [73] 107  72  86  86  39  31 112  81  50 113  34   4  13  69  25  52  22  89##  [91]  32 110  25  87  35  40 112  30  12  31 110  30  64  99  14  93  96  71## [109]  67  23  79  85  37

然后界说一个函数，获取每次的自助集：

get_bootset <- function(data){  boot_index <- sample(nrow(data)， size = nrow(data)， replace = T)  bootset <- data[boot_index，]  return(bootset)}#set.seed(123)#get_bootset(aSAH)

使用bootstrap获取1000个自助集，通过for轮回兑现：

# 每次死心王人不相同bootsets <- list()for(i in 1:1000){  bootsets[[i]] <- get_bootset(aSAH)}length(bootsets)

## [1] 1000

对每一个自助集进行1次ROC分析，通过for轮回兑现：

library(pROC)rocs <- list()for(i in 1:1000){  rocs[[i]] <- pROC::roc(bootsets[[i]][，"outcome"]， bootsets[[i]][，"s100b"]，                   quiet=T)}

画1000条ROC弧线，如故通过for轮回兑现：

# 提供一个画布plot(roc(aSAH$outcome， aSAH$s100b)，col="blue")

## Setting levels: control = Good， case = Poor

## Setting direction: controls < cases

# 画1000条ROC弧线for(i in 1:1000){  lines.roc(rocs[[i]]，col="grey")}# 画完1000条把正本的挡住了，重新画一条lines.roc(roc(aSAH$outcome， aSAH$s100b)，col="blue")

## Setting levels: control = Good， case = Poor## Setting direction: controls < cases

图片

然后是盘算1000个AUC的置信区间，和tidy的表率相同的。

盘算1000个AUC：

aucs <- list()for(i in 1:1000){  aucs[[i]] <- auc(pROC::roc(bootsets[[i]][，"outcome"]，bootsets[[i]][，"s100b"]，                   quiet=T))}aucs <- unlist(aucs)

盘算信得过区间：

sample_mean <- mean(aucs)sample_mean

## [1] 0.7312995

sample_size <- length(aucs)standard_d <- sd(aucs)se <- standard_d/sqrt(sample_size)se

## [1] 0.001569356

95%的信得过区间，参考讲义省略这个知乎的阐发注解[1]

conf_low <- sample_mean - 1.96 * seconf_low

## [1] 0.7282235

conf_high <- sample_mean + 1.96 * seconf_high

## [1] 0.7343754

这种表率由于我莫得在每次重抽样时设定种子数，导致死心是不可重叠的哈，每次王人不太相同~

boot

boot是挑升作念重抽样的经典R包，在《R话语实战》一书中有详备先容。

通过这个包也可以盘算bootstrap AUC的置信区间，然则这种表率只可盘算盘算，不可画ROC弧线。

library(boot)library(pROC)

界说一个函数，提真金不怕火AUC：

# boot的使用步地很奇怪get_auc <- function(data， ind， outcome， predictor){  d = data[ind，] #这句必须加  au <- as.numeric(auc(pROC::roc(d[，outcome]， d[，predictor]，quiet=T)))  au}get_auc(aSAH， outcome="outcome"，predictor="s100b")

## [1] 0.7313686

提供给boot使用即可：

set.seed(123)ba <- boot(aSAH， get_auc， R = 1000，           outcome="outcome"，predictor="s100b")ba

## ## ORDINARY NONPARAMETRIC BOOTSTRAP## ## ## Call:## boot(data = aSAH， statistic = get_auc， R = 1000， outcome = "outcome"， ##     predictor = "s100b")## ## ## Bootstrap Statistics :##      original       bias    std. error## t1* 0.7313686 0.0001084232  0.05365581

死心给出了原始的AUC，以及1000次bootstrap得到的AUC的圭臬误。

可以对这个死心画个图望望这1000个AUC的散播：

plot(ba)

图片

获取这1000个AUC的置信区间，默许会给出95%的置信区间，并包含4种盘算表率的死心：

boot.ci(ba)

## Warning in boot.ci(ba): bootstrap variances needed for studentized intervals

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS## Based on 1000 bootstrap replicates## ## CALL : ## boot.ci(boot.out = ba)## ## Intervals : ## Level      Normal              Basic         ## 95%   ( 0.6261，  0.8364 )   ( 0.6314，  0.8479 )  ## ## Level     Percentile            BCa          ## 95%   ( 0.6148，  0.8313 )   ( 0.6048，  0.8228 )  ## Calculations and Intervals on Original Scale

4种盘算表率的置信区间王人有了。

OVER！

参考府上[1]

置信区间盘算: https://zhuanlan.zhihu.com/p/259232881家用电视机，

本站仅提供存储处事，整个内容均由用户发布，如发现存害或侵权内容，请点击举报。

• 家用电视机
• sensitivity