Cadeia Leve Livre Sérica e Mortalidade por Todas as Causas

Trabalho Final — Análise de Sobrevivência

Autor

Arthur Pontes Motta

Data de Publicação

29 de junho de 2026

Carregar pacotes

# --- Análise de sobrevivência ---
library(survival)    # estimadores KM, Nelson-Aalen, modelo de Cox, dados flchain
library(flexsurv)    # modelos paramétricos AFT: Weibull, log-normal, gama generalizada, etc.
library(muhaz)       # estimação do hazard por kernel de Nadaraya-Watson
library(bshazard)    # estimação do hazard por B-splines (preferível por ausência de efeito de borda)

# --- Seleção de variáveis ---
library(MASS)           # stepAIC: seleção de covariáveis por critério AIC
library(My.stepwise)    # My.stepwise.coxph: seleção bidirecional com p-valores para o Cox

# --- Visualização de curvas de sobrevivência ---
library(ggsurvfit)   # curvas KM e Nelson-Aalen com ggplot2 (add_risktable, add_pvalue, etc.)
library(survminer)   # diagnósticos do Cox (ggcoxzph, ggcoxdiagnostics, ggforest)

# --- Tabelas e relatório ---
library(gtsummary)   # tabelas descritivas (tbl_summary) e de modelos (tbl_regression)
library(gt)          # tabelas gt com formatação avançada (LaTeX, cores, footnotes)

# --- Manipulação e visualização geral ---
library(tidyverse)   # dplyr, ggplot2, tidyr, purrr, stringr, forcats, readr
library(patchwork)   # composição de múltiplos gráficos ggplot2 em painéis
library(broom)       # tidy(), glance(), augment() para extrair resultados de modelos
library(corrplot)    # matriz de correlação visual (corrplot)
library(scales)      # formatação de eixos (percent_format, label_number, etc.)

theme_set(theme_minimal(base_size = 12))
theme_gtsummary_compact()

Introdução

A cadeia leve livre sérica (serum free light chain, FLC) é uma proteína produzida pelas células plasmáticas do sistema imunológico. Em indivíduos saudáveis, os níveis de FLC são mantidos dentro de uma faixa estreita. Elevações anormais, mesmo na ausência de doença hematológica diagnosticada, têm sido associadas a piores desfechos clínicos, incluindo maior mortalidade (Dispenzieri et al., 2012).

Este trabalho analisa dados de um estudo de coorte de base populacional conduzido no Condado de Olmsted, Minnesota (EUA), disponíveis no pacote survival do R sob o nome flchain. O estudo recrutou aproximadamente dois terços dos residentes do condado com 50 anos ou mais a partir de 1995, com o objetivo de determinar a prevalência da gamopatia monoclonal de significado indeterminado (MGUS) e sua relação com a mortalidade.

Pergunta principal de interesse: A concentração sérica de cadeia leve livre está associada ao tempo de sobrevivência de adultos com 50 anos ou mais, após controle por idade, sexo, creatinina sérica e diagnóstico de MGUS?

O conjunto de dados é adequado para análise de sobrevivência porque contém uma variável de tempo de acompanhamento (futime, em dias), um indicador de evento (death), mecanismo de censura à direita (sobreviventes ao final do seguimento foram censurados) e múltiplas covariáveis clínicas e demográficas que permitem modelagem multivariada.

Formalmente, define-se a função de sobrevivência como \[S(t) = P(T > t), \quad t \geq 0,\] onde \(T\) é o tempo até o evento (óbito). A função de risco instantâneo (hazard function) é \[h(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t + \Delta t \mid T \geq t)}{\Delta t} = -\frac{d}{dt}\log S(t),\] e a relação entre as duas é \[S(t) = \exp\!\left(-\int_0^t h(u)\,du\right) = e^{-H(t)},\] onde \(H(t) = \int_0^t h(u)\,du\) é o risco acumulado (Colosimo; Giolo, 2006; Klein; Moeschberger, 2003).

Descrição dos Dados e Análise Exploratória

Fonte e Variáveis

Os dados provêm de uma amostra estratificada aleatória de 50% dos participantes do estudo original (Kyle et al., 2006). O conjunto contém 7.874 indivíduos e as seguintes variáveis:

Variável	Tipo	Descrição
`futime`	Tempo	Dias desde a coleta até o óbito ou censura
`death`	Evento	1 = óbito, 0 = censurado
`age`	Contínua	Idade em anos no momento da coleta
`sex`	Categórica	F = feminino, M = masculino
`kappa`	Contínua	FLC kappa sérica (mg/dL)
`lambda`	Contínua	FLC lambda sérica (mg/dL)
`flc.grp`	Ordinal	Grupo de FLC total (1–10), usado na análise original
`creatinine`	Contínua	Creatinina sérica (mg/dL)
`mgus`	Binária	1 = diagnóstico prévio de MGUS

A variável flc.grp representa decis do valor de FLC total (kappa + lambda) e será utilizada como covariável categórica principal, conforme a análise original de Dispenzieri et al. (2012). Agrupamos os decis em quatro categorias para facilitar a interpretação: baixo (1–2), médio-baixo (3–5), médio-alto (6–8) e alto (9–10).

Preparação dos dados

data(flchain, package = "survival")

# Remover os 3 indivíduos com futime = 0
# (coleta realizada no dia do óbito — tempo zero impossível no KM)
flchain_clean <- flchain |>
  filter(futime > 0) |>
  mutate(
    # Tempo em anos (facilita interpretação)
    anos       = futime / 365.25,
    # Grupo de FLC em 4 categorias (labels sem travessão especial para evitar encoding)
    flc_cat    = case_when(
      flc.grp %in% 1:2   ~ "1-2 (Baixo)",
      flc.grp %in% 3:5   ~ "3-5 (Medio-baixo)",
      flc.grp %in% 6:8   ~ "6-8 (Medio-alto)",
      flc.grp %in% 9:10  ~ "9-10 (Alto)"
    ) |> factor(levels = c("1-2 (Baixo)", "3-5 (Medio-baixo)",
                            "6-8 (Medio-alto)", "9-10 (Alto)")),
    # Grupo de idade
    age_cat    = cut(age,
                     breaks = c(49, 59, 69, 79, 110),
                     labels = c("50–59", "60–69", "70–79", "80+"),
                     right  = TRUE),
    # Sexo como fator legível
    sexo       = factor(sex, levels = c("F","M"),
                        labels = c("Feminino","Masculino")),
    # MGUS como fator
    mgus_f     = factor(mgus, levels = c(0,1),
                        labels = c("Não","Sim")),
    # FLC total
    flc_total  = kappa + lambda
  )

cat(sprintf("Amostra final: %d indivíduos\n", nrow(flchain_clean)))

Amostra final: 7871 indivíduos

Preparação dos dados

cat(sprintf("Óbitos: %d (%.1f%%)\n",
            sum(flchain_clean$death),
            100 * mean(flchain_clean$death)))

Óbitos: 2166 (27.5%)

Preparação dos dados

cat(sprintf("Censurados: %d (%.1f%%)\n",
            sum(flchain_clean$death == 0),
            100 * mean(flchain_clean$death == 0)))

Censurados: 5705 (72.5%)

Preparação dos dados

cat(sprintf("Seguimento mediano: %.1f anos (máx: %.1f)\n",
            median(flchain_clean$anos),
            max(flchain_clean$anos)))

Seguimento mediano: 11.8 anos (máx: 14.3)

Análise Descritiva

Ocultar código

