Author

Caio Vallio

Published

September 21, 2025

Abstract
A survival analysis of a synthetic dataset

1 Introdução

Neste documento será realizado um estudo de sobrevivência com o dataset df_sint, para:

  • Estimar a distribuição de tempo até o evento de interesse (evento / censura);
  • Comparar funções de sobrevivência entre grupos (tratamento, sexo, etc.);
  • Ajustar modelos multivariados para identificar efeitos de covariáveis;
  • Verificar pressupostos (principalmente o da proporcionalidade de riscos);
  • Avaliar modelos alternativos caso algum pressuposto seja violado;
  • Oferecer interpretações estatísticas rigorosas dos resultados.

2 Carregamento dos dados

Code
df_sint <- read.csv("data/df_sint.csv", stringsAsFactors = FALSE)

# Visualizar estrutura dos dados
skimr::skim(df_sint)
Data summary
Name df_sint
Number of rows 200
Number of columns 8
_______________________
Column type frequency:
character 4
numeric 4
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
dt_inicio 0 1.00 10 10 0 173 0
dt_evento 96 0.52 10 10 0 101 0
sexo 0 1.00 1 1 0 2 0
tratamento 0 1.00 1 1 0 2 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
id 0 1 100.50 57.88 1 50.75 100.5 150.25 200 ▇▇▇▇▇
status 0 1 0.52 0.50 0 0.00 1.0 1.00 1 ▇▁▁▁▇
idade 0 1 59.38 9.92 31 52.75 60.0 65.25 92 ▁▆▇▃▁
tempo 0 1 247.29 195.14 0 91.00 199.5 364.50 751 ▇▃▃▂▂

Interpretação esperada desta etapa

  • Verificar que todos os registros têm tempo >= 0.
  • Verificar contagem de eventos vs censuras.
  • Verificar estrutura de covariáveis: idade (distribuição), categorias de sexo, tratamento, possíveis outliers.

3 Análise descritiva

Code
df_sint %>%
    summarise(
        n = n(),
        eventos = sum(status),
        censurados = sum(1 - status),
        prop_eventos = mean(status),
        media_tempo = mean(tempo),
        mediana_tempo = median(tempo),
        min_tempo = min(tempo),
        max_tempo = max(tempo)
    )
Code
df_sint %>%
    group_by(tratamento) %>%
    summarise(
        n = n(),
        eventos = sum(status),
        prop_eventos = mean(status),
        mediana_tempo = median(tempo)
    )
Code
df_sint %>%
    group_by(sexo) %>%
    summarise(
        n = n(),
        eventos = sum(status),
        prop_eventos = mean(status),
        mediana_tempo = median(tempo)
    )

Interpretação esperada

  • Qual a proporção de sujeitos que tiveram evento vs censurados (“status = 1” vs “0”).
  • Mediana de tempo de seguimento geral, e por grupo.
  • Verificar se distribuição de idade tem variabilidade, possíveis outliers.
  • Verificar número de indivíduos por grupo de tratamento e sexo para garantir poder suficiente.

4 Estimativa de Kaplan–Meier (KM)

4.1 Curva de sobrevivência geral

Code
surv_obj <- Surv(time = df_sint$tempo, event = df_sint$status)
km_fit <- survfit(surv_obj ~ 1, data = df_sint)

ggsurvplot(
    km_fit,
    conf.int = TRUE,
    risk.table = TRUE,
    xlab = "Tempo (dias)",
    ylab = "Probabilidade de sobrevivência",
    title = "Curva KM geral para df_sint"
)

Interpretação esperada

  • A curva inicia em 1 (100%) em tempo zero.
  • Em tempos intermediários, queda na sobrevivência conforme ocorrem eventos.
  • A curva de confiança mostra incerteza, que geralmente aumenta em tempos mais avançados (porque menos sujeitos em risco).
  • Tabela de risco mostra quantos indivíduos ainda sob observação (“em risco”) ao longo do tempo.

4.2 Curvas de sobrevivência por grupo (tratamento, sexo)

Code
km_trat <- survfit(surv_obj ~ tratamento, data = df_sint)
ggsurvplot(
    km_trat,
    conf.int = TRUE,
    pval = TRUE,
    risk.table = TRUE,
    legend.title = "Tratamento",
    xlab = "Tempo (dias)",
    ylab = "Sobrevivência",
    title = "KM por Tratamento"
)

