В описанных ниже примерах использовались две модели: одна с явным разделяющим порогом, вторая – без него. Эти две модели, конечно, вряд ли описывают реальные объекты, но позволяют приблизительно оценить свойства получаемых при дихотомизации количественного предиктора оценок. Сама по себе дихотомизация предикторов (равно как и исходов) негативно воспринимается многими известными статистиками [Harrell, Senn] и в реальных исследованиях ее нужно избегать.

nsim <- 1000

NB! Есть два существенных замечания по поводу результатов представленных ниже симуляций. 1) я для простоты целенаправленно проводил их так, чтобы избежать значительного дисбаланса классов, поэтому результаты больше имеют отношение к исследованиям типа случай-контроль; 2) использованные методы определения оптимальных порогов не всегда приводят к однозначным резульататам, поэтому при наличии двух возможных порогов (для одного метода) я выбирал тот, который меньше отличается от теоретического порога.

Все симуляции проводились с размером выборки \(N = 100\), количество выборок = 1000.

Модель #1.

В данной модели переменная отклик (\(Y\)) имеет условную функцию рапределения следующего вида:

\[ P(Y=y | X=x) = \begin{cases} 0.8^y(1-0.8)^{1-y}I_{\{x \ge 100\}}\\0.2^y(1-0.2)^{1-y}I_{\{x < 100\}} \end{cases} \]

где \(Y \in \{0, 1\}\) и \(X \sim N(100, 20^{2})\), т.е. \(Y\) принимает значение 1 с вероятностью 0.8, если \(X \ge 100\) (и 0 иначе) и значение 1 с вероятность 0.2 если \(X<100\). Независимая переменная \(X\) имеет нормальное распределение с матожиданием \(\mu = 100\) и среднеквадратичным отклонением \(\sigma =20\).

Несложно провести расчет отношения шансов:

\[OR = \frac{Odds_{(x \ge 100)}}{Odds_{(x < 100)}} = \frac{0.8/(1-0.8)}{0.2/(1-0.2)} = 16\]

Кроме того при испрользовании порога \(X = 100\) в качестве порога, можно оцйенить диагностические свойства.

  1. вероятность истинноотрицательного результата

\[P(Y = 0|X < 100) = 0.2^0(1-0.2)^1 \times 0.5 = 0.4\]

  1. вероятность ложноотрицательного результата

\[P(Y = 0|X \ge 100) = 0.8^0(1-0.8)^1 \times 0.5 = 0.1\]

  1. вероятность истинноположительного результата

\[P(Y = 1|X \ge 100) = 0.8^1(1-0.8)^0 \times 0.5 = 0.4\]

  1. вероятность ложноположительного результата

\[P(Y = 1|X < 100) = 0.2^1(1-0.2)^0 \times 0.5 = 0.1\]

  1. точность (Accuracy)

\[P(Y = 0|X < 100) + P(Y = 1|X \ge 100) = 0.8\]

n <- 100000
x <- rnorm(n, 100, 20)
y <- sapply(x, function(x) ifelse(x >= 100, rbinom(1, 1, 0.8), rbinom(1, 1, 0.2)))
df <- data.frame(X = x, Xge100 = as.numeric(x >= 100), Y = y)

Дополним расчет, симуляцией (классы сбалансированы: \(N_{(Y=1)} \approx N_{(Y=0)}\) при \(N_{(X\ge 100)} + N_{(X<100)} \to \infty\)). Для начала посмотрим на распрелеление и сравним полученные выше результаты с результатами при \(N_{(X\ge 100)} + N_{(X<100)} = 10^{5}\).

plot(jitter(df$Y, factor = 0.5) ~ df$X, yaxt="n", 
     col = rgb(0.15, 0.17, 0.18, alpha = 0.05), 
     pch = 20, cex = 0.2, xlab = 'X', ylab = 'Y')
axis(2, at = c(0, 1), labels = c('0', '1'), las = 2)
abline(v = 100, col = 'red')

Y <- factor(df$Y, labels = c('Y=0', 'Y=1'))
confusion <- (table(df$Xge100, Y)/n)
cbind( c('X<100', 'X≥100'),  confusion) %>% 
    kable(row.names = F, align = 'c', caption = 'Матрица ошибок(Confusion matrix)') %>% 
    kable_styling(full_width = F, bootstrap_options = 'bordered') %>% 
    column_spec(1, bold = T)
