Artigo 2

Autor

Caio Vallio

Data de Publicação

6 de abril de 2025

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

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

2 Importação dos dados

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

3 Ajuste do dataset

Código
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

Código
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
Código
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()
    )

Código
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.

Código
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

Código
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
Código
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()
    )

Código
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
Código
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

Código
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

Código
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")

Código
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

Código
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

Código
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)
Código
modelo1 |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p < .001).
Código
modelo1 |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.900).
Código
modelo1 |>
    performance::check_residuals()
OK: Simulated residuals appear as uniformly distributed (p = 0.074).
Código
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
Código
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
Código
modelo_cor1 |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p < .001).
Código
modelo_cor1 |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.927).
Código
modelo_cor1$residuals |>
    performance::check_normality()
Warning: Non-normality of raw detected (p < .001).
Código
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

Código
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

Código
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")

Código
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

Código
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

Código
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)
Código
modelo1.1 |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p < .001).
Código
modelo1.1 |>
    performance::check_heteroscedasticity()
Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.037).
Código
modelo1.1 |>
    performance::check_residuals()
Warning: Non-uniformity of simulated residuals detected (p = 0.018).
Código
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

Código
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

Código
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")

Código
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

Código
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

Código
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)
Código
modelo1.2 |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p < .001).
Código
modelo1.2 |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.862).
Código
modelo1.2 |>
    performance::check_residuals()
Warning: Non-uniformity of simulated residuals detected (p < .001).
Código
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

Código
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
Código
modelo_cor |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p < .001).
Código
modelo_cor |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.592).
Código
modelo_cor$residuals |>
    performance::check_normality()
Warning: Non-normality of raw detected (p < .001).
Código
modelo_cor$residuals |>
    performance::check_normality() |>
    plot()

4.6 MDP

4.6.1 Estatísticas descritivas

Código
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

Código
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")

Código
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

Código
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

Código
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)
Código
modelo2 |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p < .001).
Código
modelo2 |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.363).
Código
modelo2 |>
    performance::check_residuals()
OK: Simulated residuals appear as uniformly distributed (p = 0.062).
Código
modelo2 |>
    performance::check_normality()
OK: residuals appear as normally distributed (p = 0.055).
Código
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

Código
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
Código
modelo_cor02 |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p < .001).
Código
modelo_cor02 |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.502).
Código
modelo_cor02$residuals |>
    performance::check_normality()
Warning: Non-normality of raw detected (p < .001).
Código
modelo_cor02$residuals |>
    performance::check_normality() |>
    plot()

5 MDP e NPRS

5.1 Dor antes do movimento

Código
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)
Código
modelo4 |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p < .001).
Código
modelo4 |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.369).
Código
modelo4 |>
    performance::check_residuals()
OK: Simulated residuals appear as uniformly distributed (p = 0.067).
Código
modelo4 |>
    performance::check_normality()
OK: residuals appear as normally distributed (p = 0.056).
Código
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

Código
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
Código
modelo_cor02 |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p < .001).
Código
modelo_cor02 |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.496).
Código
modelo_cor02$residuals |>
    performance::check_normality()
Warning: Non-normality of raw detected (p < .001).
Código
modelo_cor02$residuals |>
    performance::check_normality() |>
    plot()

5.2 Dor e movimento interagindo em funcao da sessão

Código
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)
Código
modelo5 |>
    performance::check_autocorrelation()
Warning: Autocorrelated residuals detected (p < .001).
Código
modelo5 |>
    performance::check_heteroscedasticity()
OK: Error variance appears to be homoscedastic (p = 0.378).
Código
modelo5 |>
    performance::check_normality()
OK: residuals appear as normally distributed (p = 0.099).
Código
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 antes do movimento influencie a qualidade do movimento de forma diferente ao longo do tempo.

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.