Scoliosis Model

Autor

Caio Vallio

Data de Publicação

24 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/dataset_escoliose_00.xlsx",
    sheet = "dados",
    na = ""
)

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

# escoliometro - classificacao do local onde a escoliose esta mais severa
df <- df |>
    mutate(
        escoliometro_maior_10_graus = case_when(
            (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",
            (escoliometro_cervical > 10) ~ "cervical",
            (escoliometro_torarica > 10) ~ "toracica",
            (escoliometro_lombar > 10) ~ "lombar",
            TRUE ~ "nenhum"
        )
    )

# lenke - classificacao da curvatura da escoliose
df <- df |>
    mutate(
        lenke_curvatura = case_when(
            lenke %in% c(1, 2, 5) ~ "curva em c",
            lenke %in% c(3, 4, 6) ~ "curva em s",
            TRUE ~ NA_character_
        )
    )

# risser - classificacao do inicio e final do crescimento
df <- df |>
    mutate(
        risser_agrupado = case_when(
            risser %in% c(0, 1, 2) ~ "inicio do crescimento",
            risser %in% c(3, 4) ~ "final do crescimento",
            TRUE ~ NA_character_
        )
    )

# cifose - classificacao da cifose
df <- df |>
    mutate(
        cifose_categ = case_when(
            cifose_toracica < 10 ~ "hipocifose",
            cifose_toracica >= 10 & cifose_toracica < 40 ~ "normalidade",
            cifose_toracica >= 40 ~ "hipercifose",
            TRUE ~ NA_character_
        )
    )

# lordose - classificacao da lordose
df <- df |>
    mutate(
        lordose_categ = case_when(
            lordose_lombar < 45 ~ "hipolordose",
            lordose_lombar >= 45 & lordose_lombar < 50 ~ "normalidade",
            lordose_lombar >= 50 ~ "hiperlordose",
            TRUE ~ NA_character_
        )
    )

# imc - calculo do imc
df <- df |>
    mutate(
        imc = peso / (altura)^2
    )

# definicao do desfecho - alvo
# cobb inicial - alvo
df <- df |>
    mutate(
        cobb_inicial_maior = pmax(
            cobb_toracico_proximal,
            cobb_toracica,
            cobb_lombar,
            na.rm = TRUE
        )
    ) |>
    mutate(
        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_
        )
    )

# categorizacao do desfecho - alvo
df <- df |>
    mutate(delta = maior_curva_6_meses - cobb_inicial_maior) |>
    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),
        flexibilidade = as.factor(flexibilidade),
        regiao_cobb_inicial = as.factor(regiao_cobb_inicial),
        escoliometro_maior_10_graus = as.factor(escoliometro_maior_10_graus),
        lenke_curvatura = as.factor(lenke_curvatura),
        risser_agrupado = as.factor(risser_agrupado),
        cifose_categ = as.factor(cifose_categ),
        lordose_categ = as.factor(lordose_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: 30 × 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 lenke_curvatura                     0 0.0%       
22 risser_agrupado                     0 0.0%       
23 cifose_categ                        0 0.0%       
24 lordose_categ                       0 0.0%       
25 imc                                 0 0.0%       
26 cobb_inicial_maior                  0 0.0%       
27 regiao_cobb_inicial                 0 0.0%       
28 delta                               0 0.0%       
29 delta_cat                           0 0.0%       
30 delta_cat_f                         0 0.0%       
Código
# drop na in flexibilidade
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.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",
    "lenke",
    "risser",
    "flexibilidade",
    "regiao_cobb_inicial",
    "escoliometro_maior_10_graus",
    "lenke_curvatura",
    "risser_agrupado",
    "cifose_categ",
    "lordose_categ"
)
map_dfr(categoricas, function(v) {
    tibble(
        variavel = v,
        niveis = paste(levels(df[[v]]), collapse = " | ")
    )
}) |>
    print(n = Inf)