flchain_clean |>
  select(anos, sexo, age, age_cat, flc_cat, creatinine, mgus_f, death) |>
  tbl_summary(
    by      = death,
    label   = list(
      anos       ~ "Tempo de seguimento (anos)",
      sexo       ~ "Sexo",
      age        ~ "Idade (anos)",
      age_cat    ~ "Faixa etária",
      flc_cat    ~ "Grupo de FLC",
      creatinine ~ "Creatinina (mg/dL)",
      mgus_f     ~ "MGUS"
    ),
    statistic = list(
      all_continuous()  ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits    = list(all_continuous() ~ 1),
    missing   = "no"
  ) |>
  add_overall() |>
  add_p() |>
  modify_header(
    label     ~ "**Variável**",
    stat_0    ~ "**Total**\nN = {N}",
    stat_1    ~ "**Censurado**\nn = {n}",
    stat_2    ~ "**Óbito**\nn = {n}"
  ) |>
  modify_caption("**Tabela 1.** Características da amostra por desfecho [@kyle2006prevalence]") |>
  bold_labels()

Características da amostra por desfecho
Variável	Total N = 7871¹	Censurado n = 5705¹	Óbito n = 2166¹	p-value²
Tempo de seguimento (anos)	11.8 (7.8, 13.1)	12.6 (11.0, 13.3)	5.9 (2.6, 9.0)	<0.001
Sexo				0.082
Feminino	4,347 (55%)	3,185 (56%)	1,162 (54%)
Masculino	3,524 (45%)	2,520 (44%)	1,004 (46%)
Idade (anos)	63.0 (55.0, 72.0)	59.0 (54.0, 66.0)	74.0 (67.0, 81.0)	<0.001
Faixa etária				<0.001
50–59	3,157 (40%)	2,899 (51%)	258 (12%)
60–69	2,329 (30%)	1,848 (32%)	481 (22%)
70–79	1,623 (21%)	831 (15%)	792 (37%)
80+	762 (9.7%)	127 (2.2%)	635 (29%)
Grupo de FLC				<0.001
1-2 (Baixo)	1,580 (20%)	1,344 (24%)	236 (11%)
3-5 (Medio-baixo)	2,397 (30%)	1,945 (34%)	452 (21%)
6-8 (Medio-alto)	2,327 (30%)	1,651 (29%)	676 (31%)
9-10 (Alto)	1,567 (20%)	765 (13%)	802 (37%)
Creatinina (mg/dL)	1.0 (0.9, 1.2)	1.0 (0.9, 1.2)	1.1 (0.9, 1.3)	<0.001
MGUS				<0.001
Não	7,756 (99%)	5,606 (98%)	2,150 (99%)
Sim	115 (1.5%)	99 (1.7%)	16 (0.7%)
¹ Median (Q1, Q3); n (%)
² Wilcoxon rank sum test; Pearson’s Chi-squared test

A Tabela 1 revela diferenças substanciais e clinicamente coerentes entre os grupos. Os indivíduos que foram a óbito eram substancialmente mais velhos (mediana de 75 anos vs. 62 anos nos censurados), diferença de 13 anos que, por si só, já seria suficiente para explicar grande parte da divergência no desfecho, dada a forte dependência etária da mortalidade. Essa diferença etária entre grupos aponta para a necessidade de ajuste multivariado: análises brutas que ignorassem a idade estariam confundindo o efeito da FLC com o efeito do envelhecimento.

A distribuição por grupo de FLC é marcadamente assimétrica entre os desfechos: enquanto 24% dos censurados pertenciam ao grupo de FLC baixo, apenas 11% dos que foram a óbito estavam nesse grupo. No extremo oposto, 37% dos óbitos pertenciam ao grupo de FLC alto, contra 13% dos censurados, razão de quase 3:1. Essa separação já na análise descritiva univariada antecipa a associação que será confirmada formalmente nos modelos.

A creatinina sérica, marcador de função renal, também difere entre os grupos (mediana 1,1 vs. 1,0 mg/dL), com distribuição consistente com o esperado, pois disfunção renal é reconhecidamente associada a maior mortalidade cardiovascular e por todas as causas. É importante notar que creatinina e FLC estão correlacionadas: rins com menor capacidade de filtração eliminam menos FLC, elevando seus níveis séricos. Esse mecanismo de confundimento reforça a necessidade de ajuste conjunto pelas duas variáveis na modelagem.

A proporção de MGUS é baixa (1,5%), mas estatisticamente diferente entre os grupos (\(p < 0{,}001\)). Surpreendentemente, a proporção de MGUS é ligeiramente maior entre os censurados (1,7%) do que entre os que foram a óbito (0,7%). Isso pode parecer contraintuitivo, dado que MGUS é uma condição pré-maligna; contudo, pode refletir viés de sobrevivência (indivíduos com MGUS diagnosticado podem receber acompanhamento médico mais intensivo), o pequeno tamanho do subgrupo, ou o fato de que, no modelo multivariado ajustado por FLC, o efeito independente do MGUS não é significativo.

Distribuição do Tempo de Sobrevivência

Ocultar código

flchain_clean |>
  mutate(Desfecho = factor(death, labels = c("Censurado", "Óbito"))) |>
  ggplot(aes(x = anos, fill = Desfecho)) +
  geom_histogram(bins = 50, alpha = 0.7, color = "white", linewidth = 0.2) +
  facet_wrap(~Desfecho, scales = "free_y") +
  scale_fill_manual(values = c("#2E86AB", "#E94F37")) +
  labs(
    title = "Distribuição do Tempo de Seguimento",
    x     = "Tempo (anos)",
    y     = "Frequência"
  ) +
  theme(legend.position = "none")

Distribuição do tempo de seguimento por desfecho

Os histogramas revelam perfis completamente distintos entre censurados e óbitos. Os censurados concentram-se em tempos longos, com acúmulo visível entre 11 e 14 anos, correspondendo ao final do período de acompanhamento, quando o estudo foi encerrado e os sobreviventes foram censurados administrativamente. Esse padrão é consistente com censura à direita não informativa: os indivíduos foram censurados por fim do estudo, não por terem deixado de ser observáveis por razões relacionadas ao risco.

Os óbitos, por sua vez, distribuem-se de forma mais uniforme ao longo de todo o período, com leve concentração nos primeiros anos e decrescimento gradual à medida que a coorte inicial vai sendo esgotada. A ausência de um pico muito pronunciado nos primeiros meses descarta a presença de evento agudo inicial importante, como ocorreria em coortes pós-hospitalização, confirmando que a coorte de base populacional tem perfil de risco distribuído.

Distribuição das Variáveis Contínuas

FLC Total, Kappa, Lambda e Creatinina

A distribuição das variáveis de laboratório é fortemente assimétrica à direita, característica típica de biomarcadores séricos cuja distribuição populacional é limitada por zero e apresenta cauda superior longa devido a indivíduos com produção aumentada de proteínas. Kappa e lambda apresentam outliers extremos (valores de até 20 mg/dL, contra mediana de aproximadamente 1,3 mg/dL), correspondendo a casos de MGUS ou condições inflamatórias severas. A escala logarítmica aproxima as distribuições da normalidade em todas as variáveis, o que justifica o uso de log(creatinina) na modelagem e a interpretação dos coeficientes como efeitos multiplicativos sobre a taxa de risco (Colosimo; Giolo, 2006).

Ocultar código

flchain_clean |>
  mutate(flc_total = kappa + lambda) |>
  select(kappa, lambda, flc_total, creatinine) |>
  pivot_longer(everything(), names_to = "variavel", values_to = "valor") |>
  filter(!is.na(valor)) |>
  mutate(variavel = factor(variavel,
    levels = c("kappa","lambda","flc_total","creatinine"),
    labels = c("Kappa","Lambda","FLC Total","Creatinina"))) |>
  ggplot(aes(x = valor, fill = variavel)) +
  geom_histogram(bins = 60, color = "white", linewidth = 0.15, alpha = 0.85) +
  facet_wrap(~variavel, scales = "free", ncol = 4) +
  scale_fill_manual(values = c("#2E86AB","#3BB273","#F4A261","#E94F37")) +
  scale_x_continuous(trans = "log1p",
                     labels = scales::label_number(accuracy = 0.1)) +
  labs(
    title = "Distribuição das Variáveis Laboratoriais (escala log+1)",
    x     = "Valor (escala log+1)",
    y     = "Frequência"
  ) +
  theme(legend.position = "none")

Distribuição das variáveis laboratoriais (escala log+1)

Distribuição da FLC por Grupo e por Desfecho

Ocultar código

p1 <- flchain_clean |>
  ggplot(aes(x = flc_cat, y = log(flc_total), fill = flc_cat)) +
  geom_boxplot(alpha = 0.75, outlier.size = 0.5, outlier.alpha = 0.3) +
  scale_fill_manual(values = c("#3BB273","#2E86AB","#F4A261","#E94F37")) +
  labs(title = "FLC Total por Grupo", x = NULL, y = "log(FLC Total)") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 20, hjust = 1))

