Relatório de Análise da Incapacidade funcional

Autor

Caio Sain Vallio e Vitor Sain Vallio

Data de Publicação

27 de junho de 2025

Código 
# Libraries
library(tidyverse)
library(readxl)
library(janitor)
library(gtsummary)
library(ggplot2)
library(qqplotr)
library(mice)
library(car)

# Read data
df0 <- janitor::clean_names(read_excel("Dados_Ms_Tami.xlsx"))

# Alterando a natureza dos dados
df1 <- df0 |>
  mutate(
    sexo         = as.factor(sexo),
    motivo       = as.factor(motivo),
    comorb       = as.factor(comorb),
    queixas_msk  = as.factor(queixas_msk),
    absenteismo  = as.factor(absenteismo),
    eq5_d        = as.numeric(eq5_d)
  )

# Método padrão PMM (Predictive Mean Matching) para dados contínuos
imputacao <- mice(df1, m = 5, method = "pmm", seed = 123)

# Extração do banco imputado (1 dos 5, ou média posterior com pooling depois)
df2 <- complete(imputacao, 1)

1 Banco de dados

Código 
skimr::skim(df2)
Data summary
Name df2
Number of rows 46
Number of columns 16
_______________________
Column type frequency:
factor 5
numeric 11
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
sexo 0 1 FALSE 2 1: 33, 0: 13
motivo 0 1 FALSE 2 1: 29, 0: 17
comorb 0 1 FALSE 2 1: 36, 0: 10
queixas_msk 0 1 FALSE 2 1: 33, 0: 13
absenteismo 0 1 FALSE 2 1: 32, 0: 14

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
idade 0 1 61.96 11.09 40.00 55.00 59.00 70.00 86.00 ▂▇▇▃▃
altura_m 0 1 1.65 0.11 1.42 1.59 1.65 1.70 1.92 ▂▃▇▂▂
peso_kg 0 1 73.18 12.30 40.00 67.00 75.00 82.75 94.00 ▁▃▅▇▅
imc 0 1 26.75 3.85 17.69 24.98 26.81 28.65 37.02 ▂▅▇▂▁
tempo_de_espera 0 1 5.76 7.08 0.00 1.00 3.50 7.75 35.00 ▇▂▁▁▁
t_absente 0 1 10.80 11.22 0.00 0.00 9.00 20.00 43.00 ▇▂▅▁▁
eva 0 1 6.59 3.16 0.00 5.00 8.00 9.00 10.00 ▅▂▃▇▇
b_pcs 0 1 24.50 14.80 0.00 14.75 23.50 32.00 59.00 ▅▇▇▃▂
csi 0 1 54.54 20.02 5.00 41.00 55.50 66.00 97.00 ▂▅▇▆▂
incapacidade 0 1 44.01 22.04 9.85 26.18 44.40 53.42 92.00 ▆▅▇▂▂
eq5_d 0 1 0.51 0.25 0.07 0.32 0.50 0.71 1.00 ▆▇▇▇▂

2 Análise descritiva

2.1 Variáveis demográficas

Código 
df2 %>% 
    gtsummary::tbl_summary(
        include = c(
            sexo,
            idade,
            altura_m,
            peso_kg,
            imc,
            comorb,
            tempo_de_espera,
            queixas_msk,
            absenteismo,
            t_absente
        ),
        statistic = list(
            gtsummary::all_continuous() ~ "{mean} ({sd})",
            gtsummary::all_categorical() ~ "{n} ({p}%)"
        )
    )
Characteristic N = 461
sexo
    0 13 (28%)
    1 33 (72%)
idade 62 (11)
altura_m 1.65 (0.11)
peso_kg 73 (12)
imc 26.8 (3.8)
comorb
    0 10 (22%)
    1 36 (78%)
tempo_de_espera 5.8 (7.1)
queixas_msk
    0 13 (28%)
    1 33 (72%)
absenteismo
    0 14 (30%)
    1 32 (70%)
t_absente 11 (11)
1 n (%); Mean (SD)

2.2 Variável depentente

