df_fkg <- dados %>%
  filter(!is.na(group), is.finite(fentanyl_per_kg)) %>%
  select(group, fentanyl_val = fentanyl, fentanyl_per_kg)

wins_by_group <- function(x, g, trim = 0.05) {
  for (lv in levels(g)) {
    idx <- which(g == lv & is.finite(x))
    q   <- quantile(x[idx], c(trim, 1 - trim), na.rm = TRUE)
    x[idx] <- pmin(pmax(x[idx], q[1]), q[2])
  }
  x
}

# 1) t de Welch winsorizado (µg/kg)
x_wins <- wins_by_group(df_fkg$fentanyl_per_kg, df_fkg$group, trim = 0.05)
t.test(x_wins ~ df_fkg$group, var.equal = FALSE)

# 2) Teste de Yuen (médias aparadas a 20%)
try(WRS2::yuen(fentanyl_per_kg ~ group, data = df_fkg, tr = 0.2), silent = TRUE)

# 3) Teste de permutação (independência)
set.seed(123)
coin::independence_test(fentanyl_per_kg ~ group, data = df_fkg,
                        distribution = coin::approximate(B = 5000))

# 4) Dose total — Mann–Whitney
wilcox.test(fentanyl_val ~ group, data = df_fkg, exact = FALSE, conf.int = TRUE)

# 5) Bootstrap BCa da diferença de médias (CLO − BMRA), estratificado por grupo
boot_diff <- function(d, i) {
  d <- d[i, ]
  mean(d$fentanyl_val[d$group == "BMRA-CLO"]) - mean(d$fentanyl_val[d$group == "BMRA"])
}
set.seed(123)
br <- boot::boot(df_fkg, boot_diff, R = 2000, strata = df_fkg$group)
boot::boot.ci(br, type = "bca")

# 6) Casos influentes — distância de Cook (leave-one-out)
fit_lm <- lm(fentanyl_per_kg ~ group, data = df_fkg)
inf_df <- df_fkg %>%
  mutate(id = row_number(), cooks = cooks.distance(fit_lm), influente = cooks > 4/n())

ggplot(inf_df, aes(id, cooks, color = influente)) +
  geom_segment(aes(xend = id, y = 0, yend = cooks), color = "gray") +
  geom_point(size = 3) +
  geom_hline(yintercept = 4/nrow(inf_df), linetype = "dashed", color = "red") +
  scale_color_manual(values = c("TRUE" = "#E74C3C", "FALSE" = "#2C3E50")) +
  labs(title = "Distância de Cook (leave-one-out)",
       subtitle = "Pontos acima da linha vermelha distorcem isoladamente a média do grupo",
       x = "Paciente (índice)", y = "Distância de Cook") +
  theme(legend.position = "none")

# -------------------------------------------------------------------------
# Pacotes adicionais para análises robustas
# -------------------------------------------------------------------------
# sandwich/lmtest: erros-padrão robustos HC3 para modelos lineares.
# broom: organização de saídas de modelos.
# purrr: programação funcional para bootstrap/AUC.
# ordinal: modelo ordinal misto para EVN, usado como análise de sensibilidade.
# logistf: regressão logística penalizada de Firth; opcional, útil quando há poucos eventos
#          ou separação completa/quase completa. Se não instalar, o script continua.
# -------------------------------------------------------------------------

pkgs_robustos <- c("sandwich", "lmtest", "broom", "purrr", "ordinal")
for (p in pkgs_robustos) {
  if (!requireNamespace(p, quietly = TRUE)) {
    install.packages(p, repos = "https://cloud.r-project.org")
  }
}

# Firth é opcional. Em alguns computadores pode demorar mais para instalar.
# O código tentará usar se estiver disponível; caso contrário, usa GLM comum com cautela.
if (!requireNamespace("logistf", quietly = TRUE)) {
  message("Pacote 'logistf' não instalado. Regressões de Firth serão omitidas, mas as análises principais continuarão.")
}

# -------------------------------------------------------------------------
# Função: tabela de coeficientes com erro-padrão robusto HC3
# -------------------------------------------------------------------------
# HC3 é recomendado em amostras pequenas por ser mais conservador que HC0/HC1
# diante de heterocedasticidade e observações de alta alavancagem.
# -------------------------------------------------------------------------

