Scoliosis Model

Autor

Caio Vallio

Data de Publicação

7 de abril de 2026

Resumo

Contexto: A escoliose idiopática do adolescente pode apresentar resposta clínica variável a intervenções conservadoras. Objetivo: Explorar associações entre características basais e evolução radiográfica em 6 meses, utilizando abordagem estatística otimizada. Métodos: Estudo observacional com análise exploratória e inferencial. Implementamos três tipos de modelos: (1) modelo linear único para o desfecho contínuo \(\Delta\) (mudança na maior curva), (2) modelos logísticos para desfechos binários de melhora clínica (\(\Delta \le -5°\)) e progressão clínica (\(\Delta \ge +5°\)), e (3) modelos CART para identificação de perfis de risco. Realizamos diagnósticos completos de adequação para todos os modelos, reportando coeficientes/odds ratios com IC95%.

0.1 Como reproduzir

  • O arquivo de dados é lido de data/dataset_escoliose_01.xlsx (aba dados).
  • As análises dependem de pacotes R listados no chunk de setup; caso algum pacote esteja ausente, o documento interrompe a execução com uma mensagem explícita.

0.2 Objetivo e delineamento

Este é um relatório exploratório e inferencial: o objetivo é estimar associações ajustadas entre variáveis basais e a evolução radiográfica em 6 meses.

0.3 Estrutura dos Modelos

A análise utiliza uma abordagem estatística otimizada com três tipos de modelos complementares:

0.3.1 Modelo Linear (Desfecho Contínuo)

  • Objetivo: Modelar diretamente o \(\Delta\) (mudança na maior curva em graus)
  • Interpretação: Coeficientes negativos indicam associação com melhora, positivos com piora
  • Vantagem: Aproveita toda a informação contínua do desfecho, sem perda por dicotomização

0.3.2 Modelos para Desfechos Binários

  • Melhora Clínica (\(\Delta \le -5°\)): Modelos logístico e CART
  • Progressão Clínica (\(\Delta \ge +5°\)): Modelos logístico e CART
  • Objetivo: Identificar fatores associados aos desfechos clinicamente relevantes e perfis de risco

0.4 Definições de desfecho

  • Desfecho contínuo: \(\Delta\) = (maior curva em 6 meses − maior curva baseline), em graus. Valores negativos indicam melhora, positivos indicam piora.
  • Melhora clínica: \(\Delta \le -5°\) (redução de pelo menos 5 graus na maior curva).
  • Progressão clínica: \(\Delta \ge +5°\) (aumento de pelo menos 5 graus na maior curva).

0.5 Plano de análise estatística

  • Descrição da amostra: estatísticas descritivas das variáveis preditoras e desfechos.
  • Inferência (associações ajustadas):
    • Modelo linear único para \(\Delta\) (coeficientes com IC95%, interpretação direcionada).
    • Modelos logísticos específicos para melhora e progressão (odds ratios com IC95%).
    • Modelos CART para ambos os desfechos binários (regras interpretáveis e perfis de risco).
  • Adequação dos modelos: diagnósticos completos de pressupostos para todos os modelos (normalidade, heteroscedasticidade, colinearidade, outliers, adequação do ajuste).
  • Forma funcional: relações aproximadamente lineares entre preditores contínuos e desfechos.

Por se tratar de análise exploratória, os resultados devem ser interpretados com cautela, com ênfase em magnitude/direção e incerteza (IC95%).

Código 
# =============================================================================
# SETUP E PREPARAÇÃO DOS DADOS
# =============================================================================

## Carregamento e verificação de pacotes necessários
pkgs <- c(
    "tidyverse", # manipulação de dados e visualização
    "gtsummary", # tabelas estatísticas
    "readxl", # leitura de arquivos Excel
    "janitor", # limpeza de nomes de variáveis
    "performance", # diagnósticos de modelos
    "qqplotr", # gráficos Q-Q
    "PupillometryR" # gráficos violin plot
)

# Verificar se todos os pacotes estão instalados
missing_pkgs <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_pkgs) > 0) {
    stop("Pacotes ausentes: ", paste(missing_pkgs, collapse = ", "))
}

# Carregar pacotes
invisible(lapply(pkgs, library, character.only = TRUE))

# Configura<U+00E7><U+00F5>es de reprodutibilidade e tema
set.seed(123)
theme_set(theme_minimal())


## Leitura dos dados
df_raw <- readxl::read_excel(
    "data/dataset_escoliose_01.xlsx",
    sheet = "dados",
    na = "" # tratar células vazias como NA
)

## Limpeza dos nomes das variáveis (snake_case)
df <- df_raw |> janitor::clean_names()

# =============================================================================
# CRIAÇÃO DE VARIÁVEIS DERIVADAS
# =============================================================================

## Classificação do escoliômetro: identificar regi<U+00F5>es com deformidade > 10°
# Esta variável indica onde a gibosidade é mais pronunciada
df <- df |>
    mutate(
        escoliometro_maior_10_graus = case_when(
            # Múltiplas regi<U+00F5>es afetadas
            (escoliometro_cervical > 10) & (escoliometro_torarica > 10) & (escoliometro_lombar > 10) ~
                "cervical_toracica_lombar",
            (escoliometro_cervical > 10) & (escoliometro_torarica > 10) ~
                "cervical_toracica",
            (escoliometro_cervical > 10) & (escoliometro_lombar > 10) ~
                "cervical_lombar",
            (escoliometro_torarica > 10) & (escoliometro_lombar > 10) ~
                "toracica_lombar",

            # Região única afetada
            (escoliometro_cervical > 10) ~ "cervical",
            (escoliometro_torarica > 10) ~ "toracica",
            (escoliometro_lombar > 10) ~ "lombar",

            # Nenhuma região > 10° (referência)
            TRUE ~ "normal"
        )
    )

# Variáveis categóricas não são mais utilizadas nos modelos finais
# (mantendo apenas as variáveis numéricas conforme plano de análise)

## Cálculo do Índice de Massa Corporal (IMC)
df <- df |>
    mutate(imc = peso / (altura)^2) # altura deve estar em metros

# =============================================================================
# DEFINIÇÃO DOS DESFECHOS PRIMÁRIOS
# =============================================================================

## Identificação da maior curva no baseline
# Compara as três regi<U+00F5>es e identifica o maior ângulo de Cobb inicial
df <- df |>
    mutate(
        # Maior curva entre todas as regi<U+00F5>es medidas
        cobb_inicial_maior = pmax(
            cobb_toracico_proximal,
            cobb_toracica,
            cobb_lombar,
            na.rm = TRUE
        ),

        # Identificar qual região apresenta a maior curva (para referência)
        regiao_cobb_inicial = case_when(
            cobb_inicial_maior == cobb_toracico_proximal ~ "toracico_proximal",
            cobb_inicial_maior == cobb_toracica ~ "toracica",
            cobb_inicial_maior == cobb_lombar ~ "lombar",
            TRUE ~ NA_character_
        )
    )


# =============================================================================
# CRIAÇÃO DAS VARIÁVEIS DE DESFECHO
# =============================================================================

## Variáveis de desfecho para os dois modelos principais
df <- df |>
    # DESFECHO CONTÍNUO: Mudan<U+00E7>a na maior curva (em graus)
    # Valores negativos = melhora, valores positivos = piora
    mutate(delta = maior_curva_6_meses - cobb_inicial_maior) |>
    # MODELO 1: MELHORA CLÍNICA (delta ≤ -5°)
    mutate(
        delta_cat = if_else(delta <= -5, 1, 0),
        delta_cat_f = factor(
            delta_cat,
            levels = c(0, 1),
            labels = c("Sem melhora (MCID)", "Melhora (MCID)")
        )
    ) |>
    # MODELO 2: PROGRESS<U+00C3>O CLÍNICA (delta ≥ +5°)
    mutate(
        delta_progression = if_else(delta >= 5, 1, 0),
        delta_progression_f = factor(
            delta_progression,
            levels = c(0, 1),
            labels = c("Não progrediu", "Progrediu")
        )
    )

