media_idade   <- mean(dados$idade,   na.rm = TRUE)
mediana_idade <- median(dados$idade, na.rm = TRUE)
dp_idade      <- sd(dados$idade,     na.rm = TRUE)
bw_idade      <- diff(range(dados$idade, na.rm = TRUE)) / 20

ggplot(dados, aes(x = idade)) +
  geom_histogram(aes(fill = after_stat(count)),
                 binwidth = bw_idade, boundary = 0, closed = "left", alpha = 0.95) +
  # densidade reescalada para contagem (× binwidth) para casar com as barras
  stat_density(aes(y = after_stat(count) * bw_idade),
               geom = "line", linewidth = 1.1, na.rm = TRUE) +
  geom_vline(xintercept = media_idade,   linewidth = 0.8) +
  geom_vline(xintercept = mediana_idade, linetype = "dashed", linewidth = 0.8) +
  geom_rug(sides = "b", alpha = 0.25, linewidth = 0.3) +
  scale_fill_gradient(low = "#E9EDF5", high = "#4FA2DA") +
  labs(title = "Distribuição da idade",
       subtitle = sprintf("Média = %.1f ± %.1f | Mediana = %.1f",
                          media_idade, dp_idade, mediana_idade),
       x = "Idade (anos)", y = "Frequência") +
  guides(fill = "none")

colmap <- c("NRS RPA r"="SRPA_Repouso","NRS RPA m"="SRPA_Movimento",
            "NRS 24h r"="24h_Repouso", "NRS 24h m"="24h_Movimento",
            "NRS 48h r"="48h_Repouso", "NRS 48h m"="48h_Movimento")
present <- intersect(names(colmap), names(dados))

df_long <- dados %>%
  select(group, all_of(present)) %>%
  rename(!!!setNames(present, colmap[present])) %>%
  pivot_longer(-group, names_to = "key", values_to = "nrs") %>%
  separate(key, into = c("momento","cond"), sep = "_") %>%
  mutate(nrs = as.numeric(nrs),
         momento = factor(momento, levels = c("SRPA","24h","48h")),
         cond    = factor(cond,    levels = c("Repouso","Movimento"))) %>%
  filter(!is.na(nrs))

painel_test <- function(d) {
  ct <- table(d$group)
  if (length(ct) < 2 || any(ct < 2)) return(NA_real_)
  norm <- tapply(d$nrs, d$group, function(v) if (length(v) >= 3) shapiro.test(v)$p.value >= .05 else FALSE)
  if (all(unlist(norm))) t.test(nrs ~ group, d)$p.value
  else suppressWarnings(wilcox.test(nrs ~ group, d, exact = FALSE)$p.value)
}

tests <- df_long %>%
  group_by(momento, cond) %>%
  group_modify(~ tibble(p = painel_test(.x))) %>%
  ungroup() %>%
  mutate(p_label = ifelse(is.na(p), "",
                          ifelse(p < .001, "p < 0,001",
                                 paste0("p = ", nv(p, 3)))))
y_pos  <- df_long %>% group_by(momento, cond) %>%
  summarise(y_p = min(10, max(nrs) + 0.6), .groups = "drop")
tests  <- left_join(tests, y_pos, by = c("momento","cond"))

pd <- position_dodge(width = 0.8)
ggplot(df_long, aes(momento, nrs, fill = group)) +
  geom_boxplot(width = 0.65, position = pd, outlier.shape = 8, outlier.size = 2, alpha = 0.92) +
  stat_summary(aes(group = group), fun = mean, geom = "point",
               shape = 18, size = 3, color = "black", position = pd) +
  stat_summary(aes(label = sprintf("%.2f", after_stat(y)), group = group),
               fun = mean, geom = "text", fontface = "bold", vjust = -0.6, size = 3.2, position = pd) +
  geom_text(data = tests, aes(momento, y_p, label = p_label),
            fontface = "bold", size = 3.2, inherit.aes = FALSE) +
  facet_wrap(~ cond, nrow = 1) +
  scale_fill_manual(values = pal_grupo, name = "Grupo") +
  labs(title = "EVN de dor (0–10) em repouso e movimento — SRPA, 24h e 48h",
       x = "Tempo de avaliação", y = "EVN (0–10)") +
  coord_cartesian(ylim = c(0, 10)) +
  theme(legend.position = "right")