tabela_hc3 <- function(modelo, titulo = NULL) {
  vc <- sandwich::vcovHC(modelo, type = "HC3")
  ct <- lmtest::coeftest(modelo, vcov. = vc, df = df.residual(modelo))
  out <- as.data.frame(unclass(ct)) %>%
    tibble::rownames_to_column("Termo") %>%
    dplyr::rename(
      Estimativa = Estimate,
      EP_robusto = `Std. Error`,
      Estatistica = `t value`,
      p = `Pr(>|t|)`
    ) %>%
    dplyr::mutate(
      IC_low  = Estimativa - qt(0.975, df.residual(modelo)) * EP_robusto,
      IC_high = Estimativa + qt(0.975, df.residual(modelo)) * EP_robusto,
      `Estimativa` = nv(Estimativa, 3),
      `EP robusto HC3` = nv(EP_robusto, 3),
      `IC 95% robusto` = sprintf("[%s; %s]", nv(IC_low, 3), nv(IC_high, 3)),
      `Valor-p robusto` = ifelse(p < 0.001, "<0,001", nv(p, 3))
    ) %>%
    dplyr::select(Termo, Estimativa, `EP robusto HC3`, `IC 95% robusto`, `Valor-p robusto`)

  knitr::kable(out, align = c("l", "c", "c", "c", "c"), caption = titulo)
}

# -------------------------------------------------------------------------
# Função: extrair coeficiente do grupo BMRA-CLO em modelos lineares
# -------------------------------------------------------------------------
# Por convenção, groupBMRA-CLO representa a diferença BMRA-CLO − BMRA.
# Valor negativo para fentanil indica menor consumo no grupo clonidina.
# -------------------------------------------------------------------------

extrair_beta_grupo <- function(modelo, termo = "groupBMRA-CLO") {
  cf <- coef(modelo)
  if (!termo %in% names(cf)) return(NA_real_)
  unname(cf[termo])
}

# -------------------------------------------------------------------------
# Bootstrap estratificado por grupo para coeficiente do grupo em regressão linear
# -------------------------------------------------------------------------
# Reamostra pacientes dentro de cada grupo, preservando o tamanho de cada braço.
# Isso é mais adequado para ECR porque mantém a estrutura de alocação.
# -------------------------------------------------------------------------

bootstrap_beta_lm <- function(data, formula, termo = "groupBMRA-CLO", R = 5000,
                              strata = "group", seed = 123) {
  set.seed(seed)
  data <- data %>% dplyr::filter(stats::complete.cases(model.frame(formula, data = .)))
  data[[strata]] <- factor(data[[strata]], levels = c("BMRA", "BMRA-CLO"))

  beta_obs <- extrair_beta_grupo(lm(formula, data = data), termo = termo)
  grupos <- split(data, data[[strata]])

  betas <- replicate(R, {
    amostra <- dplyr::bind_rows(lapply(grupos, function(g) g[sample(seq_len(nrow(g)), replace = TRUE), , drop = FALSE]))
    amostra[[strata]] <- factor(amostra[[strata]], levels = c("BMRA", "BMRA-CLO"))
    tryCatch(extrair_beta_grupo(lm(formula, data = amostra), termo = termo), error = function(e) NA_real_)
  })

  betas <- betas[is.finite(betas)]
  tibble::tibble(
    Estimativa = beta_obs,
    IC_low = as.numeric(quantile(betas, 0.025, na.rm = TRUE)),
    IC_high = as.numeric(quantile(betas, 0.975, na.rm = TRUE)),
    R_validos = length(betas)
  )
}

# -------------------------------------------------------------------------
# Teste de permutação para coeficiente do grupo em regressão linear
# -------------------------------------------------------------------------
# Em um ensaio randomizado, permutar os rótulos dos grupos simula a hipótese nula
# compatível com a própria lógica da randomização.
# -------------------------------------------------------------------------