Матрица ошибок(Confusion matrix)
Y=0 Y=1
X<100 0.40281 0.09892
X≥100 0.10019 0.39808

Теперь посмотрим на ROC-кривую, применим два подхода к определению порога.

roc_fit <- roc(df$Y, df$X, ci = F, show.thres = T, auc = T)
plot(roc_fit, ylab="TPR", xlab="TNR")

AUC = 0.800482

Два метода определения порогов:

  1. J-статистика Йодена (максимизирует расстояние до дианонали ROC-кривой)

\[max[Чувствительность + Специфичность]\]

  1. Closest TopLeft (минимизирует расстояние до левого верхнего угла ROC-кривой)

\[min[(1-Чувствительность)^2 + (1-Специфичность)^2]\]

rbind(
    coords(roc_fit, "best", best.method="youden", 
       ret=c("threshold", "sens", "spec",'accuracy')),
    coords(roc_fit, "best", best.method="closest.topleft", 
       ret=c("threshold", "sens", "spec",'accuracy'))
) %>% t() %>% as.data.frame() -> thresh
fit1 <- glm(factor(df$Y) ~ factor(as.numeric(df$X>=thresh[1,1])), family = 'binomial')
fit2 <- glm(factor(df$Y) ~ factor(as.numeric(df$X>=thresh[1,2])), family = 'binomial')
colnames(thresh) <- c("Youden's J", "TopLeft")
rbind(
    thresh[1,],
    c(exp(coef(fit1))[2], exp(coef(fit2))[2]),
    thresh[2:4,]) -> thresh
rownames(thresh) <- c('Порог','Отношение шансов','Чувствительность','Специфичность','Точность')
thresh %>%
     kable(row.names = T, digits = 3) %>%
     kable_styling(full_width = F)
Youden’s J TopLeft
Порог 99.997 99.997
Отношение шансов 16.184 16.184
Чувствительность 0.801 0.801
Специфичность 0.801 0.801
Точность 0.801 0.801

Таким образом, оба метода определения порога позволяют сделать “хорошие” оценки интересующих показателей при большой выборке и сбалансированных классах.

n <- 100
result <- matrix(data = NA_real_, nrow = nsim, ncol = 12)
colnames(result) <- c('P_hat','AUC','Yoden','Sens_Yoden','Spec_Yoden','Accuracy_Yoden',
                      'OR_Yoden','TopLeft','Sens_TL','Spec_TL','Accuracy_TL','OR_TL')

for (i in 1:nsim) {
    x <- rnorm(n, 100, 20)
    y <- sapply(x, function(x) ifelse(x >= 100, rbinom(1, 1, 0.8), rbinom(1, 1, 0.2)))
    roc_fit <- roc(y, x, ci = F, show.thres = T, auc = T)
    p_hat <- mean(y)
    AUC <- as.numeric(auc(roc_fit))

    Yoden <- data.frame(coords(roc_fit, "best", best.method="youden", 
                ret=c("threshold", "sens", "spec",'accuracy'), transpose = F))
    TopLeft <- data.frame(coords(roc_fit, "best", best.method="closest.topleft", 
                  ret=c("threshold", "sens", "spec",'accuracy'), transpose = F))
    Yoden <- Yoden[which.min(abs(100 - Yoden[,1])),]
    TopLeft <- TopLeft[which.min(abs(100 - Yoden[,1])),]
    OR_Yoden <- exp(coef(glm(factor(y) ~ factor(as.numeric(x>=Yoden[1,1])), family = 'binomial')))[2]
    OR_TL <- exp(coef(glm(factor(y) ~ factor(as.numeric(x>=TopLeft[1,1])), family = 'binomial')))[2]
    result[i,] <- unlist(c(p_hat, AUC, Yoden[1,], OR_Yoden, TopLeft[1,], OR_TL))
}
res1 <- apply(result, 2, mean_sd)

Теперь проведем симуляцию при размере выборки \(N = 100\) для того, чтобы оценить свойства распределений интересующих нас статистик.

NB! Залитые области на графиках эмпирических функций плотности соответствуют 95% интервалам наибольшей плотности (Highest density interval, HDI – 95% наиболее часто встречающихся (наиболее вероятных) значений). Пунктирными линиями обозначены матожидания соотвествующих статистик, сплошными – наблюдаемые средние по 1000 выборкам.

