Scoliosis Model

Autor

Caio Vallio

Data de Publicação

18 de dezembro de 2025

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 melhora radiográfica em 6 meses. Métodos: Estudo observacional com análise exploratória e inferencial. O desfecho contínuo foi \(\Delta\) (maior curva em 6 meses − baseline, em graus). O desfecho binário (MCID) foi melhora definida como \(\Delta \le -5°\). Ajustamos um modelo linear para \(\Delta\) e um modelo logístico para o MCID, reportando estimativas com IC95% e verificações de adequação dos modelos.

0.1 Como reproduzir

  • O arquivo de dados é lido de data/modelagem final.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 Definições de desfecho

  • Desfecho contínuo: \(\Delta\) = (maior curva em 6 meses − maior curva baseline), em graus. Valores mais negativos indicam maior melhora.
  • Desfecho binário (MCID): melhora clínica definida como \(\Delta \le -5°\) (redução de pelo menos 5 graus).

0.4 Plano de análise estatística

  • Descrição da amostra: estatísticas descritivas das variáveis basais e do desfecho.
  • Inferência (associações ajustadas):
    • Modelo linear para \(\Delta\) (efeitos como betas, IC95%).
    • Modelo logístico para MCID (efeitos como odds ratios, IC95%).
  • Adequação dos modelos: diagnósticos de assunções (colinearidade, resíduos, heteroscedasticidade quando aplicável) e medidas de qualidade de ajuste (por exemplo, \(R^2\)).
  • Forma funcional: para preditores contínuos, assume-se relação aproximadamente linear.

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 (reprodutibilidade)
pkgs <- c(
    "tidyverse", "gtsummary",
    "readxl", "janitor",
    "performance", "qqplotr", "PupillometryR"
)
missing_pkgs <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_pkgs) > 0) {
    stop("Pacotes ausentes: ", paste(missing_pkgs, collapse = ", "))
}
invisible(lapply(pkgs, library, character.only = TRUE))

set.seed(123)
theme_set(theme_minimal())


# read data
df_raw <- readxl::read_excel("data/modelagem final.xlsx", sheet = "dados")

# clean names
df <- df_raw |> janitor::clean_names()

# clean escoliometro
df <- df |>
    mutate(
        escoliometro = pmax(
            escoliomertro_cervical,
            escoliometro_torarica,
            escoliometro_lombar,
            na.rm = TRUE
        )
    ) |>
    mutate(
        regiao = case_when(
            escoliometro == escoliomertro_cervical ~ "cervical",
            escoliometro == escoliometro_torarica ~ "toracica",
            escoliometro == escoliometro_lombar ~ "lombar",
            TRUE ~ NA_character_
        )
    ) |>
    select(
        -escoliomertro_cervical,
        -escoliometro_torarica,
        -escoliometro_lombar
    )

# definicao do desfecho (ver secao acima)
df <- df |>
    mutate(delta = maior_curva_6_meses - maior_curva_baseline) |>
    mutate(delta_cat = if_else(delta <= -5, 1, 0)) |>
    mutate(
        delta_cat_f = factor(
            delta_cat,
            levels = c(0, 1),
            labels = c("Sem melhora (MCID)", "Melhora (MCID)")
        )
    )

# casting de variaveis (garantir tipos adequados no ajuste)
df <- df |>
    mutate(
        idade = as.double(idade),
        lenke = as.factor(lenke),
        risser = as.factor(risser),
        sexo = as.factor(sexo),
        regiao = as.factor(regiao),
        cifose_categ = as.factor(cifose_categ)
    )

0.5 Dados faltantes e tamanho amostral

Código 
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))

missing_tbl |>
    mutate(pct_missing = scales::percent(pct_missing, accuracy = 0.1)) |>
    print(n = 50)