# A tibble: 10 × 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 regiao_cobb_inicial         lombar | toracica | toracico_proximal       
 6 escoliometro_maior_10_graus lombar | nenhum | toracica | toracica_lombar
 7 lenke_curvatura             curva em c | curva em s                     
 8 risser_agrupado             final do crescimento | inicio do crescimento
 9 cifose_categ                hipercifose | hipocifose | normalidade      
10 lordose_categ               hiperlordose | hipolordose | normalidade    




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(
            # numericas
            idade,
            altura,
            peso,
            imc,
            cifose_toracica,
            lordose_lombar,
            dif_colete,
            correcao_colete,

            # categoricas
            sexo,
            lenke,
            risser,
            flexibilidade,
            regiao_cobb_inicial,
            escoliometro_maior_10_graus,
            lenke_curvatura,
            risser_agrupado,
            cifose_categ,
            lordose_categ,
        ),
        type = list(
            idade ~ "continuous"
        ),
        statistic = list(
            gtsummary::all_continuous() ~ "{mean} ({sd})",
            gtsummary::all_categorical() ~ "{n} ({p}%)"
        )
    )
Characteristic N = 2521
idade 12.96 (1.58)
altura 1.58 (0.08)
peso 49 (27)
imc 19.76 (10.87)
cifose_toracica 22 (11)
lordose_lombar 54 (12)
dif_colete -17.4 (6.4)
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%)
regiao_cobb_inicial
    lombar 137 (54%)
    toracica 112 (44%)
    toracico_proximal 3 (1.2%)
escoliometro_maior_10_graus
    lombar 88 (35%)
    nenhum 52 (21%)
    toracica 94 (37%)
    toracica_lombar 18 (7.1%)
lenke_curvatura
    curva em c 56 (22%)
    curva em s 196 (78%)
risser_agrupado
    final do crescimento 90 (36%)
    inicio do crescimento 162 (64%)
cifose_categ
    hipercifose 20 (7.9%)
    hipocifose 25 (9.9%)
    normalidade 207 (82%)
lordose_categ
    hiperlordose 169 (67%)
    hipolordose 50 (20%)
    normalidade 33 (13%)
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(
            cobb_inicial_maior,
            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 = 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%)
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 Combinação de variáveis

2.1.1 Variáveis originais

Incluindo:

  • idade
  • imc
  • cifose_toracica
  • lordose_lombar
  • dif_colete
  • sexo
  • lenke
  • risser
  • flexibilidade
  • escoliometro_maior_10_graus

2.1.1.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 ~
        # numericas
        idade +
        # altura +
        # peso +
        imc +
        cifose_toracica +
        lordose_lombar +
        dif_colete +
        # correcao_colete +

        # categoricas
        sexo +
        lenke +
        risser +
        flexibilidade +
        # regiao_cobb_inicial +
        escoliometro_maior_10_graus
    # lenke_curvatura +
    # risser_agrupado +
    # cifose_categ +
    # lordose_categ,
    ,
    data = df
)
2.1.1.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.15 -0.65, 0.34 0.5
imc -0.03 -0.09, 0.02 0.2
cifose_toracica -0.03 -0.10, 0.03 0.3
lordose_lombar 0.01 -0.05, 0.06 0.9
dif_colete 0.29 0.19, 0.39 <0.001
sexo


    feminino
    masculino -0.55 -2.6, 1.5 0.6
lenke


    1
    2 2.0 -1.3, 5.3 0.2
    3 5.2 2.4, 8.0 <0.001
    4 5.6 2.5, 8.7 <0.001
    5 1.3 -2.0, 4.6 0.4
    6 3.5 0.71, 6.4 0.014
risser


    0
    1 -2.5 -4.4, -0.48 0.015
    2 -0.99 -2.8, 0.84 0.3
    3 -0.22 -2.2, 1.7 0.8
    4 0.02 -2.4, 2.4 >0.9
flexibilidade


    flexivel
    rigido 1.7 0.44, 3.0 0.008
escoliometro_maior_10_graus


    lombar
    nenhum 0.20 -1.6, 1.9 0.8
    toracica 0.80 -1.2, 2.7 0.4
    toracica_lombar 2.6 0.07, 5.2 0.044