permutacao_beta_lm <- function(data, formula, termo = "groupBMRA-CLO", B = 10000,
                               seed = 123) {
  set.seed(seed)
  data <- data %>% dplyr::filter(stats::complete.cases(model.frame(formula, data = .)))
  data$group <- factor(data$group, levels = c("BMRA", "BMRA-CLO"))

  beta_obs <- extrair_beta_grupo(lm(formula, data = data), termo = termo)

  betas_perm <- replicate(B, {
    dp <- data
    dp$group <- factor(sample(dp$group), levels = c("BMRA", "BMRA-CLO"))
    tryCatch(extrair_beta_grupo(lm(formula, data = dp), termo = termo), error = function(e) NA_real_)
  })

  betas_perm <- betas_perm[is.finite(betas_perm)]
  p_perm <- mean(abs(betas_perm) >= abs(beta_obs))

  tibble::tibble(
    Estimativa_observada = beta_obs,
    `p por permutação` = p_perm,
    B_validas = length(betas_perm)
  )
}

# -------------------------------------------------------------------------
# Função geral para apresentar resultados bootstrap/permutação
# -------------------------------------------------------------------------

fmt_boot <- function(x, unidade = "") {
  x %>%
    dplyr::mutate(
      Estimativa = paste0(nv(Estimativa, 3), unidade),
      `IC 95% bootstrap` = sprintf("[%s; %s]%s", nv(IC_low, 3), nv(IC_high, 3), unidade)
    ) %>%
    dplyr::select(Estimativa, `IC 95% bootstrap`, R_validos)
}

# -------------------------------------------------------------------------
# Três níveis de ajuste para o desfecho primário
# -------------------------------------------------------------------------
# 1) Modelo bruto: estima o efeito simples da randomização.
# 2) Modelo mínimo: adiciona IMC, por ser clinicamente relevante para dose/volume corporal.
# 3) Modelo ajustado principal: idade, IMC, sexo e ASA, conforme análise principal.
# -------------------------------------------------------------------------

df_fent_rob <- dados %>%
  dplyr::filter(!is.na(group), is.finite(fentanyl_per_kg), !is.na(idade), !is.na(imc), !is.na(sexo), !is.na(asa)) %>%
  dplyr::mutate(group = factor(group, levels = c("BMRA", "BMRA-CLO")))

fit_fent_bruto   <- lm(fentanyl_per_kg ~ group, data = df_fent_rob)
fit_fent_minimo  <- lm(fentanyl_per_kg ~ group + imc, data = df_fent_rob)
fit_fent_ajustado <- lm(fentanyl_per_kg ~ group + idade + imc + sexo + asa, data = df_fent_rob)

# Tabela comparativa dos modelos com estimativa clássica do coeficiente de grupo
modelos_fent <- list(
  "Bruto" = fit_fent_bruto,
  "Mínimo: grupo + IMC" = fit_fent_minimo,
  "Ajustado: grupo + idade + IMC + sexo + ASA" = fit_fent_ajustado
)

res_modelos_fent <- purrr::imap_dfr(modelos_fent, function(m, nm) {
  ci <- confint(m)
  beta <- extrair_beta_grupo(m)
  tibble::tibble(
    Modelo = nm,
    `β grupo BMRA-CLO − BMRA` = beta,
    `IC 95% low` = ci["groupBMRA-CLO", 1],
    `IC 95% high` = ci["groupBMRA-CLO", 2],
    `p clássico` = summary(m)$coefficients["groupBMRA-CLO", "Pr(>|t|)"],
    `R² ajustado` = summary(m)$adj.r.squared
  )
}) %>%
  dplyr::mutate(
    `β grupo BMRA-CLO − BMRA` = nv(`β grupo BMRA-CLO − BMRA`, 3),
    `IC 95%` = sprintf("[%s; %s]", nv(`IC 95% low`, 3), nv(`IC 95% high`, 3)),
    `p clássico` = ifelse(`p clássico` < 0.001, "<0,001", nv(`p clássico`, 3)),
    `R² ajustado` = nv(`R² ajustado`, 3)
  ) %>%
  dplyr::select(Modelo, `β grupo BMRA-CLO − BMRA`, `IC 95%`, `p clássico`, `R² ajustado`)

knitr::kable(res_modelos_fent, align = c("l", "c", "c", "c", "c"),
             caption = "Fentanil µg/kg — estabilidade do efeito do grupo em modelos com diferentes níveis de ajuste")

# Tabela robusta HC3 do modelo principal
cat("\n\n**Modelo ajustado principal com erro-padrão robusto HC3**\n\n")

tabela_hc3(fit_fent_ajustado, titulo = "ANCOVA robusta HC3 — fentanil µg/kg")