# A tibble: 21 × 3
   variavel             n_missing pct_missing
   <chr>                    <int> <chr>      
 1 id                           0 0.0%       
 2 idade                        0 0.0%       
 3 sexo                         0 0.0%       
 4 altura                       0 0.0%       
 5 peso                         0 0.0%       
 6 imc                          0 0.0%       
 7 regiao                       0 0.0%       
 8 maior_curva_baseline         0 0.0%       
 9 maior_curva_6_meses          0 0.0%       
10 delta                        0 0.0%       
11 lenke                        0 0.0%       
12 risser                       0 0.0%       
13 cifose_toracica              0 0.0%       
14 lordose_lombar               0 0.0%       
15 cifose_categ                 0 0.0%       
16 inclinacao                   0 0.0%       
17 dif_colete                   0 0.0%       
18 correcao_colete              0 0.0%       
19 escoliometro                 0 0.0%       
20 delta_cat                    0 0.0%       
21 delta_cat_f                  0 0.0%

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.6 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 
categoricas <- c("sexo", "regiao", "lenke", "risser", "cifose_categ")
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 regiao       lombar | toracica                
3 lenke        1 | 2 | 3 | 4 | 5 | 6            
4 risser       0 | 1 | 2 | 3 | 4                
5 cifose_categ hipercifose | hipocifose | normal




1 Análise descritiva

1.1 Variáveis de entrada do modelo

Variáveis coletadas na linha de base.

Código 
df |>
    gtsummary::tbl_summary(
        include = c(
            idade,
            altura,
            peso,
            imc,
            cifose_toracica,
            lordose_lombar,
            dif_colete,
            correcao_colete,
            escoliometro,
            inclinacao,
            cifose_categ,
            sexo,
            regiao,
            lenke,
            risser
        ),
        type = list(
            idade ~ "continuous"
        ),
        statistic = list(
            gtsummary::all_continuous() ~ "{mean} ({sd})",
            gtsummary::all_categorical() ~ "{n} ({p}%)"
        )
    )
Characteristic N = 1521
idade 13.07 (1.48)
altura 1.58 (0.08)
peso 48 (9)
imc 19.05 (2.91)
cifose_toracica 22 (12)
lordose_lombar 54 (12)
dif_colete -17.6 (6.2)
correcao_colete 49 (19)
escoliometro 13.8 (4.3)
inclinacao
    flexivel 96 (63%)
    rigido 56 (37%)
cifose_categ
    hipercifose 8 (5.3%)
    hipocifose 17 (11%)
    normal 127 (84%)
sexo
    feminino 130 (86%)
    masculino 22 (14%)
regiao
    lombar 80 (53%)
    toracica 72 (47%)
lenke
    1 8 (5.3%)
    2 11 (7.2%)
    3 32 (21%)
    4 13 (8.6%)
    5 16 (11%)
    6 72 (47%)
risser
    0 32 (21%)
    1 34 (22%)
    2 45 (30%)
    3 24 (16%)
    4 17 (11%)
1 Mean (SD); n (%)

1.2 Desfecho principal

  • Maior curva na linha de base e em 6 meses.
  • Delta: maior curva 6 meses - maior curva baseline.
  • Delta categórica: delta melhor em 5 graus ou mais (MCID).
Código 
df |>
    gtsummary::tbl_summary(
        include = c(
            maior_curva_baseline,
            maior_curva_6_meses,
            delta,
            delta_cat_f
        ),
        type = list(
            delta_cat_f ~ "categorical"
        ),
        statistic = list(
            gtsummary::all_continuous() ~ "{mean} ({sd})",
            gtsummary::all_categorical() ~ "{n} ({p}%)"
        )
    )
Characteristic N = 1521
maior_curva_baseline 36.6 (5.9)
maior_curva_6_meses 31 (9)
delta -5.4 (6.1)
delta_cat_f
    Sem melhora (MCID) 68 (45%)
    Melhora (MCID) 84 (55%)
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),
        # fill = "steelblue",
        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 de regressão linear