# Taxa de óbito por grupo (bruta, sem ajuste)
taxa_obito <- flchain_clean |>
  group_by(flc_cat) |>
  summarise(
    n          = n(),
    obitos     = sum(death),
    taxa       = obitos / n,
    ic_inf     = taxa - 1.96 * sqrt(taxa * (1 - taxa) / n),
    ic_sup     = taxa + 1.96 * sqrt(taxa * (1 - taxa) / n)
  )

p2 <- taxa_obito |>
  ggplot(aes(x = flc_cat, y = taxa, fill = flc_cat)) +
  geom_col(alpha = 0.85, width = 0.6) +
  geom_errorbar(aes(ymin = ic_inf, ymax = ic_sup), width = 0.2, linewidth = 0.7) +
  scale_fill_manual(values = c("#3BB273","#2E86AB","#F4A261","#E94F37")) +
  scale_y_continuous(labels = scales::percent_format(1), limits = c(0, 0.7)) +
  labs(title = "Taxa Bruta de Óbito por Grupo", x = NULL, y = "% óbitos") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 20, hjust = 1))

p1 + p2

FLC total (log) por grupo e desfecho — a proporção de óbitos cresce com o grupo

O painel esquerdo confirma que os grupos de flc_cat capturam de fato níveis progressivamente crescentes de FLC total em escala logarítmica, com medianas bem separadas e variabilidade intragrupo relativamente homogênea, indicando que a categorização em decis foi eficiente na criação de grupos internamente coesos.

O painel direito é o resultado mais impactante da análise exploratória: a taxa bruta de óbito sobe de forma monotônica e pronunciada do grupo Baixo (aproximadamente 15%) para o Alto (aproximadamente 51%), com intervalos de confiança que não se sobrepõem entre grupos adjacentes. Isso estabelece uma relação dose a resposta já na análise univariada, antes de qualquer ajuste por confundidores. Essa relação monotônica é importante porque sugere que o efeito da FLC não é limiar, ou seja, presente apenas acima de um valor crítico, mas gradual. Essa característica tem implicações clínicas diretas: mesmo elevações moderadas de FLC conferem risco adicional mensurável.

Análise de Valores Ausentes e Perfil de Recrutamento

Ocultar código

# Recrutamento por ano
p_ano <- flchain_clean |>
  count(sample.yr) |>
  ggplot(aes(x = factor(sample.yr), y = n)) +
  geom_col(fill = "#2E86AB", alpha = 0.85, width = 0.7) +
  geom_text(aes(label = n), vjust = -0.4, size = 3.2) +
  labs(
    title = "Recrutamento por Ano de Coleta",
    x     = "Ano",
    y     = "N de indivíduos"
  )

# Missing em creatinina por grupo de FLC
p_miss <- flchain_clean |>
  mutate(creat_miss = is.na(creatinine)) |>
  group_by(flc_cat) |>
  summarise(pct_miss = mean(creat_miss)) |>
  ggplot(aes(x = flc_cat, y = pct_miss, fill = flc_cat)) +
  geom_col(alpha = 0.85, width = 0.6) +
  scale_fill_manual(values = c("#3BB273","#2E86AB","#F4A261","#E94F37")) +
  scale_y_continuous(labels = scales::percent_format(1)) +
  labs(
    title = "% de Creatinina Ausente por Grupo de FLC",
    x     = NULL,
    y     = "% missing"
  ) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 20, hjust = 1))

p_ano + p_miss

Perfil de recrutamento por ano e valores ausentes em creatinina

Os valores ausentes de creatinina concentram-se sistematicamente nos grupos de FLC mais baixo, chegando a cerca de 25 a 30% de ausência no grupo Baixo, contra valores menores nos grupos mais elevados. Esse padrão não é aleatório e tem implicação direta na análise: ao excluir os aproximadamente 1.350 indivíduos sem creatinina, estamos removendo preferencialmente pessoas de menor risco. Isso significa que a amostra modelada (n = 6.521) sobrerrepresenta indivíduos de maior risco relativo, o que pode levar a uma subestimação da sobrevivência média e a uma superestimação do efeito dos grupos de FLC mais elevados. A distribuição temporal do recrutamento, com 80% nos três primeiros anos, é esperada para estudos de coorte que recrutam de forma intensiva no início e dependem de visitas ambulatoriais subsequentes para os demais.

Relação entre Covariáveis

Razão Kappa/Lambda e MGUS

O diagnóstico de MGUS é clinicamente associado a uma razão kappa/lambda fora do intervalo de referência (0,26–1,65, segundo a Mayo Clinic). O gráfico confirma esse padrão nos dados.

Ocultar código

flchain_clean |>
  mutate(
    flc_ratio = kappa / lambda,
    mgus_lab  = factor(mgus, labels = c("Sem MGUS","Com MGUS"))
  ) |>
  ggplot(aes(x = mgus_lab, y = log(flc_ratio), fill = mgus_lab)) +
  geom_violin(alpha = 0.5, draw_quantiles = c(0.25, 0.5, 0.75)) +
  geom_hline(yintercept = log(c(0.26, 1.65)),
             linetype = "dashed", color = "gray30", linewidth = 0.7) +
  annotate("text", x = 2.45, y = log(0.26) + 0.08,
           label = "Ref. inf. (0,26)", size = 3, color = "gray30") +
  annotate("text", x = 2.45, y = log(1.65) + 0.08,
           label = "Ref. sup. (1,65)", size = 3, color = "gray30") +
  scale_fill_manual(values = c("#2E86AB","#E94F37")) +
  labs(
    title    = "Razão Kappa/Lambda por Status de MGUS",
    subtitle = "Linhas tracejadas = intervalo de referência clínica (Mayo Clinic: 0,26–1,65)",
    x        = NULL,
    y        = "log(Kappa / Lambda)"
  ) +
  theme(legend.position = "none")

Razão kappa/lambda (log) por status de MGUS com linhas de referência clínica

Correlação entre Variáveis Contínuas

Ocultar código

library(corrplot)

