Artigo 2

Resumo

Verificar se a intensidade da dor durante o movimento de Step Down Lateral está mais associada ao movimento de pacientes com dor femoropatelar quando comparada a dor autorrelatada de 15 dias.

1 Bibliotecas

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

library(emmeans)
library(ggplot2)
library(dplyr)

• 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.
• lme4 (Bates et al. 2015) para análise de modelos lineares mistos.
• lmerTest (Kuznetsova, Brockhoff, e Christensen 2017) para análise de modelos lineares mistos.

2 Importação dos dados

df <- read_excel("data/Tabela 2 Estatística - Otavio.xlsx")

3 Ajuste do dataset

df_long <- df |>
    pivot_longer(
        cols = c(
            SDL_01_MDP, NPRS_SDL01,
            SDL_02_MDP, NPRS_SDL02,
            SDL_03_MDP, NPRS_SDL03,
            SDL_04_MDP, NPRS_SDL04
        ),
        names_to = "variable",
        values_to = "value"
    ) |>
    mutate(
        type = case_when(
            str_detect(variable, "MDP") ~ "MDP",
            str_detect(variable, "NPRS") ~ "NPRS"
        ),
        session = case_when(
            str_detect(variable, "01") ~ 1,
            str_detect(variable, "02") ~ 2,
            str_detect(variable, "03") ~ 3,
            str_detect(variable, "04") ~ 4
        )
    ) |>
    select(-variable) |>
    pivot_wider(
        id_cols = c(ID, session),
        names_from = type,
        values_from = value
    ) |>
    arrange(ID, session)


df_long02 <- df |>
    # Primeiro criar registros para EVA_15dias como sessão 0
    mutate(
        session = 0,
        NPRS = EVA_15dias,
    ) |>
    select(ID, session, NPRS) |>
    # Depois processar as outras sessões normalmente
    bind_rows(
        df |>
            pivot_longer(
                cols = c(
                    NPRS_SDL01,
                    NPRS_SDL02,
                    NPRS_SDL03,
                    NPRS_SDL04
                ),
                names_to = "variable",
                values_to = "value"
            ) |>
            mutate(
                type = case_when(
                    str_detect(variable, "NPRS") ~ "NPRS"
                ),
                session = case_when(
                    str_detect(variable, "01") ~ 1,
                    str_detect(variable, "02") ~ 2,
                    str_detect(variable, "03") ~ 3,
                    str_detect(variable, "04") ~ 4
                )
            ) |>
            select(-variable) |>
            pivot_wider(
                id_cols = c(ID, session),
                names_from = type,
                values_from = value
            )
    ) |>
    arrange(ID, session)


df_long03 <- df |>
    # Primeiro criar registros para NPRS_Inicial como sessão 0
    mutate(
        session = 0,
        NPRS = NPRS_Inicial,
    ) |>
    select(ID, session, NPRS) |>
    # Depois processar as outras sessões normalmente
    bind_rows(
        df |>
            pivot_longer(
                cols = c(
                    NPRS_SDL01,
                    NPRS_SDL02,
                    NPRS_SDL03,
                    NPRS_SDL04
                ),
                names_to = "variable",
                values_to = "value"
            ) |>
            mutate(
                type = case_when(
                    str_detect(variable, "NPRS") ~ "NPRS"
                ),
                session = case_when(
                    str_detect(variable, "01") ~ 1,
                    str_detect(variable, "02") ~ 2,
                    str_detect(variable, "03") ~ 3,
                    str_detect(variable, "04") ~ 4
                )
            ) |>
            select(-variable) |>
            pivot_wider(
                id_cols = c(ID, session),
                names_from = type,
                values_from = value
            )
    ) |>
    arrange(ID, session)

4 Análise descritiva

4.1 Escala de dor dos últimos 15 dias

4.1.1 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.55
Mediana	6.00
DesvioPadrão	1.53
Mínimo	3.00
Máximo	9.00

4.1.1.1 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 da EVA_15dias",
        x = "", y = "EVA_15dias"
    ) +
    geom_point(
        aes(
            # color = EVA_15dias
        ),
        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 = 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()

4.1.1.2 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$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.943	0.004
Kolmogorov-Smirnov	0.162	0.062
Anderson-Darling	1.418	0.001

4.2 Escala de dor antes do movimento

4.2.1 Estatísticas descritivas

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

Média	1.53
Mediana	0.50
DesvioPadrão	1.99
Mínimo	0.00
Máximo	7.00

4.2.1.1 Distribuição

df |>
    ggplot(aes(x = "", y = NPRS_Inicial)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribuição da NPRS_Inicial",
        x = "", y = "NPRS_Inicial"
    ) +
    geom_point(
        aes(
            # color = NPRS_Inicial
        ),
        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 = NPRS_Inicial)) +
    qqplotr::stat_qq_band(alpha = 0.2) +
    qqplotr::stat_qq_line() +
    qqplotr::stat_qq_point() +
    labs(
        title = "QQ-Plot de NPRS_Inicial",
        x = "Quantis teóricos",
        y = "Quantis amostrais"
    ) +
    theme_minimal()

4.2.1.2 Testes de normalidade

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

# 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.771	0
Kolmogorov-Smirnov	0.279	0
Anderson-Darling	6.176	0

4.3 Escala de dor durante o movimento

4.3.1 Estatísticas descritivas

df_long |>
    group_by(session) |>
    summarise(
        Média = round(mean(NPRS), 2),
        Mediana = round(median(NPRS), 2),
        DesvioPadrão = round(sd(NPRS), 2),
        Mínimo = round(min(NPRS), 2),
        Máximo = round(max(NPRS), 2)
    ) |>
    gt::gt()

session	Média	Mediana	DesvioPadrão	Mínimo	Máximo
1	2.32	2	2.24	0	8
2	2.80	2	2.32	0	8
3	3.29	3	2.33	0	9
4	3.59	3	2.40	0	10

4.3.2 Distribuição

df_long |>
    ggplot() +
    aes(
        x = factor(session, levels = sort(unique(session), decreasing = TRUE)),
        y = NPRS
    ) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2), alpha = .2
    ) +
    geom_point(
        position = position_jitter(w = .15), alpha = .2
    ) +
    geom_boxplot(width = .25, alpha = .2, outlier.shape = NA) +
    coord_flip() +
    theme_minimal() +
    theme(
        legend.position = "none"
    ) +
    labs(x = "Sessão")

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