Este modelo estima associações ajustadas com o desfecho contínuo \(\Delta\) (maior curva em 6 meses − baseline, em graus).

Código 
model_lin <- lm(
    delta ~
        idade +
        altura +
        peso +
        # imc +
        cifose_toracica +
        lordose_lombar +
        dif_colete +
        correcao_colete +
        escoliometro +
        inclinacao +
        cifose_categ +
        sexo +
        regiao +
        lenke +
        risser,
    data = df
)

2.1.1 Parâmetros do modelo

Código 
model_lin |>
    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.54 -1.3, 0.20 0.2
altura -3.0 -19, 13 0.7
peso 0.02 -0.10, 0.13 0.8
cifose_toracica -0.03 -0.14, 0.07 0.5
lordose_lombar 0.01 -0.08, 0.09 0.9
dif_colete 0.34 0.03, 0.65 0.032
correcao_colete 0.01 -0.09, 0.11 0.8
escoliometro 0.00 -0.21, 0.21 >0.9
inclinacao


    flexivel — —
    rigido 2.6 0.63, 4.5 0.010
cifose_categ


    hipercifose — —
    hipocifose -3.7 -9.5, 2.1 0.2
    normal -4.9 -9.3, -0.49 0.030
sexo


    feminino — —
    masculino 0.81 -2.0, 3.6 0.6
regiao


    lombar — —
    toracica 0.24 -2.3, 2.8 0.9
lenke


    1 — —
    2 1.0 -3.8, 5.9 0.7
    3 6.3 2.2, 10 0.003
    4 5.0 0.46, 9.6 0.031
    5 -1.6 -6.4, 3.2 0.5
    6 2.9 -1.1, 6.9 0.2
risser


    0 — —
    1 -2.0 -4.6, 0.57 0.12
    2 -1.3 -3.8, 1.1 0.3
    3 0.23 -2.7, 3.1 0.9
    4 0.58 -2.9, 4.1 0.7
Abbreviation: CI = Confidence Interval

2.1.2 Verificações de adequação do modelo

2.1.2.1 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_lin) |>
    arrange(desc(VIF)) |>
    select(Term, VIF, VIF_CI_low, VIF_CI_high) |>
    performance::print_html()
Term VIF VIF_CI_low VIF_CI_high
dif_colete 5.74 4.53 7.35
correcao_colete 5.58 4.41 7.15
lenke 4.49 3.57 5.72
regiao 2.57 2.11 3.23
cifose_toracica 2.41 1.99 3.02
altura 2.34 1.93 2.93
risser 2.30 1.90 2.88
cifose_categ 2.18 1.81 2.73
idade 1.88 1.58 2.34
peso 1.68 1.43 2.08
lordose_lombar 1.65 1.40 2.03
sexo 1.52 1.31 1.87
inclinacao 1.40 1.22 1.72
escoliometro 1.26 1.12 1.56
Código 
performance::check_collinearity(model_lin) |> plot()

2.1.2.2 Heteroscedasticidade

Verificação de heteroscedasticidade entre as variáveis independentes.

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

2.1.2.3 Outliers

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

2.1.2.4 Valores ajustados (diagnóstico)

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

2.1.2.5 Resíduos simulados

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

2.1.3 Qualidade do ajuste do modelo

2.1.3.1 Coeficiente de determinação

\(R^2\).

Código 
# Qualidade do ajuste
performance::r2(model_lin)
# R2 for Linear Regression
       R2: 0.436
  adj. R2: 0.340

2.1.3.2 Tamanho amostral efetivo (casos completos)

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




2.2 Modelo de regressão logística

Este modelo estima associações ajustadas com a chance de melhora (MCID), definida como \(\Delta \le -5°\).

