Artigo 1

Resumo

Verificar se existe associação entre cinesiofobia com domínios de função (AKPS, MDP), fatores cognitivo-comportamentais (PSEQ-10, CPSE, Catastrofização) e dor (EVA) em pacientes com dor femoropatelar.

0.1 Bibliotecas

library(tidyverse)
library(readxl)
library(gtsummary)
library(gt)
library(qqplotr)
library(nortest)
library(parameters)
library(performance)

• tidyverse (Wickham, Averick, et al. 2019) para manipulação de dados.
  
• readxl (Wickham, Bryan, et al. 2019) para leitura de arquivos Excel.
  
• gtsummary (Sjoberg et al. 2021) para criação de tabelas e estatísticas descritivas.
• gt (Iannone et al. 2025) para criação de tabelas.
• qqplotr (Almeida, Loy, e Hofmann 2018) para criação deQQ-plots.
• nortest (Gross, Ligges, e Ligges 2015) para testes de normalidade.
• parameters (Lüdecke et al. 2020) para análise de modelos.
• performance (Lüdecke et al. 2021) para análise de modelos.

0.2 Importação dos dados

df <- read_excel("data/Tabela 1 Estatistica - Otavio.xlsx")

# Convertendo a coluna Altura para numérica
df <- df |>
    mutate(Altura = as.numeric(Altura))

Dataset com 66 participantes e 12 variáveis.

0.3 Análise descritiva

0.3.1 Variáveis demográficas

df |>
    gtsummary::tbl_summary(
        include = c(
            Sexo,
            Idade,
            Altura,
            Peso
        ),
        statistic = list(
            gtsummary::all_continuous() ~ "{mean} ({sd})",
            gtsummary::all_categorical() ~ "{n} ({p}%)"
        )
    )

Characteristic	N = 661
Sexo
Feminino	52 (79%)
Masculino	14 (21%)
Idade	25.6 (5.9)
Altura	1.66 (0.09)
Peso	64 (11)
1 n (%); Mean (SD)

0.3.2 Variável dependente

0.3.2.1 Cinesiofobia

0.3.2.1.1 Estatísticas descritivas

df |>
    select(Cinesiofobia) |>
    summarise(
        Média = round(mean(Cinesiofobia), 2),
        Mediana = round(median(Cinesiofobia), 2),
        DesvioPadrão = round(sd(Cinesiofobia), 2),
        Mínimo = round(min(Cinesiofobia), 2),
        Máximo = round(max(Cinesiofobia), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)

Média	38.82
Mediana	38.00
DesvioPadrão	5.54
Mínimo	20.00
Máximo	52.00

0.3.2.1.2 Distribuição

df |>
    ggplot(aes(x = "", y = Cinesiofobia)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição da Cinesiofobia",
        x = "",
        y = "Cinesiofobia"
    ) +
    geom_point(
        aes(
            # color = Cinesiofobia
        ),
        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()
    )

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

0.3.2.1.3 Testes de normalidade

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.

shapiro_test <- shapiro.test(df$Cinesiofobia)
ks_test <- ks.test(
    df$Cinesiofobia,
    "pnorm",
    mean = mean(df$Cinesiofobia),
    sd = sd(df$Cinesiofobia)
)
ad_test <- nortest::ad.test(df$Cinesiofobia)

# 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.962	0.040
Kolmogorov-Smirnov	0.150	0.104
Anderson-Darling	1.006	0.011

0.3.3 Variáveis independentes

0.3.3.1 Domínio de função

0.3.3.1.1 AKPS

Estatísticas descritivas

df |>
    select(AKPS) |>
    summarise(
        Média = round(mean(AKPS), 2),
        Mediana = round(median(AKPS), 2),
        DesvioPadrão = round(sd(AKPS), 2),
        Mínimo = round(min(AKPS), 2),
        Máximo = round(max(AKPS), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)

Média	23.94
Mediana	24.00
DesvioPadrão	4.50
Mínimo	16.00
Máximo	33.00

Distribuição

df |>
    ggplot(aes(x = "", y = AKPS)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição de AKPS",
        x = "",
        y = "AKPS"
    ) +
    geom_point(
        aes(),
        position = position_jitter(w = .15),
    ) +
    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()
    )

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

Testes de normalidade

shapiro_test <- shapiro.test(df$AKPS)
ks_test <- ks.test(df$AKPS, "pnorm", mean = mean(df$AKPS), sd = sd(df$AKPS))
ad_test <- nortest::ad.test(df$AKPS)

# 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.972	0.133
Kolmogorov-Smirnov	0.086	0.720
Anderson-Darling	0.469	0.240

0.3.3.2 Domínio de movimento

0.3.3.2.1 MDP

Estatísticas descritivas

df |>
    select(MDP) |>
    summarise(
        Média = round(mean(MDP), 2),
        Mediana = round(median(MDP), 2),
        DesvioPadrão = round(sd(MDP), 2),
        Mínimo = round(min(MDP), 2),
        Máximo = round(max(MDP), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)

Média	15.86
Mediana	14.44
DesvioPadrão	5.05
Mínimo	8.95
Máximo	31.70

Distribuição

df |>
    ggplot(aes(x = "", y = MDP)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição de MDP",
        x = "",
        y = "MDP"
    ) +
    geom_point(
        aes(),
        position = position_jitter(w = .15),
    ) +
    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()
    )

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

Testes de normalidade

shapiro_test <- shapiro.test(df$MDP)
ks_test <- ks.test(df$MDP, "pnorm", mean = mean(df$MDP), sd = sd(df$MDP))
ad_test <- nortest::ad.test(df$MDP)

# 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.879	0.000
Kolmogorov-Smirnov	0.154	0.077
Anderson-Darling	2.300	0.000

0.3.3.3 Domínio cognitivo-comportamental

0.3.3.3.1 PSEQ-10

Estatísticas descritivas

df |>
    select(PSEQ_10) |>
    summarise(
        Média = round(mean(PSEQ_10), 2),
        Mediana = round(median(PSEQ_10), 2),
        DesvioPadrão = round(sd(PSEQ_10), 2),
        Mínimo = round(min(PSEQ_10), 2),
        Máximo = round(max(PSEQ_10), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)

Média	52.88
Mediana	55.00
DesvioPadrão	6.45
Mínimo	35.00
Máximo	60.00

Distribuição

df |>
    ggplot(aes(x = "", y = PSEQ_10)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição de PSEQ_10",
        x = "",
        y = "PSEQ_10"
    ) +
    geom_point(
        aes(),
        position = position_jitter(w = .15),
    ) +
    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()
    )

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

Testes de normalidade

shapiro_test <- shapiro.test(df$PSEQ_10)
ks_test <- ks.test(
    df$PSEQ_10,
    "pnorm",
    mean = mean(df$PSEQ_10),
    sd = sd(df$PSEQ_10)
)
ad_test <- nortest::ad.test(df$PSEQ_10)

# 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.893	0.00
Kolmogorov-Smirnov	0.159	0.07
Anderson-Darling	2.121	0.00

0.3.3.3.2 CPSE

Estatísticas descritivas

df |>
    select(CPSE) |>
    summarise(
        Média = round(mean(CPSE), 2),
        Mediana = round(median(CPSE), 2),
        DesvioPadrão = round(sd(CPSE), 2),
        Mínimo = round(min(CPSE), 2),
        Máximo = round(max(CPSE), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)

Média	93.83
Mediana	96.50
DesvioPadrão	10.97
Mínimo	60.00
Máximo	110.00

Distribuição

df |>
    ggplot(aes(x = "", y = CPSE)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição de CPSE",
        x = "",
        y = "CPSE"
    ) +
    geom_point(
        aes(),
        position = position_jitter(w = .15),
    ) +
    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()
    )

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

Testes de normalidade

shapiro_test <- shapiro.test(df$CPSE)
ks_test <- ks.test(df$CPSE, "pnorm", mean = mean(df$CPSE), sd = sd(df$CPSE))
ad_test <- nortest::ad.test(df$CPSE)

# 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.944	0.005
Kolmogorov-Smirnov	0.136	0.176
Anderson-Darling	1.167	0.004

0.3.3.3.3 Catastrofizacao

Estatísticas descritivas

df |>
    select(Catastrofizacao) |>
    summarise(
        Média = round(mean(Catastrofizacao), 2),
        Mediana = round(median(Catastrofizacao), 2),
        DesvioPadrão = round(sd(Catastrofizacao), 2),
        Mínimo = round(min(Catastrofizacao), 2),
        Máximo = round(max(Catastrofizacao), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)

Média	28.71
Mediana	28.50
DesvioPadrão	9.25
Mínimo	13.00
Máximo	52.00

Distribuição

df |>
    ggplot(aes(x = "", y = Catastrofizacao)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição de Catastrofizacao",
        x = "",
        y = "Catastrofizacao"
    ) +
    geom_point(
        aes(),
        position = position_jitter(w = .15),
    ) +
    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()
    )

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

Testes de normalidade

shapiro_test <- shapiro.test(df$Catastrofizacao)
ks_test <- ks.test(
    df$Catastrofizacao,
    "pnorm",
    mean = mean(df$Catastrofizacao),
    sd = sd(df$Catastrofizacao)
)
ad_test <- nortest::ad.test(df$Catastrofizacao)

# 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.160
Kolmogorov-Smirnov	0.073	0.870
Anderson-Darling	0.432	0.296

0.3.3.4 Domínio de dor

0.3.3.4.1 EVA 15 dias

Estatísticas descritivas

df |>
    select(EVA_15dias) |>
    summarise(
        Média = round(mean(EVA_15dias), 2),
        Mediana = round(median(EVA_15dias), 2),
        DesvioPadrão = round(sd(EVA_15dias), 2),
        Mínimo = round(min(EVA_15dias), 2),
        Máximo = round(max(EVA_15dias), 2)
    ) |>
    t() |>
    as.data.frame() |>
    rownames_to_column(var = "Métricas") |>
    gt::gt() |>
    gt::tab_options(column_labels.hidden = TRUE)

Média	5.53
Mediana	6.00
DesvioPadrão	1.56
Mínimo	2.00
Máximo	9.00

Distribuição

df |>
    ggplot(aes(x = "", y = EVA_15dias)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição de EVA_15dias",
        x = "",
        y = "EVA_15dias"
    ) +
    geom_point(
        aes(),
        position = position_jitter(w = .15),
    ) +
    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()
    )

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