4.3.3 Testes de normalidade

run_normality_tests <- function(data) {
    shapiro_test <- shapiro.test(data)
    ks_test <- ks.test(data, "pnorm", mean = mean(data), sd = sd(data))
    ad_test <- nortest::ad.test(data)

    return(data.frame(
        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)
        )
    ))
}

normality_results <- df_long %>%
    group_by(session) %>%
    summarise(
        test_results = list(run_normality_tests(NPRS)),
        .groups = "drop"
    )

results_table <- data.frame()

for (i in 1:nrow(normality_results)) {
    session_data <- normality_results$test_results[[i]]
    session_data$Teste <- c(
        "Shapiro-Wilk",
        "Kolmogorov-Smirnov",
        "Anderson-Darling"
    )
    session_data$Sessão <- normality_results$session[i]
    results_table <- rbind(results_table, session_data)
}

results_table %>%
    select(Sessão, Teste, Estatística, `Valor-p` = `Valor.p`) %>%
    arrange(Teste) %>%
    group_by(Teste) %>%
    gt::gt() %>%
    gt::tab_header(title = "Testes de Normalidade por Sessão") %>%
    gt::cols_label(
        Teste = "Teste de Normalidade",
        Sessão = "Sessão",
        Estatística = "Estatística",
        `Valor-p` = "Valor-p"
    ) %>%
    gt::tab_spanner(
        label = "Resultados",
        columns = c(Estatística, `Valor-p`)
    )

Testes de Normalidade por Sessão
Sessão	Resultados
	Estatística	Valor-p
Anderson-Darling
1	2.499	0.000
2	1.696	0.000
3	1.017	0.010
4	0.801	0.036
Kolmogorov-Smirnov
1	0.163	0.061
2	0.161	0.066
3	0.134	0.184
4	0.127	0.235
Shapiro-Wilk
1	0.883	0.000
2	0.916	0.000
3	0.948	0.008
4	0.958	0.024

4.3.4 Análise de variância

modelo1 <- lme4::lmer(NPRS ~ factor(session) + (1 | ID), data = df_long)

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

# Fixed Effects

Parameter	Coefficient	SE	95% CI	t(258)	p
(Intercept)	2.32	0.29	(1.76, 2.88)	8.11	< .001
session [2]	0.48	0.12	(0.26, 0.71)	4.21	< .001
session [3]	0.97	0.12	(0.74, 1.20)	8.43	< .001
session [4]	1.27	0.12	(1.05, 1.50)	11.06	< .001

# Random Effects

Parameter	Coefficient	SE	95% CI
SD (Intercept: ID)	2.23	0.20	(1.87, 2.65)
SD (Residual)	0.66	0.03	(0.60, 0.73)

modelo1 |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p < .001).

modelo1 |>
    performance::check_heteroscedasticity()

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

modelo1 |>
    performance::check_residuals()

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

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

Efeitos fixos:

Intercepto (Sessão 1):
– Valor estimado de 2.32, representando a média da escala de dor na sessão 1.

Sessão 2:
– Coeficiente de 0.48 (p < 0.001), indicando que, em média, a dor aumenta 0.48 pontos na sessão 2 em relação à sessão 1.

Sessão 3:
– Coeficiente de 0.97 (p < 0.001), sugerindo um aumento médio de 0.97 pontos na dor comparado à sessão 1.

Sessão 4:
– Coeficiente de 1.27 (p < 0.001), indicando um incremento médio de 1.27 pontos na dor em relação à sessão 1.

Interpretação:
Os efeitos fixos mostram um aumento progressivo e estatisticamente significativo da dor ao longo das sessões. Assim, a dor relatada tende a ser maior em sessões posteriores, em comparação com a sessão inicial.

Efeitos randomicos:

Intercepto (ID):
– O desvio padrão é de 2.23, o que indica uma variabilidade substancial entre os indivíduos na dor basal (ou seja, diferenças individuais não explicadas pelos efeitos fixos).

Erro Residual:
– O desvio padrão dos resíduos é de 0.66, demonstrando que a variabilidade não explicada dentro dos indivíduos (ou erro de medida) é relativamente baixa.

Interpretação:
O modelo evidencia que há uma variação considerável nos níveis basais de dor entre os sujeitos, enquanto a variabilidade dos dados dentro de cada sujeito é bem menor. Isso reforça a adequação do modelo em capturar a estrutura hierárquica dos dados.

4.3.4.1 Ajuste de correlação serial

library(nlme)
modelo_cor1 <- nlme::lme(
    # Define os efeitos fixos (a influência das sessões na escala de dor).
    fixed = NPRS ~ factor(session),
    # Permite que cada indivíduo tenha seu intercepto, capturando a variabilidade interindividual.
    random = ~ 1 | ID,
    data = df_long,
    # Especifica uma estrutura de correlação AR(1) para as medições dentro de cada indivíduo,
    # considerando que as observações mais próximas no tempo estão mais correlacionadas.
    correlation = corAR1(form = ~ session | ID)
)

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

# Fixed Effects

Parameter	Coefficient	SE	95% CI	t(195)	p
(Intercept)	2.32	0.29	(1.75, 2.88)	8.11	< .001
session2	0.48	0.10	(0.29, 0.67)	5.04	< .001
session3	0.97	0.13	(0.72, 1.22)	7.61	< .001
session4	1.27	0.15	(0.98, 1.56)	8.68	< .001

# Random Effects

Parameter	Coefficient
SD (Intercept: ID)	2.04
SD (Residual)	1.11

modelo_cor1 |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p < .001).

modelo_cor1 |>
    performance::check_heteroscedasticity()

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

modelo_cor1$residuals |>
    performance::check_normality()

Warning: Non-normality of raw detected (p < .001).

modelo_cor1$residuals |>
    performance::check_normality() |>
    plot()

4.4 Análise considerando a escala de dor dos últimos 15 dias como linha de base

4.4.1 Estatísticas descritivas

df_long02 |>
    group_by(session) |>
    summarise(
        Média = round(mean(NPRS), 2),
        Mediana = round(median(NPRS), 2),
        DesvioPadrão = round(sd(NPRS), 2),
        Mínimo = round(min(NPRS), 2),
        Máximo = round(max(NPRS), 2)
    ) |>
    gt::gt()