Code
km_sexo <- survfit(surv_obj ~ sexo, data = df_sint)
ggsurvplot(
    km_sexo,
    conf.int = TRUE,
    pval = TRUE,
    risk.table = TRUE,
    legend.title = "Sexo",
    xlab = "Tempo (dias)",
    ylab = "Sobrevivência",
    title = "KM por Sexo"
)

Interpretação esperada

  • Se as curvas de tratamento A e B divergem, pode haver efeito do tratamento. Observações visuais: diferenças iniciais, diferenças tardias.
  • P-valor do teste de log-rank indica se a diferença entre curvas é estatisticamente significativa.
  • Atenção à censura desigual nos grupos: se um grupo tiver muitos censurados cedo, pode influenciar visualização.

5 Teste de log-rank

Code
survdiff(surv_obj ~ tratamento, data = df_sint)
Call:
survdiff(formula = surv_obj ~ tratamento, data = df_sint)

               N Observed Expected (O-E)^2/E (O-E)^2/V
tratamento=A 103       53     54.5    0.0404    0.0856
tratamento=B  97       51     49.5    0.0444    0.0856

 Chisq= 0.1  on 1 degrees of freedom, p= 0.8 
Code
survdiff(surv_obj ~ sexo, data = df_sint)
Call:
survdiff(formula = surv_obj ~ sexo, data = df_sint)

         N Observed Expected (O-E)^2/E (O-E)^2/V
sexo=F 105       58     57.9  0.000200  0.000454
sexo=M  95       46     46.1  0.000252  0.000454

 Chisq= 0  on 1 degrees of freedom, p= 1 

Interpretação esperada

  • O teste de log-rank testa a hipótese nula de que as curvas de sobrevivência dos grupos são iguais ao longo do tempo.
  • Se p < 0.05, rejeita-se essa hipótese.
  • Mas o teste assume proporcionalidade de riscos entre os grupos; se violado, os resultados podem ser enganosos. (Ver seção de pressupostos abaixo).

5.1 Modelo de Cox Proporcional de Riscos

Code
cox_mod <- coxph(surv_obj ~ idade + sexo + tratamento, data = df_sint)
summary(cox_mod)
Call:
coxph(formula = surv_obj ~ idade + sexo + tratamento, data = df_sint)

  n= 200, number of events= 104 

                  coef  exp(coef)   se(coef)      z Pr(>|z|)
idade        0.0038932  1.0039008  0.0103307  0.377    0.706
sexoM       -0.0006181  0.9993821  0.1979864 -0.003    0.998
tratamentoB  0.0452310  1.0462695  0.1994253  0.227    0.821

            exp(coef) exp(-coef) lower .95 upper .95
idade          1.0039     0.9961    0.9838     1.024
sexoM          0.9994     1.0006    0.6780     1.473
tratamentoB    1.0463     0.9558    0.7078     1.547

Concordance= 0.485  (se = 0.034 )
Likelihood ratio test= 0.23  on 3 df,   p=1
Wald test            = 0.23  on 3 df,   p=1
Score (logrank) test = 0.23  on 3 df,   p=1

Interpretação dos resultados do summary(cox_mod)

Output O que representa
Coeficiente β para cada covariável Log do hazard ratio para uma unidade da covariável (ou comparação de categoria)
exp(coef) Hazard Ratio (HR): efeito multiplicativo no risco
erro padrão (se(coef)) Incerteza do estimador da β
z e valor-p Teste de hipótese H₀: β = 0 (ou HR = 1)
Intervalo de confiança de 95% para HR Se não incluir 1 ⇒ efeito “significativo” ao nível de 5%
Estatísticas de verossimilhança, AIC Para comparar modelos (se fizer mais de um)

6 Pressupostos do modelo de Cox e testes diagnósticos

Aqui são os principais pressupostos, como testá-los, e o que fazer se forem violados.

6.1 Pressupostos importantes

1. Riscos proporcionais (Proportional Hazards, PH)

O hazard ratio entre níveis de covariáveis (ex: tratamento A vs B) deve permanecer constante ao longo do tempo.

2. Censura não informativa (independent / random censoring)

A censura deve ocorrer de forma que não esteja relacionada com a probabilidade de ter o evento.

3. Independência entre observações

Cada indivíduo deve ser independente; se houver estrutura de agrupamento pode ser necessário usar modelos de frailty ou incluir variáveis de cluster.

4. Linearidade para covariáveis contínuas no termo do log-hazard