2.2.1 Incapacidade

Código 
df2 |>
    select(incapacidade) |>
    summarise(
        Média = round(mean(incapacidade), 2),
        Mediana = round(median(incapacidade), 2),
        DesvioPadrão = round(sd(incapacidade), 2),
        Mínimo = round(min(incapacidade), 2),
        Máximo = round(max(incapacidade), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)
Média 44.01
Mediana 44.40
DesvioPadrão 22.04
Mínimo 9.85
Máximo 92.00
Código 
df2 |>
    ggplot(aes(x = "", y = incapacidade)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição da Incapacidade",
        x = "", y = "Incapacidade"
    ) +
    geom_point(
        aes(
            # color = incapacidade
        ),
        position = position_jitter(w = .15),
        # size = 2
    ) +
    geom_boxplot(
        width = .25,
        # fill = "steelblue",
        alpha = 0.7,
        outlier.shape = NA
    ) +
    coord_flip() +
    theme_minimal() +
    theme(
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
    )

Código 
df2 |>
    ggplot(aes(sample = incapacidade)) +
    qqplotr::stat_qq_band(alpha = 0.2) +
    qqplotr::stat_qq_line() +
    qqplotr::stat_qq_point() +
    labs(
        title = "QQ-Plot de Incapacidade",
        x = "Quantis teóricos",
        y = "Quantis amostrais"
    ) +
    theme_minimal()

Teste de Shapiro-Wilk, Kolmogorov-Smirnov e Anderson-Darling consideram a hipótese nula de que a distribuição é normal. Ou seja, se o valor-p for maior que 0.05, não há evidências para rejeitar a hipótese nula. Em outras palavras, se o valor-p for maior que 0.05, a distribuição é normal.

Código 
shapiro_test <- shapiro.test(df2$incapacidade)
ks_test <- ks.test(df2$incapacidade, "pnorm",
    mean = mean(df2$incapacidade),
    sd = sd(df2$incapacidade)
)
ad_test <- nortest::ad.test(df2$incapacidade)

# Criar tabela com resultados dos testes
data.frame(
    Teste = c(
        "Shapiro-Wilk",
        "Kolmogorov-Smirnov",
        "Anderson-Darling"
    ),
    Estatística = c(
        round(shapiro_test$statistic, 3),
        round(ks_test$statistic, 3),
        round(ad_test$statistic, 3)
    ),
    `Valor-p` = c(
        round(shapiro_test$p.value, 3),
        round(ks_test$p.value, 3),
        round(ad_test$p.value, 3)
    )
) |>
    gt::gt()
Teste Estatística Valor.p
Shapiro-Wilk 0.948 0.040
Kolmogorov-Smirnov 0.133 0.392
Anderson-Darling 0.683 0.070

2.3 Variáveis indepententes

2.3.1 EVA

Código 
df2 |>
    select(eva) |>
    summarise(
        Média = round(mean(eva), 2),
        Mediana = round(median(eva), 2),
        DesvioPadrão = round(sd(eva), 2),
        Mínimo = round(min(eva), 2),
        Máximo = round(max(eva), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)
Média 6.59
Mediana 8.00
DesvioPadrão 3.16
Mínimo 0.00
Máximo 10.00
Código 
df2 |>
    ggplot(aes(x = "", y = eva)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição da eva",
        x = "", y = "eva"
    ) +
    geom_point(
        aes(
            # color = incapacidade
        ),
        position = position_jitter(w = .15),
        # size = 2
    ) +
    geom_boxplot(
        width = .25,
        # fill = "steelblue",
        alpha = 0.7,
        outlier.shape = NA
    ) +
    coord_flip() +
    theme_minimal() +
    theme(
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
    )

Código 
df2 |>
    ggplot(aes(sample = eva)) +
    qqplotr::stat_qq_band(alpha = 0.2) +
    qqplotr::stat_qq_line() +
    qqplotr::stat_qq_point() +
    labs(
        title = "QQ-Plot de eva",
        x = "Quantis teóricos",
        y = "Quantis amostrais"
    ) +
    theme_minimal()

Código 
shapiro_test <- shapiro.test(df2$eva)
ks_test <- ks.test(df2$eva, "pnorm",
    mean = mean(df2$eva),
    sd = sd(df2$eva)
)
ad_test <- nortest::ad.test(df2$eva)

# Criar tabela com resultados dos testes
data.frame(
    Teste = c(
        "Shapiro-Wilk",
        "Kolmogorov-Smirnov",
        "Anderson-Darling"
    ),
    Estatística = c(
        round(shapiro_test$statistic, 3),
        round(ks_test$statistic, 3),
        round(ad_test$statistic, 3)
    ),
    `Valor-p` = c(
        round(shapiro_test$p.value, 3),
        round(ks_test$p.value, 3),
        round(ad_test$p.value, 3)
    )
) |>
    gt::gt()
Teste Estatística Valor.p
Shapiro-Wilk 0.860 0.000
Kolmogorov-Smirnov 0.238 0.011
Anderson-Darling 2.260 0.000

2.3.2 B-PCS

Código 
df2 |>
    select(b_pcs) |>
    summarise(
        Média = round(mean(b_pcs), 2),
        Mediana = round(median(b_pcs), 2),
        DesvioPadrão = round(sd(b_pcs), 2),
        Mínimo = round(min(b_pcs), 2),
        Máximo = round(max(b_pcs), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)
Média 24.5
Mediana 23.5
DesvioPadrão 14.8
Mínimo 0.0
Máximo 59.0
Código 
df2 |>
    ggplot(aes(x = "", y = b_pcs)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição da b_pcs",
        x = "", y = "b_pcs"
    ) +
    geom_point(
        aes(
            # color = incapacidade
        ),
        position = position_jitter(w = .15),
        # size = 2
    ) +
    geom_boxplot(
        width = .25,
        # fill = "steelblue",
        alpha = 0.7,
        outlier.shape = NA
    ) +
    coord_flip() +
    theme_minimal() +
    theme(
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
    )

Código 
df2 |>
    ggplot(aes(sample = b_pcs)) +
    qqplotr::stat_qq_band(alpha = 0.2) +
    qqplotr::stat_qq_line() +
    qqplotr::stat_qq_point() +
    labs(
        title = "QQ-Plot de b_pcs",
        x = "Quantis teóricos",
        y = "Quantis amostrais"
    ) +
    theme_minimal()

Código 
shapiro_test <- shapiro.test(df2$b_pcs)
ks_test <- ks.test(df2$b_pcs, "pnorm",
    mean = mean(df2$b_pcs),
    sd = sd(df2$b_pcs)
)
ad_test <- nortest::ad.test(df2$b_pcs)

# Criar tabela com resultados dos testes
data.frame(
    Teste = c(
        "Shapiro-Wilk",
        "Kolmogorov-Smirnov",
        "Anderson-Darling"
    ),
    Estatística = c(
        round(shapiro_test$statistic, 3),
        round(ks_test$statistic, 3),
        round(ad_test$statistic, 3)
    ),
    `Valor-p` = c(
        round(shapiro_test$p.value, 3),
        round(ks_test$p.value, 3),
        round(ad_test$p.value, 3)
    )
) |>
    gt::gt()
Teste Estatística Valor.p
Shapiro-Wilk 0.975 0.416
Kolmogorov-Smirnov 0.077 0.947
Anderson-Darling 0.293 0.588

2.3.3 CSI

Código 
df2 |>
    select(csi) |>
    summarise(
        Média = round(mean(csi), 2),
        Mediana = round(median(csi), 2),
        DesvioPadrão = round(sd(csi), 2),
        Mínimo = round(min(csi), 2),
        Máximo = round(max(csi), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)
Média 54.54
Mediana 55.50
DesvioPadrão 20.02
Mínimo 5.00
Máximo 97.00
Código 
df2 |>
    ggplot(aes(x = "", y = csi)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição da csi",
        x = "", y = "csi"
    ) +
    geom_point(
        aes(
            # color = incapacidade
        ),
        position = position_jitter(w = .15),
        # size = 2
    ) +
    geom_boxplot(
        width = .25,
        # fill = "steelblue",
        alpha = 0.7,
        outlier.shape = NA
    ) +
    coord_flip() +
    theme_minimal() +
    theme(
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
    )

Código 
df2 |>
    ggplot(aes(sample = csi)) +
    qqplotr::stat_qq_band(alpha = 0.2) +
    qqplotr::stat_qq_line() +
    qqplotr::stat_qq_point() +
    labs(
        title = "QQ-Plot de csi",
        x = "Quantis teóricos",
        y = "Quantis amostrais"
    ) +
    theme_minimal()

Código 
shapiro_test <- shapiro.test(df2$csi)
ks_test <- ks.test(df2$csi, "pnorm",
    mean = mean(df2$csi),
    sd = sd(df2$csi)
)
ad_test <- nortest::ad.test(df2$csi)

# Criar tabela com resultados dos testes
data.frame(
    Teste = c(
        "Shapiro-Wilk",
        "Kolmogorov-Smirnov",
        "Anderson-Darling"
    ),
    Estatística = c(
        round(shapiro_test$statistic, 3),
        round(ks_test$statistic, 3),
        round(ad_test$statistic, 3)
    ),
    `Valor-p` = c(
        round(shapiro_test$p.value, 3),
        round(ks_test$p.value, 3),
        round(ad_test$p.value, 3)
    )
) |>
    gt::gt()
Teste Estatística Valor.p
Shapiro-Wilk 0.986 0.864
Kolmogorov-Smirnov 0.076 0.953
Anderson-Darling 0.206 0.862

2.3.4 EQ5-D

Código 
df2 |>
    select(eq5_d) |>
    summarise(
        Média = round(mean(eq5_d), 2),
        Mediana = round(median(eq5_d), 2),
        DesvioPadrão = round(sd(eq5_d), 2),
        Mínimo = round(min(eq5_d), 2),
        Máximo = round(max(eq5_d), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)
Média 0.51
Mediana 0.50
DesvioPadrão 0.25
Mínimo 0.07
Máximo 1.00
Código 
df2 |>
    ggplot(aes(x = "", y = eq5_d)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição da eq5_d",
        x = "", y = "eq5_d"
    ) +
    geom_point(
        aes(
            # color = incapacidade
        ),
        position = position_jitter(w = .15),
        # size = 2
    ) +
    geom_boxplot(
        width = .25,
        # fill = "steelblue",
        alpha = 0.7,
        outlier.shape = NA
    ) +
    coord_flip() +
    theme_minimal() +
    theme(
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
    )

Código 
df2 |>
    ggplot(aes(sample = eq5_d)) +
    qqplotr::stat_qq_band(alpha = 0.2) +
    qqplotr::stat_qq_line() +
    qqplotr::stat_qq_point() +
    labs(
        title = "QQ-Plot de eq5_d",
        x = "Quantis teóricos",
        y = "Quantis amostrais"
    ) +
    theme_minimal()

Código 
shapiro_test <- shapiro.test(df2$eq5_d)
ks_test <- ks.test(df2$eq5_d, "pnorm",
    mean = mean(df2$eq5_d),
    sd = sd(df2$eq5_d)
)
ad_test <- nortest::ad.test(df2$eq5_d)

# Criar tabela com resultados dos testes
data.frame(
    Teste = c(
        "Shapiro-Wilk",
        "Kolmogorov-Smirnov",
        "Anderson-Darling"
    ),
    Estatística = c(
        round(shapiro_test$statistic, 3),
        round(ks_test$statistic, 3),
        round(ad_test$statistic, 3)
    ),
    `Valor-p` = c(
        round(shapiro_test$p.value, 3),
        round(ks_test$p.value, 3),
        round(ad_test$p.value, 3)
    )
) |>
    gt::gt()
Teste Estatística Valor.p
Shapiro-Wilk 0.973 0.370
Kolmogorov-Smirnov 0.117 0.553
Anderson-Darling 0.338 0.489

2.4 Matriz de correlação

Código 
df2 |>
    select(incapacidade, csi, b_pcs, eva, eq5_d) |>
    GGally::ggpairs(method = "spearman") +
    theme_classic()

3 Análise bivariada

3.1 Incapacidade e EVA

Código 
mdl_eva <- glm(
    incapacidade ~ eva,
    data = df2, family = gaussian(link = "identity")
)

mdl_eva |>
    parameters::model_parameters() |>
    parameters::print_md()
Parameter Coefficient SE 95% CI t(44) p
(Intercept) 70.00 6.33 (57.60, 82.39) 11.07 < .001
eva -3.95 0.87 (-5.65, -2.25) -4.55 < .001

Para o aumento de 1 ponto na escala EVA, o corre uma redução média de 3,95 pontos na incapacidade (IC 95%: −5,65 a −2,25; p < 0,001).

Código 
mdl_nulo <- glm(
    incapacidade ~ 1,
    data = df2, family = gaussian(link = "identity")
)

performance::test_performance(mdl_nulo, mdl_eva) |>
    parameters::print_md()
Name Model BF df df_diff Chi2 p
mdl_nulo glm 2
mdl_eva glm > 1000 3 1.00 17.72 < .001

Models were detected as nested (in terms of fixed parameters) and are compared in sequential order.

Código 
performance::performance(mdl_eva) |>
    parameters::print_md()
Indices of model performance
AIC AICc BIC R2 RMSE Sigma
402.37 402.94 407.86 0.32 17.98 18.39
Código 
mdl_eva |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.918).
Código 
mdl_eva |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.699).
Código 
mdl_eva |>
    performance::check_residuals()
OK: Simulated residuals appear as uniformly distributed (p = 0.907).
Código 
mdl_eva |>
    performance::check_normality() |>
    plot()

3.2 Incapacidade e B-PCS

Código 
mdl_pcs <- glm(
    incapacidade ~ b_pcs,
    data = df2, family = gaussian(link = "identity")
)

mdl_pcs |>
    parameters::model_parameters() |>
    parameters::print_md()
Parameter Coefficient SE 95% CI t(44) p
(Intercept) 58.15 5.91 (46.57, 69.73) 9.84 < .001
b pcs -0.58 0.21 (-0.98, -0.17) -2.79 0.005

Para o aumento de 1 ponto na escala B-PCS, o corre uma redução média de 0,58 pontos na incapacidade (IC 95%: −0,98 a −0,17; p = 0,005).

Código 
performance::test_performance(mdl_nulo, mdl_pcs) |>
    parameters::print_md()
Name Model BF df df_diff Chi2 p
mdl_nulo glm 2
mdl_pcs glm 6.22 3 1.00 7.49 0.006

Models were detected as nested (in terms of fixed parameters) and are compared in sequential order.

Código 
performance::performance(mdl_pcs) |>
    parameters::print_md()
Indices of model performance
AIC AICc BIC R2 RMSE Sigma
412.61 413.18 418.10 0.15 20.10 20.55
Código 
mdl_pcs |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.926).
Código 
mdl_pcs |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.512).
Código 
mdl_pcs |>
    performance::check_residuals()
OK: Simulated residuals appear as uniformly distributed (p = 0.522).
Código 
mdl_pcs |>
    performance::check_normality() |>
    plot()

3.3 Incapacidade e CSI

Código 
mdl_csi <- glm(
    incapacidade ~ csi,
    data = df2, family = gaussian(link = "identity")
)

mdl_csi |>
    parameters::model_parameters() |>
    parameters::print_md()
Parameter Coefficient SE 95% CI t(44) p
(Intercept) 72.45 8.48 (55.82, 89.07) 8.54 < .001
csi -0.52 0.15 (-0.81, -0.23) -3.57 < .001

Para o aumento de 1 ponto na escala CSI, o corre uma redução média de 0,52 pontos na incapacidade (IC 95%: −0,81 a −0,23; p < 0,001).

Código 
performance::test_performance(mdl_nulo, mdl_csi) |>
    parameters::print_md()
Name Model BF df df_diff Chi2 p
mdl_nulo glm 2
mdl_csi glm 50.87 3 1.00 11.69 < .001

Models were detected as nested (in terms of fixed parameters) and are compared in sequential order.

Código 
performance::performance(mdl_csi) |>
    parameters::print_md()
Indices of model performance
AIC AICc BIC R2 RMSE Sigma
408.41 408.98 413.89 0.22 19.20 19.63
Código 
mdl_csi |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.936).
Código 
mdl_csi |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.192).
Código 
mdl_csi |>
    performance::check_residuals()
