Relatório de Análise da Incapacidade funcional

Autor

Caio Sain Vallio e Vitor Sain Vallio

Data de Publicação

3 de outubro 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 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.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 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.4 Matriz de correlação

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

3 Análise bivariada

3.1 EQ5-D e EVA

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

mdl_eva |>
    parameters::model_parameters() |>
    parameters::print_md()
Parameter Coefficient SE 95% CI t(44) p
(Intercept) 0.84 0.06 (0.72, 0.97) 13.11 < .001
eva -0.05 8.83e-03 (-0.07, -0.03) -5.76 < .001

Para o aumento de 1 ponto na escala EVA, o corre uma redução média de 0,05 pontos na EQ5D (IC 95%: −0,07 a −0,03; p < 0,001).

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

performance::test_performance(mdl_nulo, mdl_eva) |>
    parameters::print_md()
Name Model BF Omega2 p (Omega2) LR p (LR)
mdl_nulo glm
mdl_eva glm > 1000 0.37 > .999 25.87 0.889

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
-19.71 -19.14 -14.22 0.43 0.18 0.19
Código 
mdl_eva |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.286).
Código 
mdl_eva |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.218).
Código 
mdl_eva |>
    performance::check_normality() |>
    plot()

3.2 EQ5-D e B-PCS

Código 
mdl_pcs <- glm(
    eq5_d ~ 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) 0.67 0.07 (0.54, 0.80) 10.24 < .001
b pcs -6.57e-03 2.29e-03 (-0.01, -2.08e-03) -2.87 0.004

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

Código 
performance::test_performance(mdl_nulo, mdl_pcs) |>
    parameters::print_md()
Name Model BF Omega2 p (Omega2) LR p (LR)
mdl_nulo glm
mdl_pcs glm 7.62 0.24 0.500 7.89 0.868

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
-1.73 -1.16 3.75 0.16 0.22 0.23
Código 
mdl_pcs |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.068).
Código 
mdl_pcs |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.965).
Código 
mdl_pcs |>
    performance::check_normality() |>
    plot()

3.3 EQ5-D e CSI

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

mdl_csi |>
    parameters::model_parameters() |>
    parameters::print_md()
Parameter Coefficient SE 95% CI t(44) p
(Intercept) 0.83 0.09 (0.64, 1.01) 8.77 < .001
csi -5.82e-03 1.62e-03 (-9.00e-03, -2.64e-03) -3.59 < .001

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

Código 
performance::test_performance(mdl_nulo, mdl_csi) |>
    parameters::print_md()
Name Model BF Omega2 p (Omega2) LR p (LR)
mdl_nulo glm
mdl_csi glm 53.80 0.23 > .999 11.80 0.844

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
-5.64 -5.07 -0.16 0.23 0.21 0.22
Código 
mdl_csi |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p = 0.018).
Código 
mdl_csi |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.320).
Código 
mdl_csi |>
    performance::check_normality() |>
    plot()

3.4 EQ5-D e Incapacidade

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

mdl_inc |>
    parameters::model_parameters() |>
    parameters::print_md()
Parameter Coefficient SE 95% CI t(44) p
(Intercept) 0.14 0.05 (0.03, 0.24) 2.57 0.010
incapacidade 8.42e-03 1.09e-03 (6.27e-03, 0.01) 7.69 < .001

Para o aumento de 1 ponto na escala Incapacidade, o corre um incremento médio de 0,00842 pontos na incapacidade (IC 95%: 0,00627 a 0,01 ; p < 0,001).

Código 
performance::test_performance(mdl_nulo, mdl_csi) |>
    parameters::print_md()
Name Model BF Omega2 p (Omega2) LR p (LR)
mdl_nulo glm
mdl_csi glm 53.80 0.23 > .999 11.80 0.844

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

Código 
performance::performance(mdl_inc) |>
    parameters::print_md()
Indices of model performance
AIC AICc BIC R2 RMSE Sigma
-33.06 -32.49 -27.57 0.57 0.16 0.16
Código 
mdl_inc |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.986).
Código 
mdl_inc |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.087).
Código 
mdl_inc |>
    performance::check_normality() |>
    plot()

4 Análise multivariada

4.1 EQ5-D e todas as variáveis