session	Média	Mediana	DesvioPadrão	Mínimo	Máximo
0	5.55	6	1.53	3	9
1	2.32	2	2.24	0	8
2	2.80	2	2.32	0	8
3	3.29	3	2.33	0	9
4	3.59	3	2.40	0	10

4.4.2 Distribuição

df_long02 |>
    ggplot() +
    aes(
        x = factor(session, levels = sort(unique(session), decreasing = TRUE)),
        y = NPRS
    ) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2), alpha = .2
    ) +
    geom_point(
        position = position_jitter(w = .15), alpha = .2
    ) +
    geom_boxplot(width = .25, alpha = .2, outlier.shape = NA) +
    coord_flip() +
    theme_minimal() +
    theme(
        legend.position = "none"
    ) +
    labs(x = "Sessão")

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

4.4.3 Testes de normalidade

run_normality_tests <- function(data) {
    shapiro_test <- shapiro.test(data)
    ks_test <- ks.test(data, "pnorm", mean = mean(data), sd = sd(data))
    ad_test <- nortest::ad.test(data)

    return(data.frame(
        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)
        )
    ))
}

normality_results <- df_long02 %>%
    group_by(session) %>%
    summarise(
        test_results = list(run_normality_tests(NPRS)),
        .groups = "drop"
    )

results_table <- data.frame()

for (i in 1:nrow(normality_results)) {
    session_data <- normality_results$test_results[[i]]
    session_data$Teste <- c(
        "Shapiro-Wilk",
        "Kolmogorov-Smirnov",
        "Anderson-Darling"
    )
    session_data$Sessão <- normality_results$session[i]
    results_table <- rbind(results_table, session_data)
}

results_table %>%
    select(Sessão, Teste, Estatística, `Valor-p` = `Valor.p`) %>%
    arrange(Teste) %>%
    group_by(Teste) %>%
    gt::gt() %>%
    gt::tab_header(title = "Testes de Normalidade por Sessão") %>%
    gt::cols_label(
        Teste = "Teste de Normalidade",
        Sessão = "Sessão",
        Estatística = "Estatística",
        `Valor-p` = "Valor-p"
    ) %>%
    gt::tab_spanner(
        label = "Resultados",
        columns = c(Estatística, `Valor-p`)
    )

Testes de Normalidade por Sessão
Sessão	Resultados
	Estatística	Valor-p
Anderson-Darling
0	1.418	0.001
1	2.499	0.000
2	1.696	0.000
3	1.017	0.010
4	0.801	0.036
Kolmogorov-Smirnov
0	0.162	0.062
1	0.163	0.061
2	0.161	0.066
3	0.134	0.184
4	0.127	0.235
Shapiro-Wilk
0	0.943	0.004
1	0.883	0.000
2	0.916	0.000
3	0.948	0.008
4	0.958	0.024

4.4.4 Análise de variância

modelo1.1 <- lme4::lmer(NPRS ~ factor(session) + (1 | ID), data = df_long02)

modelo1.1 |>
    parameters::model_parameters() |>
    parameters::print_md()

# Fixed Effects

Parameter	Coefficient	SE	95% CI	t(323)	p
(Intercept)	5.55	0.27	(5.02, 6.08)	20.60	< .001
session [1]	-3.23	0.19	(-3.61, -2.85)	-16.68	< .001
session [2]	-2.74	0.19	(-3.12, -2.36)	-14.17	< .001
session [3]	-2.26	0.19	(-2.64, -1.88)	-11.67	< .001
session [4]	-1.95	0.19	(-2.34, -1.57)	-10.10	< .001

# Random Effects

Parameter	Coefficient	SE	95% CI
SD (Intercept: ID)	1.88	0.18	(1.57, 2.26)
SD (Residual)	1.11	0.05	(1.02, 1.21)

modelo1.1 |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p < .001).

modelo1.1 |>
    performance::check_heteroscedasticity()

Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.037).

modelo1.1 |>
    performance::check_residuals()

Warning: Non-uniformity of simulated residuals detected (p = 0.018).

modelo1.1 |>
    performance::check_normality() |>
    plot()

Efeitos fixos:

Intercepto (Sessão 0):
– Coeficiente de 5.55 (p < 0.001), representando a média da escala de dor nos últimos 15 dias.

Sessão 1:
– Coeficiente de -3.23 (p < 0.001): indica que, em média, a dor na sessão 1 é 3.23 pontos menor que a linha de base, correspondendo a uma média de 2.32, conforme as estatísticas descritivas.

Sessão 2:
– Coeficiente de -2.74 (p < 0.001): a dor aumenta em relação à sessão 1, sendo 2.74 pontos menor que a linha de base (média de 2.80).

Sessão 3:
– Coeficiente de -2.26 (p < 0.001): mostra um aumento adicional, com média de 3.29 pontos.

Sessão 4:
– Coeficiente de -1.95 (p < 0.001): representa uma média de 3.59 pontos, ou seja, a diferença em relação à linha de base torna-se menos acentuada.

Interpretação:
Os efeitos fixos indicam que, comparada à linha de base (sessão 0 com média 5.55), a dor diminui significativamente nas sessões seguintes. Além disso, os coeficientes menos negativos em sessões posteriores sugerem um aumento gradual dos níveis de dor entre as sessões 1 e 4, embora estes permaneçam abaixo da linha de base.

Efeitos randomicos:

Interceptos entre Indivíduos:
– O desvio padrão dos interceptos é de 1.88 (IC: [1.57, 2.26]), indicando variabilidade significativa entre os sujeitos quanto aos níveis basais de dor.

Erro Residual:
– O desvio padrão residual é de 1.11 (IC: [1.02, 1.21]), refletindo a variabilidade não explicada dentro dos indivíduos.

Interpretação:
Há diferenças notáveis entre os indivíduos (variação dos interceptos) e uma variabilidade interna moderada, o que justifica a inclusão dos efeitos randômicos no modelo.

4.5 Análise considerando a escala de dor inicial como linha de base

4.5.1 Estatísticas descritivas

df_long03 |>
    group_by(session) |>
    summarise(
        Média = round(mean(NPRS), 2),
        Mediana = round(median(NPRS), 2),
        DesvioPadrão = round(sd(NPRS), 2),
        Mínimo = round(min(NPRS), 2),
        Máximo = round(max(NPRS), 2)
    ) |>
    gt::gt()