Testes de normalidade

shapiro_test <- shapiro.test(df$EVA_15dias)
ks_test <- ks.test(
    df$EVA_15dias,
    "pnorm",
    mean = mean(df$EVA_15dias),
    sd = sd(df$EVA_15dias)
)
ad_test <- nortest::ad.test(df$EVA_15dias)

# 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.956	0.019
Kolmogorov-Smirnov	0.164	0.058
Anderson-Darling	1.301	0.002

0.3.3.5 Matriz de correlação

df |>
    select(AKPS, MDP, PSEQ_10, CPSE, Catastrofizacao, EVA_15dias) |>
    GGally::ggpairs(method = "spearman") +
    theme_classic()

0.4 Associação entre variáveis

0.4.1 Cinesiofobia e domínio de função

mdl_nulo <- glm(
    Cinesiofobia ~ 1,
    data = df,
    family = gaussian(link = "identity")
)

mdl_funcao <- glm(
    Cinesiofobia ~ AKPS,
    data = df,
    family = gaussian(link = "identity")
)

performance::test_performance(mdl_nulo, mdl_funcao) |>
    parameters::print_md()

Name	Model	BF	df	df_diff	Chi2	p
mdl_nulo	glm		2
mdl_funcao	glm	0.940	3	1	4.07	0.044

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

performance::performance(mdl_funcao) |>
    parameters::print_md()

Indices of model performance

AIC	AICc	BIC	R2	RMSE	Sigma
414.3	414.7	420.9	0.06	5.33	5.42

mdl_funcao |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.008).

mdl_funcao |>
    performance::check_heteroscedasticity()

OK: Error variance appears to be homoscedastic (p = 0.694).

mdl_funcao |>
    performance::check_residuals()

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

mdl_funcao |>
    performance::check_normality() |>
    plot()

mdl_funcao |>
    parameters::model_parameters() |>
    parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(64)	p
(Intercept)	31.61	3.64	(24.48, 38.74)	8.69	< .001
AKPS	0.30	0.15	(8.44e-03, 0.59)	2.02	0.044

Aumentos de 1 unidade em AKPS estão associados a um acréscimo médio de 0.30 pontos na cinesiofobia.