Рсапределения числа успехов (\(Y = 1\)) и оценок AUC.

par(mfrow = c(1,2), mar = c(5,5,1,1))
plot(density(result[, 'P_hat']), main = NA, xlab = 'Доля Y=1', ylab = 'Плотность')
polygon(density(result[, 'P_hat'])$x[density(result[, 'P_hat'])$y/max(density(result[, 'P_hat'])$y)>=0.025],
        density(result[, 'P_hat'])$y[density(result[, 'P_hat'])$y/max(density(result[, 'P_hat'])$y)>=0.025],
        col = rgb(0.15, 0.15, 0.15, alpha = 0.05), border = NA)
abline(v=0.5, col=rgb(0.15, 0.15, 0.15), lwd = 2)
abline(v=mean(result[, 'P_hat']), col='red', lty = 2, lwd = 2)
plot(density(result[, 'AUC']), main = NA, xlab = 'AUC', ylab = 'Плотность')
polygon(density(result[, 'AUC'])$x[density(result[, 'AUC'])$y/max(density(result[, 'AUC'])$y)>=0.025],
        density(result[, 'AUC'])$y[density(result[, 'AUC'])$y/max(density(result[, 'AUC'])$y)>=0.025],
        col = rgb(0.15, 0.15, 0.15, alpha = 0.05), border = NA)
abline(v=0.8, col=rgb(0.15, 0.15, 0.15), lwd = 2)
abline(v=mean(result[, 'AUC']), col='red', lty = 2, lwd = 2)

Оценки порогов

par(mar = c(5,5,1,1))
plot(result[,'Yoden'], result[,'TopLeft'], col = rgb(0.15, 0.17, 0.18, alpha = 0.2), 
     xlim = c(90,110), ylim = c(90,110), cex = 0.8, xlab = 'Yoden\'s J', ylab = 'TopLeft')
segments(100, 80, 100, 100, lwd = 2)
segments(mean(result[,'Yoden']), 80, mean(result[,'Yoden']), mean(result[,'TopLeft']),  
         col = rgb(1,1,0.5), lty = 2, lwd = 2)
segments(80, 100, 100, 100, lwd = 2)
segments(80, mean(result[,'TopLeft']), mean(result[,'Yoden']), mean(result[,'TopLeft']),
         col = rgb(0,1,0.7), lty = 2, lwd = 2)
rug(result[,'Yoden'], lwd = 0.5, side = 1, col = rgb(1,1,0.5, alpha = 0.2), tck = 0.05, pos = 90)
rug(result[,'TopLeft'], lwd = 0.5, side = 2, col = rgb(0,1,0.7, alpha = 0.2), tck = 0.05, pos = 90)

Оценки характеристик диагностического теста для соответсвующих порогов.

par(mfrow = c(2,2), mar = c(5,5,1,1))

plot(density(result[, 'Sens_TL']), main = NA, xlab = 'Чувствительность', ylab = 'Плотность',
     col = rgb(0,1,0.7))
lines(range(density(result[, 'Sens_TL'])$x[density(result[, 'Sens_TL'])$y/max(density(result[, 'Sens_TL'])$y)>=0.025]), rep(min(density(result[, 'Sens_TL'])$y[density(result[, 'Sens_TL'])$y/max(density(result[, 'Sens_TL'])$y)>=0.025]),2),
      col = rgb(0,1,0.7), lty = 1, lwd = 2)
lines(density(result[, 'Sens_Yoden']), col = rgb(1,1,0.5))
lines(range(density(result[, 'Sens_Yoden'])$x[density(result[, 'Sens_Yoden'])$y/max(density(result[, 'Sens_Yoden'])$y)>=0.025]), rep(min(density(result[, 'Sens_Yoden'])$y[density(result[, 'Sens_Yoden'])$y/max(density(result[, 'Sens_Yoden'])$y)>=0.025]),2),
      col = rgb(1,1,0.5), lty = 2, lwd = 2)
abline(v=mean(result[, 'Sens_TL']), col = rgb(0,1,0.7), lty = 2, lwd = 2)
abline(v=mean(result[, 'Sens_Yoden']), col = rgb(1,1,0.5), lty = 2, lwd = 2)