session	Média	Mediana	DesvioPadrão	Mínimo	Máximo
0	1.53	0.5	1.99	0	7
1	2.32	2.0	2.24	0	8
2	2.80	2.0	2.32	0	8
3	3.29	3.0	2.33	0	9
4	3.59	3.0	2.40	0	10

4.5.2 Distribuição

df_long03 |>
    ggplot() +
    aes(
        x = factor(session, levels = sort(unique(session), decreasing = TRUE)),
        y = NPRS
    ) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2), alpha = .2
    ) +
    geom_point(
        position = position_jitter(w = .15), alpha = .2
    ) +
    geom_boxplot(width = .25, alpha = .2, outlier.shape = NA) +
    coord_flip() +
    theme_minimal() +
    theme(
        legend.position = "none"
    ) +
    labs(x = "Sessão")

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

4.5.3 Testes de normalidade

run_normality_tests <- function(data) {
    shapiro_test <- shapiro.test(data)
    ks_test <- ks.test(data, "pnorm", mean = mean(data), sd = sd(data))
    ad_test <- nortest::ad.test(data)

    return(data.frame(
        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)
        )
    ))
}

normality_results <- df_long03 %>%
    group_by(session) %>%
    summarise(
        test_results = list(run_normality_tests(NPRS)),
        .groups = "drop"
    )

results_table <- data.frame()

for (i in 1:nrow(normality_results)) {
    session_data <- normality_results$test_results[[i]]
    session_data$Teste <- c(
        "Shapiro-Wilk",
        "Kolmogorov-Smirnov",
        "Anderson-Darling"
    )
    session_data$Sessão <- normality_results$session[i]
    results_table <- rbind(results_table, session_data)
}

results_table %>%
    select(Sessão, Teste, Estatística, `Valor-p` = `Valor.p`) %>%
    arrange(Teste) %>%
    group_by(Teste) %>%
    gt::gt() %>%
    gt::tab_header(title = "Testes de Normalidade por Sessão") %>%
    gt::cols_label(
        Teste = "Teste de Normalidade",
        Sessão = "Sessão",
        Estatística = "Estatística",
        `Valor-p` = "Valor-p"
    ) %>%
    gt::tab_spanner(
        label = "Resultados",
        columns = c(Estatística, `Valor-p`)
    )

Testes de Normalidade por Sessão
Sessão	Resultados
	Estatística	Valor-p
Anderson-Darling
0	6.176	0.000
1	2.499	0.000
2	1.696	0.000
3	1.017	0.010
4	0.801	0.036
Kolmogorov-Smirnov
0	0.279	0.000
1	0.163	0.061
2	0.161	0.066
3	0.134	0.184
4	0.127	0.235
Shapiro-Wilk
0	0.771	0.000
1	0.883	0.000
2	0.916	0.000
3	0.948	0.008
4	0.958	0.024

4.5.4 Análise de variância

modelo1.2 <- lme4::lmer(NPRS ~ factor(session) + (1 | ID), data = df_long03)

modelo1.2 |>
    parameters::model_parameters() |>
    parameters::print_md()

# Fixed Effects

Parameter	Coefficient	SE	95% CI	t(323)	p
(Intercept)	1.53	0.28	(0.98, 2.08)	5.50	< .001
session [1]	0.79	0.15	(0.49, 1.09)	5.21	< .001
session [2]	1.27	0.15	(0.98, 1.57)	8.42	< .001
session [3]	1.76	0.15	(1.46, 2.06)	11.62	< .001
session [4]	2.06	0.15	(1.76, 2.36)	13.63	< .001

# Random Effects

Parameter	Coefficient	SE	95% CI
SD (Intercept: ID)	2.09	0.19	(1.75, 2.49)
SD (Residual)	0.87	0.04	(0.80, 0.95)

modelo1.2 |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p < .001).

modelo1.2 |>
    performance::check_heteroscedasticity()

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

modelo1.2 |>
    performance::check_residuals()

Warning: Non-uniformity of simulated residuals detected (p < .001).

modelo1.2 |>
    performance::check_normality() |>
    plot()

Efeitos fixos: Intercepto (Sessão 0):
Coeficiente = 1.53 (95% CI: [0.98, 2.08], p < 0.001) Interpretação: Representa a média da dor na sessão 0 (linha de base).

Sessão 1:
Coeficiente = 0.79 (95% CI: [0.49, 1.09], p < 0.001) Interpretação: Em média, a dor aumenta 0.79 pontos na sessão 1 em relação à linha de base.

Sessão 2:
Coeficiente = 1.27 (95% CI: [0.98, 1.57], p < 0.001) Interpretação: A dor aumenta 1.27 pontos na sessão 2 em comparação à sessão 0.

Sessão 3:
Coeficiente = 1.76 (95% CI: [1.46, 2.06], p < 0.001) Interpretação: A dor aumenta 1.76 pontos na sessão 3 em relação à linha de base.

Sessão 4:
Coeficiente = 2.06 (95% CI: [1.76, 2.36], p < 0.001) Interpretação: A dor aumenta 2.06 pontos na sessão 4 em relação à linha de base.

Resumo dos Efeitos Fixos:
Os coeficientes positivos e crescentes demonstram que, em comparação à linha de base, a dor aumenta de forma progressiva e significativa ao longo das sessões.

Efeitos randomicos:

Intercepto (ID):
Desvio padrão = 2.09 (95% CI: [1.75, 2.49]) Interpretação: Há variabilidade significativa entre os indivíduos em seus níveis basais de dor.

Erro Residual:
Desvio padrão = 0.87 (95% CI: [0.80, 0.95]) Interpretação: A variabilidade não explicada dentro dos indivíduos é relativamente baixa.

Resumo dos Efeitos Randômicos:
A inclusão dos interceptos aleatórios permite captar as diferenças individuais nos níveis iniciais de dor, reforçando que os sujeitos variam significativamente entre si, enquanto a variabilidade intraindividual é menor.

4.5.5 Incorparando a correlação serial

