Scoliosis Model

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_00.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%).

# =============================================================================
# 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_00.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

# =============================================================================
# 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 escoliometro_cervical             262 97.8%      
 2 cobb_toracico_proximal            219 81.7%      
 3 cobb_lombar                        36 13.4%      
 4 escoliometro_lombar                27 10.1%      
 5 cobb_toracica                      23 8.6%       
 6 flexibilidade                      16 6.0%       
 7 escoliometro_torarica              13 4.9%       
 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 lenke                               0 0.0%       
14 risser                              0 0.0%       
15 cifose_toracica                     0 0.0%       
16 lordose_lombar                      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%       

# 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.

# 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.

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 = 2521
idade	12.96 (1.58)
imc	19.76 (10.87)
cifose_toracica	22 (11)
lordose_lombar	54 (12)
correcao_colete	48 (19)
sexo
feminino	225 (89%)
masculino	27 (11%)
lenke
1	14 (5.6%)
2	18 (7.1%)
3	56 (22%)
4	23 (9.1%)
5	24 (9.5%)
6	117 (46%)
risser
0	54 (21%)
1	43 (17%)
2	65 (26%)
3	49 (19%)
4	41 (16%)
flexibilidade
flexivel	158 (63%)
rigido	94 (37%)
escoliometro_maior_10_graus
normal	52 (21%)
lombar	88 (35%)
toracica	94 (37%)
toracica_lombar	18 (7.1%)
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).

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 = 2521
cobb_inicial_maior	36.7 (5.6)
maior_curva_6_meses	32 (8)
delta	-4.5 (5.5)
delta_cat_f
Sem melhora (MCID)	135 (54%)
Melhora (MCID)	117 (46%)
delta_progression_f
Não progrediu	243 (96%)
Progrediu	9 (3.6%)
1 Mean (SD); n (%)

1.3 Distribuição do desfecho principal

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()
    )

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)

O modelo linear analisa diretamente a mudança na maior curva (\(\Delta\)), fornecendo interpretação intuitiva dos coeficientes:
- Coeficientes negativos: associação com melhora clínica
- Coeficientes positivos: associação com progressão clínica
- Magnitude: quantifica o efeito em graus de mudança na curva

# =============================================================================
# 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.1 Parâmetros do modelo

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.18	-0.67, 0.32	0.5
imc	-0.03	-0.08, 0.03	0.3
cifose_toracica	-0.03	-0.09, 0.03	0.3
lordose_lombar	0.00	-0.05, 0.06	0.9
correcao_colete	-0.09	-0.12, -0.05	<0.001
sexo
feminino	—	—
masculino	-0.38	-2.4, 1.7	0.7
lenke
1	—	—
2	1.7	-1.6, 5.1	0.3
3	5.4	2.6, 8.3	<0.001
4	5.6	2.5, 8.8	<0.001
5	1.8	-1.5, 5.2	0.3
6	3.6	0.70, 6.5	0.015
risser
0	—	—
1	-2.8	-4.8, -0.79	0.006
2	-1.1	-2.9, 0.78	0.3
3	-0.43	-2.4, 1.5	0.7
4	0.27	-2.2, 2.7	0.8
flexibilidade
flexivel	—	—
rigido	1.5	0.25, 2.8	0.020
escoliometro_maior_10_graus
normal	—	—
lombar	-0.74	-2.5, 1.0	0.4
toracica	-0.09	-1.9, 1.8	>0.9
toracica_lombar	1.8	-0.82, 4.4	0.2
Abbreviation: CI = Confidence Interval

2.1.2 Análises de pressupostos

2.1.2.1 Heteroscedasticidade

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

performance::check_heteroscedasticity(model_delta)

Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.019).

performance::check_heteroscedasticity(model_delta) |> plot()

2.1.2.2 Normalidade dos resíduos

performance::check_residuals(model_delta)

OK: Simulated residuals appear as uniformly distributed (p = 0.767).

performance::check_residuals(model_delta) |> plot()

2.1.2.3 Autocorrelação dos resíduos

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

performance::check_autocorrelation(model_delta)

Warning: Autocorrelated residuals detected (p = 0.008).

2.1.2.4 Verificações preditivas posteriores

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

performance::check_predictions(model_delta) |> plot()