plot(density(result[, 'Spec_TL']), main = NA, xlab = 'Специфичность', ylab = 'Плотность',
     col = rgb(0,1,0.7))
lines(range(density(result[, 'Spec_TL'])$x[density(result[, 'Spec_TL'])$y/max(density(result[, 'Spec_TL'])$y)>=0.025]), rep(min(density(result[, 'Spec_TL'])$y[density(result[, 'Spec_TL'])$y/max(density(result[, 'Spec_TL'])$y)>=0.025]),2),
      col = rgb(0,1,0.7), lty = 1, lwd = 2)
lines(density(result[, 'Spec_Yoden']), col = rgb(1,1,0.5))
lines(range(density(result[, 'Spec_Yoden'])$x[density(result[, 'Spec_Yoden'])$y/max(density(result[, 'Spec_Yoden'])$y)>=0.025]), rep(min(density(result[, 'Spec_Yoden'])$y[density(result[, 'Spec_Yoden'])$y/max(density(result[, 'Spec_Yoden'])$y)>=0.025]),2),
      col = rgb(1,1,0.5), lty = 2, lwd = 2)
abline(v=mean(result[, 'Spec_TL']), col = rgb(0,1,0.7), lty = 2, lwd = 2)
abline(v=mean(result[, 'Spec_Yoden']), col = rgb(1,1,0.5), lty = 2, lwd = 2)

plot(density(result[, 'Accuracy_TL']), main = NA, xlab = 'Точность', ylab = 'Плотность',
     col = rgb(0,1,0.7))
lines(range(density(result[, 'Accuracy_TL'])$x[density(result[, 'Accuracy_TL'])$y/max(density(result[, 'Accuracy_TL'])$y)>=0.025]), rep(min(density(result[, 'Accuracy_TL'])$y[density(result[, 'Accuracy_TL'])$y/max(density(result[, 'Accuracy_TL'])$y)>=0.025]),2),
      col = rgb(0,1,0.7), lty = 1, lwd = 2)
lines(density(result[, 'Accuracy_Yoden']), col = rgb(1,1,0.5))
lines(range(density(result[, 'Accuracy_Yoden'])$x[density(result[, 'Accuracy_Yoden'])$y/max(density(result[, 'Accuracy_Yoden'])$y)>=0.025]), rep(min(density(result[, 'Accuracy_Yoden'])$y[density(result[, 'Accuracy_Yoden'])$y/max(density(result[, 'Accuracy_Yoden'])$y)>=0.025]),2),
      col = rgb(1,1,0.5), lty = 2, lwd = 2)
abline(v=mean(result[, 'Accuracy_TL']), col = rgb(0,1,0.7), lty = 2, lwd = 2)
abline(v=mean(result[, 'Accuracy_Yoden']), col = rgb(1,1,0.5), lty = 2, lwd = 2)

plot(density(result[, 'OR_TL']), main = NA, xlab = 'Отношение шансов', ylab = 'Плотность',
     col = rgb(0,1,0.7))
lines(range(density(result[, 'OR_TL'])$x[density(result[, 'OR_TL'])$y/max(density(result[, 'OR_TL'])$y)>=0.025]), rep(min(density(result[, 'OR_TL'])$y[density(result[, 'OR_TL'])$y/max(density(result[, 'OR_TL'])$y)>=0.025]),2),
      col = rgb(0,1,0.7), lty = 1, lwd = 2)
lines(density(result[, 'OR_Yoden']), col = rgb(1,1,0.5))
lines(range(density(result[, 'OR_Yoden'])$x[density(result[, 'OR_Yoden'])$y/max(density(result[, 'OR_Yoden'])$y)>=0.025]), rep(min(density(result[, 'OR_Yoden'])$y[density(result[, 'OR_Yoden'])$y/max(density(result[, 'OR_Yoden'])$y)>=0.025]),2),
      col = rgb(1,1,0.5), lty = 2, lwd = 2)
abline(v=mean(result[, 'OR_TL']), col = rgb(0,1,0.7), lty = 2, lwd = 2)
abline(v=mean(result[, 'OR_Yoden']), col = rgb(1,1,0.5), lty = 2, lwd = 2)

Модель #2.