library(nlme)
modelo_cor <- nlme::lme(
    # Define os efeitos fixos (a influência das sessões na escala de dor).
    fixed = NPRS ~ factor(session),
    # Permite que cada indivíduo tenha seu intercepto, capturando a variabilidade interindividual.
    random = ~ 1 | ID,
    data = df_long03,
    # Especifica uma estrutura de correlação AR(1) para as medições dentro de cada indivíduo,
    # considerando que as observações mais próximas no tempo estão mais correlacionadas.
    correlation = corAR1(form = ~ session | ID)
)


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

# Fixed Effects

Parameter	Coefficient	SE	95% CI	t(260)	p
(Intercept)	1.53	0.27	(0.99, 2.07)	5.57	< .001
session1	0.79	0.12	(0.55, 1.03)	6.48	< .001
session2	1.27	0.16	(0.95, 1.59)	7.82	< .001
session3	1.76	0.19	(1.38, 2.13)	9.28	< .001
session4	2.06	0.21	(1.65, 2.47)	9.90	< .001

# Random Effects

Parameter	Coefficient
SD (Intercept: ID)	1.61
SD (Residual)	1.54

modelo_cor |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p = 0.004).

modelo_cor |>
    performance::check_heteroscedasticity()

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

modelo_cor$residuals |>
    performance::check_normality()

Warning: Non-normality of raw detected (p < .001).

modelo_cor$residuals |>
    performance::check_normality() |>
    plot()

4.6 MDP

4.6.1 Estatísticas descritivas

df_long |>
    group_by(session) |>
    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)
    ) |>
    gt::gt()

session	Média	Mediana	DesvioPadrão	Mínimo	Máximo
1	11.89	12.00	2.40	7.39	17.41
2	11.68	11.49	2.38	7.86	17.83
3	11.82	11.29	2.56	8.09	18.67
4	12.01	11.52	2.42	8.09	16.41

4.6.2 Distribuição

df_long |>
    ggplot() +
    aes(
        x = factor(session, levels = sort(unique(session), decreasing = TRUE)),
        y = MDP
    ) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2), alpha = .2
    ) +
    geom_point(
        position = position_jitter(w = .15), alpha = .2
    ) +
    geom_boxplot(width = .25, alpha = .2, outlier.shape = NA) +
    coord_flip() +
    theme_minimal() +
    theme(
        legend.position = "none"
    ) +
    labs(x = "Sessão")

df_long |>
    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() +
    facet_wrap(~session)

4.6.3 Testes de normalidade

run_normality_tests <- function(data) {
    shapiro_test <- shapiro.test(data)
    ks_test <- ks.test(data, "pnorm", mean = mean(data), sd = sd(data))
    ad_test <- nortest::ad.test(data)

    return(data.frame(
        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)
        )
    ))
}

normality_results <- df_long %>%
    group_by(session) %>%
    summarise(
        test_results = list(run_normality_tests(MDP)),
        .groups = "drop"
    )

results_table <- data.frame()

for (i in 1:nrow(normality_results)) {
    session_data <- normality_results$test_results[[i]]
    session_data$Teste <- c(
        "Shapiro-Wilk",
        "Kolmogorov-Smirnov",
        "Anderson-Darling"
    )
    session_data$Sessão <- normality_results$session[i]
    results_table <- rbind(results_table, session_data)
}

results_table %>%
    select(Sessão, Teste, Estatística, `Valor-p` = `Valor.p`) %>%
    arrange(Teste) %>%
    group_by(Teste) %>%
    gt::gt() %>%
    gt::tab_header(title = "Testes de Normalidade por Sessão") %>%
    gt::cols_label(
        Teste = "Teste de Normalidade",
        Sessão = "Sessão",
        Estatística = "Estatística",
        `Valor-p` = "Valor-p"
    ) %>%
    gt::tab_spanner(
        label = "Resultados",
        columns = c(Estatística, `Valor-p`)
    )

Testes de Normalidade por Sessão
Sessão	Resultados
	Estatística	Valor-p
Anderson-Darling
1	0.610	0.109
2	0.519	0.181
3	0.987	0.012
4	0.999	0.012
Kolmogorov-Smirnov
1	0.097	0.532
2	0.083	0.719
3	0.120	0.277
4	0.096	0.545
Shapiro-Wilk
1	0.968	0.084
2	0.967	0.075
3	0.950	0.010
4	0.947	0.007

4.6.4 Análise de variância

modelo2 <- lme4::lmer(MDP ~ factor(session) + (1 | ID), data = df_long)

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

# Fixed Effects

Parameter	Coefficient	SE	95% CI	t(258)	p
(Intercept)	11.89	0.30	(11.30, 12.48)	39.60	< .001
session [2]	-0.21	0.17	(-0.54, 0.12)	-1.26	0.210
session [3]	-0.07	0.17	(-0.40, 0.25)	-0.45	0.656
session [4]	0.12	0.17	(-0.21, 0.45)	0.72	0.475

# Random Effects

Parameter	Coefficient	SE	95% CI
SD (Intercept: ID)	2.24	0.21	(1.88, 2.69)
SD (Residual)	0.95	0.05	(0.86, 1.05)

modelo2 |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p < .001).

modelo2 |>
    performance::check_heteroscedasticity()

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

modelo2 |>
    performance::check_residuals()

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

modelo2 |>
    performance::check_normality()

OK: residuals appear as normally distributed (p = 0.055).

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

Efeitos Fixos:

Intercepto: 11.89 (p < 0.001), indicando a média de MDP na sessão de referência (sessão 1).

Sessões:

Sessão 2: coeficiente de -0.21 (p = 0.210)

Sessão 3: coeficiente de -0.07 (p = 0.656)

Sessão 4: coeficiente de 0.12 (p = 0.475) → Nenhuma das comparações entre sessões apresenta efeito estatisticamente significativo. Isso sugere que a qualidade do movimento, medida pelo MDP, não difere significativamente entre as sessões.

Efeitos Randômicos:

O desvio padrão dos interceptos por ID é de 2.24, o que indica uma variabilidade considerável entre os indivíduos em seus níveis basais de MDP.

A variabilidade residual é de 0.95, mostrando que a variabilidade não explicada dentro dos indivíduos é relativamente pequena.

4.6.5 Incorparando a correlação serial

modelo_cor02 <- nlme::lme(
    # Define os efeitos fixos (a influência das sessões na escala de dor).
    fixed = MDP ~ factor(session),
    # Permite que cada indivíduo tenha seu intercepto, capturando a variabilidade interindividual.
    random = ~ 1 | ID,
    data = df_long,
    # Especifica uma estrutura de correlação AR(1) para as medições dentro de cada indivíduo,
    # considerando que as observações mais próximas no tempo estão mais correlacionadas.
    correlation = corAR1(form = ~ session | ID)
)

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

