diagnostic-accuracy

# Diagnostic Accuracy Analysis in R ## Overview Comprehensive diagnostic test accuracy analysis covering ROC curve analysis, optimal cutpoint determination, sensitivity and specificity estimation, likelihood ratios, decision curve analysis, inter-rater reliability measures, and diagnostic meta-analysis. ## Basic Diagnostic Measures ### 2x2 Table Analysis ```r library(epiR) # Create 2x2 table # Format: [TP, FN; FP, TN] diag_table <- matrix(c(85, 15, # Disease+ (TP, FN) 20, 180), # Disease- (FP, TN) nrow = 2, byrow = TRUE, dimnames = list( Test = c("Positive", "Negative"), Disease = c("Present", "Absent") )) # Calculate diagnostic measures results <- epi.tests(as.table(diag_table), method = "exact") print(results) # Extract specific measures results$detail # All measures with CIs # Includes: Se, Sp, PPV, NPV, LR+, LR-, DOR, accuracy, prevalence ``` ### Manual Calculations ```r # From confusion matrix elements TP <- 85; FN <- 15; FP <- 20; TN <- 180 # Sensitivity (True Positive Rate) sensitivity <- TP / (TP + FN) # Specificity (True Negative Rate) specificity <- TN / (TN + FP) # Positive Predictive Value ppv <- TP / (TP + FP) # Negative Predictive Value npv <- TN / (TN + FN) # Positive Likelihood Ratio lr_pos <- sensitivity / (1 - specificity) # Negative Likelihood Ratio lr_neg <- (1 - sensitivity) / specificity # Diagnostic Odds Ratio dor <- (TP * TN) / (FP * FN) # Youden's Index (J) youden <- sensitivity + specificity - 1 # Accuracy accuracy <- (TP + TN) / (TP + TN + FP + FN) # Wilson confidence intervals library(binom) se_ci <- binom.confint(TP, TP + FN, method = "wilson") sp_ci <- binom.confint(TN, TN + FP, method = "wilson") ``` ## ROC Curve Analysis ### Basic ROC Curve ```r library(pROC) # Create ROC object roc_obj <- roc( response = df$disease, # Binary outcome (0/1 or factor) predictor = df$biomarker, # Continuous test value levels = c(0, 1), # Control, case direction = "<" # Lower values = control ) # Summary print(roc_obj) # AUC with confidence interval auc(roc_obj) ci.auc(roc_obj, method = "delong") # DeLong method ci.auc(roc_obj, method = "bootstrap", boot.n = 2000) # Bootstrap ``` ### ROC Curve Plotting ```r library(pROC) # Basic plot plot(roc_obj, print.auc = TRUE, print.thres = TRUE) # ggplot2-based ROC curve ggroc(roc_obj) + geom_abline(intercept = 1, slope = 1, linetype = "dashed", color = "gray") + annotate("text", x = 0.25, y = 0.25, label = paste0("AUC = ", round(auc(roc_obj), 3))) + theme_bw() + labs(x = "Specificity", y = "Sensitivity", title = "ROC Curve for Biomarker X") # Multiple ROC curves roc1 <- roc(df$disease, df$biomarker1) roc2 <- roc(df$disease, df$biomarker2) roc3 <- roc(df$disease, df$biomarker3) ggroc(list(Biomarker1 = roc1, Biomarker2 = roc2, Biomarker3 = roc3)) + geom_abline(intercept = 1, slope = 1, linetype = "dashed") + theme_bw() + scale_color_brewer(palette = "Set1") ``` ### Comparing ROC Curves ```r library(pROC) # DeLong test for comparing AUCs roc_test <- roc.test(roc1, roc2, method = "delong") print(roc_test) # Bootstrap comparison roc_test_boot <- roc.test(roc1, roc2, method = "bootstrap", boot.n = 2000) # Venkatraman test (for entire curves, not just AUC) roc_test_venk <- roc.test(roc1, roc2, method = "venkatraman") # Compare multiple markers comparison <- data.frame( Marker = c("Biomarker1", "Biomarker2", "Biomarker3"), AUC = c(auc(roc1), auc(roc2), auc(roc3)), CI_Lower = c(ci.auc(roc1)[1], ci.auc(roc2)[1], ci.auc(roc3)[1]), CI_Upper = c(ci.auc(roc1)[3], ci.auc(roc2)[3], ci.auc(roc3)[3]) ) ``` ### Partial AUC ```r library(pROC) # Partial AUC for high specificity region (Sp > 0.9) pauc_spec <- auc(roc_obj, partial.auc = c(1, 0.9), partial.auc.focus = "specificity") # Partial AUC for high sensitivity region (Se > 0.9) pauc_sens <- auc(roc_obj, partial.auc = c(0.9, 1), partial.auc.focus = "sensitivity") # Standardized partial AUC (McClish correction) pauc_std <- auc(roc_obj, partial.auc = c(1, 0.9), partial.auc.focus = "specificity", partial.auc.correct = TRUE) # CI for partial AUC ci.auc(roc_obj, partial.auc = c(1, 0.9), partial.auc.focus = "specificity") ``` ## Optimal Cutpoint Selection ### Using cutpointr Package ```r library(cutpointr) # Youden's Index (maximize Se + Sp - 1) cp_youden <- cutpointr( data = df, x = biomarker, class = disease, method = maximize_metric, metric = youden ) summary(cp_youden) plot(cp_youden) # Maximize sensitivity with specificity >= 0.9 cp_constrained <- cutpointr( df, biomarker, disease, method = maximize_metric, metric = sensitivity, tol_metric = specificity, tol_threshold = 0.9 ) # Cost-based optimization cp_cost <- cutpointr( df, biomarker, disease, method = minimize_metric, metric = misclassification_cost, cost_fp = 1, # Cost of false positive cost_fn = 5 # Cost of false negative (5x higher) ) # Multiple optimal cutpoints cp_multi <- multi_cutpointr( df, biomarker, disease, method = maximize_metric, metric = youden, boot_cut = 1000 ) ``` ### Using OptimalCutpoints Package ```r library(OptimalCutpoints) # Multiple methods simultaneously opt_cut <- optimal.cutpoints( X = "biomarker", status = "disease", methods = c("Youden", "MaxSpSe", "MaxProdSpSe", "ROC01", "MinValueSp", "MinValueSe", "MaxEfficiency"), data = df, tag.healthy = 0 ) summary(opt_cut) # Cost-benefit method opt_cost <- optimal.cutpoints( X = "biomarker", status = "disease", methods = "CB", data = df, tag.healthy = 0, costs.ratio = 5, # FN cost / FP cost prevalence = 0.10 # Disease prevalence ) ``` ### Manual Cutpoint Selection ```r library(pROC) # Youden's optimal threshold coords_youden <- coords(roc_obj, "best", best.method = "youden") print(coords_youden) # Closest to (0,1) corner coords_closest <- coords(roc_obj, "best", best.method = "closest.topleft") # At specific sensitivity/specificity coords_se90 <- coords(roc_obj, x = 0.90, input = "sensitivity", ret = c("threshold", "sensitivity", "specificity")) # All coordinates all_coords <- coords(roc_obj, x = "all", ret = c("threshold", "sensitivity", "specificity", "ppv", "npv", "accuracy")) ``` ## Decision Curve Analysis ### Using dcurves Package ```r library(dcurves) # Fit prediction models model1 <- glm(cancer ~ age + psa, data = df, family = binomial) model2 <- glm(cancer ~ age + psa + dre, data = df, family = binomial) df$pred1 <- predict(model1, type = "response") df$pred2 <- predict(model2, type = "response") # Decision curve analysis dca_result <- dca( cancer ~ pred1 + pred2, data = df, thresholds = seq(0, 0.5, by = 0.01), label = list(pred1 = "PSA Model", pred2 = "PSA + DRE Model") ) # Plot decision curves plot(dca_result) # Customized plot plot(dca_result, smooth = TRUE) + ggplot2::coord_cartesian(ylim = c(-0.05, 0.2)) + ggplot2::labs(x = "Threshold Probability", y = "Net Benefit", title = "Decision Curve Analysis") ``` ### Net Benefit Calculation ```r library(dcurves) # Extract net benefit at specific threshold nb_data <- as_tibble(dca_result) # Net interventions avoided net_intervention_avoided(dca_result) # Standardized net benefit standardized_net_benefit <- function(nb, prevalence, threshold) { max_nb <- prevalence - (1 - prevalence) * threshold / (1 - threshold) nb / max_nb } ``` ### Clinical Utility Visualization ```r library(dcurves) # Clinical impact plot dca_result |> plot(type = "clinical_impact") # Net benefit with confidence intervals (bootstrap) dca_boot <- dca( cancer ~ pred1, data = df, thresholds = seq(0, 0.5, by = 0.05) ) ``` ## Inter-Rater Reliability ### Cohen's Kappa (Two Raters) ```r library(irr) # Two raters, categorical data ratings <- data.frame( rater1 = c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1), rater2 = c(1, 2, 3, 1, 2, 2, 1, 3, 3, 2) ) # Unweighted kappa (nominal categories) kappa_unweighted <- kappa2(ratings, weight = "unweighted") print(kappa_unweighted) # Linear weighted kappa (ordinal categories) kappa_linear <- kappa2(ratings, weight = "equal") # Quadratic weighted kappa kappa_quadratic <- kappa2(ratings, weight = "squared") # Interpretation: # < 0.20: Poor # 0.21-0.40: Fair # 0.41-0.60: Moderate # 0.61-0.80: Substantial # 0.81-1.00: Almost perfect ``` ### Fleiss' Kappa (Multiple Raters) ```r library(irr) # Multiple raters (each row = subject, each column = rater) ratings_multi <- matrix(c( 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 2 ), nrow = 5, byrow = TRUE) # Fleiss' kappa fleiss_k <- kappam.fleiss(ratings_multi) print(fleiss_k) # Light's kappa (average of all pairwise kappas) light_k <- kappam.light(ratings_multi) ``` ### Intraclass Correlation Coefficient (ICC) ```r library(irr) # Continuous measurements measurements <- data.frame( rater1 = c(2.5, 3.1, 4.2, 2.8, 3.5), rater2 = c(2.4, 3.3, 4.0, 2.9, 3.4), rater3 = c(2.6, 3.0, 4.1, 2.7, 3.6) ) # ICC types: # ICC(1,1): Single rater, absolute agreement # ICC(2,1): Single rater, consistency # ICC(3,1): Single rater, consistency (fixed raters) # ICC(1,k): Average of k raters, absolute agreement # ICC(2,k): Average of k raters, consistency # ICC(3,k): Average of k raters, consistency (fixed raters) # Two-way random effects, single measures, absolute agreement icc_result <- icc(measurements, model = "twoway", type = "agreement", unit = "single") print(icc_result) # Two-way mixed effects, average measures, consistency icc_avg <- icc(measurements, model = "twoway", type = "consistency", unit = "average") ``` ### Agreement for Continuous Data ```r library(BlandAltmanLeh) # Bland-Altman analysis ba <- bland.altman.stats( method1 = df$measurement1, method2 = df$measurement2 ) # Bland-Altman plot bland.altman.plot( method1 = df$measurement1, method2 = df$measurement2, main = "Bland-Altman Plot", xlab = "Mean of Methods",