# Conversão de variáveis para tipos adequados nos modelos
df <- df |>
    mutate(
        # Variáveis numéricas
        idade = as.double(idade),

        # Variáveis categóricas (fatores com níveis específicos)
        lenke = as.factor(lenke),
        risser = as.factor(risser),
        sexo = as.factor(sexo),
        flexibilidade = as.factor(flexibilidade),

        # Escoliômetro: definir "normal" como categoria de referência
        escoliometro_maior_10_graus = forcats::fct_relevel(
            as.factor(escoliometro_maior_10_graus),
            "normal"
        )
    )

0.6 Dados faltantes e tamanho amostral

Código 
# =============================================================================
# ANÁLISE DE DADOS FALTANTES
# =============================================================================

# Calcular quantidade e percentual de dados faltantes por variável
missing_tbl <- df |>
    summarise(across(everything(), ~ sum(is.na(.)))) |>
    pivot_longer(
        everything(),
        names_to = "variavel",
        values_to = "n_missing"
    ) |>
    mutate(pct_missing = n_missing / nrow(df)) |>
    arrange(desc(pct_missing)) # ordenar por maior percentual faltante

# Exibir tabela de dados faltantes (formato percentual)
missing_tbl |>
    mutate(pct_missing = scales::percent(pct_missing, accuracy = 0.1)) |>
    print(n = 50)
# A tibble: 28 × 3
   variavel                    n_missing pct_missing
   <chr>                           <int> <chr>      
 1 cobb_toracico_proximal            506 81.5%      
 2 escoliometro_cervical             266 42.8%      
 3 cobb_lombar                        38 6.1%       
 4 escoliometro_lombar                32 5.2%       
 5 cobb_toracica                      23 3.7%       
 6 escoliometro_torarica              17 2.7%       
 7 lenke                               3 0.5%       
 8 id                                  0 0.0%       
 9 idade                               0 0.0%       
10 sexo                                0 0.0%       
11 altura                              0 0.0%       
12 peso                                0 0.0%       
13 risser                              0 0.0%       
14 cifose_toracica                     0 0.0%       
15 lordose_lombar                      0 0.0%       
16 flexibilidade                       0 0.0%       
17 dif_colete                          0 0.0%       
18 correcao_colete                     0 0.0%       
19 maior_curva_6_meses                 0 0.0%       
20 escoliometro_maior_10_graus         0 0.0%       
21 imc                                 0 0.0%       
22 cobb_inicial_maior                  0 0.0%       
23 regiao_cobb_inicial                 0 0.0%       
24 delta                               0 0.0%       
25 delta_cat                           0 0.0%       
26 delta_cat_f                         0 0.0%       
27 delta_progression                   0 0.0%       
28 delta_progression_f                 0 0.0%
Código 
# Tratamento de dados faltantes
# Remover observa<U+00E7><U+00F5>es com flexibilidade faltante (variável key do modelo)
df <- df |> drop_na(flexibilidade)

Observação: os modelos abaixo usam análise por caso completo (listwise deletion), que é o comportamento padrão do glm()/lm() quando há dados faltantes nas variáveis do modelo.

0.7 Codificação de variáveis categóricas (referências)

As variáveis categóricas entram como fatores. A categoria de referência (baseline) é o primeiro nível do fator (conforme levels()), a menos que seja explicitamente reordenado.

Código 
# Lista de variáveis categóricas utilizadas nos modelos finais
categoricas <- c(
    "sexo",
    "lenke",
    "risser",
    "flexibilidade",
    "escoliometro_maior_10_graus"
)
map_dfr(categoricas, function(v) {
    tibble(
        variavel = v,
        niveis = paste(levels(df[[v]]), collapse = " | ")
    )
}) |>
    print(n = Inf)
# A tibble: 5 × 2
  variavel                    niveis                                      
  <chr>                       <chr>                                       
1 sexo                        feminino | masculino                        
2 lenke                       1 | 2 | 3 | 4 | 5 | 6                       
3 risser                      0 | 1 | 2 | 3 | 4                           
4 flexibilidade               flexivel | rigido                           
5 escoliometro_maior_10_graus normal | lombar | toracica | toracica_lombar




1 Análise descritiva

1.1 Variáveis preditoras dos modelos

Variáveis coletadas na linha de base e utilizadas nos modelos finais.

Código 
df |>
    gtsummary::tbl_summary(
        include = c(
            # Variáveis numéricas
            idade,
            imc,
            cifose_toracica,
            lordose_lombar,
            correcao_colete,

            # Variáveis categóricas
            sexo,
            lenke,
            risser,
            flexibilidade,
            escoliometro_maior_10_graus,
        ),
        type = list(
            idade ~ "continuous"
        ),
        statistic = list(
            gtsummary::all_continuous() ~ "{mean} ({sd})",
            gtsummary::all_categorical() ~ "{n} ({p}%)"
        )
    )
Characteristic N = 6211
idade 13.19 (1.70)
imc 19.21 (2.92)
cifose_toracica 23 (11)
lordose_lombar 55 (12)
correcao_colete 47 (20)
sexo
    feminino 556 (90%)
    masculino 65 (10%)
lenke
    1 37 (6.0%)
    2 41 (6.6%)
    3 137 (22%)
    4 60 (9.7%)
    5 55 (8.9%)
    6 288 (47%)
    Unknown 3
risser
    0 132 (21%)
    1 95 (15%)
    2 124 (20%)
    3 127 (20%)
    4 143 (23%)
flexibilidade
    flexivel 393 (63%)
    rigido 228 (37%)
escoliometro_maior_10_graus
    normal 147 (24%)
    lombar 204 (33%)
    toracica 224 (36%)
    toracica_lombar 46 (7.4%)
1 Mean (SD); n (%)

1.2 Desfechos principais

  • Maior curva na linha de base e em 6 meses.
  • Delta: maior curva 6 meses - maior curva baseline.
  • Melhora clínica: delta ≤ -5 graus (redução de pelo menos 5 graus).
  • Progressão clínica: delta ≥ +5 graus (aumento de pelo menos 5 graus).
Código 
df |>
    gtsummary::tbl_summary(
        include = c(
            cobb_inicial_maior,
            maior_curva_6_meses,
            delta,
            delta_cat_f,
            delta_progression_f
        ),
        type = list(
            delta_cat_f ~ "categorical",
            delta_progression_f ~ "categorical"
        ),
        statistic = list(
            gtsummary::all_continuous() ~ "{mean} ({sd})",
            gtsummary::all_categorical() ~ "{n} ({p}%)"
        )
    )
Characteristic N = 6211
cobb_inicial_maior 36.6 (5.6)
maior_curva_6_meses 32 (8)
delta -4.6 (5.4)
delta_cat_f
    Sem melhora (MCID) 301 (48%)
    Melhora (MCID) 320 (52%)
delta_progression_f
    Não progrediu 596 (96%)
    Progrediu 25 (4.0%)
1 Mean (SD); n (%)

1.3 Distribuição do desfecho principal

Código 
df |>
    ggplot(aes(x = "", y = delta)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        alpha = 0.7
    ) +
    labs(
        title = "Distribution of the outcome (delta)",
        x = "",
        y = "delta (degrees)"
    ) +
    geom_point(position = position_jitter(w = .15)) +
    geom_boxplot(
        width = .25,
        alpha = 0.7,
        outlier.shape = NA
    ) +
    coord_flip() +
    theme(
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
    )

Código 
df |>
    ggplot(aes(sample = delta)) +
    qqplotr::stat_qq_band(alpha = 0.2) +
    qqplotr::stat_qq_line() +
    qqplotr::stat_qq_point() +
    labs(
        title = "Outcome QQ-plot (delta)",
        x = "Theoretical quantiles",
        y = "Sample quantiles"
    )




2 Modelagem

2.1 Modelo Linear - Desfecho Contínuo (Delta)

2.1.1 Interpretação dos coeficientes

Os coeficientes do modelo linear representam a mudança esperada no \(\Delta\) (em graus) para cada unidade de aumento na variável preditora:

  • Coeficientes negativos: indicam associação com melhora clínica (redução na curvatura)
  • Coeficientes positivos: indicam associação com progressão clínica (aumento na curvatura)
  • Magnitude: quantifica o efeito em graus de mudança na maior curva
  • Intervalos de confiança: estimam a precisão e significância estatística dos efeitos
Código 
# =============================================================================
# MODELO LINEAR - DESFECHO CONTÍNUO (DELTA)
# =============================================================================