dados_long <- dados %>%
  transmute(group, SRPA = qor_srpa, `24h` = qor_24h, `48h` = qor_48h) %>%
  pivot_longer(c("SRPA","24h","48h"), names_to = "tempo", values_to = "qor15") %>%
  filter(!is.na(qor15)) %>%
  mutate(tempo = factor(tempo, levels = c("SRPA","24h","48h")))

sw <- dados_long %>% group_by(tempo, group) %>% shapiro_test(qor15)

testes <- dados_long %>%
  group_split(tempo) %>%
  lapply(function(d) {
    t_atual <- as.character(unique(d$tempo))
    sw_t    <- sw %>% filter(tempo == t_atual)
    n_min   <- min(table(d$group))
    if (any(sw_t$p < .05) || n_min < 8) {
      tt <- wilcox_test(d, qor15 ~ group); tt$test_used <- "Wilcoxon"
    } else {
      tt <- t_test(d, qor15 ~ group, var.equal = FALSE); tt$test_used <- "Welch"
    }
    tt$tempo <- t_atual; tt
  }) %>% bind_rows()

tests_lab <- testes %>%
  mutate(p_label = ifelse(p < .001, "p < 0,001", paste0("p = ", nv(p, 3))))
y_pos <- dados_long %>% group_by(tempo) %>%
  summarise(y_p = max(qor15, na.rm = TRUE) + 4, .groups = "drop")
tests_lab <- left_join(tests_lab, y_pos, by = "tempo") %>%
  mutate(tempo = factor(tempo, levels = c("SRPA","24h","48h")))

means_df <- dados_long %>% group_by(tempo, group) %>%
  summarise(y = mean(qor15, na.rm = TRUE), .groups = "drop")

pd <- position_dodge(width = 0.8)
ggplot(dados_long, aes(tempo, qor15, fill = group)) +
  geom_boxplot(width = 0.65, position = pd, outlier.shape = 8, outlier.size = 2, alpha = 0.90) +
  geom_point(data = means_df, aes(tempo, y, group = group),
             position = pd, shape = 18, size = 3, color = "black") +
  shadowtext::geom_shadowtext(
    data = means_df, aes(tempo, y, label = sprintf("%.1f", y), group = group),
    position = pd, fontface = "bold", size = 3.4, vjust = -0.6,
    colour = "black", bg.colour = "white", bg.r = 0.12) +
  geom_text(data = tests_lab, aes(tempo, y_p, label = p_label),
            fontface = "bold", size = 3.6, inherit.aes = FALSE) +
  scale_fill_manual(values = pal_grupo, name = "Grupo") +
  labs(title = "Qualidade de recuperação (QoR-15)", x = "Tempo de avaliação", y = "QoR-15") +
  coord_cartesian(clip = "off") +
  theme(legend.position = "right", plot.margin = margin(10, 20, 10, 10))

dl <- dados %>%
  transmute(id, group, SRPA = qor_srpa, `24h` = qor_24h, `48h` = qor_48h) %>%
  pivot_longer(c("SRPA","24h","48h"), names_to = "tempo", values_to = "qor15") %>%
  filter(!is.na(qor15), !is.na(group)) %>%
  mutate(tempo = factor(tempo, levels = c("SRPA","24h","48h")))

fit    <- lmer(qor15 ~ group * tempo + (1 | id), data = dl)
emm    <- emmeans(fit, ~ group | tempo)
emm_df <- as.data.frame(emm)

# pairwise devolve (BMRA − BMRA-CLO); inverter para (CLO − BMRA)
contr <- as.data.frame(confint(contrast(emm, method = "pairwise")))
contr$estimate <- -contr$estimate
tmp <- contr$lower.CL; contr$lower.CL <- -contr$upper.CL; contr$upper.CL <- -tmp