OK: Simulated residuals appear as uniformly distributed (p = 0.250).
Código 
mdl_csi |>
    performance::check_normality() |>
    plot()

3.4 Incapacidade e EQ5-D

Código 
mdl_eq5 <- glm(
    incapacidade ~ eq5_d,
    data = df2, family = gaussian(link = "identity")
)

mdl_eq5 |>
    parameters::model_parameters() |>
    parameters::print_md()
Parameter Coefficient SE 95% CI t(44) p
(Intercept) 9.36 4.99 (-0.41, 19.14) 1.88 0.061
eq5 d 68.13 8.86 (50.78, 85.49) 7.69 < .001

Para o aumento de 1 ponto na escala EQ-5D, o corre um incremento médio de 68,13 pontos na incapacidade (IC 95%: 50,78 a 85,49; p < 0,001).

Código 
performance::test_performance(mdl_nulo, mdl_csi) |>
    parameters::print_md()
Name Model BF df df_diff Chi2 p
mdl_nulo glm 2
mdl_csi glm 50.87 3 1.00 11.69 < .001

Models were detected as nested (in terms of fixed parameters) and are compared in sequential order.

Código 
performance::performance(mdl_eq5) |>
    parameters::print_md()
Indices of model performance
AIC AICc BIC R2 RMSE Sigma
380.88 381.45 386.37 0.57 14.24 14.56
Código 
mdl_eq5 |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.104).
Código 
mdl_eq5 |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.840).
Código 
mdl_eq5 |>
    performance::check_residuals()