flchain_clean |>
  mutate(
    flc_total = kappa + lambda,
    flc_ratio = kappa / lambda,
    log_creat = log(creatinine)
  ) |>
  select(age, kappa, lambda, flc_total, flc_ratio, log_creat, anos) |>
  filter(complete.cases(pick(everything()))) |>
  cor(method = "spearman") |>
  corrplot(
    method      = "color",
    type        = "upper",
    tl.col      = "black",
    tl.cex      = 0.85,
    addCoef.col = "black",
    number.cex  = 0.72,
    col         = colorRampPalette(c("#E94F37","white","#2E86AB"))(200),
    title       = "Correlação de Spearman — Variáveis Contínuas",
    mar         = c(0, 0, 1.5, 0)
  )

Matriz de correlação entre variáveis contínuas (Spearman)

Idade e FLC por Sexo

Ocultar código

flchain_clean |>
  mutate(flc_total = kappa + lambda) |>
  ggplot(aes(x = age, y = log(flc_total), color = sexo)) +
  geom_point(alpha = 0.07, size = 0.6) +
  geom_smooth(method = "loess", se = TRUE, linewidth = 1.2) +
  scale_color_manual(values = c("#E94F37","#2E86AB")) +
  scale_x_continuous(breaks = seq(50, 100, by = 10)) +
  labs(
    title    = "FLC Total (log) vs. Idade por Sexo",
    subtitle = "Linha suavizada LOESS com IC 95% — homens e mulheres têm trajetórias similares",
    x        = "Idade (anos)",
    y        = "log(FLC Total)",
    color    = "Sexo"
  ) +
  theme(legend.position = "bottom")

Creatinina por Sexo e Faixa Etária

Ocultar código

flchain_clean |>
  filter(!is.na(creatinine)) |>
  ggplot(aes(x = age_cat, y = log(creatinine), fill = sexo)) +
  geom_boxplot(alpha = 0.75, outlier.size = 0.5, outlier.alpha = 0.3,
               position = position_dodge(0.8)) +
  scale_fill_manual(values = c("#E94F37","#2E86AB")) +
  labs(
    title    = "Creatinina (log) por Faixa Etária e Sexo",
    subtitle = "Homens apresentam creatinina consistentemente maior em todas as faixas",
    x        = "Faixa etária",
    y        = "log(Creatinina)",
    fill     = "Sexo"
  ) +
  theme(legend.position = "bottom")

Creatinina (log) por sexo e faixa etária

Homens apresentam creatinina consistentemente maior que mulheres em todas as faixas etárias, o que reflete a maior massa muscular masculina, pois a creatinina sérica é produto da degradação da creatina muscular. Esse dimorfismo sexual na creatinina é relevante para a modelagem: ao ajustar por log(creatinina) sem interação com sexo, assume-se que o efeito da creatinina sobre a mortalidade é o mesmo para ambos os sexos, hipótese plausível biologicamente, mas que poderia ser testada formalmente. Adicionalmente, a creatinina aumenta com a idade em ambos os sexos, refletindo o declínio da taxa de filtração glomerular com o envelhecimento.

Causas de Óbito

Ocultar código

# Simplificar categorias de chapter (mesmo critério do repositório de referência)
capitulos_principais <- c("Circulatory","Neoplasms","Respiratory","Mental","Nervous")

flchain_clean |>
  filter(death == 1, !is.na(chapter)) |>
  mutate(
    causa = if_else(chapter %in% capitulos_principais,
                    as.character(chapter), "Outras"),
    causa = factor(causa, levels = c(capitulos_principais, "Outras")),
    causa_pt = fct_recode(causa,
      "Circulatório"  = "Circulatory",
      "Neoplasias"    = "Neoplasms",
      "Respiratório"  = "Respiratory",
      "Mental"        = "Mental",
      "Nervoso"       = "Nervous",
      "Outras"        = "Outras"
    )
  ) |>
  count(flc_cat, causa_pt) |>
  group_by(flc_cat) |>
  mutate(pct = n / sum(n)) |>
  ggplot(aes(x = flc_cat, y = pct, fill = causa_pt)) +
  geom_col(position = "stack", alpha = 0.9, width = 0.7) +
  scale_fill_brewer(palette = "Set2") +
  scale_y_continuous(labels = scales::percent_format(1)) +
  labs(
    title    = "Causas de Óbito por Grupo de FLC",
    subtitle = "Proporção relativa entre os que foram a óbito em cada grupo",
    x        = "Grupo de FLC",
    y        = "Proporção",
    fill     = "Causa (CID-9)"
  ) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 15, hjust = 1))

Distribuição das causas de óbito por grupo de FLC

A composição de causas de óbito por grupo de FLC traz uma perspectiva importante para a interpretação dos resultados. Em todos os grupos, as doenças circulatórias (cardiovasculares) são a causa mais frequente de morte, seguidas pelas neoplasias, padrão esperado para uma coorte de adultos com 50 anos ou mais nos EUA. Notavelmente, a proporção relativa de neoplasias tende a ser ligeiramente maior nos grupos de FLC mais elevado, o que faz sentido biologicamente: FLC elevada é marcador de desregulação imunológica que pode preceder condições malignas hematológicas. A proporção de causas respiratórias e neurológicas é relativamente estável entre os grupos.

Esse padrão tem implicação metodológica importante: o desfecho analisado é mortalidade por todas as causas, e a FLC parece estar associada a múltiplas causas de morte e não apenas a uma causa específica. Isso sugere que a FLC pode ser um marcador de saúde geral e de envelhecimento biológico acelerado, e não apenas um marcador de risco para uma doença específica.

Curvas de Kaplan-Meier

Curva Global

O estimador de Kaplan-Meier (Therneau, 2024) da função de sobrevivência é definido como:

\[\hat{S}(t) = \prod_{j:\,t_j \leq t} \left(1 - \frac{d_j}{n_j}\right),\]

onde o produto é sobre todos os tempos de evento \(t_j \leq t\), com \(d_j\) óbitos e \(n_j\) indivíduos em risco. O estimador é não-paramétrico e lida naturalmente com censura à direita.

Ocultar código

survfit2(Surv(anos, death) ~ 1, data = flchain_clean) |>
  ggsurvfit(linewidth = 1.1, color = "#2E86AB") +
  add_confidence_interval(alpha = 0.15, fill = "#2E86AB") +
  add_risktable(
    risktable_stats = c("n.risk", "cum.event"),
    stats_label     = c("Em risco", "Óbitos acum.")
  ) +
  add_quantile(y_value = 0.5, linetype = "dashed", color = "gray50") +
  scale_ggsurvfit() +
  labs(
    title    = "Curva de Sobrevivência Global",
    subtitle = "Condado de Olmsted, Minnesota — Adultos ≥ 50 anos",
    x        = "Tempo (anos)",
    y        = expression(hat(S)(t))
  )

Curva de Kaplan-Meier global com IC 95% de Hall-Wellner

A curva de Kaplan-Meier (Therneau, 2024) global exibe o comportamento típico de uma coorte de base populacional de adultos: queda lenta nos primeiros anos, pois os indivíduos mais jovens da coorte (com 50 a 59 anos) dominam numericamente e têm baixo risco de mortalidade a curto prazo, e aceleração da queda após 8 a 10 anos, quando os indivíduos mais velhos que ainda estavam sob risco passam a contribuir mais fortemente. A tabela de risco mostra que a amostra de 7.871 indivíduos vai se esgotando progressivamente, o que naturalmente alarga o intervalo de confiança nas caudas da curva. Como mais de 50% da amostra sobreviveu ao seguimento, a mediana global não é estimável, dado biologicamente plausível, pois a coorte inclui indivíduos com 50 anos no início do estudo, muitos dos quais têm expectativa de vida superior a 14 anos.

Por Grupo de FLC

Ocultar código