# -------------------------------------------------------------------------
# Bootstrap e permutação do coeficiente de grupo no modelo ajustado
# -------------------------------------------------------------------------

boot_fent <- bootstrap_beta_lm(
  data = df_fent_rob,
  formula = fentanyl_per_kg ~ group + idade + imc + sexo + asa,
  R = 5000,
  seed = 2026
)

perm_fent <- permutacao_beta_lm(
  data = df_fent_rob,
  formula = fentanyl_per_kg ~ group + idade + imc + sexo + asa,
  B = 10000,
  seed = 2026
)

knitr::kable(fmt_boot(boot_fent, unidade = " µg/kg"), align = c("c", "c", "c"),
             caption = "Bootstrap estratificado — efeito ajustado do grupo sobre fentanil µg/kg")

perm_tab <- perm_fent %>%
  dplyr::mutate(
    Estimativa_observada = paste0(nv(Estimativa_observada, 3), " µg/kg"),
    `p por permutação` = ifelse(`p por permutação` < 0.001, "<0,001", nv(`p por permutação`, 3))
  )
knitr::kable(perm_tab, align = c("c", "c", "c"),
             caption = "Teste de permutação — efeito ajustado do grupo sobre fentanil µg/kg")

# -------------------------------------------------------------------------
# Análise leave-one-out
# -------------------------------------------------------------------------
# Retira um paciente por vez e recalcula o coeficiente do grupo.
# Se o efeito permanece negativo e de magnitude semelhante, a conclusão não
# depende de um único paciente extremo.
# -------------------------------------------------------------------------

loo_fent <- purrr::map_dfr(seq_len(nrow(df_fent_rob)), function(i) {
  dloo <- df_fent_rob[-i, , drop = FALSE]
  fit <- lm(fentanyl_per_kg ~ group + idade + imc + sexo + asa, data = dloo)
  tibble::tibble(
    paciente_retirado = as.character(df_fent_rob$id[i]),
    beta_grupo = extrair_beta_grupo(fit),
    cooks_no_modelo_original = cooks.distance(fit_fent_ajustado)[i]
  )
})

res_loo <- tibble::tibble(
  `β original` = extrair_beta_grupo(fit_fent_ajustado),
  `Menor β leave-one-out` = min(loo_fent$beta_grupo, na.rm = TRUE),
  `Maior β leave-one-out` = max(loo_fent$beta_grupo, na.rm = TRUE),
  `Maior Cook no modelo original` = max(loo_fent$cooks_no_modelo_original, na.rm = TRUE),
  `Limite Cook 4/n` = 4 / nrow(df_fent_rob)
) %>%
  dplyr::mutate(dplyr::across(dplyr::everything(), ~ nv(.x, 3)))

knitr::kable(res_loo, align = c("c", "c", "c", "c", "c"),
             caption = "Leave-one-out e distância de Cook — fentanil µg/kg")

ggplot(loo_fent, aes(x = seq_along(beta_grupo), y = beta_grupo)) +
  geom_hline(yintercept = extrair_beta_grupo(fit_fent_ajustado), linetype = "dashed") +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_point(size = 2.5) +
  labs(
    title = "Leave-one-out — estabilidade do efeito da clonidina sobre fentanil µg/kg",
    subtitle = "Cada ponto representa o coeficiente do grupo após retirar um paciente",
    x = "Paciente retirado, em ordem do banco",
    y = "β grupo BMRA-CLO − BMRA"
  )

# -------------------------------------------------------------------------
# Área sob a curva da EVN: carga total de dor em 48 horas
# -------------------------------------------------------------------------
# AUC pelo método trapezoidal usando tempos: SRPA = 0h, 24h e 48h.
# A AUC tem unidade EVN*h. Uma AUC menor indica menor carga total de dor.
# Também calculamos AUC média = AUC/48, interpretável como EVN média ponderada
# pelo tempo ao longo das primeiras 48 horas.
# -------------------------------------------------------------------------

calc_auc <- function(y, tempo = c(0, 24, 48)) {
  ok <- is.finite(y) & is.finite(tempo)
  y <- y[ok]
  tempo <- tempo[ok]
  if (length(y) < 2) return(NA_real_)
  ord <- order(tempo)
  y <- y[ord]
  tempo <- tempo[ord]
  sum(diff(tempo) * (head(y, -1) + tail(y, -1)) / 2)
}