Abbreviation: CI = Confidence Interval
2.1.1.1.2 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
lenke 2.59 2.18 3.14
escoliometro_maior_10_graus 2.50 2.11 3.02
risser 2.04 1.75 2.45
idade 1.83 1.58 2.19
cifose_toracica 1.48 1.31 1.76
lordose_lombar 1.34 1.20 1.59
dif_colete 1.23 1.11 1.47
sexo 1.18 1.08 1.42
flexibilidade 1.13 1.05 1.39
imc 1.07 1.01 1.42
Código
performance::check_collinearity(model_lin) |> plot()

2.1.1.1.3 Heteroscedasticidade

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

Código
performance::check_heteroscedasticity(model_lin)
Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.033).
Código
performance::check_heteroscedasticity(model_lin) |> plot()

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

2.1.1.1.5 Valores ajustados (diagnóstico)
Código
performance::check_predictions(model_lin) |> plot()

2.1.1.1.6 Resíduos simulados
Código
performance::check_residuals(model_lin)
OK: Simulated residuals appear as uniformly distributed (p = 0.803).
Código
performance::check_residuals(model_lin) |> plot()

2.1.1.1.7 Coeficiente de determinação

\(R^2\).

Código
# Qualidade do ajuste
performance::r2(model_lin)
# R2 for Linear Regression
       R2: 0.343
  adj. R2: 0.289
2.1.1.1.8 Tamanho amostral efetivo (casos completos)
Código
tibble(n_modelo = nobs(model_lin))




2.1.1.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 ~
        # numericas
        idade +
        # altura +
        # peso +
        imc +
        cifose_toracica +
        lordose_lombar +
        dif_colete +
        # correcao_colete +

        # categoricas
        sexo +
        lenke +
        risser +
        flexibilidade +
        # regiao_cobb_inicial +
        escoliometro_maior_10_graus
    # lenke_curvatura +
    # risser_agrupado +
    # cifose_categ +
    # lordose_categ,
    ,
    data = df,
    family = "binomial"
)
2.1.1.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 0.96 0.74, 1.24 0.7
imc 1.03 1.00, 1.12 0.3
cifose_toracica 1.03 1.0, 1.06 0.11
lordose_lombar 0.99 0.96, 1.02 0.5
dif_colete 0.86 0.81, 0.91 <0.001
sexo


    feminino
    masculino 2.23 0.78, 6.55 0.14
lenke


    1
    2 0.18 0.03, 1.02 0.061
    3 0.08 0.01, 0.34 0.001
    4 0.07 0.01, 0.38 0.003
    5 0.33 0.05, 1.76 0.2
    6 0.19 0.03, 0.79 0.033
risser


    0
    1 2.92 1.07, 8.22 0.039
    2 1.41 0.55, 3.69 0.5
    3 0.66 0.24, 1.83 0.4
    4 1.81 0.53, 6.30 0.3
flexibilidade


    flexivel
    rigido 0.44 0.22, 0.84 0.014
escoliometro_maior_10_graus


    lombar
    nenhum 1.27 0.54, 3.01 0.6
    toracica 1.34 0.50, 3.67 0.6
    toracica_lombar 1.07 0.29, 3.99 >0.9
Abbreviations: CI = Confidence Interval, OR = Odds Ratio
2.1.1.2.2 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 2.52 2.13 3.05
escoliometro_maior_10_graus 2.46 2.08 2.97
risser 2.11 1.80 2.54
idade 1.92 1.65 2.30
cifose_toracica 1.49 1.31 1.78
lordose_lombar 1.29 1.16 1.54
sexo 1.24 1.12 1.49
dif_colete 1.21 1.09 1.44
flexibilidade 1.13 1.05 1.39
imc 1.06 1.01 1.53
Código
performance::check_collinearity(model_bin) |> plot()

2.1.1.2.3 Resíduos binarizados
Código
performance::binned_residuals(model_bin) |> plot()

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