survfit2(Surv(anos, death) ~ flc_cat, data = flchain_clean) |>
  ggsurvfit(linewidth = 1.1) +
  add_confidence_interval(alpha = 0.10) +
  add_risktable(risktable_stats = "n.risk", stats_label = "Em risco") +
  add_pvalue(caption = "Log-rank: {p.value}") +
  scale_color_manual(values = c("#3BB273","#2E86AB","#F4A261","#E94F37")) +
  scale_fill_manual(values  = c("#3BB273","#2E86AB","#F4A261","#E94F37")) +
  scale_ggsurvfit() +
  labs(
    title = "Sobrevivência por Grupo de FLC",
    x     = "Tempo (anos)",
    y     = expression(hat(S)(t))
  ) +
  theme(legend.position = "bottom")

Curvas de Kaplan-Meier por grupo de FLC (decis agrupados)

As curvas de Kaplan-Meier estratificadas por grupo de FLC constituem o resultado não paramétrico central deste trabalho. Três aspectos merecem destaque.

Primeiro, a separação entre os grupos é visível já no primeiro ano de seguimento e se aprofunda progressivamente ao longo de todo o período de 14 anos. Esse padrão de separação contínua, em oposição ao cruzamento de curvas, é consistente com a suposição de riscos proporcionais do modelo de Cox, embora o teste formal de Schoenfeld, discutido adiante, revele evidência estatística de alguma variação temporal.

Segundo, a magnitude clínica da diferença é expressiva. Ao final de 11 anos, a probabilidade de sobrevivência estimada no grupo Baixo é aproximadamente 87%, enquanto no grupo Alto é 49%. Essa diferença de 38 pontos percentuais em uma variável que pode ser medida em um único exame de sangue tem implicação clínica direta: a FLC poderia ser incorporada como marcador prognóstico em avaliações geriátricas de rotina.

Terceiro, o fato de que os grupos Baixo, Médio-baixo e Médio-alto não atingem \(\hat{S}(t) = 0{,}50\) durante o acompanhamento não é uma limitação metodológica, mas sim um resultado positivo: mais da metade dos indivíduos nesses grupos sobreviveram ao período completo de 14 anos. Isso reforça que a população de referência (FLC baixa) tem prognóstico relativamente favorável.

Por Sexo e por Faixa Etária

Ocultar código

p_sex <- survfit2(Surv(anos, death) ~ sexo, data = flchain_clean) |>
  ggsurvfit(linewidth = 1.1) +
  add_confidence_interval(alpha = 0.12) +
  add_pvalue(caption = "Log-rank: {p.value}") +
  scale_color_manual(values = c("#E94F37","#2E86AB")) +
  scale_fill_manual(values  = c("#E94F37","#2E86AB")) +
  scale_ggsurvfit() +
  labs(title = "Por Sexo", x = "Tempo (anos)", y = expression(hat(S)(t))) +
  theme(legend.position = "bottom")

p_age <- survfit2(Surv(anos, death) ~ age_cat, data = flchain_clean) |>
  ggsurvfit(linewidth = 1.1) +
  add_confidence_interval(alpha = 0.10) +
  add_pvalue(caption = "Log-rank: {p.value}") +
  scale_color_manual(values = c("#3BB273","#2E86AB","#F4A261","#E94F37")) +
  scale_fill_manual(values  = c("#3BB273","#2E86AB","#F4A261","#E94F37")) +
  scale_ggsurvfit() +
  labs(title = "Por Faixa Etária", x = "Tempo (anos)", y = "") +
  theme(legend.position = "bottom")

p_sex + p_age

Curvas de Kaplan-Meier por sexo (esq.) e faixa etária (dir.)

O gráfico por sexo confirma a sobremortalidade masculina bem estabelecida na literatura epidemiológica. As curvas se separam progressivamente ao longo do seguimento, sugerindo que o risco relativo dos homens aumenta com o tempo. Em 11 anos, homens têm probabilidade de sobrevivência de aproximadamente 70% vs. 77% para mulheres. É notável que essa diferença seja relativamente modesta em comparação com o efeito da faixa etária.

O gráfico por faixa etária é o mais impactante visualmente: as curvas das quatro faixas são completamente separadas e sem sobreposição dos intervalos de confiança. O grupo 80+ apresenta probabilidade de sobrevivência em 5 anos de aproximadamente 50%, enquanto o grupo 50 a 59 anos mantém probabilidade acima de 95% nesse mesmo horizonte, razão de risco implícita na ordem de 10:1. Isso justifica plenamente o ajuste por idade em qualquer análise de mortalidade nessa coorte.

Tempo Mediano de Sobrevivência por Grupo

Ocultar código

# A mediana (S(t) = 0.5) é indefinida nos grupos com FLC mais baixo porque
# a curva KM não cai abaixo de 50% durante os ~11 anos de seguimento.
# Isso indica que MAIS de 50% desses indivíduos sobreviveram ao período todo.
# Mostramos o percentil 25 (S(t) = 0.75) como medida alternativa.

survfit(Surv(anos, death) ~ flc_cat, data = flchain_clean) |>
  tbl_survfit(
    times        = c(2, 5, 8, 11),
    label_header = "**$\\hat{{S}}(t = {time})$**"
  ) |>
  add_n() |>
  modify_caption("**Tabela 2.** Estimativas de $\\hat{{S}}(t)$ por grupo de FLC em tempos fixos") |>
  bold_labels()

**Tabela 2.** Estimativas de \(\hat{S}(t)\) por grupo de FLC em tempos fixos
Characteristic	N	\(\hat{S}(t = 2)\)	\(\hat{S}(t = 5)\)	\(\hat{S}(t = 8)\)	\(\hat{S}(t = 11)\)
flc_cat	7,871
1-2 (Baixo)		98% (97%, 99%)	96% (95%, 97%)	91% (90%, 93%)	87% (85%, 89%)
3-5 (Medio-baixo)		97% (97%, 98%)	93% (92%, 94%)	88% (87%, 89%)	82% (81%, 84%)
6-8 (Medio-alto)		96% (95%, 96%)	88% (87%, 90%)	81% (80%, 83%)	73% (71%, 75%)
9-10 (Alto)		85% (83%, 87%)	72% (70%, 75%)	60% (58%, 63%)	49% (47%, 52%)

Nota: A mediana de sobrevivência (tempo em que \(\hat{S}(t) = 0{,}50\)) é indefinida para os grupos Baixo, Médio-baixo e Médio-alto porque mais de 50% dos indivíduos nesses grupos sobreviveram durante todo o período de acompanhamento (~11 anos) — um resultado favorável. Apenas o grupo Alto apresenta mediana estimável (~11 anos). Por isso, a tabela acima apresenta \(\hat{S}(t)\) em tempos fixos como alternativa.

Estimador de Nelson-Aalen (Função de Risco Acumulada)

A função de risco acumulada \(\hat{\Lambda}(t)\) complementa o estimador de Kaplan-Meier (Klein; Moeschberger, 2003). Enquanto \(\hat{S}(t)\) mede a probabilidade de sobreviver além de \(t\), \(\hat{\Lambda}(t)\) acumula a intensidade de risco instantâneo. O estimador de Nelson-Aalen (Klein; Moeschberger, 2003, cap. 4) é definido como:

\[\hat{\Lambda}(t) = \sum_{j:\,t_j \leq t} \frac{d_j}{n_j},\]