nv1  <- function(x) formatC(x, digits = 1, format = "f", decimal.mark = ",")
sgn  <- function(x) ifelse(x >= 0, paste0("+", nv1(x)), paste0("−", nv1(abs(x))))
delta_labs <- contr %>%
  transmute(tempo, delta_lab = sprintf("Δ=%s [%s; %s]", sgn(estimate), nv1(lower.CL), nv1(upper.CL)))

mcid_w <- 0.46
base_mcid <- emm_df %>% filter(group == "BMRA") %>%
  mutate(xn = as.numeric(tempo), xmin = xn - mcid_w, xmax = xn + mcid_w,
         ymin = emmean - 8, ymax = emmean + 8)

y_top <- emm_df %>% group_by(tempo) %>%
  summarise(y = max(upper.CL, na.rm = TRUE) + 3, .groups = "drop") %>%
  left_join(delta_labs, by = "tempo")
y_lo <- max(0, floor(min(c(emm_df$lower.CL, base_mcid$ymin), na.rm = TRUE) - 3))
y_hi <- ceiling(max(c(emm_df$upper.CL, base_mcid$ymax, y_top$y), na.rm = TRUE) + 2)

ggplot(emm_df, aes(tempo, emmean, color = group, group = group)) +
  geom_rect(data = base_mcid, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            inherit.aes = FALSE, fill = "grey60", alpha = 0.14) +
  geom_line(linewidth = 1) + geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .08, linewidth = 0.8) +
  geom_text(data = y_top, aes(tempo, y, label = delta_lab),
            inherit.aes = FALSE, size = 4, fontface = "bold") +
  scale_color_manual(values = pal_grupo, name = "Grupo") +
  labs(title = "QoR-15 — médias marginais estimadas (EMM), IC95% e faixa de MCID (±8)",
       x = "Tempo de avaliação", y = "QoR-15 (EMM)") +
  scale_y_continuous(limits = c(y_lo, y_hi), expand = expansion(mult = c(.02, .08))) +
  theme(legend.position = "right")

fit_ancova <- lm(fentanyl_per_kg ~ group + idade + imc + sexo + asa, data = dados)

# Tabela de coeficientes
coef_tab <- parameters::model_parameters(fit_ancova, ci = 0.95) %>%
  as.data.frame() %>%
  transmute(
    Termo       = Parameter,
    Coeficiente = nv(Coefficient, 3),
    `IC 95%`    = sprintf("[%s; %s]", nv(CI_low, 3), nv(CI_high, 3)),
    `Valor-p`   = ifelse(p < 0.001, "&lt;0,001", nv(p, 3))
  )
knitr::kable(coef_tab, align = c("l","c","c","c"),
             caption = "ANCOVA — coeficientes do modelo (fentanil µg/kg)")

# Diferença média ajustada (CLO − BMRA)
emm_anc       <- emmeans(fit_ancova, ~ group)
contr_anc_obj <- contrast(emm_anc, method = "pairwise")
diff_anc      <- as.data.frame(summary(contr_anc_obj, infer = c(TRUE, TRUE)))
diff_anc$estimate <- -diff_anc$estimate
tmp <- diff_anc$lower.CL; diff_anc$lower.CL <- -diff_anc$upper.CL; diff_anc$upper.CL <- -tmp

ajust <- tibble::tibble(
  Contraste                         = "BMRA-CLO − BMRA",
  `Diferença ajustada (µg/kg)`      = nv(diff_anc$estimate, 3),
  `IC 95%`                          = sprintf("[%s; %s]", nv(diff_anc$lower.CL, 3), nv(diff_anc$upper.CL, 3)),
  `Valor-p`                         = ifelse(diff_anc$p.value < 0.001, "&lt;0,001", nv(diff_anc$p.value, 3))
)
knitr::kable(ajust, align = c("l","c","c","c"),
             caption = "Diferença média ajustada por ANCOVA (com IC 95% e p)")

# η² parcial (tamanho de efeito)
effectsize::eta_squared(fit_ancova, partial = TRUE)