# Modelo único para o desfecho contínuo (delta em graus)
# Interpretação direta: coeficientes negativos = melhora, positivos = piora
# Utiliza toda a informação contínua sem perda por dicotomização

model_delta <- lm(
    delta ~
        # Variáveis numéricas
        idade + # idade do paciente
        imc + # índice de massa corporal
        cifose_toracica + # ângulo da cifose torácica
        lordose_lombar + # ângulo da lordose lombar
        correcao_colete + # percentual de correção com colete

        # Variáveis categóricas
        sexo + # sexo do paciente
        lenke + # classificação de Lenke (numérica como factor)
        risser + # sinal de Risser (numérica como factor)
        flexibilidade + # flexibilidade da curva
        escoliometro_maior_10_graus, # região com gibosidade > 10°
    data = df
)

2.1.2 Parâmetros do modelo

Código 
model_delta |>
    gtsummary::tbl_regression(conf.level = 0.95) |>
    # gtsummary::add_global_p() |>
    gtsummary::bold_p() |>
    gtsummary::bold_labels() |>
    gtsummary::italicize_levels()
Characteristic Beta 95% CI p-value
idade 0.06 -0.21, 0.32 0.7
imc 0.02 -0.10, 0.13 0.8
cifose_toracica -0.04 -0.07, 0.00 0.039
lordose_lombar -0.01 -0.04, 0.02 0.6
correcao_colete -0.12 -0.14, -0.11 <0.001
sexo


    feminino — —
    masculino -0.63 -1.8, 0.55 0.3
lenke


    1 — —
    2 0.22 -1.7, 2.1 0.8
    3 3.7 2.2, 5.3 <0.001
    4 4.2 2.4, 5.9 <0.001
    5 1.1 -0.74, 3.0 0.2
    6 2.3 0.76, 3.9 0.004
risser


    0 — —
    1 -2.6 -3.7, -1.4 <0.001
    2 -1.9 -3.1, -0.83 <0.001
    3 -1.3 -2.4, -0.16 0.025
    4 -0.89 -2.2, 0.40 0.2
flexibilidade


    flexivel — —
    rigido 0.91 0.18, 1.6 0.015
escoliometro_maior_10_graus


    normal — —
    lombar -0.40 -1.4, 0.56 0.4
    toracica -0.03 -1.1, 0.98 >0.9
    toracica_lombar 2.7 1.3, 4.1 <0.001
Abbreviation: CI = Confidence Interval

2.1.3 Análises de pressupostos

2.1.3.1 Heteroscedasticidade

Verificação de heteroscedasticidade dos resíduos (teste de Breusch-Pagan).

Código 
performance::check_heteroscedasticity(model_delta)
OK: Error variance appears to be homoscedastic (p = 0.095).
Código 
performance::check_heteroscedasticity(model_delta) |> plot()

2.1.3.2 Normalidade dos resíduos

Código 
performance::check_residuals(model_delta)
OK: Simulated residuals appear as uniformly distributed (p = 0.487).
Código 
performance::check_residuals(model_delta) |> plot()

2.1.3.3 Autocorrelação dos resíduos

Independência dos resíduos (Durbin-Watson test).

Código 
performance::check_autocorrelation(model_delta)
OK: Residuals appear to be independent and not autocorrelated (p = 0.090).

2.1.3.4 Verificações preditivas posteriores

Simulação de dados replicados sob o modelo ajustado e posterior comparação com os dados observados.

Código 
performance::check_predictions(model_delta) |> plot()

2.1.3.5 Colinearidade

Verificação de colinearidade entre as variáveis independentes (Variance Inflation Factor). Valores maiores que 10 indicam colinearidade entre as variáveis.

Código 
performance::check_collinearity(model_delta) |>
    arrange(desc(VIF)) |>
    select(Term, VIF, VIF_CI_low, VIF_CI_high) |>
    performance::print_html()
Term VIF VIF_CI_low VIF_CI_high
lenke 2.26 2.02 2.56
escoliometro_maior_10_graus 2.17 1.94 2.46
risser 1.97 1.77 2.22
idade 1.80 1.63 2.03
cifose_toracica 1.45 1.33 1.62
lordose_lombar 1.37 1.26 1.52
sexo 1.17 1.09 1.31
correcao_colete 1.11 1.05 1.26
flexibilidade 1.11 1.05 1.25
imc 1.10 1.04 1.25
Código 
performance::check_collinearity(model_delta) |> plot()

2.1.3.6 Outliers

Código 
performance::check_outliers(model_delta)
OK: No outliers detected.
- Based on the following method and threshold: cook (0.9).
- For variable: (Whole model)
Código 
performance::check_outliers(model_delta) |> plot()

2.1.3.7 Coeficiente de determinação

\(R^2\) - proporção da variância explicada pelo modelo.

Código 
performance::r2(model_delta)
# R2 for Linear Regression
       R2: 0.400
  adj. R2: 0.381

2.1.3.8 Tamanho amostral efetivo (casos completos)

Código 
tibble(n_modelo = nobs(model_delta))




2.2 Desfecho Binário: Melhora Clínica

Análise da melhora clínica, definida como redução ≥ 5 graus na maior curva (\(\Delta \le -5°\)).

2.2.1 Modelo Logístico

Estima associações ajustadas com a chance de melhora clínica, definida como \(\Delta \le -5°\).

Código 
# =============================================================================
# MODELO LOGÍSTICO - MELHORA CLÍNICA
# =============================================================================

# Modelo para chance de melhora clínica (delta ≤ -5°)
# Utiliza as mesmas variáveis do modelo linear

model_melhora <- glm(
    delta_cat ~
        # Variáveis numéricas
        idade + # idade do paciente
        imc + # índice de massa corporal
        cifose_toracica + # ângulo da cifose torácica
        lordose_lombar + # ângulo da lordose lombar
        correcao_colete + # percentual de correção com colete

        # Variáveis categóricas
        sexo + # sexo do paciente
        lenke + # classificação de Lenke
        risser + # sinal de Risser
        flexibilidade + # flexibilidade da curva
        escoliometro_maior_10_graus, # região com gibosidade > 10°
    data = df,
    family = "binomial" # distribuição binomial para regressão logística
)

2.2.1.1 Parâmetros do modelo

Parâmetros do modelo de regressão logística. Resposta em Odds Ratio.

Código 
model_melhora |>
    gtsummary::tbl_regression(exponentiate = TRUE, conf.level = 0.95) |>
    # gtsummary::add_global_p() |>
    gtsummary::bold_p() |>
    gtsummary::bold_labels() |>
    gtsummary::italicize_levels()
Characteristic OR 95% CI p-value
idade 0.93 0.78, 1.10 0.4
imc 1.05 0.98, 1.13 0.15
cifose_toracica 1.02 1.00, 1.04 0.049
lordose_lombar 0.99 0.97, 1.01 0.3
correcao_colete 1.10 1.08, 1.12 <0.001
sexo


    feminino — —
    masculino 3.05 1.44, 6.64 0.004
lenke


    1 — —
    2 0.41 0.12, 1.38 0.2
    3 0.17 0.06, 0.44 <0.001
    4 0.13 0.04, 0.40 <0.001
    5 0.37 0.11, 1.18 0.10
    6 0.20 0.07, 0.53 0.002
risser


    0 — —
    1 6.16 2.90, 13.5 <0.001
    2 2.34 1.16, 4.77 0.019
    3 1.19 0.59, 2.41 0.6
    4 1.96 0.88, 4.42 0.10
flexibilidade


    flexivel — —
    rigido 0.61 0.39, 0.96 0.034
escoliometro_maior_10_graus


    normal — —
    lombar 1.03 0.58, 1.82 >0.9
    toracica 0.96 0.50, 1.84 >0.9
    toracica_lombar 0.31 0.12, 0.76 0.011
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

2.2.1.2 Resíduos binarizados

Código 
performance::binned_residuals(model_melhora) |> plot()

2.2.1.3 Autocorrelação dos resíduos

Independência dos resíduos (Durbin-Watson test).

Código 
performance::check_autocorrelation(model_melhora)
OK: Residuals appear to be independent and not autocorrelated (p = 0.070).

2.2.1.4 Verificações preditivas posteriores

Simulação de dados replicados sob o modelo ajustado e posterior comparação com os dados observados.