dor_auc <- dados %>%
  dplyr::mutate(group = factor(group, levels = c("BMRA", "BMRA-CLO"))) %>%
  dplyr::transmute(
    id, group, idade, imc, sexo, asa,
    auc_evn_repouso = purrr::pmap_dbl(
      list(`NRS RPA r`, `NRS 24h r`, `NRS 48h r`),
      function(a, b, c) calc_auc(c(a, b, c))
    ),
    auc_evn_movimento = purrr::pmap_dbl(
      list(`NRS RPA m`, `NRS 24h m`, `NRS 48h m`),
      function(a, b, c) calc_auc(c(a, b, c))
    )
  ) %>%
  dplyr::mutate(
    evn_media_48h_repouso = auc_evn_repouso / 48,
    evn_media_48h_movimento = auc_evn_movimento / 48
  )

comparar_auc <- function(var_auc, label) {
  form <- as.formula(paste0(var_auc, " ~ group + idade + imc + sexo + asa"))
  d <- dor_auc %>% dplyr::filter(stats::complete.cases(model.frame(form, data = .)))
  fit <- lm(form, data = d)
  vc <- sandwich::vcovHC(fit, type = "HC3")
  ct <- lmtest::coeftest(fit, vcov. = vc, df = df.residual(fit))
  beta <- ct["groupBMRA-CLO", "Estimate"]
  se <- ct["groupBMRA-CLO", "Std. Error"]
  p <- ct["groupBMRA-CLO", "Pr(>|t|)"]
  ic_low <- beta - qt(0.975, df.residual(fit)) * se
  ic_high <- beta + qt(0.975, df.residual(fit)) * se

  boot <- bootstrap_beta_lm(d, form, R = 3000, seed = 2026)
  perm <- permutacao_beta_lm(d, form, B = 5000, seed = 2026)

  tibble::tibble(
    Desfecho = label,
    `β ajustado HC3` = beta,
    `IC 95% HC3` = sprintf("[%s; %s]", nv(ic_low, 2), nv(ic_high, 2)),
    `p HC3` = p,
    `IC 95% bootstrap` = sprintf("[%s; %s]", nv(boot$IC_low, 2), nv(boot$IC_high, 2)),
    `p permutação` = perm$`p por permutação`
  )
}

auc_tab <- dplyr::bind_rows(
  comparar_auc("auc_evn_repouso", "AUC EVN repouso 0–48h"),
  comparar_auc("auc_evn_movimento", "AUC EVN movimento 0–48h"),
  comparar_auc("evn_media_48h_repouso", "EVN média ponderada repouso"),
  comparar_auc("evn_media_48h_movimento", "EVN média ponderada movimento")
) %>%
  dplyr::mutate(
    `β ajustado HC3` = nv(`β ajustado HC3`, 2),
    `p HC3` = ifelse(`p HC3` < 0.001, "<0,001", nv(`p HC3`, 3)),
    `p permutação` = ifelse(`p permutação` < 0.001, "<0,001", nv(`p permutação`, 3))
  )

knitr::kable(auc_tab, align = c("l", "c", "c", "c", "c", "c"),
             caption = "Dor pós-operatória — AUC e EVN média ponderada em 0–48h")

# Gráfico exploratório da carga total de dor em movimento

ggplot(dor_auc, aes(x = group, y = auc_evn_movimento, fill = group)) +
  geom_boxplot(width = 0.55, alpha = 0.75, outlier.shape = NA) +
  geom_jitter(width = 0.08, alpha = 0.75, size = 2) +
  scale_fill_manual(values = pal_grupo) +
  labs(
    title = "Carga total de dor em movimento — AUC EVN 0–48h",
    subtitle = "Menor AUC indica menor carga acumulada de dor no período pós-operatório",
    x = NULL,
    y = "AUC EVN movimento (EVN*h)"
  ) +
  theme(legend.position = "none")

# -------------------------------------------------------------------------
# Modelo ordinal misto para EVN em movimento
# -------------------------------------------------------------------------
# A EVN é uma escala discreta de 0 a 10. O modelo linear misto é útil e
# interpretável, mas o modelo ordinal respeita melhor a natureza ordinal do escore.
# Esta análise deve ser interpretada como sensibilidade, não como substituta do LMM.
# -------------------------------------------------------------------------