Código 
model_bin <- glm(
    delta_cat ~
        idade +
        altura +
        peso +
        # imc +
        cifose_toracica +
        lordose_lombar +
        dif_colete +
        correcao_colete +
        escoliometro +
        inclinacao +
        cifose_categ +
        sexo +
        regiao +
        lenke +
        risser,
    data = df, family = "binomial"
)

2.2.1 Parâmetros do modelo

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

Código 
model_bin |>
    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 1.19 0.79, 1.85 0.4
altura 0.65 0.00, 3,725 >0.9
peso 1.01 0.94, 1.07 0.8
cifose_toracica 1.01 0.95, 1.07 0.7
lordose_lombar 0.98 0.93, 1.02 0.3
dif_colete 0.85 0.72, 1.01 0.064
correcao_colete 1.02 0.96, 1.08 0.5
escoliometro 1.12 1.00, 1.28 0.074
inclinacao


    flexivel — —
    rigido 0.23 0.07, 0.64 0.007
cifose_categ


    hipercifose — —
    hipocifose 4.92 0.14, 185 0.4
    normal 16.7 1.33, 326 0.042
sexo


    feminino — —
    masculino 1.46 0.29, 7.93 0.6
regiao


    lombar — —
    toracica 1.32 0.33, 5.54 0.7
lenke


    1 — —
    2 0.42 0.02, 7.63 0.5
    3 0.05 0.00, 0.41 0.009
    4 0.17 0.01, 1.57 0.13
    5 2.21 0.15, 30.2 0.6
    6 0.36 0.04, 2.60 0.3
risser


    0 — —
    1 3.94 1.00, 16.7 0.054
    2 1.93 0.51, 7.54 0.3
    3 0.69 0.13, 3.55 0.7
    4 1.60 0.26, 10.5 0.6
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

2.2.2 Verificações de adequação do modelo

2.2.2.1 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_bin) |>
    arrange(desc(VIF)) |>
    select(Term, VIF, VIF_CI_low, VIF_CI_high) |>
    performance::print_html()
Term VIF VIF_CI_low VIF_CI_high
lenke 5.66 4.48 7.26
dif_colete 3.85 3.09 4.89
correcao_colete 3.84 3.08 4.87
risser 2.93 2.38 3.70
altura 2.71 2.21 3.41
regiao 2.67 2.18 3.36
cifose_categ 2.30 1.90 2.88
cifose_toracica 2.28 1.88 2.85
idade 1.97 1.65 2.45
sexo 1.77 1.50 2.20
peso 1.75 1.48 2.17
lordose_lombar 1.57 1.34 1.93
inclinacao 1.46 1.27 1.80
escoliometro 1.43 1.24 1.76
Código 
performance::check_collinearity(model_bin) |> plot()

2.2.2.2 Resíduos binarizados

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

2.2.2.3 Outliers

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

2.2.2.4 Valores ajustados (diagnóstico)

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

2.2.2.5 Resíduos simulados

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

2.2.3 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 
performance::performance_hosmer(model_bin)
# Hosmer-Lemeshow Goodness-of-Fit Test

  Chi-squared: 7.009
           df: 8    
      p-value: 0.536

2.2.3.1 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_bin)
# R2 for Logistic Regression
  Tjur's R2: 0.448

2.2.3.2 Tamanho amostral efetivo (casos completos)

Código 
mf_bin <- model.frame(model_bin)
bin_n <- nrow(mf_bin)
bin_events <- sum(mf_bin[[1]] == 1, na.rm = TRUE)
bin_nonevents <- sum(mf_bin[[1]] == 0, na.rm = TRUE)
tibble(
    n_modelo = bin_n,
    eventos_mcid = bin_events,
    nao_eventos = bin_nonevents
)




2.3 CART (árvore de decisão) para MCID (delta_cat)

2.3.1 Racional teórico