Данная модель несколько более приближена к реальности, в ней отсутствует четкая разделяющая граница между классами, зависящая от значения переменной \(X\). При дихотомизации на уровне \(X = 100\), \(p = 0.5\). Причем классы сбалансированы при данном пороге (\(N_{(Y=1)} \approx N_{(Y=0)}\) при \(N_{(X\ge 100)} + N_{(X<100)} \to \infty\)).

В целом модель имеет следующий вид:

\[log(\frac{p_i}{1-p_i}) = \beta_0 + \beta_1\times x_i + \epsilon_i\]

где \(\beta_1 = log(OR) = \frac{1}{2}\) при изменении количественного предиктора на 1, а \(e_i \sim N(0, 0.1)\).

Проверим данные утверждения на большой выборке (\(N= 100\)).

n <- 100000
proba <- seq(0.1, 0.9, length.out = n)
y <- rbinom(n, 1, proba)
logits <- log(proba/(1-proba))
x <- 100 + 2 * logits + rnorm(n, 0, 0.1)
df <- data.frame(X = x, Xge100 = as.numeric(x >= 100), Y = y)
fit <- glm(Y ~ X, df, family = 'binomial')
df$pred <- predict(fit, type = 'response')

plot(jitter(df$Y, factor = 0.5) ~ df$X, yaxt="n", 
     col = rgb(0.158, 0.168, 0.184, alpha = 0.05), 
     pch = 20, cex = 0.2, xlab = 'X', ylab = 'Y')
axis(2, at = c(0, 0.5, 1), labels = c('0','0.5', '1'), las = 2)
lines(df$X, df$pred, col = 'blue')
segments(100, -0.5, 100, 0.5, col = 'red')
segments(94, 0.5, 100, 0.5, col = 'red')

Соотношение сигнал/шум в модели

plot(df$pred, proba, type = 'l', col = rgb(0.158, 0.168, 0.184, alpha = 0.3), 
     xlab = 'Предсказанная вероятность Y=1', ylab = 'Вероятность Y=1')

summary(fit)$coef %>% 
    as.table() -> coeffs
colnames(coeffs) <- c('Оценка', 'Стандартная\nошибка', 'z-статистика', 'p-значение')
rownames(coeffs) <- c('$\\beta_0$', '$\\beta_1$')
coeffs %>% 
    kable(caption = 'Коэффициенты регрессионного уравнения', escape = F) %>% kable_styling(full_width = F)
Коэффициенты регрессионного уравнения
Оценка Стандартная ошибка z-статистика p-значение
\(\beta_0\) -50.0474400 0.3754004 -133.3175 0
\(\beta_1\) 0.5004046 0.0037531 133.3304 0 
roc_fit <- roc(y, x, ci = F, show.thres = T)
plot(roc_fit, ylab="TPR", xlab="TNR")

auc_num <- as.numeric(auc(roc_fit))

AUC = 0.7669859

rbind(
    coords(roc_fit, "best", best.method="youden", 
       ret=c("threshold", "sens", "spec",'accuracy')),
    coords(roc_fit, "best", best.method="closest.topleft", 
       ret=c("threshold", "sens", "spec",'accuracy'))
) %>% t() %>% as.data.frame() -> thresh
fit1 <- glm(factor(df$Y) ~ factor(as.numeric(df$X>=thresh[1,1])), family = 'binomial')
fit2 <- glm(factor(df$Y) ~ factor(as.numeric(df$X>=thresh[1,2])), family = 'binomial')
colnames(thresh) <- c("Youden's J", "TopLeft")
rbind(
    thresh[1,],
    c(exp(coef(fit1))[2], exp(coef(fit2))[2]),
    thresh[2:4,]) -> thresh
rownames(thresh) <- c('Порог','Отношение шансов','Чувствительность','Специфичность','Точность')
thresh %>%
     kable(row.names = T, digits = 3) %>%
     kable_styling(full_width = F)
Youden’s J TopLeft
Порог 99.934 100.012
Отношение шансов 5.493 5.483
Чувствительность 0.711 0.699
Специфичность 0.690 0.702
Точность 0.701 0.701 
n <- 100
result <- matrix(data = NA_real_, nrow = nsim, ncol = 13)
colnames(result) <- c('P_hat','AUC','OR','Yoden','Sens_Yoden','Spec_Yoden','Accuracy_Yoden',
                      'OR_Yoden','TopLeft','Sens_TL','Spec_TL','Accuracy_TL','OR_TL')

