诊断准确性试验（DTA）Meta 分析方法学指南——以影像学诊断策略为例

library(dplyr)
library(readr)
library(meta)
library(mada)
library(lme4)
library(ggplot2)

# 读取清洗后的四格表数据
df <- read_csv("analysis_dataset.csv") %>%
  filter(stream == "A" & final_decision == "Include") %>%
  select(
    study_id,
    tp, fp, fn, tn,
    reference_layer,
    strategy_group,
    index_test
  )

# 校验样本量
df <- df %>%
  mutate(
    total_2x2 = tp + fp + fn + tn,
    se_raw = tp / (tp + fn),
    sp_raw = tn / (tn + fp)
  )

# 敏感度合并（使用连续性校正后的 adj_* 列）
m_se <- metaprop(
  event = adj_tp,
  n = adj_tp + adj_fn,
  studlab = study_id,
  data = df,
  sm = "PLOGIT",
  method = "Inverse",
  method.tau = "REML",
  incr = 0,
  method.incr = "only0",
  method.ci = "NAsm",
  warn = FALSE
)

# 特异度合并
m_sp <- metaprop(
  event = adj_tn,
  n = adj_tn + adj_fp,
  studlab = study_id,
  data = df,
  sm = "PLOGIT",
  method = "Inverse",
  method.tau = "REML",
  incr = 0,
  method.incr = "only0",
  method.ci = "NAsm",
  warn = FALSE
)

# 提取结果（I² 需乘以 100 转为百分比）
pooled_se <- plogis(m_se$TE.random)
pooled_sp <- plogis(m_sp$TE.random)
i2_se <- m_se$I2 * 100
i2_sp <- m_sp$I2 * 100

# 构建双变量模型所需的数据框（mada 要求大写列名）
# 使用原始 tp/fp/fn/tn，不使用 adj_* 校正值
biv_df <- data.frame(
  TP = df$tp,
  FP = df$fp,
  FN = df$fn,
  TN = df$tn
)

# 拟合双变量模型（Reitsma et al. 方法）
biv_model <- reitsma(biv_df, correction = 0.5, correction.control = "single")

# 提取相关系数和 AUC
summary(biv_model)
Psi <- biv_model$Psi
rho <- Psi[1, 2] / sqrt(Psi[1, 1] * Psi[2, 2])
auc <- summary(biv_model)$AUC$AUC

# 判断 rho 是否到达边界
if (abs(abs(rho) - 1) < 0.01) {
  message("Bivariate rho at boundary — use logit REML for Se/Sp.")
}

# GLMM 不使用 0.5 连续性校正
glmm_model <- tryCatch(
  glmer(
    cbind(tp, fn) ~ 1 + (1 | study_id),
    family = binomial(link = "logit"),
    data = df,
    control = glmerControl(optimizer = "bobyqa")
  ),
  warning = function(w) {
    message("GLMM warning: ", w$message)
    NULL
  },
  error = function(e) {
    message("GLMM failed to converge — recording and moving on.")
    NULL
  }
)

if (is.null(glmm_model)) {
  message("GLMM did not converge — using univariate REML as primary.")
}

# 基于 mada 包绘制 SROC 曲线
# mada 提供 SROC() 函数绘制概要 ROC
SROC(biv_model,
     xlim = c(0, 1), ylim = c(0, 1),
     xlab = "1 - Specificity",
     ylab = "Sensitivity",
     main = "Summary ROC Curve")

# 若 rho 到达边界，标注为探索性
if (abs(rho) > 0.99) {
  message("Bivariate rho at boundary — SROC is exploratory only.")
}

# 敏感度森林图
forest(m_se,
       sortvar = TE,
       leftcols = c("studlab", "event", "n"),
       rightcols = c("effect", "ci"),
       xlab = "Sensitivity (logit scale)",
       cex = 0.8)

# 保存为 SVG
svg("forest_sensitivity.svg", width = 10, height = 8)
forest(m_se)
dev.off()

# 按参考标准分层
biopsy_layer <- df %>% filter(reference_layer == "Biopsy")
clinical_layer <- df %>% filter(reference_layer == "Clinical")

# 活检层合并（研究数≥4 时可做定量合并）
m_biopsy_se <- metaprop(
  event = adj_tp, n = adj_tp + adj_fn,
  studlab = study_id,
  data = biopsy_layer,
  sm = "PLOGIT", method = "Inverse",
  method.tau = "REML",
  incr = 0, method.incr = "only0",
  method.ci = "NAsm", warn = FALSE
)

# 临床诊断层（仅 2 项，描述性报告）
clinical_layer %>%
  select(study_id, se_raw, sp_raw) %>%
  print()

# 仅保留同研究内配对比较
strategy_gain <- df %>%
  filter(strategy_group %in% c("Standard", "Extended")) %>%
  select(study_id, strategy_group, tp, fp, fn, tn) %>%
  pivot_wider(
    names_from = strategy_group,
    values_from = c(tp, fp, fn, tn),
    names_sep = "_"
  ) %>%
  filter(!is.na(tp_Standard) & !is.na(tp_Extended)) %>%
  mutate(
    delta_se = (tp_Extended / (tp_Extended + fn_Extended)) -
               (tp_Standard / (tp_Standard + fn_Standard)),
    additional_tp = tp_Extended - tp_Standard,
    reduced_fn = fn_Standard - fn_Extended,
    additional_fp = fp_Extended - fp_Standard
  )

# 生成临床后果表
scenarios <- expand.grid(
  scenario = c("Core DTA Benchmark", "Biopsy Layer", "Standard Strategy", "Extended Strategy"),
  prevalence = c(0.10, 0.25, 0.50)
)

clinical_consequences <- scenarios %>%
  left_join(se_sp_estimates, by = "scenario") %>%
  mutate(
    tp = round(1000 * prevalence * se),
    fn = round(1000 * prevalence * (1 - se)),
    tn = round(1000 * (1 - prevalence) * sp),
    fp = round(1000 * (1 - prevalence) * (1 - sp))
  )

# QUADAS-2 汇总：各域内 LOW/UNCLEAR/HIGH 的比例
quadas2_summary <- df %>%
  select(study_id, starts_with("quadas_")) %>%
  pivot_longer(
    cols = starts_with("quadas_"),
    names_to = "domain",
    values_to = "rating"
  ) %>%
  group_by(domain, rating) %>%
  summarise(n = n(), .groups = "drop_last") %>%
  mutate(pct = n / sum(n) * 100)

# 堆叠条形图
ggplot(quadas2_summary, aes(x = domain, y = pct, fill = rating)) +
  geom_col(position = "stack") +
  scale_fill_manual(values = c("LOW" = "#4CAF50", "UNCLEAR" = "#FFC107", "HIGH" = "#F44336")) +
  labs(y = "Percentage (%)", fill = "Risk of Bias") +
  theme_minimal()

# 亚组分析示例（使用 study-level 数据，避免重复计数）
subgroup_factors <- c("steroid_exposure", "operator_training",
                       "probe_frequency", "region")

for (factor in subgroup_factors) {
  result <- tryCatch({
    metaprop(
      event = tp, n = tp + fn,
      studlab = study_id,
      data = study_level_df,
      sm = "PLOGIT", method = "Inverse",
      method.tau = "REML",
      incr = 0, method.incr = "only0",
      subgroup = study_level_df[[factor]]
    )
  }, error = function(e) {
    message(paste("Subgroup", factor, "did not converge:", e$message))
    NULL
  })
}