OK: Simulated residuals appear as uniformly distributed (p = 0.287).
Código 
mdl_eq5 |>
    performance::check_normality() |>
    plot()

4 Análise multivariada

4.1 Incapacidade e todas as variáveis

Código 
mdl_all <- glm(
    incapacidade ~ eva + b_pcs + csi + eq5_d 
    # + idade 
    # + sexo
    # + motivo
    # + imc
    # + comorb
    # + tempo_de_espera
    # + queixas_msk
    # + absenteismo
    # + t_absente
    ,data = df2, family = gaussian(link = "identity")
)

mdl_all |>
    parameters::model_parameters() |>
    parameters::print_md()
Parameter Coefficient SE 95% CI t(41) p
(Intercept) 30.79 14.48 (2.42, 59.17) 2.13 0.033
eva -0.83 0.91 (-2.62, 0.96) -0.91 0.364
b pcs -0.05 0.19 (-0.42, 0.31) -0.29 0.775
csi -0.14 0.14 (-0.42, 0.14) -0.97 0.331
eq5 d 54.37 12.78 (29.32, 79.43) 4.25 < .001

Para o aumento de 1 ponto na escala EQ-5D, o corre um incremento médio de 54,37 pontos na incapacidade (IC 95%: 29,32 a 79,43; p < 0,001).
Todas as outras variáveis não apresentaram associação estatisticamente significante.