onde \(d_j\) é o número de eventos e \(n_j\) o número em risco no instante \(t_j\). A relação com o estimador de Kaplan-Meier é \(\hat{S}(t) \approx e^{-\hat{\Lambda}(t)}\), sendo a igualdade exata quando se usa o estimador de Fleming-Harrington. O estimador de Nelson-Aalen é preferível ao KM para visualizar a forma da função de risco e verificar suposições de modelos paramétricos (ver Seção 3.3):

Se \(\hat{\Lambda}(t)\) for aproximadamente linear em \(t\): distribuição exponencial (\(h(t) = \lambda\), constante)
Se linear em \(\log(t)\): distribuição Weibull (\(h(t) = \lambda \gamma t^{\gamma-1}\))
Se apresentar padrão em sino após transformação log: log-normal ou log-logística

Ocultar código

# survfit com type="fh" retorna o estimador de Fleming-Harrington (Nelson-Aalen)
km_na <- survfit(Surv(anos, death) ~ flc_cat, data = flchain_clean,
                 type = "fh")

tidy(km_na) |>
  mutate(
    grupo = str_remove(strata, "flc_cat=") |>
            factor(levels = levels(flchain_clean$flc_cat))
  ) |>
  ggplot(aes(x = time, y = -log(estimate), color = grupo)) +
  geom_step(linewidth = 0.9, alpha = 0.9) +
  scale_color_manual(values = c("#3BB273","#2E86AB","#F4A261","#E94F37")) +
  scale_x_continuous(breaks = seq(0, 14, by = 2)) +
  labs(
    title    = "Função de Risco Acumulada — Estimador de Nelson-Aalen",
    subtitle = "Curvatura crescente indica risco não-constante (distribuição não-exponencial)",
    x        = "Tempo (anos)",
    y        = expression(hat(Lambda)(t)),
    color    = "Grupo de FLC"
  ) +
  theme(legend.position = "bottom")

Estimador de Nelson-Aalen da função de risco acumulada por grupo de FLC

A curvatura crescente de \(\hat{\Lambda}(t)\) em todos os grupos confirma que o risco não é constante ao longo do tempo, o que torna a distribuição exponencial inadequada para esses dados. O comportamento aproximadamente linear em escala log-log, visto no gráfico de adequação Weibull, é consistente com a distribuição Weibull com parâmetro de forma \(\gamma > 1\), indicando risco crescente com a idade. Além disso, o gráfico evidencia que a separação entre grupos começa imediatamente e se mantém proporcional ao longo de quase todo o seguimento, pois as curvas são aproximadamente paralelas em escala logarítmica. Nos anos finais (além de 10 anos), observa-se ligeira convergência das curvas dos grupos Baixo e Médio-baixo com o grupo Médio-alto, resultado consistente com o teste de Schoenfeld e que reforça a ressalva sobre a interpretação dos HRs como estritamente constantes ao longo de todo o período.

Análise de Subgrupo: MGUS

O MGUS (monoclonal gammopathy of undetermined significance) é uma condição hematológica pré-maligna presente em 1,5% da amostra (n = 115). Dado seu papel clínico como marcador de progressão para mieloma múltiplo, investigamos se o efeito da FLC sobre a sobrevivência difere entre indivíduos com e sem MGUS.

Ocultar código

p_mgus_global <- survfit2(Surv(anos, death) ~ mgus_f, data = flchain_clean) |>
  ggsurvfit(linewidth = 1.1) +
  add_confidence_interval(alpha = 0.15) +
  add_risktable(risktable_stats = "n.risk", stats_label = "Em risco") +
  add_pvalue(caption = "Log-rank: {p.value}") +
  scale_color_manual(values = c("#2E86AB","#E94F37")) +
  scale_fill_manual(values  = c("#2E86AB","#E94F37")) +
  scale_ggsurvfit() +
  labs(title = "Sobrevivência por MGUS",
       x = "Tempo (anos)", y = expression(hat(S)(t))) +
  theme(legend.position = "bottom")

p_mgus_flc <- survfit2(Surv(anos, death) ~ flc_cat,
                        data = filter(flchain_clean, mgus == 1)) |>
  ggsurvfit(linewidth = 1.1) +
  add_confidence_interval(alpha = 0.12) +
  add_pvalue(caption = "Log-rank: {p.value}") +
  scale_color_manual(values = c("#3BB273","#2E86AB","#F4A261","#E94F37")) +
  scale_fill_manual(values  = c("#3BB273","#2E86AB","#F4A261","#E94F37")) +
  scale_ggsurvfit() +
  labs(title = "Sobrevivência por FLC (apenas MGUS)",
       x = "Tempo (anos)", y = "") +
  theme(legend.position = "bottom")

p_mgus_global + p_mgus_flc

Curvas de Kaplan-Meier estratificadas por MGUS e grupo de FLC

Ocultar código

survfit(Surv(anos, death) ~ mgus_f, data = flchain_clean) |>
  tbl_survfit(
    times        = c(2, 5, 8, 11),
    label_header = "**$\\hat{{S}}(t = {time})$**"
  ) |>
  add_n() |>
  add_p() |>
  modify_caption("**Tabela 2b.** Sobrevivência estimada por status de MGUS") |>
  bold_labels()

**Tabela 2b.** Sobrevivência estimada por status de MGUS
Characteristic	N	\(\hat{S}(t = 2)\)	\(\hat{S}(t = 5)\)	\(\hat{S}(t = 8)\)	\(\hat{S}(t = 11)\)	p-value¹
mgus_f	7,871					<0.001
Não		94% (94%, 95%)	88% (87%, 89%)	81% (80%, 82%)	74% (73%, 75%)
Sim		98% (96%, 100%)	96% (92%, 99%)	93% (88%, 98%)	88% (83%, 95%)
¹ Log-rank test

Indivíduos com MGUS apresentam, surpreendentemente, maior sobrevivência estimada do que os sem MGUS na Tabela 2b (\(p < 0{,}001\)). Esse resultado contraintuitivo merece discussão cuidadosa. O MGUS é uma condição pré-maligna que, em teoria, deveria aumentar o risco de mortalidade. Três explicações são plausíveis: primeiro, viés de detecção, pois indivíduos com MGUS diagnosticado recebem acompanhamento médico mais intensivo, com monitoramento regular e detecção precoce de complicações, o que pode reduzir sua mortalidade a curto e médio prazos; segundo, confundimento por FLC, pois no modelo univariado sem ajuste por FLC o efeito do MGUS sobre a mortalidade pode estar parcialmente capturado pela própria FLC elevada, que é um critério diagnóstico do MGUS, diluindo o efeito aparente; terceiro, tamanho amostral reduzido, pois com apenas 115 indivíduos com MGUS os intervalos de confiança são muito amplos e o resultado pode ser afetado por variabilidade amostral. No modelo de Cox multivariado ajustado por FLC, o efeito do MGUS é de fato não significativo (\(p \approx 0{,}40\)), o que apoia a segunda explicação.

Dentro do subgrupo com MGUS, o gradiente de risco por grupo de FLC se mantém na direção esperada, embora com menor poder estatístico dado o tamanho reduzido da amostra (n = 115). Esse achado reforça que a FLC é um preditor de mortalidade independente do diagnóstico formal de MGUS.

Modelagem

Justificativa da Estratégia de Modelagem

A estratégia de modelagem segue três etapas sequenciais e complementares.

Primeiro, ajustamos o modelo de Cox (semiparamétrico), que não impõe distribuição ao tempo de sobrevivência e permite interpretação direta via razão de riscos (hazard ratio, HR). A seleção de covariáveis seguiu o procedimento estruturado de Collett (2003), com seleção univariada, eliminação retroativa e inclusão prospectiva, complementada pelo critério AIC via stepAIC (Venables; Ripley, 2002). Os diagnósticos incluem teste de riscos proporcionais (Grambsch; Therneau, 1994), linearidade das contínuas (Therneau; Grambsch; Fleming, 1990), influência (dfbeta/dfbetas) (Belsley; Kuh; Welsch, 1980) e outliers (deviance) (Therneau; Grambsch, 2000).