0.4.2 Cinesiofobia e domínio de movimento

mdl_movimento <- glm(
    Cinesiofobia ~ MDP,
    data = df,
    family = gaussian(link = "identity")
)

performance::test_performance(mdl_nulo, mdl_movimento) |>
    parameters::print_md()

Name	Model	BF	df	df_diff	Chi2	p
mdl_nulo	glm		2
mdl_movimento	glm	2.59	3	1	6.09	0.014

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

performance::test_performance(mdl_movimento, mdl_funcao) |>
    parameters::print_md()

Name	Model	BF
mdl_movimento	glm
mdl_funcao	glm	0.363

Each model is compared to mdl_movimento.

performance::compare_performance(
    mdl_funcao,
    mdl_movimento,
    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_movimento	glm	0.09	5.25	5.34	0.734	0.734	0.734	100.00%
mdl_funcao	glm	0.06	5.33	5.42	0.266	0.266	0.266	0.00%

mdl_movimento |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.018).

mdl_movimento |>
    performance::check_heteroscedasticity()

OK: Error variance appears to be homoscedastic (p = 0.390).

mdl_movimento |>
    performance::check_residuals()

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

mdl_movimento |>
    performance::check_normality() |>
    plot()

mdl_movimento |>
    parameters::model_parameters() |>
    parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(64)	p
(Intercept)	33.64	2.18	(29.37, 37.92)	15.43	< .001
MDP	0.33	0.13	(0.07, 0.58)	2.49	0.013

Aumentos de 1 unidade em MDP estão associados a um acréscimo médio de 0.33 pontos na cinesiofobia.

0.4.3 Cinesiofobia e domínios cognitivo-comportamentais

mdl_cognicao <- glm(
    Cinesiofobia ~ PSEQ_10 + CPSE + Catastrofizacao,
    data = df,
    family = gaussian(link = "identity")
)

performance::test_performance(mdl_nulo, mdl_cognicao) |>
    parameters::print_md()

Name	Model	BF	df	df_diff	Chi2	p
mdl_nulo	glm		2
mdl_cognicao	glm	0.795	5	3	12.11	0.007

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

performance::test_performance(mdl_cognicao, mdl_funcao, mdl_movimento) |>
    parameters::print_md()

Name	Model	BF
mdl_cognicao	glm
mdl_funcao	glm	1.18
mdl_movimento	glm	3.26

Each model is compared to mdl_cognicao.

performance::compare_performance(
    mdl_funcao,
    mdl_movimento,
    mdl_cognicao,
    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_cognicao	glm	0.17	5.02	5.18	0.668	0.597	0.184	83.33%
mdl_movimento	glm	0.09	5.25	5.34	0.244	0.296	0.599	42.02%
mdl_funcao	glm	0.06	5.33	5.42	0.088	0.107	0.217	1.35%

mdl_cognicao |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.034).

mdl_cognicao |>
    performance::check_heteroscedasticity()

OK: Error variance appears to be homoscedastic (p = 0.450).

mdl_cognicao |>
    performance::check_residuals()

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

mdl_cognicao |>
    performance::check_normality() |>
    plot()

mdl_cognicao |>
    parameters::model_parameters() |>
    parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(62)	p
(Intercept)	55.83	9.15	(37.90, 73.75)	6.10	< .001
PSEQ 10	-0.38	0.14	(-0.66, -0.10)	-2.68	0.007
CPSE	0.03	0.09	(-0.14, 0.21)	0.36	0.721
Catastrofizacao	7.94e-03	0.09	(-0.16, 0.18)	0.09	0.928

Cada unidade adicional no PSEQ 10 está associada a uma redução média de 0.38 pontos na cinesiofobia, considerando CPSE e Catastrofizacao constantes.

0.4.4 Cinesiofobia e dor

mdl_dor <- glm(
    Cinesiofobia ~ EVA_15dias,
    data = df,
    family = gaussian(link = "identity")
)

performance::test_performance(mdl_nulo, mdl_dor) |>
    parameters::print_md()

Name	Model	BF	df	df_diff	Chi2	p
mdl_nulo	glm		2
mdl_dor	glm	1.03	3	1	4.24	0.039

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

performance::test_performance(
    mdl_cognicao,
    mdl_dor,
    mdl_funcao,
    mdl_movimento
) |>
    parameters::print_md()

Name	Model	BF
mdl_cognicao	glm
mdl_dor	glm	1.29
mdl_funcao	glm	1.18
mdl_movimento	glm	3.26

Each model is compared to mdl_cognicao.

performance::compare_performance(
    mdl_funcao,
    mdl_movimento,
    mdl_cognicao,
    mdl_dor,
    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_cognicao	glm	0.17	5.02	5.18	0.609	0.534	0.148	83.33%
mdl_movimento	glm	0.09	5.25	5.34	0.222	0.265	0.484	42.02%
mdl_dor	glm	0.06	5.33	5.41	0.088	0.105	0.192	3.99%
mdl_funcao	glm	0.06	5.33	5.42	0.081	0.096	0.176	1.35%

mdl_dor |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.028).

mdl_dor |>
    performance::check_heteroscedasticity()

OK: Error variance appears to be homoscedastic (p = 0.054).

mdl_dor |>
    performance::check_residuals()

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

mdl_dor |>
    performance::check_normality() |>
    plot()

mdl_dor |>
    parameters::model_parameters() |>
    parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(64)	p
(Intercept)	33.92	2.47	(29.08, 38.76)	13.74	< .001
EVA 15dias	0.89	0.43	(0.04, 1.73)	2.06	0.039

Um aumento de 1 unidade na avaliação da dor está associado a um acréscimo médio de 0.89 pontos na cinesiofobia.

0.4.5 Cinesiofobia e todas as variáveis

mdl_all <- glm(
    Cinesiofobia ~ AKPS + MDP + PSEQ_10 + CPSE + Catastrofizacao + EVA_15dias,
    data = df,
    family = gaussian(link = "identity")
)

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	0.036	8	6	18.49	0.005

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

performance::test_performance(
    mdl_all,
    mdl_funcao,
    mdl_movimento,
    mdl_cognicao,
    mdl_dor
) |>
    parameters::print_md()

Name	Model	BF
mdl_all	glm
mdl_funcao	glm	26.13
mdl_movimento	glm	72.02
mdl_cognicao	glm	22.09
mdl_dor	glm	28.53