Efeitos das variáveis contínuas devem ser lineares no log-hazard. Se não for, considerar transformações ou splines.

5. Covariáveis constantes no tempo, ou modelagem apropriada se variáveis mudam ao longo do tempo.

6.2 Teste da proporcionalidade de riscos

Code
zph_res <- cox.zph(cox_mod)
zph_res
           chisq df     p
idade       2.73  1 0.098
sexo        2.16  1 0.142
tratamento  2.49  1 0.115
GLOBAL      6.34  3 0.096
Code
plot(zph_res)

Code
ggcoxzph(zph_res)

  • cox.zph() retorna testes estatísticos (baseado nos resíduos de Schoenfeld) para cada covariável e um teste global.
  • Um p valor < 0.05 sugere violação da proporcionalidade para aquela covariável.
  • Os gráficos de resíduos escalonados de Schoenfeld permitem visualizar se há tendência de variação do β ao longo do tempo (linha de tendência não horizontal indica suspeita de violação).

6.3 Outras verificações diagnósticas

  • Resíduos de Martingale: para verificar forma funcional de covariáveis contínuas (linearidade). Se padrão sistemático no gráfico do resíduo vs valor da covariável, pode indicar necessidade de transformação ou de usar termos não-lineares/splines.
  • Resíduos de deviance ou deviance residuals: para detectar influências ou observações atípicas.

7 Ações se algum pressuposto for violado

Code
# Exemplo 1: estratificação por covariável que viola proporcionalidade
cox_mod_strata <- coxph(surv_obj ~ idade + sexo + strata(tratamento), data = df_sint)
summary(cox_mod_strata)
Call:
coxph(formula = surv_obj ~ idade + sexo + strata(tratamento), 
    data = df_sint)

  n= 200, number of events= 104 

           coef exp(coef)  se(coef)      z Pr(>|z|)
idade  0.001963  1.001965  0.010433  0.188    0.851
sexoM -0.005011  0.995002  0.198308 -0.025    0.980

      exp(coef) exp(-coef) lower .95 upper .95
idade     1.002      0.998    0.9817     1.023
sexoM     0.995      1.005    0.6746     1.468

Concordance= 0.483  (se = 0.034 )
Likelihood ratio test= 0.04  on 2 df,   p=1
Wald test            = 0.04  on 2 df,   p=1
Score (logrank) test = 0.04  on 2 df,   p=1
Code
# Exemplo 2: incluir efeito dependente do tempo
# Suponha tratamento violou, então:
df_sint$tratB <- as.numeric(df_sint$tratamento == "B")

cox_mod_td <- coxph(
    surv_obj ~ idade + sexo + tratB + tt(tratB),
    data = df_sint,
    tt = function(x, t, ...) x * log(t + 1)
)

summary(cox_mod_td)
Call:
coxph(formula = surv_obj ~ idade + sexo + tratB + tt(tratB), 
    data = df_sint, tt = function(x, t, ...) x * log(t + 1))

  n= 200, number of events= 104 

               coef exp(coef)  se(coef)      z Pr(>|z|)
idade      0.002786  1.002790  0.010374  0.269    0.788
sexoM     -0.001373  0.998628  0.198054 -0.007    0.994
tratB     -0.675499  0.508902  0.715141 -0.945    0.345
tt(tratB)  0.155727  1.168507  0.147895  1.053    0.292

          exp(coef) exp(-coef) lower .95 upper .95
idade        1.0028     0.9972    0.9826     1.023
sexoM        0.9986     1.0014    0.6774     1.472
tratB        0.5089     1.9650    0.1253     2.067
tt(tratB)    1.1685     0.8558    0.8745     1.561

Concordance= 0.534  (se = 0.031 )
Likelihood ratio test= 1.35  on 4 df,   p=0.9
Wald test            = 1.34  on 4 df,   p=0.9
Score (logrank) test = 1.35  on 4 df,   p=0.9
Code
# alternativa: tt() serviço do survival para interações com tempo
  • Estratificação: a variável que violou é usada como strata(…), assim não se estima coeficiente para ela, mas permite que a base‐hazard varie por estrato.
  • Interação com tempo / coeficientes dependentes do tempo: modelar β_j(t) para aquela covariável.
  • Modelos paramétricos ou de intervalo de tempo: ex: modelo Weibull, Log-normal, ou dividir seguimento em pedaços de tempo onde PH pode valer.

8 Modelos paramétricos alternativos