# Pressupostos
shapiro.test(residuals(fit_ancova))   # normalidade dos resíduos

car::ncvTest(fit_ancova)              # homocedasticidade (Breusch–Pagan)

# Forest plot dos coeficientes
fp <- parameters::model_parameters(fit_ancova, ci = 0.95) %>%
  as.data.frame() %>%
  filter(Parameter != "(Intercept)") %>%
  mutate(Parameter = factor(Parameter, levels = rev(Parameter)))

ggplot(fp, aes(Coefficient, Parameter)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "#7f8c8d") +
  geom_pointrange(aes(xmin = CI_low, xmax = CI_high), size = 0.6, color = cor_controle) +
  labs(title = "ANCOVA — coeficientes (IC 95%) para fentanil (µg/kg)",
       x = "Coeficiente", y = NULL)

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

coef_tab2 <- parameters::model_parameters(fit_qor48, ci = 0.95) %>%
  as.data.frame() %>%
  transmute(
    Termo       = Parameter,
    Coeficiente = nv(Coefficient, 2),
    `IC 95%`    = sprintf("[%s; %s]", nv(CI_low, 2), nv(CI_high, 2)),
    `Valor-p`   = ifelse(p < 0.001, "&lt;0,001", nv(p, 3))
  )
knitr::kable(coef_tab2, align = c("l","c","c","c"),
             caption = "Regressão linear múltipla — QoR-15 em 48h")

# R² ajustado e diagnósticos
cat(sprintf("R² ajustado: %s\n", nv(summary(fit_qor48)$adj.r.squared, 3)))

shapiro.test(residuals(fit_qor48))

car::ncvTest(fit_qor48)

de_mov <- dados %>%
  select(id, group, `NRS RPA m`, `NRS 24h m`, `NRS 48h m`) %>%
  pivot_longer(-c(id, group), names_to = "tempo", values_to = "nrs") %>%
  mutate(tempo = factor(dplyr::recode(tempo,
                                      "NRS RPA m" = "SRPA",
                                      "NRS 24h m" = "24h",
                                      "NRS 48h m" = "48h"),
                        levels = c("SRPA","24h","48h"))) %>%
  filter(!is.na(nrs))

fit_lmm_evn <- lmer(nrs ~ group * tempo + (1 | id), data = de_mov)

# Tabela ANOVA tipo III com Satterthwaite
anova_tab <- as.data.frame(anova(fit_lmm_evn)) %>%
  tibble::rownames_to_column("Termo") %>%
  transmute(
    Termo,
    `F`       = nv(`F value`, 2),
    `GL num`  = nv(NumDF, 1),
    `GL den`  = nv(DenDF, 1),
    `Valor-p` = ifelse(`Pr(>F)` < 0.001, "&lt;0,001", nv(`Pr(>F)`, 3))
  )
knitr::kable(anova_tab, align = c("l","c","c","c","c"),
             caption = "LMM — ANOVA tipo III (Satterthwaite) para EVN em movimento")

# EMMs por tempo e contraste BMRA-CLO − BMRA em cada momento
emm_evn       <- emmeans(fit_lmm_evn, ~ group | tempo)
contr_evn_obj <- contrast(emm_evn, method = "pairwise")
contr_ev      <- as.data.frame(summary(contr_evn_obj, infer = c(TRUE, TRUE)))
contr_ev$estimate <- -contr_ev$estimate
tmp <- contr_ev$lower.CL; contr_ev$lower.CL <- -contr_ev$upper.CL; contr_ev$upper.CL <- -tmp

contr_tab <- contr_ev %>%
  transmute(
    Tempo                    = tempo,
    `Δ ajustada (CLO − BMRA)` = nv(estimate, 2),
    `IC 95%`                 = sprintf("[%s; %s]", nv(lower.CL, 2), nv(upper.CL, 2)),
    `Valor-p`                = ifelse(p.value < 0.001, "&lt;0,001", nv(p.value, 3))
  )
knitr::kable(contr_tab, align = c("l","c","c","c"),
             caption = "LMM — diferenças entre grupos em cada tempo (EVN movimento)")