Código 
performance::test_performance(mdl_nulo, mdl_all) |>
    parameters::print_md()
Name Model BF df df_diff Chi2 p
mdl_nulo glm 2
mdl_all glm > 1000 6 4.00 42.08 < .001

Models were detected as nested (in terms of fixed parameters) and are compared in sequential order.

Código 
performance::test_performance(
  mdl_all, 
  mdl_eva, 
  mdl_pcs, 
  mdl_csi, 
  mdl_eq5) |>
    parameters::print_md()
Name Model BF
mdl_all glm
mdl_eva glm 0.002
mdl_pcs glm < 0.001
mdl_csi glm < 0.001
mdl_eq5 glm 74.42

Each model is compared to mdl_all.

Código 
performance::compare_performance(
  mdl_all, 
  mdl_eva, 
  mdl_pcs, 
  mdl_csi, 
  mdl_eq5,
    rank = TRUE, verbose = FALSE
) |>
    parameters::print_md()
Comparison of Model Performance Indices
Name Model R2 RMSE Sigma AIC weights AICc weights BIC weights Performance-Score
mdl_eq5 glm 0.57 14.24 14.56 0.827 0.914 0.987 97.89%
mdl_all glm 0.60 13.80 14.62 0.173 0.086 0.013 55.11%
mdl_eva glm 0.32 17.98 18.39 1.78e-05 1.97e-05 2.13e-05 17.91%
mdl_csi glm 0.22 19.20 19.63 8.71e-07 9.62e-07 1.04e-06 7.68%
mdl_pcs glm 0.15 20.10 20.55 1.07e-07 1.18e-07 1.27e-07 0.00%
Código 
mdl_all |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.078).
Código 
mdl_all |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.718).
Código 
mdl_all |>
    performance::check_residuals()