Each model is compared to mdl_all.

performance::compare_performance(
    mdl_funcao,
    mdl_movimento,
    mdl_cognicao,
    mdl_dor,
    mdl_all,
    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.24	4.78	5.06	0.424	0.231	0.007	74.46%
mdl_cognicao	glm	0.17	5.02	5.18	0.351	0.411	0.147	65.40%
mdl_movimento	glm	0.09	5.25	5.34	0.128	0.204	0.481	35.55%
mdl_dor	glm	0.06	5.33	5.41	0.051	0.081	0.191	7.76%
mdl_funcao	glm	0.06	5.33	5.42	0.046	0.074	0.174	5.90%

mdl_all |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.024).

mdl_all |>
    performance::check_heteroscedasticity()

OK: Error variance appears to be homoscedastic (p = 0.926).

mdl_all |>
    performance::check_residuals()

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

mdl_all |>
    performance::check_normality() |>
    plot()

mdl_all |>
    parameters::model_parameters() |>
    parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(59)	p
(Intercept)	49.15	9.76	(30.03, 68.27)	5.04	< .001
AKPS	0.16	0.16	(-0.15, 0.46)	0.99	0.322
MDP	0.48	0.30	(-0.11, 1.07)	1.59	0.112
PSEQ 10	-0.39	0.14	(-0.66, -0.11)	-2.76	0.006
CPSE	0.05	0.09	(-0.12, 0.22)	0.55	0.579
Catastrofizacao	-5.53e-03	0.09	(-0.18, 0.17)	-0.06	0.950
EVA 15dias	-1.02	0.99	(-2.96, 0.92)	-1.03	0.302

Apenas o PSEQ 10 é significativo na predição da cinesiofobia.

A variável AKPS perdeu a associação estatisticamente significante.

As variáveis catastrofizacao e EVA 15 dias mudaram a direção da associação, sugerindo multicolinearidade.

0.4.5.1 Análise de multicolinearidade

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
AKPS	1.27	1.08	1.88	1.13	0.79	0.53	0.92
MDP	5.83	4.05	8.63	2.41	0.17	0.12	0.25
PSEQ_10	2.07	1.58	2.96	1.44	0.48	0.34	0.63
CPSE	2.39	1.79	3.45	1.55	0.42	0.29	0.56
Catastrofizacao	1.70	1.34	2.42	1.30	0.59	0.41	0.74
EVA_15dias	6.04	4.20	8.96	2.46	0.17	0.11	0.24

As variáveis EVA 15 dias e MDP apresentam multicolinearidade.

0.4.5.2 Removendo variáveis com multicolinearidade

Removendo EVA 15 dias

mdl_sem_eva <- glm(
    Cinesiofobia ~ AKPS + MDP + PSEQ_10 + CPSE + Catastrofizacao,
    data = df,
    family = gaussian(link = "identity")
)

performance::test_performance(mdl_all, mdl_sem_eva) |>
    parameters::print_md()

Name	Model	BF	df	df_diff	Chi2	p
mdl_all	glm		8
mdl_sem_eva	glm	4.50	7	-1	1.18	0.277

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

performance::compare_performance(
    mdl_all,
    mdl_sem_eva,
    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.24	4.78	5.06	0.399	0.330	0.182	50.00%
mdl_sem_eva	glm	0.23	4.83	5.06	0.601	0.670	0.818	50.00%

mdl_sem_eva |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.008).

mdl_sem_eva |>
    performance::check_heteroscedasticity()

OK: Error variance appears to be homoscedastic (p = 0.924).

mdl_sem_eva |>
    performance::check_residuals()

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

mdl_sem_eva |>
    performance::check_normality() |>
    plot()

mdl_sem_eva |>
    parameters::model_parameters() |>
    parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(60)	p
(Intercept)	47.77	9.67	(28.82, 66.72)	4.94	< .001
AKPS	0.15	0.16	(-0.16, 0.46)	0.95	0.341
MDP	0.21	0.14	(-0.08, 0.49)	1.43	0.153
PSEQ 10	-0.38	0.14	(-0.65, -0.10)	-2.71	0.007
CPSE	0.05	0.09	(-0.12, 0.22)	0.58	0.559
Catastrofizacao	-0.02	0.09	(-0.19, 0.15)	-0.24	0.807

Removendo a variável EVA 15 dias, somente o PSEQ 10 é significativo.

Removendo MDP

mdl_sem_mdp <- glm(
    Cinesiofobia ~ AKPS + PSEQ_10 + CPSE + Catastrofizacao + EVA_15dias,
    data = df,
    family = gaussian(link = "identity")
)

performance::test_performance(mdl_all, mdl_sem_eva, mdl_sem_mdp) |>
    parameters::print_md()

Name	Model	BF
mdl_all	glm
mdl_sem_eva	glm	4.50
mdl_sem_mdp	glm	2.04

Each model is compared to mdl_all.

performance::compare_performance(
    mdl_all,
    mdl_sem_eva,
    mdl_sem_mdp,
    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_sem_eva	glm	0.23	4.83	5.06	0.473	0.514	0.597	85.21%
mdl_all	glm	0.24	4.78	5.06	0.314	0.253	0.133	57.66%
mdl_sem_mdp	glm	0.21	4.88	5.12	0.214	0.233	0.270	4.93%

mdl_sem_mdp |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.018).

mdl_sem_mdp |>
    performance::check_heteroscedasticity()

OK: Error variance appears to be homoscedastic (p = 0.699).

mdl_sem_mdp |>
    performance::check_residuals()

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

mdl_sem_mdp |>
    performance::check_normality() |>
    plot()

mdl_sem_mdp |>
    parameters::model_parameters() |>
    parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(60)	p
(Intercept)	48.69	9.88	(29.33, 68.05)	4.93	< .001
AKPS	0.20	0.16	(-0.11, 0.51)	1.25	0.210
PSEQ 10	-0.37	0.14	(-0.65, -0.10)	-2.65	0.008
CPSE	0.04	0.09	(-0.13, 0.22)	0.47	0.641
Catastrofizacao	-0.02	0.09	(-0.20, 0.15)	-0.27	0.784
EVA 15dias	0.36	0.48	(-0.58, 1.30)	0.75	0.452

Removendo a variável MDP, somente o PSEQ 10 é significativo.

Removendo EVA 15 dias e MDP