for (i in 1:nsim) {
    proba <- seq(0.1, 0.9, length.out = n)
    y <- rbinom(n, 1, proba)
    logits <- log(proba/(1-proba))
    x <- 100 + 2 * logits + rnorm(n, 0, 0.1)
    roc_fit <- roc(y, x, ci = F, show.thres = T, auc = T)
    p_hat <- mean(y)
    AUC <- as.numeric(auc(roc_fit))
    OR <- exp(coef(glm(factor(y) ~ x, family = 'binomial')))[2]

    Yoden <- data.frame(coords(roc_fit, "best", best.method="youden", 
                ret=c("threshold", "sens", "spec",'accuracy'), transpose = F))
    TopLeft <- data.frame(coords(roc_fit, "best", best.method="closest.topleft", 
                  ret=c("threshold", "sens", "spec",'accuracy'), transpose = F))
    Yoden <- Yoden[which.min(abs(100 - Yoden[,1])),]
    TopLeft <- TopLeft[which.min(abs(100 - Yoden[,1])),]
    OR_Yoden <- exp(coef(glm(factor(y) ~ factor(as.numeric(x>=Yoden[1,1])), family = 'binomial')))[2]
    OR_TL <- exp(coef(glm(factor(y) ~ factor(as.numeric(x>=TopLeft[1,1])), family = 'binomial')))[2]
    result[i,] <- unlist(c(p_hat, AUC, OR, Yoden[1,], OR_Yoden, TopLeft[1,], OR_TL))
}
res2 <- apply(result, 2, mean_sd)[-3]

Распределения числа успехов (\(Y = 1\)), оценок OR для количественного предиктора и оценок AUC.

par(mfrow = c(1,3), mar = c(5,5,1,1))
plot(density(result[, 'P_hat']), main = NA, xlab = 'Доля Y=1', ylab = 'Плотность')
polygon(density(result[, 'P_hat'])$x[density(result[, 'P_hat'])$y/max(density(result[, 'P_hat'])$y)>=0.025],
        density(result[, 'P_hat'])$y[density(result[, 'P_hat'])$y/max(density(result[, 'P_hat'])$y)>=0.025],
        col = rgb(0.15, 0.15, 0.15, alpha = 0.05), border = NA)
abline(v=0.5, col=rgb(0.15, 0.15, 0.15), lwd = 2)
abline(v=mean(result[, 'P_hat']), col='red', lty = 2, lwd = 2)
plot(density(result[, 'OR']), main = NA, xlab = 'OR', ylab = 'Плотность')
polygon(density(result[, 'OR'])$x[density(result[, 'OR'])$y/max(density(result[, 'OR'])$y)>=0.025],
        density(result[, 'OR'])$y[density(result[, 'OR'])$y/max(density(result[, 'OR'])$y)>=0.025],
        col = rgb(0.15, 0.15, 0.15, alpha = 0.05), border = NA)
abline(v=exp(1/2), col=rgb(0.15, 0.15, 0.15), lwd = 2)
abline(v=mean(result[, 'OR']), col='red', lty = 2, lwd = 2)
plot(density(result[, 'AUC']), main = NA, xlab = 'AUC', ylab = 'Плотность')
polygon(density(result[, 'AUC'])$x[density(result[, 'AUC'])$y/max(density(result[, 'AUC'])$y)>=0.025],
        density(result[, 'AUC'])$y[density(result[, 'AUC'])$y/max(density(result[, 'AUC'])$y)>=0.025],
        col = rgb(0.15, 0.15, 0.15, alpha = 0.05), border = NA)
abline(v=auc_num, col=rgb(0.15, 0.15, 0.15), lwd = 2)
abline(v=mean(result[, 'AUC']), col='red', lty = 2, lwd = 2)

Оценки порогов

par(mfrow = c(1,2), mar = c(5,5,1,1))
plot(result[,'Yoden'], result[,'TopLeft'], col = rgb(0.15, 0.17, 0.18, alpha = 0.2), 
     xlim = c(97,103), ylim = c(97,103), cex = 0.8, xlab = 'Yoden\'s J', ylab = 'TopLeft')