OK: Simulated residuals appear as uniformly distributed (p = 0.459).
Código 
mdl_all |>
    performance::check_normality() |>
    plot()

4.2 Análise de multicolinearidade

Código 
mdl_all |>
    performance::check_collinearity() |>
    parameters::print_md()
Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance Tolerance_CI_low Tolerance_CI_high
eva 1.76 1.34 2.70 1.33 0.57 0.37 0.75
b_pcs 1.62 1.26 2.49 1.27 0.62 0.40 0.80
csi 1.76 1.34 2.70 1.33 0.57 0.37 0.75
eq5_d 2.07 1.52 3.18 1.44 0.48 0.31 0.66

• Não há colinearidade significativa entre eva, b_pcs, csi e eq5_d.
• As variáveis podem ser mantidas no modelo sem risco de distorção nos coeficientes por dependência estatística entre elas.

4.3 Stepwise

Código 
mdl_step <- step(mdl_all, direction = "backward")
Start:  AIC=384.01
incapacidade ~ eva + b_pcs + csi + eq5_d

        Df Deviance    AIC
- b_pcs  1   8778.2 382.11
- eva    1   8936.8 382.93
- csi    1   8962.6 383.06
<none>       8760.7 384.01
- eq5_d  1  12626.2 398.83

Step:  AIC=382.11
incapacidade ~ eva + csi + eq5_d

        Df Deviance    AIC