de_mov_ord <- de_mov %>%
  dplyr::mutate(
    nrs_round = round(nrs),
    nrs_round = pmax(0, pmin(10, nrs_round)),
    nrs_ord = ordered(nrs_round, levels = 0:10),
    group = factor(group, levels = c("BMRA", "BMRA-CLO")),
    tempo = factor(tempo, levels = c("SRPA", "24h", "48h"))
  ) %>%
  dplyr::filter(!is.na(nrs_ord), !is.na(group), !is.na(tempo), !is.na(id))

fit_ord <- tryCatch(
  ordinal::clmm(nrs_ord ~ group * tempo + (1 | id), data = de_mov_ord, Hess = TRUE),
  error = function(e) e
)

if (inherits(fit_ord, "error")) {
  cat("O modelo ordinal misto não convergiu nesta amostra. Mensagem do R:\n")
  cat(fit_ord$message)
} else {
  ord_coef <- coef(summary(fit_ord)) %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Termo") %>%
    dplyr::filter(grepl("group|tempo", Termo)) %>%
    dplyr::mutate(
      `Odds ratio` = exp(Estimate),
      `IC 95% low` = exp(Estimate - 1.96 * `Std. Error`),
      `IC 95% high` = exp(Estimate + 1.96 * `Std. Error`),
      `Valor-p` = 2 * pnorm(abs(`z value`), lower.tail = FALSE)
    ) %>%
    dplyr::transmute(
      Termo,
      `OR acumulativo` = nv(`Odds ratio`, 2),
      `IC 95%` = sprintf("[%s; %s]", nv(`IC 95% low`, 2), nv(`IC 95% high`, 2)),
      `Valor-p` = ifelse(`Valor-p` < 0.001, "<0,001", nv(`Valor-p`, 3))
    )

  knitr::kable(ord_coef, align = c("l", "c", "c", "c"),
               caption = "Modelo ordinal misto — EVN em movimento")
}

# -------------------------------------------------------------------------
# QoR-15 em 48h: ANCOVA robusta, bootstrap e permutação
# -------------------------------------------------------------------------
# Como o QoR-15 pode apresentar efeito teto, a análise linear é complementada
# por bootstrap/permutação e por análise categórica de recuperação excelente.
# -------------------------------------------------------------------------

df_qor_rob <- dados %>%
  dplyr::filter(!is.na(group), is.finite(qor_48h), is.finite(qor_srpa), !is.na(idade), !is.na(imc), !is.na(sexo)) %>%
  dplyr::mutate(group = factor(group, levels = c("BMRA", "BMRA-CLO")))

fit_qor48_rob <- lm(qor_48h ~ group + qor_srpa + idade + sexo + imc, data = df_qor_rob)

cat("\n\n**QoR-15 em 48h — modelo ajustado com erro-padrão robusto HC3**\n\n")

tabela_hc3(fit_qor48_rob, titulo = "QoR-15 48h — regressão robusta HC3")

boot_qor <- bootstrap_beta_lm(
  data = df_qor_rob,
  formula = qor_48h ~ group + qor_srpa + idade + sexo + imc,
  R = 5000,
  seed = 2026
)

perm_qor <- permutacao_beta_lm(
  data = df_qor_rob,
  formula = qor_48h ~ group + qor_srpa + idade + sexo + imc,
  B = 10000,
  seed = 2026
)

knitr::kable(fmt_boot(boot_qor, unidade = " pontos"), align = c("c", "c", "c"),
             caption = "Bootstrap estratificado — efeito ajustado do grupo sobre QoR-15 48h")

perm_qor_tab <- perm_qor %>%
  dplyr::mutate(
    Estimativa_observada = paste0(nv(Estimativa_observada, 2), " pontos"),
    `p por permutação` = ifelse(`p por permutação` < 0.001, "<0,001", nv(`p por permutação`, 3))
  )
knitr::kable(perm_qor_tab, align = c("c", "c", "c"),
             caption = "Teste de permutação — efeito ajustado do grupo sobre QoR-15 48h")