segments(100, 80, 100, 100, lwd = 2)
segments(mean(result[,'Yoden']), 80, mean(result[,'Yoden']), mean(result[,'TopLeft']),  
         col = rgb(1,1,0.5), lty = 2, lwd = 2)
segments(80, 100, 100, 100, lwd = 2)
segments(80, mean(result[,'TopLeft']), mean(result[,'Yoden']), mean(result[,'TopLeft']),
         col = rgb(0,1,0.7), lty = 2, lwd = 2)
rug(result[,'Yoden'], lwd = 0.5, side = 1, col = rgb(1,1,0.5, alpha = 0.2), tck = 0.05, pos = 97)
rug(result[,'TopLeft'], lwd = 0.5, side = 2, col = rgb(0,1,0.7, alpha = 0.2), tck = 0.05, pos = 97)
plot(result[, 'P_hat']-0.0025, result[,'Yoden'], col = rgb(1,1,0.5, alpha = 0.5), 
     pch = 20, cex = 0.5, xlab = 'Доля Y=1', ylab = 'Оценки порогов')
points(result[, 'P_hat']+0.0025, result[,'TopLeft'], col = rgb(0,1,0.7, alpha = 0.5), 
       pch = 20, cex = 0.5)

Обобщенные характеристики диагностического теста для соответсвующих порогов.

par(mfrow = c(2,2), mar = c(5,5,1,1))

plot(density(result[, 'Sens_TL']), main = NA, xlab = 'Чувствительность', ylab = 'Плотность',
     col = rgb(0,1,0.7))
lines(range(density(result[, 'Sens_TL'])$x[density(result[, 'Sens_TL'])$y/max(density(result[, 'Sens_TL'])$y)>=0.025]), rep(min(density(result[, 'Sens_TL'])$y[density(result[, 'Sens_TL'])$y/max(density(result[, 'Sens_TL'])$y)>=0.025]),2),
      col = rgb(0,1,0.7), lty = 1, lwd = 2)
lines(density(result[, 'Sens_Yoden']), col = rgb(1,1,0.5))
lines(range(density(result[, 'Sens_Yoden'])$x[density(result[, 'Sens_Yoden'])$y/max(density(result[, 'Sens_Yoden'])$y)>=0.025]), rep(min(density(result[, 'Sens_Yoden'])$y[density(result[, 'Sens_Yoden'])$y/max(density(result[, 'Sens_Yoden'])$y)>=0.025]),2),
      col = rgb(1,1,0.5), lty = 2, lwd = 2)
abline(v=mean(result[, 'Sens_TL']), col = rgb(0,1,0.7), lty = 2, lwd = 2)
abline(v=mean(result[, 'Sens_Yoden']), col = rgb(1,1,0.5), lty = 2, lwd = 2)

plot(density(result[, 'Spec_TL']), main = NA, xlab = 'Специфичность', ylab = 'Плотность',
     col = rgb(0,1,0.7))
lines(range(density(result[, 'Spec_TL'])$x[density(result[, 'Spec_TL'])$y/max(density(result[, 'Spec_TL'])$y)>=0.025]), rep(min(density(result[, 'Spec_TL'])$y[density(result[, 'Spec_TL'])$y/max(density(result[, 'Spec_TL'])$y)>=0.025]),2),
      col = rgb(0,1,0.7), lty = 1, lwd = 2)
lines(density(result[, 'Spec_Yoden']), col = rgb(1,1,0.5))
lines(range(density(result[, 'Spec_Yoden'])$x[density(result[, 'Spec_Yoden'])$y/max(density(result[, 'Spec_Yoden'])$y)>=0.025]), rep(min(density(result[, 'Spec_Yoden'])$y[density(result[, 'Spec_Yoden'])$y/max(density(result[, 'Spec_Yoden'])$y)>=0.025]),2),
      col = rgb(1,1,0.5), lty = 2, lwd = 2)
abline(v=mean(result[, 'Spec_TL']), col = rgb(0,1,0.7), lty = 2, lwd = 2)
abline(v=mean(result[, 'Spec_Yoden']), col = rgb(1,1,0.5), lty = 2, lwd = 2)

plot(density(result[, 'Accuracy_TL']), main = NA, xlab = 'Точность', ylab = 'Плотность',
     col = rgb(0,1,0.7))