2.1.1.2.5 Valores ajustados (diagnóstico)
Código
performance::check_predictions(model_bin) |> plot()

2.1.1.2.6 Resíduos simulados
Código
performance::check_residuals(model_bin)
OK: Simulated residuals appear as uniformly distributed (p = 0.173).
Código
performance::check_residuals(model_bin) |> plot()

2.1.1.2.7 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: 11.480
           df:  8    
      p-value:  0.176
2.1.1.2.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_bin)
# R2 for Logistic Regression
  Tjur's R2: 0.283
2.1.1.2.9 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.1.1.3 CART (árvore de decisão) para MCID (delta_cat)

2.1.1.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.1.1.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,
        # numericas
        idade,
        # altura,
        # peso,
        imc,
        cifose_toracica,
        lordose_lombar,
        dif_colete,
        # correcao_colete,

        # categoricas
        sexo,
        lenke,
        risser,
        flexibilidade,
        # regiao_cobb_inicial,
        escoliometro_maior_10_graus,
        # lenke_curvatura,
        # risser_agrupado,
        # cifose_categ,
        # lordose_categ,
    )

# Verificar niveis (evento = "melhora", segundo nivel)
levels(df_cart$delta_cat_cart)
[1] "nao_melhora" "melhora"    
2.1.1.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.1.1.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.1.1.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(-5, -2)),
    tree_depth(range = c(4L, 10L)),
    min_n(range = c(5L, 20L)),
    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.1.1.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.1.1.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.1.1.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.1.1.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.1.1.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.1.1.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.1.1.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 96 19.8% dif_colete>=-19.5 E lenke=2,3,4,6 E cifose_toracica< 24.5
Grupo 2 nao_melhora 23 17.4% dif_colete>=-19.5 E lenke=2,3,4,6 E cifose_toracica>=24.5 E lordose_lombar>=54.5 E escoliometro_maior_10_graus=lombar,toracica
Grupo 4 nao_melhora 8 37.5% dif_colete>=-19.5 E lenke=2,3,4,6 E cifose_toracica>=24.5 E lordose_lombar< 54.5 E escoliometro_maior_10_graus=nenhum,toracica
Grupo 6 nao_melhora 8 12.5% dif_colete>=-19.5 E lenke=1,5 E dif_colete>=-15.5 E flexibilidade=rigido
Grupo 9 melhora 75 77.3% dif_colete< -19.5
Grupo 7 melhora 12 66.7% dif_colete>=-19.5 E lenke=1,5 E dif_colete>=-15.5 E flexibilidade=flexivel
Grupo 3 melhora 11 54.5% dif_colete>=-19.5 E lenke=2,3,4,6 E cifose_toracica>=24.5 E lordose_lombar>=54.5 E escoliometro_maior_10_graus=nenhum
Grupo 8 melhora 10 100.0% dif_colete>=-19.5 E lenke=1,5 E dif_colete< -15.5
Grupo 5 melhora 9 88.9% dif_colete>=-19.5 E lenke=2,3,4,6 E cifose_toracica>=24.5 E lordose_lombar< 54.5 E escoliometro_maior_10_graus=lombar,toracica_lombar
2.1.1.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.1.1.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).






2.1.2 Variáveis discretizadas

Incluindo:

  • idade
  • imc
  • cifose_toracica_categ
  • lordose_lombar_categ
  • dif_colete
  • sexo
  • lenke_curvatura
  • risser_agrupado
  • flexibilidade
  • escoliometro_maior_10_graus

2.1.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 ~
        # numericas
        idade +
        # altura +
        # peso +
        imc +
        # cifose_toracica +
        # lordose_lombar +
        dif_colete +
        # correcao_colete +

        # categoricas
        sexo +
        # lenke +
        # risser +
        flexibilidade +
        # regiao_cobb_inicial +
        escoliometro_maior_10_graus +
        lenke_curvatura +
        risser_agrupado +
        cifose_categ +
        lordose_categ,
    data = df
)
2.1.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.21 -0.65, 0.22 0.3
imc -0.03 -0.09, 0.02 0.2
dif_colete 0.31 0.21, 0.41 <0.001
sexo


    feminino
    masculino -0.71 -2.7, 1.3 0.5