# Fixed Effects

Parameter	Coefficient	SE	95% CI	t(195)	p
(Intercept)	11.89	0.30	(11.30, 12.48)	39.62	< .001
session2	-0.21	0.16	(-0.53, 0.11)	-1.29	0.198
session3	-0.07	0.17	(-0.41, 0.26)	-0.44	0.664
session4	0.12	0.17	(-0.22, 0.46)	0.69	0.488

# Random Effects

Parameter	Coefficient
SD (Intercept: ID)	2.23
SD (Residual)	0.98

modelo_cor02 |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p < .001).

modelo_cor02 |>
    performance::check_heteroscedasticity()

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

modelo_cor02$residuals |>
    performance::check_normality()

Warning: Non-normality of raw detected (p < .001).

modelo_cor02$residuals |>
    performance::check_normality() |>
    plot()

5 MDP e NPRS

5.1 Dor durante o movimento

modelo4 <- lme4::lmer(MDP ~ NPRS + factor(session) + (1 | ID), data = df_long)

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

# Fixed Effects

Parameter	Coefficient	SE	95% CI	t(257)	p
(Intercept)	11.58	0.35	(10.90, 12.27)	33.12	< .001
NPRS	0.13	0.08	(-0.02, 0.29)	1.66	0.098
session [2]	-0.27	0.17	(-0.61, 0.06)	-1.60	0.111
session [3]	-0.20	0.18	(-0.56, 0.16)	-1.10	0.271
session [4]	-0.05	0.19	(-0.43, 0.33)	-0.25	0.800

# Random Effects

Parameter	Coefficient	SE	95% CI
SD (Intercept: ID)	2.22	0.20	(1.85, 2.66)
SD (Residual)	0.95	0.05	(0.86, 1.05)

modelo4 |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p < .001).

modelo4 |>
    performance::check_heteroscedasticity()

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

modelo4 |>
    performance::check_residuals()

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

modelo4 |>
    performance::check_normality()

OK: residuals appear as normally distributed (p = 0.056).

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

Efeitos Fixos:

Intercepto: O valor de 11.58 indica a média de MDP na sessão de referência (sessão 1) para um valor de NPRS igual a zero.

NPRS: O coeficiente de 0.13 (p = 0.098) indica que, para cada ponto adicional em NPRS, há um aumento estimado de 0.13 pontos em MDP. Contudo, esse efeito é marginal e não atinge significância estatística ao nível de 5%.

Sessões: Os efeitos das sessões (sessão 2: -0.27, sessão 3: -0.20, sessão 4: -0.05) não são estatisticamente significativos, sugerindo que, após ajustar pela dor, a qualidade do movimento não difere significativamente entre as sessões.

Efeitos Randômicos:

Há uma variabilidade considerável entre os indivíduos (desvio padrão dos interceptos ≈ 2.22), enquanto a variabilidade residual é menor (≈ 0.95). Isso indica que as diferenças basais entre os indivíduos são importantes para a variação observada no MDP.

5.1.1 Incorparando a correlação serial

modelo_cor02 <- nlme::lme(
    # Define os efeitos fixos (a influência das sessões na escala de dor).
    fixed = MDP ~ NPRS + factor(session),
    # Permite que cada indivíduo tenha seu intercepto, capturando a variabilidade interindividual.
    random = ~ 1 | ID,
    data = df_long,
    # Especifica uma estrutura de correlação AR(1) para as medições dentro de cada indivíduo,
    # considerando que as observações mais próximas no tempo estão mais correlacionadas.
    correlation = corAR1(form = ~ session | ID)
)

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

# Fixed Effects

Parameter	Coefficient	SE	95% CI	t(194)	p
(Intercept)	11.59	0.35	(10.90, 12.28)	33.05	< .001
NPRS	0.13	0.08	(-0.03, 0.29)	1.63	0.106
session2	-0.27	0.17	(-0.60, 0.06)	-1.64	0.103
session3	-0.20	0.19	(-0.57, 0.17)	-1.07	0.284
session4	-0.05	0.20	(-0.44, 0.34)	-0.24	0.811

# Random Effects

Parameter	Coefficient
SD (Intercept: ID)	2.20
SD (Residual)	0.98

modelo_cor02 |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p < .001).

modelo_cor02 |>
    performance::check_heteroscedasticity()

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

modelo_cor02$residuals |>
    performance::check_normality()

Warning: Non-normality of raw detected (p < .001).

modelo_cor02$residuals |>
    performance::check_normality() |>
    plot()

5.2 Dor e movimento interagindo em funcao da sessão

modelo5 <- lme4::lmer(MDP ~ NPRS * factor(session) + (1 | ID), data = df_long)

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

# Fixed Effects

Parameter	Coefficient	SE	95% CI	t(254)	p
(Intercept)	11.60	0.37	(10.87, 12.33)	31.45	< .001
NPRS	0.13	0.09	(-0.06, 0.31)	1.34	0.181
session [2]	-0.13	0.26	(-0.63, 0.38)	-0.50	0.619
session [3]	-0.22	0.28	(-0.77, 0.33)	-0.79	0.429
session [4]	-0.20	0.29	(-0.76, 0.37)	-0.68	0.494
NPRS × session [2]	-0.05	0.07	(-0.20, 0.10)	-0.69	0.493
NPRS × session [3]	7.23e-03	0.07	(-0.14, 0.15)	0.10	0.923
NPRS × session [4]	0.04	0.07	(-0.10, 0.19)	0.58	0.559

# Random Effects

Parameter	Coefficient	SE	95% CI
SD (Intercept: ID)	2.22	0.20	(1.85, 2.66)
SD (Residual)	0.96	0.05	(0.87, 1.06)

modelo5 |>
    performance::check_autocorrelation()

Warning: Autocorrelated residuals detected (p < .001).

modelo5 |>
    performance::check_heteroscedasticity()

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

modelo5 |>
    performance::check_normality()

OK: residuals appear as normally distributed (p = 0.099).

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

Efeitos Fixos

Intercepto (11.60):
Representa a média de MDP na sessão de referência (sessão 1) quando NPRS é zero.