Código 
performance::check_predictions(model_melhora) |> plot()

2.2.1.5 Qualidade do ajuste do modelo

Teste de Hosmer-Lemeshow para avaliar a qualidade do ajuste de modelos binomiais. P-valor < 0.05 indica que o modelo não ajusta bem os dados.

Código 
# Teste de Hosmer-Lemeshow pode falhar com poucos eventos de melhora
tryCatch(
    {
        performance::performance_hosmer(model_melhora)
    },
    error = function(e) {
        cat("Aviso: Teste de Hosmer-Lemeshow não pode ser executado.\n")
        cat("Possível causa: poucos eventos de melhora ou falta de variação nas predi<U+00E7><U+00F5>es.\n")
        cat("Erro:", e$message, "\n")
    }
)
# Hosmer-Lemeshow Goodness-of-Fit Test

  Chi-squared: 18.657
           df:  8    
      p-value:  0.017

2.2.1.6 Colinearidade

Verificação de colinearidade entre as variáveis independentes (Variance Inflation Factor). Valores maiores que 10 indicam colinearidade entre as variáveis.

Código 
performance::check_collinearity(model_melhora) |>
    arrange(desc(VIF)) |>
    select(Term, VIF, VIF_CI_low, VIF_CI_high) |>
    performance::print_html()
Term VIF VIF_CI_low VIF_CI_high
lenke 2.45 2.18 2.78
escoliometro_maior_10_graus 2.33 2.08 2.64
risser 2.22 1.98 2.51
idade 1.90 1.71 2.14
cifose_toracica 1.45 1.33 1.62
lordose_lombar 1.35 1.24 1.50
sexo 1.24 1.15 1.38
correcao_colete 1.15 1.08 1.29
imc 1.13 1.06 1.27
flexibilidade 1.12 1.06 1.27
Código 
performance::check_collinearity(model_melhora) |> plot()

2.2.1.7 Outliers

Código 
performance::check_outliers(model_melhora)
OK: No outliers detected.
- Based on the following method and threshold: cook (0.9).
- For variable: (Whole model)
Código 
performance::check_outliers(model_melhora) |> plot()

2.2.1.8 Coeficiente de determinação (Tjur)

Coeficiente de determinação de Tjur (\(R^2_{Tjur}\)). Teste específico para modelos binomiais.

Código 
performance::r2(model_melhora)
# R2 for Logistic Regression
  Tjur's R2: 0.403

2.2.1.9 Tamanho amostral efetivo (casos completos)

Código 
mf_melhora <- model.frame(model_melhora)
n_melhora <- nrow(mf_melhora)
eventos_melhora <- sum(mf_melhora[[1]] == 1, na.rm = TRUE)
sem_melhora <- sum(mf_melhora[[1]] == 0, na.rm = TRUE)
tibble(
    n_modelo = n_melhora,
    eventos_melhora = eventos_melhora,
    sem_melhora = sem_melhora
)




2.2.2 Modelo CART (Árvore de Decisão)

Objetivo: Identificar perfis de pacientes com maior probabilidade de melhora clínica através de árvores de decisão.

Disclaimer: Resultados de modelos CART requerem validação externa. Análise baseada em validação cruzada.

2.2.2.1 Setup e preparação dos dados

Código 
# =============================================================================
# MODELO 1: CART - PREPARAÇÃO DOS DADOS
# =============================================================================

# Carregamento de bibliotecas específicas para machine learning
library(tidymodels) # framework tidyverse para ML
library(rpart.plot) # visualização de árvores
library(vip) # importância de variáveis

set.seed(123) # reprodutibilidade

# Preparação do dataset para CART
# Recodificação do desfecho para labels descritivos
df_cart <- df |>
    mutate(
        delta_cat_cart = factor(
            delta_cat,
            levels = c(0, 1),
            labels = c("nao_melhora", "melhora")
        )
    ) |>
    # Seleção das mesmas variáveis dos modelos parametricos
    select(
        delta_cat_cart, # desfecho recodificado

        # Variáveis preditoras numéricas
        idade,
        imc,
        cifose_toracica,
        lordose_lombar,
        correcao_colete,

        # Variáveis preditoras categóricas
        sexo,
        lenke,
        risser,
        flexibilidade,
        escoliometro_maior_10_graus
    )

# Verificação: evento de interesse = "melhora" (segundo nível)
levels(df_cart$delta_cat_cart)
[1] "nao_melhora" "melhora"

2.2.2.2 Recipe (pre-processamento)

O pre-processamento é feito dentro do resampling (evita leakage):

  • step_zv(): remove variáveis com variância zero
  • step_unknown() + step_novel(): robustez a níveis novos/raros durante CV
  • step_impute_median() + step_impute_mode(): imputação de missings
Código 
cart_rec <- recipe(delta_cat_cart ~ ., data = df_cart) |>
    step_zv(all_predictors()) |>
    step_unknown(all_nominal_predictors()) |>
    step_novel(all_nominal_predictors()) |>
    step_impute_median(all_numeric_predictors()) |>
    step_impute_mode(all_nominal_predictors())

cart_rec

2.2.2.3 Modelo e hiperparâmetros

Hiperparâmetros do modelo CART:

  • cost_complexity (cp): controla poda/complexidade do modelo.
  • tree_depth: limita profundidade (3-5 níveis).
  • min_n: mínimo de observações por nó terminal (15-40).
Código 
cart_spec <- decision_tree(
    cost_complexity = tune(),
    tree_depth = tune(),
    min_n = tune()
) |>
    set_engine("rpart") |>
    set_mode("classification")

cart_wf <- workflow() |>
    add_recipe(cart_rec) |>
    add_model(cart_spec)

cart_wf
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: decision_tree()

── Preprocessor ────────────────────────────────────────────────────────────────
5 Recipe Steps

• step_zv()
• step_unknown()
• step_novel()
• step_impute_median()
• step_impute_mode()

── Model ───────────────────────────────────────────────────────────────────────
Decision Tree Model Specification (classification)

Main Arguments:
  cost_complexity = tune()
  tree_depth = tune()
  min_n = tune()

Computational engine: rpart

2.2.2.4 Validação cruzada e tuning

Usamos CV estratificada repetida (10-fold, 5 repetições) para estimar desempenho e selecionar hiperparâmetros.

Código 
set.seed(123)
folds <- vfold_cv(df_cart, v = 10, repeats = 5, strata = delta_cat_cart)

# Metricas de interesse
metricas <- metric_set(roc_auc, accuracy, sens, spec)

# Grid de hiperparametros
set.seed(123)
cart_grid <- grid_latin_hypercube(
    cost_complexity(range = c(-4, -2)), # menos poda extrema
    tree_depth(range = c(3L, 5L)), # árvores menos profundas
    min_n(range = c(15L, 40L)), # mínimo mais alto por nó
    size = 30
)

# Tuning
set.seed(123)
cart_tuned <- tune_grid(
    cart_wf,
    resamples = folds,
    grid = cart_grid,
    metrics = metricas,
    control = control_grid(save_pred = TRUE)
)

2.2.2.5 Melhores hiperparâmetros

Selecionamos pelo maior AUC ROC medio na CV:

Código 
# Top 10 combinacoes por AUC
cart_tuned |>
    collect_metrics() |>
    filter(.metric == "roc_auc") |>
    arrange(desc(mean)) |>
    slice_head(n = 10)
Código 
# Melhor combinacao
best_params <- cart_tuned |>
    select_best(metric = "roc_auc")

best_params

2.2.2.6 Desempenho final (CV)

Avaliamos o modelo finalizado com os melhores hiperparâmetros:

Código 
cart_final_wf <- finalize_workflow(cart_wf, best_params)

set.seed(123)
cart_res <- fit_resamples(
    cart_final_wf,
    resamples = folds,
    metrics = metricas,
    control = control_resamples(save_pred = TRUE)
)

# Metricas agregadas (media e SD)
cart_metrics <- cart_res |>
    collect_metrics() |>
    select(.metric, mean, std_err, n) |>
    mutate(
        ci_lower = mean - 1.96 * std_err,
        ci_upper = mean + 1.96 * std_err
    )

cart_metrics

2.2.2.7 Métricas por fold (dispersão)

Visualização da incerteza nas estimativas de desempenho:

Código 
cart_res |>
    collect_metrics(summarize = FALSE) |>
    ggplot(aes(x = .metric, y = .estimate, fill = .metric)) +
    geom_boxplot(alpha = 0.7, show.legend = FALSE) +
    geom_jitter(width = 0.1, alpha = 0.3, size = 1) +
    scale_fill_brewer(palette = "Set2") +
    labs(
        title = "Distribution of metrics by fold (CV repeated)",
        subtitle = "10-fold CV x 5 repetitions = 50 estimates per metric",
        x = "Metric",
        y = "Value"
    ) +
    theme_minimal(base_size = 12) +
    theme(
        plot.title = element_text(face = "bold"),
        axis.text.x = element_text(angle = 0)
    )

2.2.2.8 Matriz de confusão agregada (normalizada)

Código 
# Predicoes agregadas de todos os folds
preds_all <- cart_res |>
    collect_predictions()

# Matriz de confusao
conf_mat_data <- preds_all |>
    conf_mat(truth = delta_cat_cart, estimate = .pred_class)

# Normalizar por linha (proporcao por classe verdadeira)
conf_mat_tbl <- conf_mat_data$table |>
    as.data.frame() |>
    group_by(Truth) |>
    mutate(
        prop = Freq / sum(Freq),
        label = scales::percent(prop, accuracy = 0.1)
    ) |>
    ungroup()

# Visualizacao heatmap normalizado
ggplot(conf_mat_tbl, aes(x = Prediction, y = Truth, fill = prop)) +
    geom_tile(color = "white", linewidth = 1) +
    geom_text(aes(label = label), size = 5, fontface = "bold") +
    scale_fill_gradient(
        low = "white",
        high = "steelblue",
        labels = scales::percent,
        name = "Proporcao"
    ) +
    scale_y_discrete(limits = rev) +
    labs(
        title = "Confucion Matrix",
        subtitle = "Proportion by true class (rows sum to 100%)",
        x = "Prediction",
        y = "True"
    ) +
    theme_minimal(base_size = 12) +
    theme(
        plot.title = element_text(face = "bold"),
        panel.grid = element_blank()
    )

2.2.2.9 Árvore final

Ajustamos o modelo final no dataset completo para visualização e interpretação:

Código 
cart_fit_full <- fit(cart_final_wf, data = df_cart)
cart_rpart <- extract_fit_engine(cart_fit_full)

# Plot da arvore com rpart.plot
rpart.plot(
    cart_rpart,
    type = 4,
    extra = 104,
    under = TRUE,
    fallen.leaves = TRUE,
    roundint = FALSE,
    box.palette = "BuGn",
    shadow.col = "gray",
    main = "Decision Tree for MCID (delta <= -5 degrees)"
)

2.2.2.10 Importância de variáveis

Código 
# Grafico de importancia com vip
cart_fit_full |>
    extract_fit_parsnip() |>
    vip(
        num_features = 15,
        geom = "col",
        aesthetics = list(fill = "steelblue", alpha = 0.8)
    ) +
    labs(
        title = "Variable Importance (CART)",
        subtitle = "Based on impurity reduction (Gini)",
        x = "Importance",
        y = "Variable"
    ) +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(face = "bold"))

2.2.2.11 Análise de estabilidade (Bootstrap)

Código 
# Função para extrair variáveis de split principais de uma árvore (únicas por árvore)
extract_main_splits <- function(rpart_model, max_levels = 3) {
    frame <- rpart_model$frame
    # Filtrar apenas nós internos (não folhas)
    internal_nodes <- frame[frame$var != "<leaf>", ]

    # Extrair apenas as variáveis (não os valores de corte)
    if (nrow(internal_nodes) > 0) {
        # Limitar aos primeiros níveis (mais estáveis)
        rownames_numeric <- as.numeric(rownames(internal_nodes))
        main_splits <- internal_nodes[rownames_numeric <= (2^max_levels - 1), ]
        # Retornar apenas variáveis ÚNICAS por árvore
        return(unique(main_splits$var))
    } else {
        return(character(0))
    }
}

# Bootstrap para verificar consistência dos splits
set.seed(123)
n_bootstrap <- 100
bootstrap_splits <- list()

for (i in 1:n_bootstrap) {
    # Amostra bootstrap
    boot_sample <- df_cart[sample(nrow(df_cart), replace = TRUE), ]

    # Ajustar modelo com mesmos hiperparâmetros
    boot_model <- decision_tree(
        cost_complexity = best_params$cost_complexity,
        tree_depth = best_params$tree_depth,
        min_n = best_params$min_n
    ) |>
        set_engine("rpart") |>
        set_mode("classification") |>
        fit(delta_cat_cart ~ ., data = boot_sample)

    # Extrair splits principais
    boot_rpart <- extract_fit_engine(boot_model)
    bootstrap_splits[[i]] <- extract_main_splits(boot_rpart)
}

# Análise de estabilidade
main_vars_stability <- table(unlist(bootstrap_splits))
stability_summary <- sort(main_vars_stability / n_bootstrap, decreasing = TRUE)

cat("Estabilidade das variáveis nos primeiros splits (% de árvores que usam cada variável):\n")
Estabilidade das variáveis nos primeiros splits (% de árvores que usam cada variável):
Código 
print(round(stability_summary * 100, 1))

            correcao_colete                      risser 
                        100                          75 
                        imc             cifose_toracica 
                         43                          35 
                      lenke                       idade 
                         29                          21 
             lordose_lombar escoliometro_maior_10_graus 
                         18                          15 
                       sexo               flexibilidade 
                          8                           2

2.2.2.12 Tabela de folhas (grupos de predicao)

Esta tabela mostra os nós terminais (folhas) da árvore, com as regras que definem cada grupo, o tamanho do grupo e a proporção de melhora:

Código 
# Extrair informacoes dos nos terminais
frame <- cart_rpart$frame
leaves_idx <- which(frame$var == "<leaf>")

# Funcao para extrair regras de cada no
get_path_rules <- function(tree, node_id) {
    path <- rpart::path.rpart(tree, nodes = node_id, print.it = FALSE)
    rules <- path[[1]][-1] # remover "root"
    if (length(rules) == 0) {
        return("(raiz)")
    }
    paste(rules, collapse = " E ")
}

# Construir tabela de folhas
leaves_tbl <- tibble(
    no = as.integer(rownames(frame)[leaves_idx]),
    n = frame$n[leaves_idx],
    n_melhora = frame$n[leaves_idx] * frame$yval2[leaves_idx, 5],
    n_nao_melhora = frame$n[leaves_idx] * frame$yval2[leaves_idx, 4],
    prop_melhora = frame$yval2[leaves_idx, 5],
    predicao = ifelse(frame$yval[leaves_idx] == 2, "melhora", "nao_melhora")
) |>
    rowwise() |>
    mutate(regras = get_path_rules(cart_rpart, no)) |>
    ungroup() |>
    mutate(
        prop_melhora_pct = scales::percent(prop_melhora, accuracy = 0.1),
        grupo = paste0("Grupo ", row_number())
    ) |>
    select(grupo, predicao, n, prop_melhora_pct, regras) |>
    arrange(desc(predicao), desc(n))

leaves_tbl |>
    rename(
        Grupo = grupo,
        Predicao = predicao,
        N = n,
        `% Melhora` = prop_melhora_pct,
        Regras = regras
    ) |>
    knitr::kable(align = c("l", "l", "r", "r", "l"))