mdl_sem_mdp_eva <- glm(
    Cinesiofobia ~ AKPS + PSEQ_10 + CPSE + Catastrofizacao,
    data = df,
    family = gaussian(link = "identity")
)

performance::test_performance(
    mdl_all,
    mdl_sem_eva,
    mdl_sem_mdp,
    mdl_sem_mdp_eva
) |>
    parameters::print_md()

Name	Model	BF
mdl_all	glm
mdl_sem_eva	glm	4.50
mdl_sem_mdp	glm	2.04
mdl_sem_mdp_eva	glm	12.15

Each model is compared to mdl_all.

performance::compare_performance(
    mdl_all,
    mdl_sem_eva,
    mdl_sem_mdp,
    mdl_sem_mdp_eva,
    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_sem_eva	glm	0.23	4.83	5.06	0.331	0.322	0.229	72.57%
mdl_all	glm	0.24	4.78	5.06	0.220	0.159	0.051	57.38%
mdl_sem_mdp_eva	glm	0.20	4.91	5.10	0.299	0.374	0.617	51.80%
mdl_sem_mdp	glm	0.21	4.88	5.12	0.150	0.146	0.103	7.74%

mdl_sem_mdp_eva |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.020).

mdl_sem_mdp_eva |>
    performance::check_heteroscedasticity()

OK: Error variance appears to be homoscedastic (p = 0.581).

mdl_sem_mdp_eva |>
    performance::check_residuals()

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

mdl_sem_mdp_eva |>
    performance::check_normality() |>
    plot()

mdl_sem_mdp_eva |>
    parameters::model_parameters() |>
    parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(61)	p
(Intercept)	50.42	9.57	(31.66, 69.18)	5.27	< .001
AKPS	0.24	0.14	(-0.04, 0.53)	1.68	0.092
PSEQ 10	-0.38	0.14	(-0.65, -0.10)	-2.67	0.008
CPSE	0.03	0.09	(-0.14, 0.20)	0.35	0.724
Catastrofizacao	-0.02	0.09	(-0.19, 0.16)	-0.19	0.852

Removendo a variável EVA 15 dias e MDP, somente o PSEQ 10 é significativo.

Removendo PSEQ 10

mdl_sem_pseq <- glm(
    Cinesiofobia ~ AKPS + MDP + CPSE + Catastrofizacao + EVA_15dias,
    data = df,
    family = gaussian(link = "identity")
)

performance::test_performance(
    mdl_all,
    mdl_sem_eva,
    mdl_sem_mdp,
    mdl_sem_mdp_eva,
    mdl_sem_pseq
) |>
    parameters::print_md()

Name	Model	BF
mdl_all	glm
mdl_sem_eva	glm	4.50
mdl_sem_mdp	glm	2.04
mdl_sem_mdp_eva	glm	12.15
mdl_sem_pseq	glm	0.146

Each model is compared to mdl_all.

performance::compare_performance(
    mdl_all,
    mdl_sem_eva,
    mdl_sem_mdp,
    mdl_sem_mdp_eva,
    mdl_sem_pseq,
    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_sem_eva	glm	0.23	4.83	5.06	0.328	0.318	0.227	82.11%
mdl_sem_mdp_eva	glm	0.20	4.91	5.10	0.296	0.370	0.612	81.86%
mdl_all	glm	0.24	4.78	5.06	0.218	0.157	0.050	68.84%
mdl_sem_mdp	glm	0.21	4.88	5.12	0.148	0.144	0.103	51.01%
mdl_sem_pseq	glm	0.15	5.08	5.33	0.011	0.010	0.007	0.00%

mdl_sem_pseq |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.028).

mdl_sem_pseq |>
    performance::check_heteroscedasticity()

OK: Error variance appears to be homoscedastic (p = 0.820).

mdl_sem_pseq |>
    performance::check_residuals()

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

mdl_sem_pseq |>
    performance::check_normality() |>
    plot()

mdl_sem_pseq |>
    parameters::model_parameters() |>
    parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(60)	p
(Intercept)	40.25	9.71	(21.22, 59.27)	4.15	< .001
AKPS	0.17	0.17	(-0.16, 0.49)	1.02	0.308
MDP	0.43	0.32	(-0.19, 1.05)	1.37	0.172
CPSE	-0.09	0.08	(-0.24, 0.06)	-1.17	0.240
Catastrofizacao	0.03	0.09	(-0.15, 0.21)	0.36	0.722
EVA 15dias	-0.87	1.04	(-2.91, 1.17)	-0.84	0.403

Removendo a variável PSEQ 10, nenhuma variável é significativa.

0.4.5.3 Seleção de variáveis

Método stepwise

mdl_step <- step(mdl_all, direction = "both", trace = FALSE)

performance::test_performance(
    mdl_step,
    mdl_all,
    mdl_sem_eva,
    mdl_sem_mdp,
    mdl_sem_mdp_eva,
    mdl_sem_pseq
) |>
    parameters::print_md()

Name	Model	BF
mdl_step	glm
mdl_all	glm	< 0.001
mdl_sem_eva	glm	0.004
mdl_sem_mdp	glm	0.002
mdl_sem_mdp_eva	glm	0.011
mdl_sem_pseq	glm	< 0.001

Each model is compared to mdl_step.

performance::compare_performance(
    mdl_step,
    mdl_all,
    mdl_sem_eva,
    mdl_sem_mdp,
    mdl_sem_mdp_eva,
    mdl_sem_pseq,
    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_step	glm	0.21	4.88	5.00	0.747	0.844	0.982	88.80%
mdl_all	glm	0.24	4.78	5.06	0.055	0.024	9.07e-04	48.67%
mdl_sem_eva	glm	0.23	4.83	5.06	0.083	0.050	0.004	45.01%
mdl_sem_mdp_eva	glm	0.20	4.91	5.10	0.075	0.058	0.011	33.98%
mdl_sem_mdp	glm	0.21	4.88	5.12	0.038	0.022	0.002	33.91%
mdl_sem_pseq	glm	0.15	5.08	5.33	0.003	0.002	1.32e-04	0.00%

mdl_step |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.030).

mdl_step |>
    performance::check_heteroscedasticity()

OK: Error variance appears to be homoscedastic (p = 0.903).

mdl_step |>
    performance::check_residuals()

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

mdl_step |>
    performance::check_normality() |>
    plot()