Segundo, avaliamos graficamente a forma da função de risco via hazard suavizado (estimador B-spline, pacote bshazard (Rebora; Salim; Reilly, 2014)) e gráficos de linearização para orientar a escolha da família paramétrica.

Terceiro, ajustamos seis distribuições paramétricas (exponencial, Weibull, log-normal, log-logística, gama e gama generalizada) via flexsurv::flexsurvreg() (Jackson, 2016), comparadas por AIC, BIC e testes formais LRT (Burnham; Anderson, 2002, 2004). A gama generalizada (\(\mu, \sigma, Q\)) é o modelo guarda-chuva: log-normal (\(Q=0\)), Weibull (\(Q=1\)) e gama (\(Q=\sigma\)) são casos aninhados testáveis via \(\chi^2(1)\) (Self; Liang, 1987).

Preparação para modelagem

flchain_mod <- flchain_clean |>
  filter(!is.na(creatinine)) |>
  mutate(log_creat = log(creatinine))

n_mod    <- nrow(flchain_mod)
n_events <- sum(flchain_mod$death)
epv      <- n_events / 5   # eventos por variável (5 covariáveis candidatas)

cat(sprintf("N para modelagem : %d (%.1f%% da amostra)\n",
            n_mod, 100 * n_mod / nrow(flchain_clean)))

N para modelagem : 6521 (82.8% da amostra)

Preparação para modelagem

cat(sprintf("Eventos          : %d\n", n_events))

Eventos          : 1959

Preparação para modelagem

cat(sprintf("Covariáveis cand.: 5\n"))

Covariáveis cand.: 5

Preparação para modelagem

cat(sprintf("EPV (eventos/var): %.1f — %s\n",
            epv, ifelse(epv >= 10, "ADEQUADO (EPV >= 10)", "INSUFICIENTE")))

EPV (eventos/var): 391.8 — ADEQUADO (EPV >= 10)

Preparação para modelagem

cat(sprintf("Excluídos (NA creatinina): %d\n",
            nrow(flchain_clean) - n_mod))

Excluídos (NA creatinina): 1350

Critério EPV: Peduzzi et al. (1995a) e Peduzzi et al. (1995b) estabelecem que são necessários ao menos 10 eventos por variável candidata no modelo de Cox para estimativas não viesadas. Com 1959 eventos e 5 covariáveis, o EPV é adequado e procedimentos de seleção automática como stepAIC não introduzem viés relevante.

Análise da Função de Risco (Pré-modelagem Paramétrica)

Antes de escolher a família paramétrica, estimamos a função de risco não parametricamente para cada grupo de FLC usando o estimador B-spline do pacote bshazard (Rebora; Salim; Reilly, 2014), que lida naturalmente com censura à direita e estima o parâmetro de suavização automaticamente. A forma do hazard guia a escolha: crescente e monotônico sugere Weibull (\(\gamma > 1\)); constante sugere exponencial; forma em sino sugere log-normal ou log-logística; em banheira sugere gama generalizada.

Ocultar código

cores_flc  <- c("#3BB273","#2E86AB","#F4A261","#E94F37")
niveis_flc <- levels(flchain_mod$flc_cat)

haz_list <- lapply(niveis_flc, function(g) {
  d    <- filter(flchain_mod, flc_cat == g)
  tmax <- quantile(d$anos[d$death == 1], 0.90)
  fit  <- bshazard(Surv(anos, death) ~ 1,
                   data    = d,        # todos os dados, sem filtrar
                   verbose = FALSE)
  # cortar a grade de saida no percentil 90 dos eventos
  tibble(time   = fit$time,
         hazard = fit$hazard,
         grupo  = g) |>
    filter(time <= tmax)
})

bind_rows(haz_list) |>
  mutate(grupo = factor(grupo, levels = niveis_flc)) |>
  ggplot(aes(x = time, y = hazard, color = grupo)) +
  geom_line(linewidth = 1.1, alpha = 0.9) +
  scale_color_manual(values = cores_flc) +
  scale_x_continuous(breaks = seq(0, 14, by = 2)) +
  scale_y_continuous(limits = c(0, NA)) +
  labs(
    title    = "Função de Risco Suavizada por Grupo de FLC",
    subtitle = "Estimador B-spline (bshazard) com suavização automática — sem escolha manual de bandwidth",
    x        = "Tempo (anos)",
    y        = expression(hat(h)(t)),
    color    = "Grupo de FLC"
  ) +
  theme(legend.position = "bottom")

Função de risco suavizada por grupo de FLC (bshazard — B-splines, suavização automática)

O hazard suavizado revela padrões distintos entre os grupos. Os grupos Baixo, Médio-baixo e Médio-alto apresentam risco crescente ao longo de todo o seguimento, com separação progressiva e consistente, compatível com a hipótese de riscos proporcionais. O grupo Alto exibe um padrão diferente: hazard elevado no início (~0,13), queda até aproximadamente 4 anos (~0,06) e leve recuperação posterior. Esse comportamento é biologicamente plausível e reflete heterogeneidade não observada (frailty): os indivíduos mais frágeis do grupo Alto morrem precocemente, e os que sobrevivem aos primeiros anos constituem uma subpopulação selecionada com risco relativo menor. Esse padrão em forma de U é consistente com a gama generalizada como distribuição adequada para o grupo de maior risco, já que essa família acomoda formas de hazard não-monotônicas. A distribuição exponencial (risco constante) é descartada para todos os grupos, e o log-normal e log-logística são candidatos secundários frente à flexibilidade da gama generalizada.

Gráficos de Linearização (Diagnóstico Pré-ajuste)

Os gráficos de linearização verificam graficamente a adequação de cada família: se a distribuição for correta, os pontos seguem uma reta.

Distribuição	Eixo X	Eixo Y	Forma linear esperada
Exponencial	\(t\)	\(-\log\hat{S}(t)\)	Reta pela origem
Weibull	\(\log t\)	\(\log(-\log\hat{S}(t))\)	Retas paralelas por grupo
Log-normal	\(\log t\)	\(\Phi^{-1}(1-\hat{S}(t))\)	Reta
Log-logística	\(\log t\)	\(\text{logit}(1-\hat{S}(t))\)	Reta

Ocultar código

km_lin <- survfit(Surv(anos, death) ~ flc_cat, data = flchain_mod)

df_lin <- tidy(km_lin) |>
  filter(estimate > 0 & estimate < 1) |>
  mutate(
    grupo    = str_remove(strata, "flc_cat=") |>
               factor(levels = niveis_flc),
    logt     = log(time),
    # Exponencial: -log S(t) vs t
    y_exp    = -log(estimate),
    # Weibull: log(-log S(t)) vs log t
    y_wb     = log(-log(estimate)),
    # Log-normal: qnorm(1 - S(t)) vs log t
    y_ln     = qnorm(1 - estimate),
    # Log-logística: logit(1 - S(t)) vs log t
    y_ll     = log((1 - estimate) / estimate)
  )

make_lin_plot <- function(df, xvar, yvar, xlab, ylab, title) {
  ggplot(df, aes(x = .data[[xvar]], y = .data[[yvar]], color = grupo)) +
    geom_line(linewidth = 0.7, alpha = 0.8) +
    scale_color_manual(values = cores_flc, name = "FLC") +
    labs(title = title, x = xlab, y = ylab) +
    theme(legend.position = "bottom",
          plot.title = element_text(size = 10, face = "bold"))
}