lines(range(density(result[, 'Accuracy_TL'])$x[density(result[, 'Accuracy_TL'])$y/max(density(result[, 'Accuracy_TL'])$y)>=0.025]), rep(min(density(result[, 'Accuracy_TL'])$y[density(result[, 'Accuracy_TL'])$y/max(density(result[, 'Accuracy_TL'])$y)>=0.025]),2),
      col = rgb(0,1,0.7), lty = 1, lwd = 2)
lines(density(result[, 'Accuracy_Yoden']), col = rgb(1,1,0.5))
lines(range(density(result[, 'Accuracy_Yoden'])$x[density(result[, 'Accuracy_Yoden'])$y/max(density(result[, 'Accuracy_Yoden'])$y)>=0.025]), rep(min(density(result[, 'Accuracy_Yoden'])$y[density(result[, 'Accuracy_Yoden'])$y/max(density(result[, 'Accuracy_Yoden'])$y)>=0.025]),2),
      col = rgb(1,1,0.5), lty = 2, lwd = 2)
abline(v=mean(result[, 'Accuracy_TL']), col = rgb(0,1,0.7), lty = 2, lwd = 2)
abline(v=mean(result[, 'Accuracy_Yoden']), col = rgb(1,1,0.5), lty = 2, lwd = 2)

plot(density(result[, 'OR_TL']), main = NA, xlab = 'Отношение шансов', ylab = 'Плотность',
     col = rgb(0,1,0.7))
lines(range(density(result[, 'OR_TL'])$x[density(result[, 'OR_TL'])$y/max(density(result[, 'OR_TL'])$y)>=0.025]), rep(min(density(result[, 'OR_TL'])$y[density(result[, 'OR_TL'])$y/max(density(result[, 'OR_TL'])$y)>=0.025]),2),
      col = rgb(0,1,0.7), lty = 1, lwd = 2)
lines(density(result[, 'OR_Yoden']), col = rgb(1,1,0.5))
lines(range(density(result[, 'OR_Yoden'])$x[density(result[, 'OR_Yoden'])$y/max(density(result[, 'OR_Yoden'])$y)>=0.025]), rep(min(density(result[, 'OR_Yoden'])$y[density(result[, 'OR_Yoden'])$y/max(density(result[, 'OR_Yoden'])$y)>=0.025]),2),
      col = rgb(1,1,0.5), lty = 2, lwd = 2)
abline(v=mean(result[, 'OR_TL']), col = rgb(0,1,0.7), lty = 2, lwd = 2)
abline(v=mean(result[, 'OR_Yoden']), col = rgb(1,1,0.5), lty = 2, lwd = 2)

Сравнение интересующих характеристик для двух моделей (среднее ± ст.отклонение).

df <- as.data.frame(cbind(res1, res2))
row.names(df) <- c('Доля Y=1','AUC','Порог (Yoden)','Чувствительность (Yoden)','Специфичность (Yoden)','Точность (Yoden)', 'Отношение шансов (Yoden)', 'Порог (TopLeft)','Чувствительнсть (TopLeft)','Специфичность (TopLeft)','Точность (TopLeft)','Отношение шансов (TopLeft)')
names(df) <- c('Модель #1','Модель #2')
df %>% kable(row.names = T) %>%
     kable_styling(bootstrap_options = 'hover',full_width = F)
Модель #1 Модель #2
Доля Y=1 0.503 (0.05) 0.5 (0.044)
AUC 0.8 (0.046) 0.772 (0.045)
Порог (Yoden) 100.04 (1.405) 99.997 (0.73)
Чувствительность (Yoden) 0.805 (0.059) 0.736 (0.114)
Специфичность (Yoden) 0.808 (0.059) 0.735 (0.115)
Точность (Yoden) 0.806 (0.039) 0.735 (0.039)
Отношение шансов (Yoden) 22.941 (19.247) 10.649 (6.586)
Порог (TopLeft) 100.03 (1.187) 99.998 (0.39)
Чувствительнсть (TopLeft) 0.804 (0.054) 0.731 (0.068)
Специфичность (TopLeft) 0.808 (0.054) 0.73 (0.069)
Точность (TopLeft) 0.806 (0.039) 0.73 (0.04)
Отношение шансов (TopLeft) 22.258 (18.72) 8.614 (4.245)