flexibilidade


    flexivel
    rigido 1.8 0.53, 3.1 0.006
escoliometro_maior_10_graus


    lombar
    nenhum 0.70 -0.99, 2.4 0.4
    toracica 2.3 0.84, 3.7 0.002
    toracica_lombar 3.3 0.88, 5.8 0.008
lenke_curvatura


    curva em c
    curva em s 3.4 2.0, 4.9 <0.001
risser_agrupado


    final do crescimento
    inicio do crescimento -1.1 -2.5, 0.21 0.10
cifose_categ


    hipercifose
    hipocifose -1.0 -4.0, 2.0 0.5
    normalidade -1.3 -3.6, 0.98 0.3
lordose_categ


    hiperlordose
    hipolordose 1.1 -0.46, 2.7 0.2
    normalidade 0.77 -1.1, 2.6 0.4
Abbreviation: CI = Confidence Interval
2.1.2.1.2 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
idade 1.38 1.23 1.65
escoliometro_maior_10_graus 1.29 1.16 1.54
lordose_categ 1.25 1.12 1.49
risser_agrupado 1.24 1.12 1.49
cifose_categ 1.24 1.12 1.49
sexo 1.18 1.08 1.43
dif_colete 1.14 1.05 1.40
flexibilidade 1.12 1.04 1.39
imc 1.09 1.02 1.41
lenke_curvatura 1.09 1.02 1.41
Código
performance::check_collinearity(model_lin) |> plot()

2.1.2.1.3 Heteroscedasticidade

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

Código
performance::check_heteroscedasticity(model_lin)
Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.011).
Código
performance::check_heteroscedasticity(model_lin) |> plot()

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

2.1.2.1.5 Valores ajustados (diagnóstico)
Código
performance::check_predictions(model_lin) |> plot()

2.1.2.1.6 Resíduos simulados
Código
performance::check_residuals(model_lin)
OK: Simulated residuals appear as uniformly distributed (p = 0.882).
Código
performance::check_residuals(model_lin) |> plot()

2.1.2.1.7 Coeficiente de determinação

\(R^2\).

Código
# Qualidade do ajuste
performance::r2(model_lin)
# R2 for Linear Regression
       R2: 0.316
  adj. R2: 0.276
2.1.2.1.8 Tamanho amostral efetivo (casos completos)
Código
tibble(n_modelo = nobs(model_lin))




2.1.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 ~
        # numericas
        idade +
        # altura +
        # peso +
        imc +
        # cifose_toracica +
        # lordose_lombar +
        dif_colete +
        # correcao_colete +

        # categoricas
        sexo +
        # lenke +
        # risser +
        flexibilidade +
        # regiao_cobb_inicial +
        escoliometro_maior_10_graus +
        lenke_curvatura +
        risser_agrupado +
        cifose_categ +
        lordose_categ,
    data = df,
    family = "binomial"
)
2.1.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.09 0.88, 1.35 0.4
imc 1.03 1.00, 1.13 0.3
dif_colete 0.87 0.82, 0.92 <0.001
sexo


    feminino
    masculino 2.00 0.73, 5.60 0.2
flexibilidade


    flexivel
    rigido 0.45 0.24, 0.83 0.012
escoliometro_maior_10_graus


    lombar
    nenhum 1.03 0.45, 2.33 >0.9
    toracica 0.71 0.35, 1.46 0.4
    toracica_lombar 0.75 0.22, 2.48 0.6
lenke_curvatura


    curva em c
    curva em s 0.33 0.16, 0.67 0.002
risser_agrupado


    final do crescimento
    inicio do crescimento 1.71 0.89, 3.34 0.11
cifose_categ


    hipercifose
    hipocifose 0.86 0.19, 3.73 0.8
    normalidade 1.47 0.50, 4.58 0.5