Um modelo CART (Classification and Regression Tree) é útil aqui porque:

  • Não linearidade: identifica automaticamente limiares (pontos de corte) em preditores contínuos.
  • Interações: a sequência de divisões (splits) representa interações condicionais sem precisar pre-especificar termos de interação.
  • Interpretabilidade: gera regras do tipo “se… então…”, úteis para comunicação clínica e geração de hipóteses.

Limitações importantes: árvores isoladas tendem a ter alta variância e podem superajustar. Por isso, esta seção usa validação cruzada (CV) e tuning/poda para controlar complexidade e reporta desempenho interno (sem validação externa).

2.3.2 Setup e preparação dos dados

Código 
library(tidymodels)
library(rpart.plot)
library(vip)

set.seed(123)

# Desfecho como fator: "melhora" (evento) vs "nao_melhora"
df_cart <- df |>
    mutate(
        delta_cat_cart = factor(
            delta_cat,
            levels = c(0, 1),
            labels = c("nao_melhora", "melhora")
        )
    ) |>
    select(
        delta_cat_cart,
        idade, altura, peso,
        cifose_toracica, lordose_lombar,
        dif_colete, correcao_colete,
        escoliometro, inclinacao,
        cifose_categ, sexo, regiao, lenke, risser
    )

# Verificar niveis (evento = "melhora", segundo nivel)
levels(df_cart$delta_cat_cart)
[1] "nao_melhora" "melhora"

2.3.3 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.3.4 Modelo e hiperparâmetros

Hiperparâmetros a serem tunados:

  • cost_complexity (cp): controla poda/complexidade; valores maiores = árvores menores.
  • tree_depth: limita profundidade (proxy do grau de interações).
  • min_n: mínimo de observações em no terminal; evita regras baseadas em poucos pacientes.
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.3.5 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, -1)),
    tree_depth(range = c(1L, 6L)),
    min_n(range = c(5L, 30L)),
    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.3.6 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.3.7 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.3.8 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.3.9 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.3.10 Á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.3.11 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.3.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 1 nao_melhora 43 18.6% correcao_colete< 51.39 E inclinacao=rigido
Grupo 4 nao_melhora 18 44.4% correcao_colete< 51.39 E inclinacao=flexivel E cifose_toracica>=18.5 E lordose_lombar>=55.5
Grupo 2 nao_melhora 10 10.0% correcao_colete< 51.39 E inclinacao=flexivel E cifose_toracica< 18.5 E escoliometro< 13.5
Grupo 6 nao_melhora 7 14.3% correcao_colete>=51.39 E lenke=3,4,6 E lordose_lombar>=67
Grupo 7 melhora 36 83.3% correcao_colete>=51.39 E lenke=3,4,6 E lordose_lombar< 67
Grupo 8 melhora 20 100.0% correcao_colete>=51.39 E lenke=1,2,5
Grupo 5 melhora 13 100.0% correcao_colete< 51.39 E inclinacao=flexivel E cifose_toracica>=18.5 E lordose_lombar< 55.5
Grupo 3 melhora 5 60.0% correcao_colete< 51.39 E inclinacao=flexivel E cifose_toracica< 18.5 E escoliometro>=13.5

2.3.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.14 Interpretação

Como ler a árvore:

  • Cada nó interno representa uma condição (split) sobre uma variável.
  • Cada caminho da raiz ate um nó terminal representa uma regra de classificação.
  • A sequencia de splits indica interações condicionais: um split so importa dentro da região definida por splits anteriores.
  • Os limiares (pontos de corte) sugerem não-linearidades nos efeitos.

Como ler a importância:

  • Variáveis com maior importância contribuem mais para a redução de impureza (Gini) nos splits.
  • Isso não implica causalidade, apenas relevância preditiva no contexto do modelo.

Limitações:

  • Desempenho reportado e interno (CV), sem validação externa.
  • Árvores são instáveis: pequenas mudanças nos dados podem alterar a estrutura.
  • Usar como ferramenta exploratória complementar aos modelos inferenciais (linear/logístico).