#data load ##########
LATElong <- read_sav("Data/LATElong.sav")
hip<-LATElong %>% filter(region_type==1)
We use a bootstrapping method with 1,000 samples. To determine the optimal cutoff, we maximize the Youden index, which combines sensitivity and specificity.
set.seed(5620)
opt_cut_h <- cutpointr(data=hip, MRI, LATE_status,
direction = ">=", pos_class = 1,
neg_class = 0,
method = maximize_metric,
metric = youden,
boot_runs = 1000)
## Running bootstrap...
summary(opt_cut_h)
## Method: maximize_metric
## Predictor: MRI
## Outcome: LATE_status
## Direction: >=
## Nr. of bootstraps: 1000
##
## AUC n n_pos n_neg
## 0.2474 51 17 34
##
## optimal_cutpoint youden acc sensitivity specificity tp fn fp tn
## 1854.1 0.0294 0.6667 0.0588 0.9706 1 16 1 33
##
## Predictor summary:
## Data Min. 5% 1st Qu. Median Mean 3rd Qu. 95% Max.
## Overall 544.1 735.20 1030.800 1315.30 1280.324 1509.450 1801.60 1859.5
## 0 840.3 901.93 1196.925 1407.25 1379.791 1579.525 1789.75 1859.5
## 1 544.1 637.46 771.400 1102.00 1081.388 1316.900 1536.02 1854.1
## SD NAs
## 335.1179 0
## 289.6575 0
## 338.8659 0
##
## Bootstrap summary:
## Variable Min. 5% 1st Qu. Median Mean 3rd Qu. 95% Max. SD
## optimal_cutpoint 918.30 1315.30 1854.10 1854.10 Inf Inf Inf Inf NaN
## AUC_b 0.03 0.13 0.19 0.24 0.24 0.29 0.38 0.57 0.08
## AUC_oob 0.00 0.08 0.17 0.25 0.25 0.32 0.43 0.71 0.11
## youden_b -0.32 0.00 0.00 0.03 0.04 0.07 0.16 0.33 0.06
## youden_oob -1.00 -0.38 -0.07 0.00 -0.05 0.00 0.00 0.00 0.14
## acc_b 0.31 0.51 0.63 0.67 0.66 0.71 0.78 0.88 0.08
## acc_oob 0.00 0.26 0.59 0.65 0.63 0.71 0.80 0.94 0.14
## sensitivity_b 0.00 0.00 0.00 0.06 0.10 0.11 0.57 1.00 0.18
## sensitivity_oob 0.00 0.00 0.00 0.00 0.02 0.00 0.11 0.83 0.09
## specificity_b 0.09 0.50 0.97 1.00 0.94 1.00 1.00 1.00 0.16
## specificity_oob 0.00 0.31 0.93 1.00 0.93 1.00 1.00 1.00 0.21
## cohens_kappa_b -0.32 0.00 0.00 0.03 0.05 0.09 0.19 0.43 0.07
## cohens_kappa_oob -0.82 -0.30 -0.09 0.00 -0.05 0.00 0.00 0.00 0.11
## NAs
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
plot(opt_cut_h)
## Warning: Removed 389 rows containing non-finite outside the scale range
## (`stat_density()`).
The identified cutoff point is 1854.1, with very low sensitivity (5.88%) and very high specificity (97.06%).
same method
amig<-LATElong %>% filter(region_type==7)
opt_cut_a <- cutpointr(data=amig, MRI, LATE_status,
direction = ">=", pos_class = 1,
neg_class = 0,
method = maximize_metric,
metric = youden,
boot_runs = 1000)
## Multiple optimal cutpoints found, applying break_ties.
## Running bootstrap...
summary(opt_cut_a)
## Method: maximize_metric
## Predictor: MRI
## Outcome: LATE_status
## Direction: >=
## Nr. of bootstraps: 1000
##
## AUC n n_pos n_neg
## 0.3028 51 17 34
##
## optimal_cutpoint youden acc sensitivity specificity tp fn fp tn
## 4296.3 -0.0294 0.6275 0.0588 0.9118 1 16 3 31
##
## Predictor summary:
## Data Min. 5% 1st Qu. Median Mean 3rd Qu. 95% Max.
## Overall 2069.8 2479.50 3321.35 3516.7 3526.984 3906.200 4389.55 4623.1
## 0 2429.1 3125.57 3396.75 3692.9 3684.144 3921.575 4444.66 4623.1
## 1 2069.8 2253.48 2610.50 3199.8 3212.665 3833.600 4238.64 4392.4
## SD NAs
## 579.2936 0
## 446.4308 0
## 692.9475 0
##
## Bootstrap summary:
## Variable Min. 5% 1st Qu. Median Mean 3rd Qu. 95% Max. SD
## optimal_cutpoint 2529.90 2990.31 3926.00 4200.20 Inf Inf Inf Inf NaN
## AUC_b 0.07 0.17 0.24 0.30 0.30 0.36 0.46 0.59 0.09
## AUC_oob 0.00 0.08 0.22 0.30 0.30 0.38 0.52 0.87 0.13
## youden_b -0.41 0.00 0.00 0.05 0.06 0.10 0.22 0.37 0.08
## youden_oob -0.85 -0.47 -0.21 -0.08 -0.13 0.00 0.00 0.18 0.16
## acc_b 0.24 0.47 0.59 0.65 0.64 0.71 0.75 0.84 0.09
## acc_oob 0.09 0.26 0.47 0.58 0.55 0.67 0.76 0.90 0.15
## sensitivity_b 0.00 0.00 0.00 0.13 0.21 0.29 0.92 1.00 0.25
## sensitivity_oob 0.00 0.00 0.00 0.00 0.08 0.11 0.50 0.86 0.17
## specificity_b 0.03 0.12 0.82 0.94 0.85 1.00 1.00 1.00 0.23
## specificity_oob 0.00 0.00 0.69 0.90 0.79 1.00 1.00 1.00 0.27
## cohens_kappa_b -0.33 0.00 0.00 0.05 0.07 0.11 0.23 0.39 0.08
## cohens_kappa_oob -0.75 -0.42 -0.21 -0.10 -0.13 0.00 0.00 0.22 0.15
## NAs
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
## 0
plot(opt_cut_a)
## Warning: Removed 274 rows containing non-finite outside the scale range
## (`stat_density()`).
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:cutpointr':
##
## auc, roc
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
hip$LATE_status <- as.factor(hip$LATE_status)
amig$LATE_status <- as.factor(amig$LATE_status)
# Calculate AUC for each model
roc_MTAscore <- roc( hip$LATE_status,hip$MTAscore)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc_hipvol <- roc(hip$LATE_status,hip$MRI )
## Setting levels: control = 0, case = 1
## Setting direction: controls > cases
roc_amigvol<- roc(amig$LATE_status,amig$MRI )
## Setting levels: control = 0, case = 1
## Setting direction: controls > cases
auc_MTA <- auc(roc_MTAscore)
auc_hip <- auc(roc_hipvol)
auc_amig <- auc(roc_amigvol)
ci_MTA <- ci.auc(roc_MTAscore)
ci_hip <- ci.auc(roc_hipvol)
ci_amig <- ci.auc(roc_amigvol)
# Print results
cat("AUC for MTA score:", auc_MTA, " (95% CI:", ci_MTA[1], "-", ci_MTA[3], ")\n")
## AUC for MTA score: 0.7147177 (95% CI: 0.5499649 - 0.8794706 )
cat("AUC for hippocampal volume:", auc_hip, " (95% CI:", ci_hip[1], "-", ci_hip[3], ")\n")
## AUC for hippocampal volume: 0.7525952 (95% CI: 0.6041152 - 0.9010751 )
cat("AUC for amygdala volume:", auc_amig, " (95% CI:", ci_amig[1], "-", ci_amig[3], ")\n")
## AUC for amygdala volume: 0.6972318 (95% CI: 0.5201028 - 0.8743609 )
Although the confidence intervals provide a clear indication, I conducted DeLong’s test to compare the AUCs and assess whether they differ. The results indicate that they do not.
# Compare AUCs using DeLong's test
test_1_vs_2 <- roc.test(roc_MTAscore, roc_hipvol)
test_1_vs_3 <- roc.test(roc_MTAscore, roc_amigvol)
test_2_vs_3 <- roc.test(roc_hipvol, roc_amigvol)
# Print test results
cat("Comparison of MTA score vs hipocampal volume:\n")
## Comparison of MTA score vs hipocampal volume:
print(test_1_vs_2)
##
## Bootstrap test for two correlated ROC curves
##
## data: roc_MTAscore and roc_hipvol
## D = -0.46869, boot.n = 2000, boot.stratified = 1, p-value = 0.6393
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.7147177 0.7500000
cat("Comparison of MTA score vs amigdala volume:\n")
## Comparison of MTA score vs amigdala volume:
print(test_1_vs_3)
##
## Bootstrap test for two correlated ROC curves
##
## data: roc_MTAscore and roc_amigvol
## D = 0.095664, boot.n = 2000, boot.stratified = 1, p-value = 0.9238
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.7147177 0.7056452
cat("Comparison of hipocampal volume vs amigdala volume:\n")
## Comparison of hipocampal volume vs amigdala volume:
print(test_2_vs_3)
##
## DeLong's test for two correlated ROC curves
##
## data: roc_hipvol and roc_amigvol
## Z = 1.0179, p-value = 0.3087
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
## -0.0512403 0.1619669
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.7525952 0.6972318
library(ggplot2)
library(dplyr)
# Crear listas de especificidad y sensibilidad con longitudes iguales
roc_list <- list(
"MTA score" = roc_MTAscore,
"Hippocampal volume" = roc_hipvol,
"Amygdala volume" = roc_amigvol
)
roc_data <- bind_rows(lapply(names(roc_list), function(name) {
roc_obj <- roc_list[[name]]
data.frame(
specificity = rev(roc_obj$specificities), # Revertir para asegurar el orden correcto
sensitivity = rev(roc_obj$sensitivities),
model = name
)
}))
# Gráfico de curvas ROC con ggplot2
ggplot(roc_data, aes(x = 1 - specificity, y = sensitivity, color = model)) +
geom_line(size = 1) +
labs(title = "ROC Curves",
x = "1 - Specificity",
y = "Sensitivity",
color = "Model") +
theme_minimal() +
theme(legend.position = "bottom")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.