mdl_step |>
    parameters::model_parameters() |>
    parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(63)	p
(Intercept)	51.36	6.00	(39.61, 63.12)	8.56	< .001
MDP	0.24	0.13	(-5.70e-03, 0.49)	1.91	0.056
PSEQ 10	-0.31	0.10	(-0.50, -0.12)	-3.14	0.002

O método stepwise selecionou as variáveis MDP e PSEQ 10.

Para o acréscimo de 1 unidade em PSEQ 10, a cinesiofobia diminui em 0.31 pontos.

Para o acréscimo de 1 unidade em MDP, a cinesiofobia aumenta em 0.24 pontos, porém não é estatisticamente significante.

0.5 Conclusão

0.5.1 Modelo inicial completo

O primeiro ajuste incluiu todas as variáveis candidatas (AKPS, MDP, PSEQ-10, CPSE, catastrofização e EVA-15 dias). Embora o coeficiente de determinação ( R² = 24 %, RMSE = 4,78) indicasse algum poder explicativo, a inspeção de colinearidade revelou potenciais fragilidades:

Termo	VIF	Intervalo 95 %	Tolerância
AKPS	1,27	1,08 – 1,88	0,79
MDP	5,83	4,05 – 8,63	0,17
PSEQ-10	2,07	1,58 – 2,96	0,48
CPSE	2,39	1,79 – 3,45	0,42
Catastrofização	1,70	1,34 – 2,42	0,59
EVA-15 dias	6,04	4,20 – 8,96	0,17

Limiares usuais: VIF > 5 sinaliza multicolinearidade substantiva; tolerância < 0,2 reforça essa suspeita. Assim, MDP e EVA-15 dias apresentaram forte sobreposição de informação com outros preditores, aumentando erros-padrão e reduzindo precisão dos coeficientes – o que ajuda a explicar p-valores não significativos mesmo com efeitos aparentes.

0.5.2 Aplicação do procedimento stepwise

Para mitigar redundância e obter um modelo mais parcimonioso, empregou-se o stepwise bidirecional guiado por AIC/BIC. Esse algoritmo:

• avalia sequencialmente a contribuição conjunta de cada variável para o ajuste global;
• remove termos cujo ganho de informação não compense o aumento de complexidade;
• tende a reduzir multicolinearidade, pois exclui variáveis que competem por mesma parcela de variância explicada.

O método eliminou AKPS, CPSE, catastrofização e EVA-15 dias. A ligeira redução de R² (24 % → 21 %) e aumento marginal no RMSE (4,78 → 4,88) foram compensados por menores AIC/BIC, indicando melhor equilíbrio ajuste-complexidade.

0.5.3 Modelo final após seleção

O modelo resultante retém dois preditores:

Termo	β	IC 95 %	p
Intercepto	51,36	39,61 – 63,12	—
PSEQ-10	–0,31	–0,50 – –0,12	< 0,01
MDP	0,24	–0,006 – 0,49	0,06

• PSEQ-10 permaneceu estatisticamente significativo, corroborando a hipótese de que maior autoeficácia para dor associa-se a menor cinesiofobia.

• MDP apresentou p = 0,06 e intervalo que tangencia o zero, mas foi mantido porque:

  1. Reduziu AIC/BIC mesmo após remoção de variáveis redundantes – sinal de contribuição informativa adicional.
  2. Seu coeficiente permaneceu estável (β ≈ 0,24–0,48) entre modelos, sugerindo efeito consistente ainda que impreciso.
  3. A retirada de EVA-15 dias diminuiu a colinearidade global; espera-se que o VIF do MDP caia (< 5) no modelo reduzido, aumentando confiabilidade na estimativa.

0.5.4 Interpretação e implicações

O resultado indica que autoeficácia (PSEQ-10) é o principal fator explicativo da cinesiofobia. O MDP, embora borderline em significância, fornece indícios de que padrões de movimento possam exercer influência independente – hipótese que merece investigação com amostras maiores.

Síntese para a discussão do artigo

“No modelo inicial, a análise de multicolinearidade evidenciou VIF > 5 para MDP e EVA-15 dias, sugerindo sobreposição de informação e possível inflacionamento dos erros-padrão. A aplicação do procedimento stepwise, orientado por AIC/BIC, resultou em um modelo mais parcimonioso (R² = 21 %, RMSE = 4,88) que manteve o PSEQ-10 como preditor negativo significativo e conservou o MDP, cujo efeito permaneceu estável (β = 0,24) e contribuiu para leve melhoria nos critérios de informação. Apesar de a estimativa do MDP não alcançar significância convencional (p = 0,06), sua inclusão foi justificada pelo ganho de ajuste, pela redução simultânea da colinearidade e pela plausibilidade de captar aspectos de padrão de movimento relacionados à cinesiofobia. Assim, interpretamos o MDP como associação potencial e recomendamos estudos futuros para confirmar sua relevância.”

1 Análises adficionais

Análises sugereidas pelo revisor

\(R^2\) ajustado

mdl_step |>
    performance::r2(adjusted = TRUE)

  R2: 0.212

F-teste global do modelo

# F-teste global do modelo
mdl_nulo <- glm(
    Cinesiofobia ~ 1,
    data = df,
    family = gaussian(link = "identity")
)
anova(mdl_nulo, mdl_step, test = "F")

Diagnostic residuals

# Verificar padrão da autocorrelação
car::durbinWatsonTest(mdl_step)

 lag Autocorrelation D-W Statistic p-value
   1      -0.2563231      2.507395   0.016
 Alternative hypothesis: rho != 0

mdl_step |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.038).

# Plotar resíduos vs ordem de observação
plot(
    residuals(mdl_step),
    type = "o",
    main = "Resíduos vs Ordem das Observações",
    ylab = "Resíduos",
    xlab = "Índice"
)
abline(h = 0, col = "red", lty = 2)

# Verificar se há padrões por grupos
df$residuos <- residuals(mdl_step)
df$ordem <- 1:nrow(df)
ggplot(df, aes(x = ordem, y = residuos)) +
    geom_point() +
    geom_line() +
    geom_hline(yintercept = 0, color = "red", linetype = "dashed")

# Erros-padrão robustos à heterocedasticidade E autocorrelação
mdl_step_robust <- parameters::model_parameters(
    mdl_step,
    vcov = "HC3",
    ci_method = "wald"
)

# Comparar resultados
print("Modelo original:")

