library(haven)
LATElong <- read_sav("LATElong.sav")
lat2<-read_sav("ROC.sav")library(cutpointr)## Warning: package 'cutpointr' was built under R version 4.4.3hip<-LATElong %>% filter(region_type==7)
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.6972 51 17 34
##
## optimal_cutpoint youden acc sensitivity specificity tp fn fp tn
## 3301.2 0.5 0.8039 0.5882 0.9118 10 7 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%
## optimal_cutpoint 2568.80 2899.40 3199.80 3301.20 3223.25 3301.20 3301.20
## AUC_b 0.36 0.55 0.64 0.70 0.70 0.76 0.84
## AUC_oob 0.07 0.46 0.61 0.69 0.69 0.78 0.91
## youden_b 0.14 0.32 0.43 0.51 0.52 0.60 0.72
## youden_oob -0.24 0.06 0.29 0.41 0.40 0.52 0.72
## acc_b 0.47 0.71 0.78 0.82 0.81 0.84 0.90
## acc_oob 0.35 0.63 0.72 0.78 0.77 0.83 0.90
## sensitivity_b 0.24 0.40 0.50 0.59 0.59 0.67 0.80
## sensitivity_oob 0.00 0.17 0.38 0.50 0.49 0.60 0.80
## specificity_b 0.30 0.83 0.90 0.94 0.92 0.97 1.00
## specificity_oob 0.09 0.80 0.88 0.92 0.91 1.00 1.00
## cohens_kappa_b 0.11 0.33 0.46 0.55 0.54 0.63 0.75
## cohens_kappa_oob -0.25 0.07 0.31 0.44 0.43 0.56 0.74
## Max. SD NAs
## 3926.00 151.81 0
## 0.99 0.09 0
## 1.00 0.13 0
## 0.95 0.13 0
## 1.00 0.20 0
## 0.98 0.06 0
## 1.00 0.09 0
## 1.00 0.13 0
## 1.00 0.20 0
## 1.00 0.07 0
## 1.00 0.10 0
## 0.95 0.13 0
## 1.00 0.20 0amig<-LATElong %>% filter(region_type==1)
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.7526 51 17 34
##
## optimal_cutpoint youden acc sensitivity specificity tp fn fp tn
## 1256.5 0.2941 0.6471 0.6471 0.6471 11 6 12 22
##
## 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.
## optimal_cutpoint 739.10 918.30 1132.20 1141.60 1218.38 1371.40 1371.40 1456.50
## AUC_b 0.44 0.62 0.71 0.76 0.75 0.81 0.87 0.94
## AUC_oob 0.31 0.56 0.68 0.76 0.75 0.83 0.93 1.00
## youden_b 0.13 0.32 0.43 0.51 0.51 0.58 0.68 0.82
## youden_oob -0.38 -0.02 0.17 0.30 0.29 0.42 0.61 1.00
## acc_b 0.45 0.63 0.69 0.75 0.74 0.80 0.86 0.92
## acc_oob 0.31 0.50 0.60 0.65 0.65 0.71 0.81 1.00
## sensitivity_b 0.21 0.50 0.68 0.80 0.78 0.90 1.00 1.00
## sensitivity_oob 0.00 0.17 0.43 0.64 0.62 0.83 1.00 1.00
## specificity_b 0.26 0.48 0.61 0.74 0.73 0.84 0.97 1.00
## specificity_oob 0.00 0.38 0.54 0.69 0.67 0.82 1.00 1.00
## cohens_kappa_b 0.09 0.28 0.39 0.47 0.47 0.55 0.67 0.79
## cohens_kappa_oob -0.38 -0.02 0.15 0.26 0.26 0.38 0.55 1.00
## SD NAs
## 155.72 0
## 0.08 0
## 0.11 0
## 0.11 0
## 0.19 0
## 0.07 0
## 0.09 0
## 0.15 0
## 0.26 0
## 0.15 0
## 0.19 0
## 0.12 0
## 0.18 0library(pROC)## Warning: package 'pROC' was built under R version 4.4.3## 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, varhip$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 < casesroc_hipvol <- roc(hip$LATE_status,hip$MRI )## Setting levels: control = 0, case = 1## Setting direction: controls > casesroc_amigvol<- roc(amig$LATE_status,amig$MRI )## Setting levels: control = 0, case = 1
## Setting direction: controls > casesauc_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)
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.6972318 (95% CI: 0.5201028 - 0.8743609 )cat("AUC for amygdala volume:", auc_amig, " (95% CI:", ci_amig[1], "-", ci_amig[3], ")\n")## AUC for amygdala volume: 0.7525952 (95% CI: 0.6041152 - 0.9010751 )modelo_logit <- glm(LATE_status ~ MTAscore, data = hip, family = binomial)
hip_clean <- na.omit(hip[, c("LATE_status", "MTAscore")])
# Ajustar modelo de regresión logística
modelo_logit <- glm(LATE_status ~ MTAscore, data = hip_clean, family = binomial)
# Obtener las probabilidades predichas
hip_clean$predicted_prob <- predict(modelo_logit, type = "response")
# Calcular y graficar la curva ROC
auc(roc(hip_clean$LATE_status, hip_clean$predicted_prob))## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.7147LATElong<-LATElong %>% mutate(mri_corrected=MRI/ICV)library(cutpointr)
hip<-LATElong %>% filter(region_type==7)
set.seed(5620)
opt_cut_h <- cutpointr(data=hip, mri_corrected, 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_h)## Method: maximize_metric
## Predictor: mri_corrected
## Outcome: LATE_status
## Direction: <=
## Nr. of bootstraps: 1000
##
## AUC n n_pos n_neg
## 0.7145 51 17 34
##
## optimal_cutpoint youden acc sensitivity specificity tp fn fp tn
## 0.0023 0.2941 0.5882 0.8235 0.4706 14 3 18 16
##
## Predictor summary:
## Data Min. 5% 1st Qu. Median Mean
## Overall 0.001507977 0.001700485 0.001973293 0.002176247 0.002193918
## 0 0.001642794 0.001791450 0.002078655 0.002269997 0.002274361
## 1 0.001507977 0.001668781 0.001806545 0.002090361 0.002033033
## 3rd Qu. 95% Max. SD NAs
## 0.002399000 0.002740250 0.002853419 0.0003167225 0
## 0.002461308 0.002763442 0.002853419 0.0003088574 0
## 0.002180441 0.002378062 0.002538169 0.0002747328 0
##
## Bootstrap summary:
## Variable Min. 5% 1st Qu. Median Mean 3rd Qu. 95% Max. SD NAs
## optimal_cutpoint 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0
## AUC_b 0.45 0.58 0.67 0.72 0.72 0.77 0.83 0.93 0.07 0
## AUC_oob 0.14 0.52 0.64 0.72 0.71 0.78 0.88 1.00 0.11 0
## youden_b 0.16 0.29 0.41 0.47 0.47 0.54 0.63 0.82 0.10 0
## youden_oob -0.57 -0.09 0.13 0.27 0.26 0.40 0.55 0.83 0.20 0
## acc_b 0.39 0.57 0.65 0.70 0.70 0.75 0.82 0.94 0.08 0
## acc_oob 0.35 0.45 0.55 0.60 0.61 0.67 0.75 0.88 0.09 0
## sensitivity_b 0.22 0.48 0.81 0.89 0.85 0.95 1.00 1.00 0.16 0
## sensitivity_oob 0.00 0.14 0.56 0.75 0.69 0.88 1.00 1.00 0.28 0
## specificity_b 0.19 0.40 0.51 0.59 0.62 0.69 0.97 1.00 0.17 0
## specificity_oob 0.08 0.31 0.44 0.54 0.57 0.67 0.92 1.00 0.18 0
## cohens_kappa_b 0.11 0.23 0.33 0.41 0.41 0.48 0.59 0.84 0.11 0
## cohens_kappa_oob -0.43 -0.09 0.11 0.22 0.21 0.32 0.48 0.71 0.17 0amig<-LATElong %>% filter(region_type==1)
opt_cut_a <- cutpointr(data=amig, mri_corrected, 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_corrected
## Outcome: LATE_status
## Direction: <=
## Nr. of bootstraps: 1000
##
## AUC n n_pos n_neg
## 0.7768 51 17 34
##
## optimal_cutpoint youden acc sensitivity specificity tp fn fp tn
## 8e-04 0.3824 0.6667 0.7647 0.6176 13 4 13 21
##
## Predictor summary:
## Data Min. 5% 1st Qu. Median Mean
## Overall 0.0003964104 0.0005219036 0.0006749899 0.0008141004 0.0007911822
## 0 0.0005501749 0.0005669756 0.0006995834 0.0008792629 0.0008477504
## 1 0.0003964104 0.0004510557 0.0005718825 0.0006838900 0.0006780456
## 3rd Qu. 95% Max. SD NAs
## 0.0009413504 0.0010248821 0.0012616914 0.0001796073 0
## 0.0009682141 0.0010426641 0.0012616914 0.0001672496 0
## 0.0007600523 0.0008855203 0.0009531907 0.0001507173 0
##
## Bootstrap summary:
## Variable Min. 5% 1st Qu. Median Mean 3rd Qu. 95% Max. SD NAs
## optimal_cutpoint 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0
## AUC_b 0.53 0.66 0.73 0.78 0.78 0.83 0.88 0.94 0.07 0
## AUC_oob 0.38 0.60 0.71 0.78 0.78 0.85 0.93 1.00 0.10 0
## youden_b 0.08 0.36 0.46 0.53 0.53 0.60 0.70 0.90 0.10 0
## youden_oob -0.50 0.00 0.20 0.33 0.32 0.46 0.63 0.90 0.20 0
## acc_b 0.39 0.61 0.69 0.73 0.73 0.78 0.84 0.92 0.07 0
## acc_oob 0.29 0.50 0.59 0.65 0.64 0.71 0.79 0.92 0.09 0
## sensitivity_b 0.29 0.64 0.80 0.88 0.86 0.95 1.00 1.00 0.12 0
## sensitivity_oob 0.00 0.25 0.57 0.75 0.71 0.88 1.00 1.00 0.24 0
## specificity_b 0.21 0.45 0.58 0.67 0.67 0.76 0.89 1.00 0.13 0
## specificity_oob 0.00 0.36 0.50 0.62 0.62 0.73 0.88 1.00 0.17 0
## cohens_kappa_b 0.08 0.29 0.39 0.47 0.47 0.55 0.66 0.80 0.11 0
## cohens_kappa_oob -0.38 0.00 0.16 0.28 0.28 0.39 0.55 0.85 0.17 0library(pROC)
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 < casesroc_hipvol <- roc(hip$LATE_status,hip$mri_corrected )## Setting levels: control = 0, case = 1## Setting direction: controls > casesroc_amigvol<- roc(amig$LATE_status,amig$mri_corrected )## Setting levels: control = 0, case = 1
## Setting direction: controls > casesauc_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)
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.7145329 (95% CI: 0.5686067 - 0.860459 )cat("AUC for amygdala volume:", auc_amig, " (95% CI:", ci_amig[1], "-", ci_amig[3], ")\n")## AUC for amygdala volume: 0.7768166 (95% CI: 0.6480574 - 0.9055758 )#curve test
# 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.099051, boot.n = 2000, boot.stratified = 1, p-value = 0.9211
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.7147177 0.7056452cat("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.79907, boot.n = 2000, boot.stratified = 1, p-value = 0.4242
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.7147177 0.7721774cat("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.0178, p-value = 0.3088
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
## -0.1822194 0.0576519
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.7145329 0.7768166