Grupo Predicao N % Melhora Regras
Grupo 2 nao_melhora 195 23.6% correcao_colete< 49.42 E correcao_colete>=29.45 E risser=0,2,3,4 E lenke=2,3,4,6 E sexo=feminino
Grupo 1 nao_melhora 66 3.0% correcao_colete< 49.42 E correcao_colete< 29.45
Grupo 6 nao_melhora 11 27.3% correcao_colete< 49.42 E correcao_colete>=29.45 E risser=1 E idade< 11.5
Grupo 4 nao_melhora 10 0.0% correcao_colete< 49.42 E correcao_colete>=29.45 E risser=0,2,3,4 E lenke=1,5 E lordose_lombar< 45
Grupo 7 nao_melhora 9 44.4% correcao_colete< 49.42 E correcao_colete>=29.45 E risser=1 E idade>=11.5 E cifose_toracica< 11
Grupo 9 nao_melhora 9 22.2% correcao_colete>=49.42 E correcao_colete< 61.72 E escoliometro_maior_10_graus=toracica_lombar
Grupo 10 nao_melhora 9 11.1% correcao_colete>=49.42 E correcao_colete< 61.72 E escoliometro_maior_10_graus=normal,lombar,toracica E imc< 15.84 E cifose_toracica< 19.5
Grupo 13 melhora 115 96.5% correcao_colete>=49.42 E correcao_colete>=61.72
Grupo 12 melhora 99 80.8% correcao_colete>=49.42 E correcao_colete< 61.72 E escoliometro_maior_10_graus=normal,lombar,toracica E imc>=15.84
Grupo 5 melhora 33 72.7% correcao_colete< 49.42 E correcao_colete>=29.45 E risser=0,2,3,4 E lenke=1,5 E lordose_lombar>=45
Grupo 8 melhora 26 84.6% correcao_colete< 49.42 E correcao_colete>=29.45 E risser=1 E idade>=11.5 E cifose_toracica>=11
Grupo 3 melhora 25 64.0% correcao_colete< 49.42 E correcao_colete>=29.45 E risser=0,2,3,4 E lenke=2,3,4,6 E sexo=masculino
Grupo 11 melhora 14 64.3% correcao_colete>=49.42 E correcao_colete< 61.72 E escoliometro_maior_10_graus=normal,lombar,toracica E imc< 15.84 E cifose_toracica>=19.5

2.2.2.13 Curva ROC

Código 
# Curva ROC a partir das predicoes CV
roc_data <- preds_all |>
    roc_curve(truth = delta_cat_cart, .pred_melhora, event_level = "second")

# AUC para o titulo
auc_val <- preds_all |>
    roc_auc(truth = delta_cat_cart, .pred_melhora, event_level = "second") |>
    pull(.estimate)

ggplot(roc_data, aes(x = 1 - specificity, y = sensitivity)) +
    geom_path(linewidth = 1.2, color = "steelblue") +
    geom_abline(linetype = "dashed", color = "gray50") +
    geom_ribbon(
        aes(ymin = 0, ymax = sensitivity),
        alpha = 0.1,
        fill = "steelblue"
    ) +
    annotate(
        "text",
        x = 0.6,
        y = 0.3,
        label = paste0("AUC = ", round(auc_val, 3)),
        size = 5,
        fontface = "bold"
    ) +
    labs(
        title = "ROC Curve (CV)",
        subtitle = "Based on aggregated predictions from all folds",
        x = "1 - Specificity (False Positive Rate)",
        y = "Sensitivity (True Positive Rate)"
    ) +
    coord_equal() +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(face = "bold"))




2.3 Desfecho Binário: Progressão Clínica

Análise da progressão clínica, definida como piora ≥ 5 graus na maior curva (\(\Delta \ge +5°\)).

2.3.1 Modelo Logístico

Estima associações ajustadas com a chance de progressão clínica, definida como \(\Delta \ge +5°\).

Código 
# =============================================================================
# MODELO LOGÍSTICO - PROGRESS<U+00C3>O CLÍNICA
# =============================================================================

# Modelo para chance de progressão clínica (delta ≥ +5°)
# Utiliza as mesmas variáveis dos outros modelos

model_progressao <- glm(
    delta_progression ~
        # Variáveis numéricas
        idade + # idade do paciente
        imc + # índice de massa corporal
        cifose_toracica + # ângulo da cifose torácica
        lordose_lombar + # ângulo da lordose lombar
        correcao_colete + # percentual de correção com colete

        # Variáveis categóricas
        sexo + # sexo do paciente
        lenke + # classificação de Lenke
        risser + # sinal de Risser
        flexibilidade + # flexibilidade da curva
        escoliometro_maior_10_graus, # região com gibosidade > 10°
    data = df,
    family = "binomial" # distribuição binomial para regressão logística
)

2.3.1.1 Parâmetros do modelo

Parâmetros do modelo de regressão logística. Resposta em Odds Ratio.

Código 
model_progressao |>
    gtsummary::tbl_regression(exponentiate = TRUE, conf.level = 0.95) |>
    gtsummary::bold_p() |>
    gtsummary::bold_labels() |>
    gtsummary::italicize_levels()
Characteristic OR 95% CI p-value
idade 0.85 0.56, 1.29 0.4
imc 1.10 0.91, 1.31 0.3
cifose_toracica 0.90 0.83, 0.96 0.003
lordose_lombar 1.01 0.95, 1.06 0.8
correcao_colete 0.96 0.92, 0.99 0.033
sexo


    feminino — —
    masculino 3.53 0.53, 21.7 0.2
lenke


    1 — —
    2 1.30 0.08, 39.4 0.9
    3 1.17 0.14, 26.6 >0.9
    4 3.29 0.35, 79.0 0.4
    5 0.00 0.00, 419,627,348,273,439,586,821,317,894,482,929,123,328 >0.9
    6 0.11 0.00, 3.82 0.2
risser


    0 — —
    1 0.35 0.07, 1.43 0.2
    2 0.08 0.01, 0.42 0.006
    3 0.00 0.00, 4,840,904,359,509,076 >0.9
    4 0.24 0.03, 1.31 0.12
flexibilidade


    flexivel — —
    rigido 3.35 1.08, 11.9 0.046
escoliometro_maior_10_graus


    normal — —
    lombar 0.96 0.03, 17.2 >0.9
    toracica 0.94 0.21, 5.33 >0.9
    toracica_lombar 7.49 1.10, 62.8 0.046
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

2.3.1.2 Resíduos binarizados

Código 
performance::binned_residuals(model_progressao) |> plot()

2.3.1.3 Autocorrelação dos resíduos

Independência dos resíduos (Durbin-Watson test).

Código 
performance::check_autocorrelation(model_progressao)
OK: Residuals appear to be independent and not autocorrelated (p = 0.786).

2.3.1.4 Verificações preditivas posteriores

Simulação de dados replicados sob o modelo ajustado e posterior comparação com os dados observados.

Código 
performance::check_predictions(model_progressao) |> plot()

2.3.1.5 Qualidade do ajuste do modelo

Teste de Hosmer-Lemeshow para avaliar a qualidade do ajuste de modelos binomiais. P-valor < 0.05 indica que o modelo não ajusta bem os dados.

Código 
# Teste de Hosmer-Lemeshow pode falhar com poucos eventos de progressão
tryCatch(
    {
        performance::performance_hosmer(model_progressao)
    },
    error = function(e) {
        cat("Aviso: Teste de Hosmer-Lemeshow não pode ser executado.\n")
        cat("Possível causa: poucos eventos de progressão ou falta de variação nas predições.\n")
        cat("Erro:", e$message, "\n")
    }
)
# Hosmer-Lemeshow Goodness-of-Fit Test

  Chi-squared: 1.030
           df: 8    
      p-value: 0.998

2.3.1.6 Colinearidade

Verificação de colinearidade entre as variáveis independentes (Variance Inflation Factor). Valores maiores que 10 indicam colinearidade entre as variáveis.

Código 
performance::check_collinearity(model_progressao) |>
    arrange(desc(VIF)) |>
    select(Term, VIF, VIF_CI_low, VIF_CI_high) |>
    performance::print_html()
Term VIF VIF_CI_low VIF_CI_high
lenke 3.89 3.42 4.45
escoliometro_maior_10_graus 3.32 2.93 3.79
risser 2.45 2.18 2.78
idade 2.12 1.90 2.40
cifose_toracica 1.58 1.43 1.77
sexo 1.58 1.43 1.76
lordose_lombar 1.49 1.36 1.66
flexibilidade 1.34 1.23 1.49
imc 1.25 1.16 1.39
correcao_colete 1.23 1.14 1.38
Código 
performance::check_collinearity(model_progressao) |> plot()

2.3.1.7 Outliers

Código 
performance::check_outliers(model_progressao)
OK: No outliers detected.
- Based on the following method and threshold: cook (0.9).
- For variable: (Whole model)
Código 
performance::check_outliers(model_progressao) |> plot()

2.3.1.8 Coeficiente de determinação (Tjur)

Coeficiente de determinação de Tjur (\(R^2_{Tjur}\)). Teste específico para modelos binomiais.