[1] "Modelo original:"

model_parameters(mdl_step)

print("\nModelo com erros robustos:")

[1] "\nModelo com erros robustos:"

mdl_step_robust

standardized effect size

# ============================================================================
# 1. COEFICIENTES PADRONIZADOS (Beta Standardized)
# ============================================================================

# Usando parameters com standardize = "refit"
mdl_step_std <- parameters::model_parameters(
    mdl_step,
    standardize = "refit", # Refit com variáveis padronizadas
    vcov = "HC3", # Erros robustos
    ci_method = "wald"
)

print("Coeficientes padronizados com erros robustos:")

[1] "Coeficientes padronizados com erros robustos:"

mdl_step_std |> parameters::print_md()

Parameter	Coefficient	SE	95% CI	t(63)	p
(Intercept)	-1.50e-16	0.11	(-0.23, 0.23)	-1.33e-15	> .999
MDP	0.22	0.08	(0.05, 0.39)	2.62	0.011
PSEQ 10	-0.36	0.13	(-0.62, -0.10)	-2.76	0.008

# ============================================================================
# 2. COHEN'S f² (Tamanho de efeito local para cada preditor)
# ============================================================================

# Cohen's f² = (R²_completo - R²_sem_variável) / (1 - R²_completo)
# Interpretação: 0.02 = pequeno, 0.15 = médio, 0.35 = grande

# Calcular Cohen's f² para cada preditor
r2_full <- performance::r2(mdl_step)$R2

# Modelo sem PSEQ_10
mdl_sem_pseq10 <- glm(
    Cinesiofobia ~ MDP,
    data = df,
    family = gaussian(link = "identity")
)
r2_sem_pseq <- performance::r2(mdl_sem_pseq10)$R2
f2_pseq <- (r2_full - r2_sem_pseq) / (1 - r2_full)

# Modelo sem MDP
mdl_sem_mdp2 <- glm(
    Cinesiofobia ~ PSEQ_10,
    data = df,
    family = gaussian(link = "identity")
)
r2_sem_mdp <- performance::r2(mdl_sem_mdp2)$R2
f2_mdp <- (r2_full - r2_sem_mdp) / (1 - r2_full)

cat(sprintf("Cohen's f² PSEQ_10: %.3f\n", f2_pseq))

Cohen's f² PSEQ_10: 0.157

cat(sprintf("Cohen's f² MDP: %.3f\n", f2_mdp))

Cohen's f² MDP: 0.058

# ============================================================================
# 3. PARTIAL R² e SEMI-PARTIAL R²
# ============================================================================

# Partial R² (contribuição única de cada variável)
effectsize::eta_squared(mdl_step, partial = TRUE, ci = 0.95)

# ============================================================================
# 4. TABELA COMPLETA COM TODOS OS EFEITOS
# ============================================================================

# Criar tabela consolidada
library(dplyr)

# Obter coeficientes não padronizados com erros robustos
params_raw <- parameters::model_parameters(
    mdl_step,
    vcov = "HC3",
    ci_method = "wald"
) |>
    as.data.frame()

# Obter coeficientes padronizados com erros robustos
params_std <- parameters::model_parameters(
    mdl_step,
    standardize = "refit",
    vcov = "HC3",
    ci_method = "wald"
) |>
    as.data.frame()

# Combinar resultados
results_table <- params_raw |>
    filter(Parameter != "(Intercept)") |>
    select(Parameter, Coefficient, SE, CI_low, CI_high, t, p) |>
    mutate(
        Std_Coefficient = params_std$Coefficient[
            params_std$Parameter != "(Intercept)"
        ],
        Cohens_f2 = c(f2_mdp, f2_pseq),
        Effect_Size = case_when(
            Cohens_f2 < 0.02 ~ "Trivial",
            Cohens_f2 < 0.15 ~ "Pequeno",
            Cohens_f2 < 0.35 ~ "Médio",
            TRUE ~ "Grande"
        )
    )

print(results_table)

  Parameter Coefficient         SE      CI_low     CI_high         t
1       MDP   0.2409546 0.09206073  0.05698588  0.42492326  2.617344
2   PSEQ_10  -0.3094811 0.11232312 -0.53394093 -0.08502123 -2.755275
            p Std_Coefficient  Cohens_f2 Effect_Size
1 0.011082462       0.2193312 0.05819063     Pequeno
2 0.007658072      -0.3599827 0.15675280       Médio

# ============================================================================
# 5. VISUALIZAÇÃO DOS TAMANHOS DE EFEITO
# ============================================================================

library(ggplot2)

# Plot dos coeficientes padronizados
plot(mdl_step_std, show_labels = TRUE) +
    labs(
        title = "Coeficientes Padronizados com IC 95% (Erros Robustos)",
        x = "Coeficiente Padronizado (β)",
        y = ""
    )

# Plot dos Cohen's f²
df_f2 <- data.frame(
    Variable = c("PSEQ_10", "MDP"),
    Cohens_f2 = c(f2_pseq, f2_mdp),
    Effect = factor(
        c(
            ifelse(
                f2_pseq < 0.02,
                "Trivial",
                ifelse(
                    f2_pseq < 0.15,
                    "Pequeno",
                    ifelse(f2_pseq < 0.35, "Médio", "Grande")
                )
            ),
            ifelse(
                f2_mdp < 0.02,
                "Trivial",
                ifelse(
                    f2_mdp < 0.15,
                    "Pequeno",
                    ifelse(f2_mdp < 0.35, "Médio", "Grande")
                )
            )
        ),
        levels = c("Trivial", "Pequeno", "Médio", "Grande")
    )
)

ggplot(df_f2, aes(x = Variable, y = Cohens_f2, fill = Effect)) +
    geom_col() +
    geom_hline(
        yintercept = c(0.02, 0.15, 0.35),
        linetype = "dashed",
        alpha = 0.5
    ) +
    scale_fill_manual(
        values = c(
            "Trivial" = "#d3d3d3",
            "Pequeno" = "#90EE90",
            "Médio" = "#FFD700",
            "Grande" = "#FF6347"
        )
    ) +
    labs(
        title = "Tamanho de Efeito Local (Cohen's f²)",
        y = "Cohen's f²",
        x = "",
        caption = "Linhas pontilhadas: 0.02 (pequeno), 0.15 (médio), 0.35 (grande)"
    ) +
    theme_minimal()