Code
# Diagnóstico
# summary(df_sint$tempo)
# table(df_sint$tempo == 0, useNA = "ifany")

# Remover registros inválidos (negativos) e ajustar zeros
df_sint_fix <- df_sint |>
    filter(!is.na(tempo) & tempo >= 0) |>
    mutate(tempo_aft = ifelse(tempo <= 0, 0.5, tempo)) # meia-dia mínimo (ajuste prático)

# Verifique de novo
# summary(df_sint_fix$tempo_aft)

# AFT Weibull / Log-normal
weibull_mod <- survreg(Surv(tempo_aft, status) ~ idade + sexo + tratamento,
    data = df_sint_fix, dist = "weibull"
)

lognormal_mod <- survreg(Surv(tempo_aft, status) ~ idade + sexo + tratamento,
    data = df_sint_fix, dist = "lognormal"
)

summary(weibull_mod)

Call:
survreg(formula = Surv(tempo_aft, status) ~ idade + sexo + tratamento, 
    data = df_sint_fix, dist = "weibull")
               Value Std. Error     z      p
(Intercept)  6.53748    0.73370  8.91 <2e-16
idade       -0.00517    0.01201 -0.43  0.667
sexoM       -0.00805    0.22815 -0.04  0.972
tratamentoB -0.00645    0.22868 -0.03  0.978
Log(scale)   0.14288    0.08404  1.70  0.089

Scale= 1.15 

Weibull distribution
Loglik(model)= -743.5   Loglik(intercept only)= -743.6
    Chisq= 0.19 on 3 degrees of freedom, p= 0.98 
Number of Newton-Raphson Iterations: 5 
n= 200 
Code
summary(lognormal_mod)

Call:
survreg(formula = Surv(tempo_aft, status) ~ idade + sexo + tratamento, 
    data = df_sint_fix, dist = "lognormal")
               Value Std. Error     z       p
(Intercept)  5.57219    0.90371  6.17 7.0e-10
idade        0.00498    0.01457  0.34    0.73
sexoM       -0.16721    0.28944 -0.58    0.56
tratamentoB  0.02713    0.28889  0.09    0.93
Log(scale)   0.59961    0.07356  8.15 3.6e-16

Scale= 1.82 

Log Normal distribution
Loglik(model)= -749.8   Loglik(intercept only)= -750.1
    Chisq= 0.47 on 3 degrees of freedom, p= 0.92 
Number of Newton-Raphson Iterations: 3 
n= 200 

Interpretação

  • O survreg fornece estimativas em escala de log‐tempo etc.; os coeficientes não são hazard ratios diretamente, mas podem ser convertidos ou interpretados conforme a escala acelerada‐“time” (AFT: Accelerated Failure Time).
  • Comparar AIC ou critério similar entre modelos (Cox vs Weibull vs Lognormal) para ver qual modelo parece se ajustar melhor aos seus dados.
  • Se modelo paramétrico ajustar bem, pode oferecer estimativas diretas de quantis de tempo de sobrevivência (medianas, percentis) e prever tempos de sobrevivência para perfis de covariáveis.

Em AFT, os coeficientes estão na escala do log-tempo; não são HR direto. Converta/relate corretamente na interpretação (p. ex., tempo “acelerado” ou “desacelerado” por um fator ao exponentiar).


9 Visualizações finais e resumo de resultados

Code
ggforest(cox_mod, data = df_sint)

Code
tidy(cox_mod, exponentiate = TRUE, conf.int = TRUE)

O que reportar no texto:

  • Hazard Ratios estimados, seus intervalos de confiança, e quais covariáveis tiveram efeitos estatisticamente significativos.
  • Mediana de sobrevivência se estimada e tempos de sobrevivência em pontos de interesse.
  • Se houver violação de proporcionalidade, descrição do que foi feito (estratificação, interação, uso de modelo alternativo).
  • Limitações: censura, tamanho de amostra, observações extremas, covariáveis não medidas ou tempo‐dependentes, supostos não verificáveis (e.g. censura não informativa).

10 Conclusão

Aqui você vai sintetizar:

  • Principais achados: quais covariáveis estão associadas com maior risco, qual o efeito estimado do tratamento, idade etc.
  • Qualidade do ajuste: se pressupostos atendidos ou violados, como isso pode afetar interpretação.
  • Sugestões para análises futuras ou melhorias (por exemplo, coleta de variáveis tempo‐dependentes, usar modelos de riscos competitivos, validação externa).