p_exp <- make_lin_plot(df_lin, "time", "y_exp",
                       "t (anos)", expression(-log~hat(S)(t)),
                       "Exponencial")
p_wb  <- make_lin_plot(df_lin, "logt", "y_wb",
                       "log(t)", expression(log(-log~hat(S)(t))),
                       "Weibull")
p_ln  <- make_lin_plot(df_lin, "logt", "y_ln",
                       "log(t)", expression(Phi^{-1}*(1-hat(S)(t))),
                       "Log-normal")
p_ll  <- make_lin_plot(df_lin, "logt", "y_ll",
                       "log(t)", expression(logit*(1-hat(S)(t))),
                       "Log-logística")

(p_exp + p_wb) / (p_ln + p_ll) +
  plot_annotation(
    title    = "Gráficos de Linearização — Diagnóstico Pré-ajuste",
    subtitle = "Retas paralelas (Weibull) ou únicas (demais): linearidade indica boa adequação"
  )

Gráficos de linearização para quatro famílias paramétricas (diagnóstico pré-ajuste)

O gráfico Weibull (superior direito) apresenta retas aproximadamente lineares e paralelas entre grupos, sugerindo boa adequação da família Weibull e proporcionalidade de riscos entre grupos. O gráfico exponencial (superior esquerdo) mostra curvatura acentuada, descartando essa distribuição, pois a linearidade de \(-\log\hat{S}(t)\) em \(t\) pressupõe hazard constante, hipótese já refutada na Seção 2.9. Log-normal e log-logística apresentam curvatura moderada nas caudas, indicando ajuste imperfeito. O diagnóstico gráfico aponta o Weibull como candidato plausível, a ser confirmado formalmente via AIC e LRT na Seção 3.7.

Seleção de Covariáveis para o Modelo de Cox

Procedimento de Collett (2003)

Seguimos o procedimento em quatro passos descrito em Collett (2003, seç. 3.6.1):

Passo 1 (univariada): ajuste univariado de cada covariável; incluir no conjunto candidato aquelas com \(p \leq 0{,}20\).

Passo 2 (backward): ajuste multivariado com todas as candidatas do Passo 1; eliminar sequencialmente as não-significativas com \(p > 0{,}10\).

Passo 3 (forward das descartadas): reintroduzir individualmente as variáveis eliminadas no Passo 1; incluir se \(p \leq 0{,}10\).

Passo 4 (poda final + interações): stepwise bidirecional com \(p_{\mathrm entrada} = p_{\mathrm saída} = 0{,}10\); avaliar interações clinicamente relevantes entre as variáveis no modelo.

Passo 1 — Triagem univariada

vars_candidatas <- c("flc_cat", "age", "sexo", "log_creat", "mgus_f")

# Ajuste univariado de cada covariável
uni_results <- map_dfr(vars_candidatas, function(v) {
  f   <- as.formula(paste("Surv(anos, death) ~", v))
  fit <- coxph(f, data = flchain_mod)
  s   <- summary(fit)
  tibble(
    Covariável   = v,
    `chi2_LRT`  = round(s$logtest["test"], 3),
    gl           = s$logtest["df"],
    `p-valor`    = round(s$logtest["pvalue"], 6),
    Selecionada  = s$logtest["pvalue"] <= 0.20
  )
})

uni_results |>
  gt() |>
  tab_header(
    title    = md("**Passo 1 — Triagem Univariada**"),
    subtitle = md("Variáveis com *p* ≤ 0,20 avançam para o modelo multivariado")
  ) |>
  cols_label(`chi2_LRT` = md("$\\chi^2$ (LRT)"),
             gl          = "g.l.",
             `p-valor`   = md("*p*-valor"),
             Selecionada = md("Entra no Passo 2?")) |>
  fmt_scientific(columns = `p-valor`, decimals = 2) |>
  tab_style(
    style     = cell_fill(color = "#e8f4f8"),
    locations = cells_body(rows = Selecionada == TRUE)
  )

Covariável	\(\chi^2\) (LRT)	g.l.	p-valor	Entra no Passo 2?
Passo 1 — Triagem Univariada
Variáveis com p ≤ 0,20 avançam para o modelo multivariado
flc_cat	635.864	3	0.00	TRUE
age	2185.489	1	0.00	TRUE
sexo	0.573	1	4.49 × 10⁻¹	FALSE
log_creat	241.680	1	0.00	TRUE
mgus_f	9.910	1	1.64 × 10⁻³	TRUE

Passo 2 — Backward elimination (p > 0,10)

# Modelo com todas as candidatas selecionadas no Passo 1
vars_p1 <- uni_results |> filter(Selecionada) |> pull(Covariável)
formula_p1 <- as.formula(
  paste("Surv(anos, death) ~", paste(vars_p1, collapse = " + "))
)
cox_p2 <- coxph(formula_p1, data = flchain_mod)

# Backward via MASS::stepAIC com k = qchisq(0.90, 1) ≈ 2.706 (equivalente p = 0.10)
cox_backward <- stepAIC(cox_p2,
                        direction = "backward",
                        k         = qchisq(0.90, df = 1),
                        trace     = FALSE)

cat("Modelo após backward elimination (Passo 2):\n")

Modelo após backward elimination (Passo 2):

Passo 2 — Backward elimination (p > 0,10)

print(formula(cox_backward))

Surv(anos, death) ~ flc_cat + age + log_creat
<environment: 0x00000271a79bf9d8>

Passos 3 e 4 — Forward + stepwise bidirecional

# Passo 3: forward das descartadas
vars_descartadas <- setdiff(vars_candidatas, vars_p1)
if (length(vars_descartadas) > 0) {
  cat("Variáveis descartadas no Passo 1 (p > 0,20):", vars_descartadas, "\n")
  cat("Tentando reintrodução individual (Passo 3)...\n\n")
  for (v in vars_descartadas) {
    f_test <- update(formula(cox_backward), paste(". ~ . +", v))
    fit_test <- coxph(f_test, data = flchain_mod)
    p_lrt <- 1 - pchisq(
      2 * (fit_test$loglik[2] - cox_backward$loglik[2]), df = 1)
    cat(sprintf("  + %s: p-valor LRT = %.4f %s\n", v, p_lrt,
                ifelse(p_lrt <= 0.10, "=> INCLUIR", "=> manter fora")))
  }
} else {
  cat("Todas as variáveis avançaram do Passo 1. Passo 3 não se aplica.\n")
}

Variáveis descartadas no Passo 1 (p > 0,20): sexo 
Tentando reintrodução individual (Passo 3)...

  + sexo: p-valor LRT = 0.0001 => INCLUIR

Passos 3 e 4 — Forward + stepwise bidirecional

# Passo 4: stepwise bidirecional com critério AIC
cat("\nPasso 4 — Stepwise bidirecional (AIC):\n")


Passo 4 — Stepwise bidirecional (AIC):

Passos 3 e 4 — Forward + stepwise bidirecional

cox_step <- stepAIC(cox_backward,
                    direction = "both",
                    scope     = list(
                      lower = ~1,
                      upper = as.formula(
                        paste("~", paste(vars_candidatas, collapse = " + ")))
                    ),
                    trace = FALSE)
cat("Modelo final após Passo 4:\n")

Modelo final após Passo 4:

Passos 3 e 4 — Forward + stepwise bidirecional

print(formula(cox_step))

Surv(anos, death) ~ flc_cat + age + log_creat + sexo
<environment: 0x00000271a8947158>