# -------------------------------------------------------------------------
# Recuperação excelente: análise categórica do QoR-15
# -------------------------------------------------------------------------
# Dois limiares são apresentados:
# - QoR-15 >= 145: recuperação muito próxima do teto.
# - QoR-15 >= 135: pelo menos 90% do escore máximo de 150.
# -------------------------------------------------------------------------

df_qor_bin <- df_qor_rob %>%
  dplyr::mutate(
    qor48_excelente_145 = as.integer(qor_48h >= 145),
    qor48_excelente_135 = as.integer(qor_48h >= 135),
    melhora_qor_8pts = as.integer((qor_48h - qor_srpa) >= 8)
  )

calcular_risco_binario <- function(data, evento, label) {
  d <- data %>% dplyr::filter(!is.na(.data[[evento]]), !is.na(group))
  g1 <- d %>% dplyr::filter(group == "BMRA-CLO")
  g0 <- d %>% dplyr::filter(group == "BMRA")

  a <- sum(g1[[evento]] == 1, na.rm = TRUE)
  n1 <- nrow(g1)
  c <- sum(g0[[evento]] == 1, na.rm = TRUE)
  n0 <- nrow(g0)

  p1 <- a / n1
  p0 <- c / n0
  rd <- p1 - p0
  se_rd <- sqrt(p1 * (1 - p1) / n1 + p0 * (1 - p0) / n0)
  rd_low <- rd - 1.96 * se_rd
  rd_high <- rd + 1.96 * se_rd

  # Correção de Haldane-Anscombe para RR quando houver célula zero
  aa <- a; bb <- n1 - a; cc <- c; dd <- n0 - c
  if (any(c(aa, bb, cc, dd) == 0)) {
    aa <- aa + 0.5; bb <- bb + 0.5; cc <- cc + 0.5; dd <- dd + 0.5
  }
  rr <- (aa / (aa + bb)) / (cc / (cc + dd))
  se_log_rr <- sqrt(1 / aa - 1 / (aa + bb) + 1 / cc - 1 / (cc + dd))
  rr_low <- exp(log(rr) - 1.96 * se_log_rr)
  rr_high <- exp(log(rr) + 1.96 * se_log_rr)

  fisher_p <- fisher.test(matrix(c(a, n1 - a, c, n0 - c), nrow = 2, byrow = TRUE))$p.value

  tibble::tibble(
    Desfecho = label,
    `BMRA-CLO` = fmt_n_pct(a, n1),
    BMRA = fmt_n_pct(c, n0),
    `Diferença absoluta de risco` = rd,
    `IC 95% RD` = sprintf("[%s; %s]", nv(rd_low, 2), nv(rd_high, 2)),
    `Risco relativo` = rr,
    `IC 95% RR` = sprintf("[%s; %s]", nv(rr_low, 2), nv(rr_high, 2)),
    `p Fisher` = fisher_p
  )
}

qor_bin_tab <- dplyr::bind_rows(
  calcular_risco_binario(df_qor_bin, "qor48_excelente_145", "QoR-15 48h ≥ 145"),
  calcular_risco_binario(df_qor_bin, "qor48_excelente_135", "QoR-15 48h ≥ 135"),
  calcular_risco_binario(df_qor_bin, "melhora_qor_8pts", "Aumento QoR-15 ≥ 8 pontos desde SRPA")
) %>%
  dplyr::mutate(
    `Diferença absoluta de risco` = nv(`Diferença absoluta de risco`, 2),
    `Risco relativo` = nv(`Risco relativo`, 2),
    `p Fisher` = ifelse(`p Fisher` < 0.001, "<0,001", nv(`p Fisher`, 3))
  )

knitr::kable(qor_bin_tab, align = c("l", "c", "c", "c", "c", "c", "c", "c"),
             caption = "QoR-15 48h — recuperação excelente e melhora clinicamente relevante")

# Regressão logística penalizada de Firth, se o pacote estiver disponível.
# Útil quando há poucos eventos ou separação completa/quase completa.
if (requireNamespace("logistf", quietly = TRUE)) {
  firth_145 <- tryCatch(
    logistf::logistf(qor48_excelente_145 ~ group + qor_srpa + idade + sexo + imc, data = df_qor_bin),
    error = function(e) e
  )
  if (!inherits(firth_145, "error")) {
    cat("\n\n**Regressão logística de Firth para QoR-15 48h ≥ 145**\n\n")
    print(summary(firth_145))
  }
} else {
  cat("\n\nRegressão de Firth omitida porque o pacote 'logistf' não está instalado.\n")
}