Código 
performance::r2(model_progressao)
# R2 for Logistic Regression
  Tjur's R2: 0.345

2.3.1.9 Tamanho amostral efetivo (casos completos)

Código 
mf_progressao <- model.frame(model_progressao)
n_progressao <- nrow(mf_progressao)
eventos_progressao <- sum(mf_progressao[[1]] == 1, na.rm = TRUE)
sem_progressao <- sum(mf_progressao[[1]] == 0, na.rm = TRUE)
tibble(
    n_modelo = n_progressao,
    eventos_progressao = eventos_progressao,
    sem_progressao = sem_progressao
)

2.3.2 Modelo CART (Árvore de Decisão)

Objetivo: Identificar perfis de pacientes com maior probabilidade de progressão clínica através de árvores de decisão.

Disclaimer: Resultados de modelos CART requerem validação externa. Análise baseada em validação cruzada.

2.3.2.1 Setup e preparação dos dados

Código 
# =============================================================================
# CART - PROGRESS<U+00C3>O CLÍNICA: PREPARAÇÃO DOS DADOS
# =============================================================================

# Carregamento de bibliotecas específicas para machine learning
library(tidymodels) # framework tidyverse para ML
library(rpart.plot) # visualização de árvores
library(vip) # importância de variáveis

set.seed(123) # reprodutibilidade

# Preparação do dataset para CART
# Recodificação do desfecho para labels descritivos
df_cart_progressao <- df |>
    mutate(
        delta_progression_cart = factor(
            delta_progression,
            levels = c(0, 1),
            labels = c("nao_progrediu", "progrediu")
        )
    ) |>
    # Seleção das mesmas variáveis dos modelos parametricos
    select(
        delta_progression_cart, # desfecho recodificado

        # Variáveis preditoras numéricas
        idade,
        imc,
        cifose_toracica,
        lordose_lombar,
        correcao_colete,

        # Variáveis preditoras categóricas
        sexo,
        lenke,
        risser,
        flexibilidade,
        escoliometro_maior_10_graus
    )

# Verificação: evento de interesse = "progrediu" (segundo nível)
levels(df_cart_progressao$delta_progression_cart)
[1] "nao_progrediu" "progrediu"

2.3.2.2 Recipe (pre-processamento)

O pre-processamento é feito dentro do resampling (evita leakage):

  • step_zv(): remove variáveis com variância zero
  • step_unknown() + step_novel(): robustez a níveis novos/raros durante CV
  • step_impute_median() + step_impute_mode(): imputação de missings
Código 
cart_rec_progressao <- recipe(delta_progression_cart ~ ., data = df_cart_progressao) |>
    step_zv(all_predictors()) |>
    step_unknown(all_nominal_predictors()) |>
    step_novel(all_nominal_predictors()) |>
    step_impute_median(all_numeric_predictors()) |>
    step_impute_mode(all_nominal_predictors())

cart_rec_progressao

2.3.2.3 Modelo e hiperparâmetros

Hiperparâmetros do modelo CART:

  • cost_complexity (cp): controla poda/complexidade do modelo.
  • tree_depth: limita profundidade (3-5 níveis).
  • min_n: mínimo de observações por nó terminal (15-40).
Código 
cart_spec_progressao <- decision_tree(
    cost_complexity = tune(),
    tree_depth = tune(),
    min_n = tune()
) |>
    set_engine("rpart") |>
    set_mode("classification")

cart_wf_progressao <- workflow() |>
    add_recipe(cart_rec_progressao) |>
    add_model(cart_spec_progressao)

cart_wf_progressao
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: decision_tree()

── Preprocessor ────────────────────────────────────────────────────────────────
5 Recipe Steps

• step_zv()
• step_unknown()
• step_novel()
• step_impute_median()
• step_impute_mode()

── Model ───────────────────────────────────────────────────────────────────────
Decision Tree Model Specification (classification)

Main Arguments:
  cost_complexity = tune()
  tree_depth = tune()
  min_n = tune()

Computational engine: rpart

2.3.2.4 Validação cruzada e tuning

Usamos CV estratificada repetida (10-fold, 5 repetições) para estimar desempenho e selecionar hiperparâmetros.

Código 
set.seed(123)
folds_progressao <- vfold_cv(df_cart_progressao, v = 10, repeats = 5, strata = delta_progression_cart)

# Metricas de interesse
metricas <- metric_set(roc_auc, accuracy, sens, spec)

# Grid de hiperparametros
set.seed(123)
cart_grid_progressao <- grid_latin_hypercube(
    cost_complexity(range = c(-4, -2)), # menos poda extrema
    tree_depth(range = c(3L, 5L)), # árvores menos profundas
    min_n(range = c(15L, 40L)), # mínimo mais alto por nó
    size = 30
)

# Tuning
set.seed(123)
cart_tuned_progressao <- tune_grid(
    cart_wf_progressao,
    resamples = folds_progressao,
    grid = cart_grid_progressao,
    metrics = metricas,
    control = control_grid(save_pred = TRUE)
)

2.3.2.5 Melhores hiperparâmetros

Selecionamos pelo maior AUC ROC medio na CV:

Código 
# Top 10 combinacoes por AUC
cart_tuned_progressao |>
    collect_metrics() |>
    filter(.metric == "roc_auc") |>
    arrange(desc(mean)) |>
    slice_head(n = 10)
Código 
# Melhor combinacao
best_params_progressao <- cart_tuned_progressao |>
    select_best(metric = "roc_auc")

best_params_progressao

2.3.2.6 Desempenho final (CV)

Avaliamos o modelo finalizado com os melhores hiperparâmetros:

Código 
cart_final_wf_progressao <- finalize_workflow(cart_wf_progressao, best_params_progressao)

set.seed(123)
cart_res_progressao <- fit_resamples(
    cart_final_wf_progressao,
    resamples = folds_progressao,
    metrics = metricas,
    control = control_resamples(save_pred = TRUE)
)

# Metricas agregadas (media e SD)
cart_metrics_progressao <- cart_res_progressao |>
    collect_metrics() |>
    select(.metric, mean, std_err, n) |>
    mutate(
        ci_lower = mean - 1.96 * std_err,
        ci_upper = mean + 1.96 * std_err
    )

cart_metrics_progressao

2.3.2.7 Métricas por fold (dispersão)

Variabilidade das métricas entre os folds de validação cruzada.

Código 
cart_res_progressao |>
    collect_metrics(summarize = FALSE) |>
    ggplot(aes(x = .metric, y = .estimate, fill = .metric)) +
    geom_boxplot(alpha = 0.7, show.legend = FALSE) +
    geom_jitter(width = 0.1, alpha = 0.3, size = 1) +
    scale_fill_brewer(palette = "Set2") +
    labs(
        title = "Distribution of metrics by fold (CV repeated)",
        subtitle = "10-fold CV x 5 repetitions = 50 estimates per metric",
        x = "Metric",
        y = "Value"
    ) +
    theme_minimal(base_size = 12) +
    theme(
        plot.title = element_text(face = "bold"),
        axis.text.x = element_text(angle = 0)
    )

2.3.2.8 Matriz de confusão agregada (normalizada)

Código 
# Matriz de confusão agregada
preds_progressao <- cart_res_progressao |>
    collect_predictions()

# Predicoes agregadas de todos os folds (mantendo fonte de dados do contexto do documento)
preds_all_progressao <- cart_res_progressao |>
    collect_predictions()

# Matriz de confusao
conf_mat_data_progressao <- preds_all_progressao |>
    conf_mat(truth = delta_progression_cart, estimate = .pred_class)

# Normalizar por linha (proporcao por classe verdadeira)
conf_mat_tbl_progressao <- conf_mat_data_progressao$table |>
    as.data.frame() |>
    dplyr::group_by(Truth) |>
    dplyr::mutate(
        prop = Freq / sum(Freq),
        label = scales::percent(prop, accuracy = 0.1)
    ) |>
    dplyr::ungroup()

# Visualizacao heatmap normalizado
ggplot(conf_mat_tbl_progressao, aes(x = Prediction, y = Truth, fill = prop)) +
    geom_tile(color = "white", linewidth = 1) +
    geom_text(aes(label = label), size = 5, fontface = "bold") +
    scale_fill_gradient(
        low = "white",
        high = "steelblue",
        labels = scales::percent,
        name = "Proporcao"
    ) +
    scale_y_discrete(limits = rev) +
    labs(
        title = "Matriz de Confusão - Progressão (Normalizada por linha)",
        subtitle = "Agregada de todos os folds de validação cruzada"
    ) +
    theme_minimal(base_size = 12) +
    theme(
        plot.title = element_text(face = "bold"),
        panel.grid = element_blank()
    )

