#data load ##########
LATElong <- read_sav("Data/LATElong.sav")

Step 1: Estimate the optimal cutoff for hippocampal volume.

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%).

Step 2: Estimate the optimal cutoff for amygdala volume.

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()`).

Step 3: Determine the area under the curve (AUC) for the hippocampus, amygdala, and MTA score.

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 )

Comparison of the areas under the curve (AUC).

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

plot

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.