# ============================================================================
# 6. INTERPRETAÇÃO DOS COEFICIENTES PADRONIZADOS
# ============================================================================

1.0.0.1 INTERPRETAÇÃO DOS TAMANHOS DE EFEITO

1. COEFICIENTES PADRONIZADOS (β):.
- Representam a mudança em desvios-padrão de Y para cada.
mudança de 1 DP em X (controlando outras variáveis).
- |β| < 0.1: efeito pequeno.
- |β| < 0.3: efeito médio. - |β| ≥ 0.5: efeito grande.

2. COHEN’S f²:.
- Mede a contribuição única de cada preditor.
- f² < 0.02: efeito trivial.
- f² < 0.15: efeito pequeno.
- f² < 0.35: efeito médio.
- f² ≥ 0.35: efeito grande.

1.0.0.2 Tamanhos de Efeito Padronizados

Para melhor interpretação dos resultados, calculamos os coeficientes padronizados utilizando erros-padrão robustos (HC3) e métricas complementares de tamanho de efeito (Tabela X).

Os coeficientes padronizados indicaram que:

• PSEQ-10 (β = -0.36, IC 95% [-0.62, -0.10], t(63) = -2.76, p = 0.008): Para cada aumento de 1 desvio-padrão no PSEQ-10, a cinesiofobia diminui 0.36 desvios-padrão, mantidas constantes as demais variáveis. Este efeito é classificado como médio segundo os critérios de Cohen.

• MDP (β = 0.22, IC 95% [0.05, 0.39], t(63) = 2.62, p = 0.011): Para cada aumento de 1 desvio-padrão no MDP, a cinesiofobia aumenta 0.22 desvios-padrão, controlando-se para autoeficácia. Este efeito é classificado como pequeno a médio.

A análise de Cohen’s f² revelou que PSEQ-10 apresenta contribuição única média (f² = 0.16), enquanto MDP mostrou contribuição pequena (f² = 0.06). Os valores de eta-quadrado parcial confirmaram esses padrões, com PSEQ-10 (η²p = 0.14) explicando 14% da variância residual e MDP (η²p = 0.10) explicando 10%.

Em conjunto, estes resultados indicam que a autoeficácia para dor (PSEQ-10) é o preditor mais robusto e de maior magnitude na explicação da cinesiofobia, enquanto o padrão de movimento disfuncional (MDP) contribui de forma independente, embora com menor intensidade.

---

Tabela X. Coeficientes e tamanhos de efeito do modelo final com erros-padrão robustos

Variável	B (não padronizado)	EP	β (padronizado)	IC 95%	t(63)	p	Cohen’s f²	η²p	Interpretação
MDP	0.24	0.09	0.22	[0.05, 0.39]	2.62	0.011	0.06	0.10	Pequeno
PSEQ-10	-0.31	0.11	-0.36	[-0.62, -0.10]	-2.76	0.008	0.16	0.14	Médio

Nota. B = coeficiente não padronizado; EP = erro-padrão robusto (HC3); β = coeficiente padronizado; IC = intervalo de confiança; η²p = eta-quadrado parcial. Critérios de Cohen para f²: < 0.02 = trivial, 0.02-0.15 = pequeno, 0.15-0.35 = médio, > 0.35 = grande.

1.0.0.3 Tamanhos de Efeito Padronizados e Cohen’s f²

A análise dos coeficientes padronizados com erros robustos (HC3) revelou que PSEQ-10 (β = -0.36, IC 95% [-0.62, -0.10], p = 0.008) apresentou efeito médio sobre a cinesiofobia, enquanto MDP (β = 0.22, IC 95% [0.05, 0.39], p = 0.011) mostrou efeito pequeno.

A contribuição única de cada preditor, avaliada por Cohen’s f², confirmou que PSEQ-10 (f² = 0.16, médio) explica maior proporção da variância do que MDP (f² = 0.06, pequeno). Os valores de eta-quadrado parcial foram η²p = 0.14 para PSEQ-10 e η²p = 0.10 para MDP, indicando que a autoeficácia para dor é o principal fator associado à cinesiofobia, com padrão de movimento exercendo influência independente de menor magnitude.

Para discussão do artigo

Os tamanhos de efeito padronizados corroboram a centralidade da autoeficácia para dor na explicação da cinesiofobia. O coeficiente padronizado de -0.36 para PSEQ-10 indica que incrementos na crença de capacidade funcional associam-se a reduções clinicamente relevantes no medo do movimento. O MDP, embora com efeito menor (β = 0.22), demonstra contribuição independente (f² = 0.06), sugerindo que aspectos biomecânicos do movimento podem ter papel complementar aos fatores cognitivos na manutenção da cinesiofobia.

Referências

Almeida, Alexandre, Adam Loy, e Heike Hofmann. 2018. “ggplot2 Compatible Quantile-Quantile Plots in”.

Gross, Juergen, Uwe Ligges, e Maintainer Uwe Ligges. 2015. “Package ‘nortest’”. Five omnibus tests for testing the composite hypothesis of normality.

Iannone, Richard, Joe Cheng, Barret Schloerke, Ellis Hughes, Alexandra Lauer, JooYoung Seo, Ken Brevoort, e Olivier Roy. 2025. gt: Easily Create Presentation-Ready Display Tables. https://gt.rstudio.com.

Lüdecke, Daniel, Mattan S Ben-Shachar, Indrajeet Patil, e Dominique Makowski. 2020. “Extracting, computing and exploring the parameters of statistical models using R”. Journal of Open Source Software 5 (53): 2445.

Lüdecke, Daniel, Mattan S Ben-Shachar, Indrajeet Patil, Philip Waggoner, e Dominique Makowski. 2021. “performance: An R package for assessment, comparison and testing of statistical models”. Journal of open source software 6 (60).

Sjoberg, Daniel D, Karissa Whiting, Michael Curry, Jessica A Lavery, e Joseph Larmarange. 2021. “Reproducible summary tables with the gtsummary package”. The R journal 13 (1): 570–80.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the Tidyverse”. Journal of open source software 4 (43): 1686.

Wickham, Hadley, Jennifer Bryan, Marcin Kalicinski, Komarov Valery, Christophe Leitienne, Bob Colbert, David Hoerl, Evan Miller, e Maintainer Jennifer Bryan. 2019. “Package ‘readxl’”. Version 13: 1.