- eva    1   8960.5 381.05
- csi    1   9142.7 381.98
<none>       8778.2 382.11
- eq5_d  1  12727.0 397.19

Step:  AIC=381.05
incapacidade ~ csi + eq5_d

        Df Deviance    AIC
- csi    1   9324.0 380.88
<none>       8960.5 381.05
- eq5_d  1  16962.6 408.41

Step:  AIC=380.88
incapacidade ~ eq5_d

        Df Deviance    AIC
<none>         9324 380.88
- eq5_d  1    21869 418.10
Código 
mdl_step |>
    parameters::model_parameters() |>
    parameters::print_md()
Parameter Coefficient SE 95% CI t(44) p
(Intercept) 9.36 4.99 (-0.41, 19.14) 1.88 0.061
eq5 d 68.13 8.86 (50.78, 85.49) 7.69 < .001
Código 
performance::performance(mdl_step) |>
    parameters::print_md()
Indices of model performance
AIC AICc BIC R2 RMSE Sigma
380.88 381.45 386.37 0.57 14.24 14.56
Código 
mdl_step |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.116).
Código 
mdl_step |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.840).
Código 
mdl_step |>
    performance::check_residuals()
OK: Simulated residuals appear as uniformly distributed (p = 0.287).
Código 
mdl_step |>
    performance::check_normality() |>
    plot()