Código 
mdl_all <- glm(
    eq5_d ~ eva + b_pcs + csi + incapacidade
    # + 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) 0.52 0.13 (0.26, 0.78) 3.89 < .001
eva -0.03 8.55e-03 (-0.04, -8.33e-03) -2.93 0.003
b pcs -6.51e-04 1.91e-03 (-4.39e-03, 3.08e-03) -0.34 0.733
csi -1.36e-03 1.47e-03 (-4.25e-03, 1.52e-03) -0.93 0.353
incapacidade 5.63e-03 1.32e-03 (3.04e-03, 8.22e-03) 4.25 < .001

Para o aumento de 1 ponto na escala EVA, o corre um decrecimo médio de 0,03 pontos na EQ5-D (IC 95%: -0.04 a -0.00833; p = 0,003).

Para o aumento de 1 ponto na escala INCAPACIDADE, o corre um incremento médio de 0,00563 pontos na EQ5-D (IC 95%: -0.00304 a -0.00822; p = 0,003).

As demais variáveis (CSI e B-PCS) não apresentaram associação estatisticamente significante.

Código 
performance::test_performance(mdl_nulo, mdl_all) |>
    parameters::print_md()
Name Model BF Omega2 p (Omega2) LR p (LR)
mdl_nulo glm
mdl_all glm > 1000 0.76 0.500 50.21 > .999

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_inc) |>
    parameters::print_md()
Name Model BF Omega2 p (Omega2) LR p (LR)
mdl_all glm
mdl_eva glm 0.002 0.59 0.500 2.33 0.010
mdl_pcs glm < 0.001 0.87 0.500 3.34 < .001
mdl_csi glm < 0.001 0.75 0.500 3.27 < .001
mdl_inc glm 1.27 0.13 0.500 2.22 0.013

Each model is compared to mdl_all.

Código 
performance::compare_performance(
  mdl_all, 
  mdl_eva, 
  mdl_pcs, 
  mdl_csi, 
  mdl_inc,
    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_all glm 0.66 0.14 0.15 0.924 0.847 0.439 96.41%
mdl_inc glm 0.57 0.16 0.16 0.076 0.153 0.560 61.68%
mdl_eva glm 0.43 0.18 0.19 9.58e-05 1.94e-04 7.07e-04 25.56%
mdl_csi glm 0.23 0.21 0.22 8.45e-08 1.71e-07 6.24e-07 6.14%
mdl_pcs glm 0.16 0.22 0.23 1.20e-08 2.42e-08 8.83e-08 0.00%
Código 
mdl_all |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.754).
Código 
mdl_all |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.464).
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.48 1.18 2.30 1.22 0.67 0.43 0.85
b_pcs 1.62 1.26 2.49 1.27 0.62 0.40 0.80
csi 1.76 1.34 2.71 1.33 0.57 0.37 0.75
incapacidade 1.73 1.32 2.66 1.32 0.58 0.38 0.76

• Não há colinearidade significativa entre eva, b_pcs, csi e incapacidade.
• 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=-38.06
eq5_d ~ eva + b_pcs + csi + incapacidade

               Df Deviance     AIC
- b_pcs         1  0.90975 -39.926
- csi           1  0.92623 -39.100
<none>             0.90716 -38.057
- eva           1  1.09765 -31.289
- incapacidade  1  1.30742 -23.245

Step:  AIC=-39.93
eq5_d ~ eva + csi + incapacidade

               Df Deviance     AIC
- csi           1  0.94654 -40.102
<none>             0.90975 -39.926
- eva           1  1.10462 -32.998
- incapacidade  1  1.31899 -24.839

Step:  AIC=-40.1
eq5_d ~ eva + incapacidade

               Df Deviance     AIC
<none>             0.94654 -40.102
- eva           1  1.15224 -33.057
- incapacidade  1  1.54018 -19.708
Código 
mdl_step |>
    parameters::model_parameters() |>
    parameters::print_md()
Parameter Coefficient SE 95% CI t(43) p
(Intercept) 0.40 0.10 (0.21, 0.60) 4.04 < .001
eva -0.03 8.49e-03 (-0.04, -9.31e-03) -3.06 0.002
incapacidade 6.32e-03 1.22e-03 (3.93e-03, 8.70e-03) 5.19 < .001

Para o aumento de 1 ponto na escala EVA, o corre um decrecimo médio de 0,03 pontos na EQ5-D (IC 95%: -0.04 a -0.00931; p = 0,002).

Para o aumento de 1 ponto na escala INCAPACIDADE, o corre um incremento médio de 0,00632 pontos na EQ5-D (IC 95%: -0.00393 a -0.0087; p < 0,001).

Código 
performance::performance(mdl_step) |>
    parameters::print_md()
Indices of model performance
AIC AICc BIC R2 RMSE Sigma
-40.10 -39.13 -32.79 0.65 0.14 0.15
Código 
mdl_step |>
    performance::check_autocorrelation()
OK: Residuals appear to be independent and not autocorrelated (p = 0.812).
Código 
mdl_step |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.703).
Código 
mdl_step |>
    performance::check_normality() |>
    plot()