2.1.2.5 Colinearidade

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

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.54	2.15	3.08
escoliometro_maior_10_graus	2.42	2.05	2.93
risser	1.97	1.69	2.37
idade	1.83	1.58	2.19
cifose_toracica	1.48	1.31	1.76
lordose_lombar	1.34	1.20	1.60
sexo	1.18	1.08	1.42
flexibilidade	1.16	1.07	1.41
correcao_colete	1.16	1.07	1.41
imc	1.07	1.01	1.43

performance::check_collinearity(model_delta) |> plot()

2.1.2.6 Outliers

performance::check_outliers(model_delta)

OK: No outliers detected.
- Based on the following method and threshold: cook (0.9).
- For variable: (Whole model)

performance::check_outliers(model_delta) |> plot()

2.1.2.7 Coeficiente de determinação

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

performance::r2(model_delta)

# R2 for Linear Regression
       R2: 0.326
  adj. R2: 0.271

2.1.2.8 Tamanho amostral efetivo (casos completos)

tibble(n_modelo = nobs(model_delta))

2.1.3 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

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°\).

# =============================================================================
# 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.

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.97	0.75, 1.25	0.8
imc	1.03	1.0, 1.12	0.3
cifose_toracica	1.02	0.99, 1.06	0.14
lordose_lombar	0.99	0.96, 1.02	0.5
correcao_colete	1.05	1.03, 1.07	<0.001
sexo
feminino	—	—
masculino	2.00	0.71, 5.76	0.2
lenke
1	—	—
2	0.18	0.03, 1.00	0.059
3	0.07	0.01, 0.29	<0.001
4	0.07	0.01, 0.35	0.003
5	0.25	0.04, 1.34	0.12
6	0.18	0.03, 0.76	0.028
risser
0	—	—
1	3.47	1.28, 9.76	0.016
2	1.41	0.56, 3.62	0.5
3	0.73	0.26, 1.99	0.5
4	1.62	0.48, 5.57	0.4
flexibilidade
flexivel	—	—
rigido	0.49	0.25, 0.94	0.034
escoliometro_maior_10_graus
normal	—	—
lombar	1.03	0.44, 2.43	>0.9
toracica	1.47	0.57, 3.92	0.4
toracica_lombar	1.17	0.31, 4.38	0.8
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

2.2.1.2 Resíduos binarizados

performance::binned_residuals(model_melhora) |> plot()

2.2.1.3 Autocorrelação dos resíduos

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

performance::check_autocorrelation(model_melhora)

Warning: Autocorrelated residuals detected (p < .001).

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.

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.

# 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: 6.672
           df: 8    
      p-value: 0.572

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.

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.55	2.15	3.09
escoliometro_maior_10_graus	2.44	2.06	2.94
risser	2.05	1.75	2.47
idade	1.91	1.64	2.29
cifose_toracica	1.48	1.31	1.76
lordose_lombar	1.29	1.16	1.53
sexo	1.27	1.14	1.51
correcao_colete	1.14	1.05	1.39
flexibilidade	1.14	1.05	1.39
imc	1.06	1.01	1.50

performance::check_collinearity(model_melhora) |> plot()

2.2.1.7 Outliers

performance::check_outliers(model_melhora)

OK: No outliers detected.
- Based on the following method and threshold: cook (0.9).
- For variable: (Whole model)

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.

performance::r2(model_melhora)

# R2 for Logistic Regression
  Tjur's R2: 0.273

2.2.1.9 Tamanho amostral efetivo (casos completos)

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

# =============================================================================
# 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

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-7 níveis).
• min_n: mínimo de observações por nó terminal (15-40).

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.

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, 7L)), # á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:

# Top 10 combinacoes por AUC
cart_tuned |>
    collect_metrics() |>
    filter(.metric == "roc_auc") |>
    arrange(desc(mean)) |>
    slice_head(n = 10)

# 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:

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:

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)

# 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:

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

# 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)

# 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):

print(round(stability_summary * 100, 1))

            correcao_colete                       lenke 
                        100                          56 
            cifose_toracica                      risser 
                         52                          45 
             lordose_lombar                         imc 
                         36                          23 
              flexibilidade                       idade 
                         19                           3 
escoliometro_maior_10_graus                        sexo 
                          2                           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:

# 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 1	nao_melhora	77	20.8%	correcao_colete< 54.86 E risser=0,3
Grupo 2	nao_melhora	38	23.7%	correcao_colete< 54.86 E risser=1,2,4 E flexibilidade=rigido
Grupo 3	nao_melhora	23	30.4%	correcao_colete< 54.86 E risser=1,2,4 E flexibilidade=flexivel E lordose_lombar>=60
Grupo 5	melhora	76	76.3%	correcao_colete>=54.86
Grupo 4	melhora	38	71.1%	correcao_colete< 54.86 E risser=1,2,4 E flexibilidade=flexivel E lordose_lombar< 60

2.2.2.13 Curva ROC

# 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°\).

# =============================================================================
# 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.

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.78	0.33, 1.70	0.5
imc	1.01		0.8
cifose_toracica	0.91	0.79, 1.02	0.15
lordose_lombar	1.02	0.92, 1.14	0.7
correcao_colete	0.99	0.92, 1.05	0.8
sexo
feminino	—	—
masculino	3.20	0.15, 65.7	0.4
lenke
1	—	—
2	208,263,068	0.00,	>0.9
3	91,264,991	0.00,	>0.9
4	822,832,895	0.00,	>0.9
5	4.41	0.00, 213,643,943,039,839,572,242,701,537,181,993,185,291,627,547,844,333,032,119,234,035,110,563,300,607,674,679,843,943,868,613,007,244,978,325,609,716,324,317,714,844,126,537,089,258,414,911,448,957,405,625,451,453,990,495,814,410,271,352,860,096,973,030,788,085,493,871,822,456,540,364,800	>0.9
6	25,607,199	0.00,	>0.9
risser
0	—	—
1	0.00		>0.9
2	0.13	0.01, 1.16	0.093
3	0.00		>0.9
4	0.46	0.01, 8.41	0.6
flexibilidade
flexivel	—	—
rigido	12.0	1.23, 309	0.068
escoliometro_maior_10_graus
normal	—	—
lombar	0.00		>0.9
toracica	0.42	0.02, 12.3	0.6
toracica_lombar	9.62	0.30, 691	0.2
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

2.3.1.2 Resíduos binarizados

performance::binned_residuals(model_progressao) |> plot()

2.3.1.3 Autocorrelação dos resíduos

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

performance::check_autocorrelation(model_progressao)

OK: Residuals appear to be independent and not autocorrelated (p = 0.356).

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.

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.

# 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")
    }
)

Aviso: Teste de Hosmer-Lemeshow não pode ser executado.
Possível causa: poucos eventos de progressão ou falta de variação nas predições.
Erro: 'breaks' are not unique 

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.

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
escoliometro_maior_10_graus	3.34	2.78	4.08
lenke	3.25	2.71	3.96
risser	2.48	2.10	3.00
cifose_toracica	1.99	1.70	2.38
idade	1.96	1.68	2.35
flexibilidade	1.92	1.65	2.30
lordose_lombar	1.86	1.60	2.23
sexo	1.79	1.54	2.14
imc	1.17	1.07	1.41
correcao_colete	1.14	1.05	1.39

performance::check_collinearity(model_progressao) |> plot()

2.3.1.7 Outliers

performance::check_outliers(model_progressao)

1 outlier detected: case 186.
- Based on the following method and threshold: cook (0.9).
- For variable: (Whole model).

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.

performance::r2(model_progressao)

# R2 for Logistic Regression
  Tjur's R2: 0.332

2.3.1.9 Tamanho amostral efetivo (casos completos)

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

# =============================================================================
# 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

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-7 níveis).
• min_n: mínimo de observações por nó terminal (15-40).

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.

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, 7L)), # á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:

# Top 10 combinacoes por AUC
cart_tuned_progressao |>
    collect_metrics() |>
    filter(.metric == "roc_auc") |>
    arrange(desc(mean)) |>
    slice_head(n = 10)

# 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:

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.

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)

# 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:

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

# 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")
    }
)

Aviso: Árvore de progressão não possui splits suficientes para análise de importância.
Possível causa: poucos eventos de progressão ou parâmetros muito restritivos.

2.3.2.11 Análise de estabilidade (Bootstrap)

# 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):

print(round(stability_summary_progressao * 100, 1))

                      lenke                      risser 
                          9                           5 
            cifose_toracica escoliometro_maior_10_graus 
                          3                           3 
              flexibilidade              lordose_lombar 
                          3                           3 
                      idade             correcao_colete 
                          2                           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:

# 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 1	nao_progrediu	252	3.6%	(raiz)

2.3.2.13 Curva ROC

# 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"))