NPRS (coeficiente = 0.13, p = 0.181): A associação principal entre NPRS e MDP na sessão de referência não é estatisticamente significativa. Isso indica que, para cada ponto adicional em NPRS, o aumento estimado em MDP é de 0.13 pontos, mas esse efeito não alcança significância.

Sessões (comparadas à sessão 1):

Sessão 2: coeficiente = -0.13 (p = 0.619)

Sessão 3: coeficiente = -0.22 (p = 0.429)

Sessão 4: coeficiente = -0.20 (p = 0.494)

Esses efeitos não são significativos, sugerindo que, após ajustar pela dor, as médias de MDP entre as sessões não diferem de forma relevante.

Interações NPRS × Sessão:

NPRS × Sessão 2: coeficiente = -0.05 (p = 0.493)

NPRS × Sessão 3: coeficiente = 0.007 (p = 0.923)

NPRS × Sessão 4: coeficiente = 0.04 (p = 0.559)

Nenhuma das interações é significativa. Isso indica que o efeito da dor sobre a qualidade do movimento não varia significativamente entre as sessões.

Efeitos Randômicos

Interceptos por ID:
O desvio padrão de 2.22 indica uma variabilidade considerável entre os indivíduos nos níveis basais de MDP.

Residual:
O desvio padrão residual de 0.96 reflete a variabilidade não explicada dentro dos indivíduos.

Diagnóstico do Modelo

Heterocedasticidade e Normalidade dos Resíduos:
Os testes indicam que os resíduos são homocedásticos (p = 0.378) e aproximadamente normais (p = 0.099).

Conclusão

Associação Geral:
Não há evidência estatística de que a dor (NPRS) tenha um efeito significativo sobre a qualidade do movimento (MDP), nem que esse efeito varie entre as sessões.

Interação:
Os termos de interação não são significativos, o que sugere que a relação entre NPRS e MDP é estável ao longo do tempo.

Em resumo, esse modelo sugere que, ao ajustar por sessão, não encontramos evidência de que a dor durante o movimento influencie a qualidade do movimento de forma diferente ao longo do tempo.

6 Análise complementar - Correlações

Objetivo: Verificar se o MDP está mais associado com: EVA 15 dias; EVA inical; ou EVA durante o movimento.

6.1 existe diferença do MDP entre as sessões?

Modelo

# Modelo
mod_misto <- lmer(MDP ~ factor(session) + (1|ID), data = df_long)
mod_misto |> parameters()

ANOVA

# Análise de variância e comparações múltiplas
anova(mod_misto)

Qualidade do ajuste

# Qualidade do ajuste
r2_vals  <- performance::r2(mod_misto)   # R2 marginal/condicional
icc_vals <- performance::icc(mod_misto)  # ICC

print(r2_vals)

# R2 for Mixed Models

  Conditional R2: 0.847
     Marginal R2: 0.002

print(icc_vals)

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.847
  Unadjusted ICC: 0.845

Comparações múltiplas

# Contrastes pareados com ajuste de Tukey
emm_sess <- emmeans(mod_misto, ~ session)
contr_tukey <- contrast(emm_sess, method = "pairwise", adjust = "tukey")
contr_tukey |> parameters()

Gráfico das EMMs

# Médias marginais por sessão (EMMs) + IC95%
emm_sess <- emmeans(mod_misto, ~ session)         # por padrão, IC95%
emm_df   <- as.data.frame(emm_sess)

# Gráfico das EMMs com IC95%
p_emm <- ggplot(emm_df, aes(x = session, y = emmean)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.12) +
  geom_line(aes(group = 1), linewidth = 0.6) +
  labs(
    title = "MDP per session (EMMs with 95% CI)",
    x = "Session",
    y = "MDP (estimated marginal mean)"
  ) +
  theme_minimal(base_size = 12)

print(p_emm)

6.1.1 Interpretação dos resultados

Não há evidência de diferença de MDP entre as sessões. O efeito global de sessão não foi significativo (ANOVA tipo III, Kenward–Roger: F(3,195)=1,37; p=0,255). Todas as comparações pareadas ajustadas por Tukey foram não significativas.

Principais resultados
- Modelo: MDP ~ factor(session) + (1|ID) (REML).
- Efeitos fixos (vs. Sessão 1):
S2: −0,21 (IC95% −0,54 a 0,12; p=0,210) · S3: −0,07 (−0,40 a 0,25; p=0,656) · S4: +0,12 (−0,21 a 0,45; p=0,475).
→ Todas as CIs incluem 0; efeitos pequenos e não significativos.

• ANOVA (global): F(3,195)=1,365; p=0,255 → sem efeito de sessão.

• Comparações pareadas (Tukey):
  Maior contraste em magnitude: S2 vs S4 = −0,33 (IC95% −0,66 a 0,00; p ajust.=0,362) → não significativo após ajuste.
  Todos os demais contrastes: p ajustados ≥ 0,807.

• Variância / consistência intraindivíduo:
  SD_intercepto(ID)=2,24 (IC95% 1,88–2,69); SD_residual=0,95 (0,86–1,05).
  ICC ajustado = 0,847 (não ajustado 0,845) → a maior parte da variância é entre indivíduos, não entre sessões.

• Qualidade do ajuste (R² misto):
  R² marginal = 0,002 (efeitos fixos ≈ 0% da variância) · R² condicional = 0,847 (modelo completo ≈ 85%).
  → O fator sessão explica quase nada da variabilidade total do MDP.

Interpretação prática
- O MDP permaneceu estável entre as quatro sessões; diferenças médias são pequenas (≈ 0,1–0,3 ponto) e clinicamente modestas.
- O alto ICC indica que as diferenças entre pessoas são muito maiores do que as variações do MDP ao longo das sessões no mesmo indivíduo.

Notas de reporte
- Texto sugerido:

“Em modelo linear misto com intercepto aleatório por participante, o efeito de sessão sobre MDP não foi significativo (F(3,195)=1,37; p=0,255). As comparações pareadas de Tukey não mostraram diferenças entre sessões (todos p ajustados > 0,36). R² marginal=0,002; R² condicional=0,847; ICC=0,847, evidenciando que a variabilidade é predominantemente entre indivíduos.”

6.2 Correlações entre MDP e EVA