2.3.2.9 Árvore final

Ajustamos o modelo final no dataset completo para visualização e interpretação:

Código 
cart_fit_full_progressao <- fit(cart_final_wf_progressao, data = df_cart_progressao)
cart_rpart_progressao <- extract_fit_engine(cart_fit_full_progressao)

# Plot da arvore com rpart.plot
rpart.plot(
    cart_rpart_progressao,
    type = 4,
    extra = 104,
    under = TRUE,
    fallen.leaves = TRUE,
    roundint = FALSE,
    box.palette = "BuGn",
    shadow.col = "gray",
    main = "Decision Tree for Progression (delta >= +5 degrees)"
)

2.3.2.10 Importância de variáveis

Código 
# Verificar se a árvore tem splits antes de calcular importância
tryCatch(
    {
        # Verificar se há splits na árvore
        frame_check <- cart_rpart_progressao$frame
        internal_nodes <- sum(frame_check$var != "<leaf>")

        if (internal_nodes > 0) {
            # Grafico de importancia com vip
            cart_fit_full_progressao |>
                extract_fit_parsnip() |>
                vip(
                    num_features = min(15, internal_nodes),
                    geom = "col",
                    aesthetics = list(fill = "darkred", alpha = 0.8)
                ) +
                labs(
                    title = "Variable Importance (CART - Progressão)",
                    subtitle = "Based on impurity reduction (Gini)",
                    x = "Importance",
                    y = "Variable"
                ) +
                theme_minimal(base_size = 12) +
                theme(plot.title = element_text(face = "bold"))
        } else {
            cat("Aviso: Árvore de progressão não possui splits suficientes para análise de importância.\n")
            cat("Possível causa: poucos eventos de progressão ou parâmetros muito restritivos.\n")

            # Plot informativo
            ggplot() +
                annotate("text",
                    x = 0.5, y = 0.5,
                    label = "Árvore sem splits suficientes\npara análise de importância",
                    size = 6, hjust = 0.5
                ) +
                labs(
                    title = "Variable Importance (CART - Progressão)",
                    subtitle = "Modelo muito simples para análise de importância"
                ) +
                theme_void()
        }
    },
    error = function(e) {
        cat("Erro na análise de importância:", e$message, "\n")
        cat("Possível causa: modelo CART muito simples ou dados insuficientes.\n")
    }
)

2.3.2.11 Análise de estabilidade (Bootstrap)

Código 
# Bootstrap para verificar consistência dos splits
set.seed(123)
n_bootstrap <- 100
bootstrap_splits_progressao <- list()

for (i in 1:n_bootstrap) {
    # Amostra bootstrap
    boot_sample <- df_cart_progressao[sample(nrow(df_cart_progressao), replace = TRUE), ]

    # Ajustar modelo com mesmos hiperparâmetros
    boot_model <- decision_tree(
        cost_complexity = best_params_progressao$cost_complexity,
        tree_depth = best_params_progressao$tree_depth,
        min_n = best_params_progressao$min_n
    ) |>
        set_engine("rpart") |>
        set_mode("classification") |>
        fit(delta_progression_cart ~ ., data = boot_sample)

    # Extrair splits principais
    boot_rpart <- extract_fit_engine(boot_model)
    bootstrap_splits_progressao[[i]] <- extract_main_splits(boot_rpart)
}

# Análise de estabilidade
main_vars_stability_progressao <- table(unlist(bootstrap_splits_progressao))
stability_summary_progressao <- sort(main_vars_stability_progressao / n_bootstrap, decreasing = TRUE)

cat("Estabilidade das variáveis nos primeiros splits (% de árvores que usam cada variável):\n")
Estabilidade das variáveis nos primeiros splits (% de árvores que usam cada variável):
Código 
print(round(stability_summary_progressao * 100, 1))

                     risser                       lenke 
                         66                          58 
escoliometro_maior_10_graus             cifose_toracica 
                         54                          31 
            correcao_colete                         imc 
                         31                          25 
                      idade              lordose_lombar 
                         12                          12 
                       sexo 
                          1

2.3.2.12 Tabela de folhas (grupos de predicao)

Esta tabela mostra os nós terminais (folhas) da árvore, com as regras que definem cada grupo, o tamanho do grupo e a proporção de progressão:

Código 
# Extrair informacoes dos nos terminais

frame <- cart_rpart_progressao$frame
leaves_idx <- which(frame$var == "<leaf>")

# Funcao para extrair regras de cada no
get_path_rules <- function(tree, node_id) {
    path <- rpart::path.rpart(tree, nodes = node_id, print.it = FALSE)
    rules <- path[[1]][-1] # remover "root"
    if (length(rules) == 0) {
        return("(raiz)")
    }
    paste(rules, collapse = " E ")
}

# Construir tabela de folhas
leaves_tbl <- tibble(
    no = as.integer(rownames(frame)[leaves_idx]),
    n = frame$n[leaves_idx],
    n_melhora = frame$n[leaves_idx] * frame$yval2[leaves_idx, 5],
    n_nao_melhora = frame$n[leaves_idx] * frame$yval2[leaves_idx, 4],
    prop_melhora = frame$yval2[leaves_idx, 5],
    predicao = ifelse(frame$yval[leaves_idx] == 2, "progrediu", "nao_progrediu")
) |>
    rowwise() |>
    mutate(regras = get_path_rules(cart_rpart_progressao, no)) |>
    ungroup() |>
    mutate(
        prop_melhora_pct = scales::percent(prop_melhora, accuracy = 0.1),
        grupo = paste0("Grupo ", row_number())
    ) |>
    select(grupo, predicao, n, prop_melhora_pct, regras) |>
    arrange(desc(predicao), desc(n))

leaves_tbl |>
    rename(
        Grupo = grupo,
        Predicao = predicao,
        N = n,
        `% Melhora` = prop_melhora_pct,
        Regras = regras
    ) |>
    knitr::kable(align = c("l", "l", "r", "r", "l"))
Grupo Predicao N % Melhora Regras
Grupo 5 progrediu 11 72.7% risser=0 E escoliometro_maior_10_graus=toracica_lombar
Grupo 4 progrediu 6 66.7% risser=0 E escoliometro_maior_10_graus=normal,lombar,toracica E lenke=2,4 E imc>=22.19
Grupo 1 nao_progrediu 489 1.8% risser=1,2,3,4
Grupo 2 nao_progrediu 94 2.1% risser=0 E escoliometro_maior_10_graus=normal,lombar,toracica E lenke=1,3,5,6
Grupo 3 nao_progrediu 21 9.5% risser=0 E escoliometro_maior_10_graus=normal,lombar,toracica E lenke=2,4 E imc< 22.19

2.3.2.13 Curva ROC

Código 
# Curva ROC a partir das predicoes CV
roc_data_progressao <- preds_progressao |>
    roc_curve(truth = delta_progression_cart, .pred_progrediu, event_level = "second")

# AUC para o titulo
auc_val_progressao <- preds_progressao |>
    roc_auc(truth = delta_progression_cart, .pred_progrediu, event_level = "second") |>
    pull(.estimate)

ggplot(roc_data_progressao, aes(x = 1 - specificity, y = sensitivity)) +
    geom_path(linewidth = 1.2, color = "steelblue") +
    geom_abline(linetype = "dashed", color = "gray50") +
    geom_ribbon(
        aes(ymin = 0, ymax = sensitivity),
        alpha = 0.1,
        fill = "steelblue"
    ) +
    annotate(
        "text",
        x = 0.6,
        y = 0.3,
        label = paste0("AUC = ", round(auc_val_progressao, 3)),
        size = 5,
        fontface = "bold"
    ) +
    labs(
        title = "ROC Curve (CV)",
        subtitle = "Based on aggregated predictions from all folds",
        x = "1 - Specificity (False Positive Rate)",
        y = "Sensitivity (True Positive Rate)"
    ) +
    coord_equal() +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(face = "bold"))