# -------------------------------------------------------------------------
# Eventos adversos e segurança
# -------------------------------------------------------------------------
# Em amostras pequenas, p-valores de Fisher raramente demonstram diferença.
# Por isso, o mais informativo é reportar a diferença absoluta de risco, o risco
# relativo e os respectivos intervalos de confiança, reconhecendo a incerteza.
# -------------------------------------------------------------------------

# Lista de eventos de segurança disponíveis no banco.
# Ajuste/adicione eventos caso novos campos sejam incluídos no REDCap.
eventos_seg <- c(
  bradicardya_rpa = "Bradicardia na SRPA",
  hypotension_rpa = "Hipotensão na SRPA",
  ponv_rpa = "NVPO na SRPA",
  general_anest_convertion = "Conversão para anestesia geral",
  hospital_internation_rpa = "Internação hospitalar",
  abdwall_hematoma_rpa = "Hematoma de parede abdominal"
)

eventos_seg <- eventos_seg[names(eventos_seg) %in% names(dados)]

seg_tab <- purrr::imap_dfr(eventos_seg, function(label, var) {
  calcular_risco_binario(dados %>% dplyr::mutate(group = factor(group, levels = c("BMRA", "BMRA-CLO"))), var, label)
}) %>%
  dplyr::mutate(
    `Diferença absoluta de risco` = nv(`Diferença absoluta de risco`, 2),
    `Risco relativo` = nv(`Risco relativo`, 2),
    `p Fisher` = ifelse(`p Fisher` < 0.001, "<0,001", nv(`p Fisher`, 3))
  )

knitr::kable(seg_tab, align = c("l", "c", "c", "c", "c", "c", "c", "c"),
             caption = "Eventos adversos — risco absoluto, risco relativo e teste exato de Fisher")

# -------------------------------------------------------------------------
# Controle exploratório da multiplicidade
# -------------------------------------------------------------------------
# O desfecho primário deve ser interpretado hierarquicamente e não precisa ser
# penalizado da mesma forma que os secundários. Esta tabela aplica FDR de
# Benjamini-Hochberg a um conjunto de desfechos secundários/exploratórios.
# -------------------------------------------------------------------------

p_secundarios <- tibble::tibble(
  Desfecho = c(
    "QoR-15 48h ajustado",
    "AUC EVN repouso 0–48h",
    "AUC EVN movimento 0–48h",
    "EVN média ponderada repouso",
    "EVN média ponderada movimento"
  ),
  p = c(
    as.numeric(summary(fit_qor48_rob)$coefficients["groupBMRA-CLO", "Pr(>|t|)"]),
    as.numeric(gsub(",", ".", auc_tab$`p HC3`[auc_tab$Desfecho == "AUC EVN repouso 0–48h"])),
    as.numeric(gsub(",", ".", auc_tab$`p HC3`[auc_tab$Desfecho == "AUC EVN movimento 0–48h"])),
    as.numeric(gsub(",", ".", auc_tab$`p HC3`[auc_tab$Desfecho == "EVN média ponderada repouso"])),
    as.numeric(gsub(",", ".", auc_tab$`p HC3`[auc_tab$Desfecho == "EVN média ponderada movimento"]))
  )
) %>%
  dplyr::mutate(
    `p FDR` = p.adjust(p, method = "BH"),
    p = ifelse(p < 0.001, "<0,001", nv(p, 3)),
    `p FDR` = ifelse(`p FDR` < 0.001, "<0,001", nv(`p FDR`, 3))
  )

knitr::kable(p_secundarios, align = c("l", "c", "c"),
             caption = "Desfechos secundários/exploratórios — ajuste FDR de Benjamini-Hochberg")

# -------------------------------------------------------------------------
# Síntese final gerada a partir dos objetos das análises robustas
# -------------------------------------------------------------------------

beta_fent <- extrair_beta_grupo(fit_fent_ajustado)
p_fent_perm <- perm_fent$`p por permutação`

beta_qor <- extrair_beta_grupo(fit_qor48_rob)
p_qor_perm <- perm_qor$`p por permutação`

cat("### Interpretação operacional\n\n")