lordose_categ


    hiperlordose
    hipolordose 0.58 0.26, 1.26 0.2
    normalidade 0.72 0.29, 1.75 0.5
Abbreviations: CI = Confidence Interval, OR = Odds Ratio
2.1.2.2.2 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
idade 1.38 1.22 1.65
escoliometro_maior_10_graus 1.31 1.17 1.57
cifose_categ 1.25 1.12 1.49
risser_agrupado 1.24 1.12 1.49
sexo 1.21 1.10 1.46
lordose_categ 1.20 1.09 1.45
dif_colete 1.11 1.03 1.39
lenke_curvatura 1.11 1.03 1.39
flexibilidade 1.11 1.03 1.39
imc 1.06 1.01 1.55
Código
performance::check_collinearity(model_bin) |> plot()

2.1.2.2.3 Resíduos binarizados
Código
performance::binned_residuals(model_bin) |> plot()

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

2.1.2.2.5 Valores ajustados (diagnóstico)
Código
performance::check_predictions(model_bin) |> plot()

2.1.2.2.6 Resíduos simulados
Código
performance::check_residuals(model_bin)
OK: Simulated residuals appear as uniformly distributed (p = 0.408).
Código
performance::check_residuals(model_bin) |> plot()

2.1.2.2.7 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: 21.170
           df:  8    
      p-value:  0.007
2.1.2.2.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_bin)
# R2 for Logistic Regression
  Tjur's R2: 0.240
2.1.2.2.9 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.1.2.3 CART (árvore de decisão) para MCID (delta_cat)

2.1.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.1.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,
        # numericas
        idade,
        # altura,
        # peso,
        imc,
        # cifose_toracica,
        # lordose_lombar,
        dif_colete,
        # correcao_colete,

        # categoricas
        sexo,
        # lenke,
        # risser,
        flexibilidade,
        # regiao_cobb_inicial,
        escoliometro_maior_10_graus,
        lenke_curvatura,
        risser_agrupado,
        cifose_categ,
        lordose_categ,
    )

# Verificar niveis (evento = "melhora", segundo nivel)
levels(df_cart$delta_cat_cart)
[1] "nao_melhora" "melhora"    
2.1.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.1.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.1.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(-5, -2)),
    tree_depth(range = c(4L, 10L)),
    min_n(range = c(5L, 20L)),
    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.1.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.1.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.1.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.1.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.1.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.1.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.1.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 136 25.0% dif_colete>=-19.5 E lenke_curvatura=curva em s E imc>=14.19 E imc< 30.51
Grupo 4 nao_melhora 9 22.2% dif_colete>=-19.5 E lenke_curvatura=curva em c E dif_colete>=-16.5 E imc>=19.5
Grupo 9 nao_melhora 7 42.9% dif_colete< -19.5 E imc>=16.41 E cifose_categ=normalidade E lordose_categ=hipolordose
Grupo 7 nao_melhora 5 0.0% dif_colete< -19.5 E imc< 16.41 E idade< 12.5
Grupo 10 melhora 48 83.3% dif_colete< -19.5 E imc>=16.41 E cifose_categ=normalidade E lordose_categ=hiperlordose,normalidade
Grupo 5 melhora 19 57.9% dif_colete>=-19.5 E lenke_curvatura=curva em c E dif_colete>=-16.5 E imc< 19.5
Grupo 11 melhora 10 100.0% dif_colete< -19.5 E imc>=16.41 E cifose_categ=hipercifose,hipocifose
Grupo 6 melhora 9 88.9% dif_colete>=-19.5 E lenke_curvatura=curva em c E dif_colete< -16.5
Grupo 8 melhora 5 100.0% dif_colete< -19.5 E imc< 16.41 E idade>=12.5
Grupo 2 melhora 2 100.0% dif_colete>=-19.5 E lenke_curvatura=curva em s E imc>=14.19 E imc>=30.51
Grupo 3 melhora 2 100.0% dif_colete>=-19.5 E lenke_curvatura=curva em s E imc< 14.19
2.1.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.1.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).