tomando como a primeira sessão como referência para o MDP, existe correlação entre o MDP e a EVA em 3 situações diferentes: 15 dias, inicial e durante o movimento.

6.2.1 EVA 15 dias

# tomando como a primeira sessão como referência para o MDP, existe correlação entre o MDP e a EVA em 3 situações diferentes: 15 dias, inicial e durante o movimento.


# EVA 15 dias
cor.test(df$EVA_15dias, df$SDL_01_MDP, method = "spearman")

    Spearman's rank correlation rho

data:  df$EVA_15dias and df$SDL_01_MDP
S = 42230, p-value = 0.3435
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.1184592 

Resposta: Não existe correlação significativa entre o MDP e a EVA 15 dias.

6.2.2 EVA inicial

# EVA inicial
cor.test(df$NPRS_Inicial, df$SDL_01_MDP, method = "spearman")

    Spearman's rank correlation rho

data:  df$NPRS_Inicial and df$SDL_01_MDP
S = 40761, p-value = 0.2321
alternative hypothesis: true rho is not equal to 0
sample estimates:
     rho 
0.149122 

Resposta: Não existe correlação significativa entre o MDP e a EVA inicial.

6.2.3 EVA durante o movimento (fixando a sessão 1)

# EVA durante o movimento (Referência - sessão 1)
cor.test(df$NPRS_SDL01, df$SDL_01_MDP, method = "spearman")

    Spearman's rank correlation rho

data:  df$NPRS_SDL01 and df$SDL_01_MDP
S = 35392, p-value = 0.03414
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.2612016 

Resposta: Existe correlação significativa entre o MDP e a EVA durante o movimento.

6.3 Existe diferença da EVA entre as sessões?**

Modelo

# Modelo
mod_misto <- lmer(NPRS ~ factor(session) + (1|ID), data = df_long)
mod_misto |> parameters()

ANOVA

# Análise de variância e comparações múltiplas
anova(mod_misto)

Qualidade do ajuste

# Qualidade do ajuste
r2_vals  <- performance::r2(mod_misto)   # R2 marginal/condicional
icc_vals <- performance::icc(mod_misto)  # ICC

print(r2_vals)

# R2 for Mixed Models

  Conditional R2: 0.922
     Marginal R2: 0.042

print(icc_vals)

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.919
  Unadjusted ICC: 0.881

Comparações múltiplas

# Contrastes pareados com ajuste de Tukey
emm_sess <- emmeans(mod_misto, ~ session)
contr_tukey <- contrast(emm_sess, method = "pairwise", adjust = "tukey")
contr_tukey |> parameters()

Gráfico das EMMs

# Médias marginais por sessão (EMMs) + IC95%
emm_sess <- emmeans(mod_misto, ~ session)         # por padrão, IC95%
emm_df   <- as.data.frame(emm_sess)

# Gráfico das EMMs com IC95%
p_emm <- ggplot(emm_df, aes(x = session, y = emmean)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.12) +
  geom_line(aes(group = 1), linewidth = 0.6) +
  labs(
    title = "NPRS per session (EMMs with 95% CI)",
    x = "Session",
    y = "NPRS (estimated marginal mean)"
  ) +
  theme_minimal(base_size = 12)

print(p_emm)

6.3.1 Interpretação dos resultados

Obs. Esse resultados são os mesmos do tópico 4.3

Conclusão

Há diferença de EVA entre as sessões. O efeito global de sessão é forte e estatisticamente significativo (ANOVA tipo III, Kenward–Roger: F(3,195)=47,09; p<2×10⁻¹⁶).

Principais resultados
- Modelo: NPRS ~ factor(session) + (1|ID) (REML).

• Efeitos fixos (referência = Sessão 1):

  • S2: +0,48 (EP=0,12; IC95% 0,26–0,71; p<0,001)
    
  • S3: +0,97 (EP=0,12; IC95% 0,74–1,20; p<0,001)
    
  • S4: +1,27 (EP=0,12; IC95% 1,05–1,50; p<0,001)
    ➜ A dor aumenta de forma monotônica de S1→S4.

• Médias estimadas (EMMs, aprox.): S1=2,32 → S2=2,80 → S3=3,29 → S4=3,59.

• Comparações pareadas (Tukey):

  • Todas as comparações são significativas (p<0,001) exceto S3 vs S4, que é não significativa após ajuste (dif.=0,30; p=0,094).
    ➜ O maior salto ocorre entre S1 e S3/S4; de S3 para S4 há tendência de aumento, mas sem significância após correção.

• Componentes de variância: SD_intercepto(ID)=2,23 (IC95% 1,87–2,65); SD_residual=0,66 (0,60–0,73).

• Qualidade do ajuste / consistência intraindivíduo:

  • R² marginal = 0,042 (sessão explica ~4% da variância total)
    
  • R² condicional = 0,922 (modelo completo ~92%)
    
  • ICC ajustado = 0,919 (não ajustado 0,881) → grande parcela da variabilidade é entre indivíduos.

Interpretação prática
- A EVA aumenta, em média, ~0,5 ponto de S1→S2 e chega a ~+1,27 pontos em S4 vs S1.

• O padrão é consistente e clinicamente interpretável como piora progressiva da dor ao longo das sessões; contudo, a variação total permanece majoritariamente entre pessoas (ICC alto), não dentro do indivíduo entre sessões.

Sugestão de texto para o relatório

“Em modelo linear misto com intercepto aleatório por participante, observou-se efeito significativo de sessão sobre a EVA (F(3,195)=47,09; p<0,001). As médias estimadas aumentaram de 2,32 (S1) para 3,59 (S4). Comparações de Tukey mostraram aumentos significativos entre S1–S2, S1–S3, S1–S4 e S2–S3/S4 (p<0,001); S3–S4 não foi significativo (p=0,094). O R² marginal foi 0,042 e o condicional 0,922; o ICC ajustado foi 0,919, indicando que a variabilidade é predominantemente entre indivíduos.”

Referências

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

Bates, Douglas, Martin Mächler, Ben Bolker, e Steve Walker. 2015. “Fitting linear mixed-effects models using lme4”. Journal of statistical software 67: 1–48.

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.

Kuznetsova, Alexandra, Per B Brockhoff, e Rune HB Christensen. 2017. “lmerTest package: tests in linear mixed effects models”. Journal of statistical software 82: 1–26.

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.