Smartphone based intervention for depressive symptoms: a randomized non-inferiority clinical trial

Código para limpeza e análise de dados

Autor
Afiliação

Bruno Braga Montezano

Data de Publicação

15 de dezembro de 2023

Carregamento dos dados

No total, foram carregados 17 arquivos para a sessão. São eles:

  1. Onda 1 (coletado pelo SurveyMonkey): w1.xlsx
  2. Onda 2 (coletado pelo SurveyMonkey): w2.xlsx
  3. Onda 3 (coletado pelo SurveyMonkey): w3.xlsx
  4. Onda 4 (coletado pelo SurveyMonkey): w4.xlsx
  5. Onda de follow-up (coletado pelo SurveyMonkey): followup.xlsx
  6. Onda 1 da lista de espera (coletado pelo SurveyMonkey): w1_wl.xlsx
  7. Onda final da lista de espera (coletado pelo SurveyMonkey): final_wl.xlsx
  8. Dados de triagem (coletado pelo SurveyMonkey): screening.xlsx
  9. Dados sociodemográficos (coletado pelo SurveyMonkey): thrive_so.xlsx
  10. Dados de triagem — versão 1 (coletado pelo ODK): odk_v1.xlsx
  11. Dados de triagem — versão 2 (coletado pelo ODK): odk_v2.xlsx
  12. Dados de triagem — versão 3 (coletado pelo ODK): odk_v3.xlsx
  13. Dados de triagem — versão 4 (coletado pelo ODK): odk_v4.xlsx
  14. Dados do aplicativo Thrive (formato largo): thrive_wide_data.csv
  15. Dados do aplicativo Thrive (formato longo): thrive_long_data.csv
  16. Dados do projeto COVIDPSY: covidpsy_data.sav
  17. Tabela com grupos de tratamento: grupos_thrive.xlsx

Os dados em formato xlsx foram carregados através da função read_xlsx do pacote readxl. Arquivos separados por vírgula (csv) foram carregados com a função read_csv do pacote readr e o arquivo em savreferente ao projeto COVIDPSY foi carregado por meio da função read_sav do pacote haven.

Código
# Ler dados do SurveyMonkey
w1 <- readxl::read_xlsx("./data/sm_clean_data/w1.xlsx")
w2 <- readxl::read_xlsx("./data/sm_clean_data/w2.xlsx")
w3 <- readxl::read_xlsx("./data/sm_clean_data/w3.xlsx")
w4 <- readxl::read_xlsx("./data/sm_clean_data/w4.xlsx")
fu <- readxl::read_xlsx("./data/sm_clean_data/followup.xlsx")
w1_wl <- readxl::read_xlsx("./data/sm_clean_data/w1_wl.xlsx")
final_wl <- readxl::read_xlsx("./data/sm_clean_data/final_wl.xlsx")
screening <- readxl::read_xlsx("./data/sm_clean_data/screening.xlsx")
socio <- readxl::read_xlsx("./data/sm_clean_data/thrive_so.xlsx")

# Ler dados do ODK
odk_1 <- readxl::read_xlsx("./data/odk_data/odk_v1.xlsx")
odk_2 <- readxl::read_xlsx("./data/odk_data/odk_v2.xlsx")
odk_3 <- readxl::read_xlsx("./data/odk_data/odk_v3.xlsx")
odk_4 <- readxl::read_xlsx("./data/odk_data/odk_v4.xlsx")

# Ler dados do aplicativo Thrive
thrive_wide <- readr::read_csv("./data/thrive_data/thrive_wide_data.csv")
thrive_long <- readr::read_csv("./data/thrive_data/thrive_long_data.csv")

# Ler dados do estudo COVIDPSY
covidpsy <- haven::read_sav("./data/covidpsy_data/covidpsy_data.sav")

# Ler tabela de grupos validados e faltas na TCCG
group_table <- readxl::read_xlsx("./data/misc_data/grupos_thrive.xlsx")
faltas_table <- readxl::read_xlsx("./data/misc_data/faltas_tccg.xlsx")

União dos dados

Logo na sequência, os dados das quatro ondas de avaliação foram unidos através da função bind_rows do pacote dplyr. Para identificação das waves, criou-se uma nova coluna chamada time (com valores de 1 a 4).

Código
# Juntar tabelas das ondas 1, 2, 3 e 4
all_waves <- dplyr::bind_rows(w1, w2, w3, w4, .id = "time")

Limpeza de emails e duplicatas

O próximo passo foi realizar uma limpeza nas observações. A variável de email foi processada através da aplicação da função str_to_lower do pacote stringr para transformar os emails para letra minúscula e a função str_squish foi usada para remover espaços no início e final dos emails, assim como remover espaços duplicados. Removemos registros de teste do banco de dados também.

Após, foi realizado um processo para solucionar o problema de múltiplos registros de resposta de um sujeito em uma mesma wave (onda de avaliação). Para tal, a base foi agrupada por tempo (time) e email (o identificador único), e foram mantidos apenas o último registro para cada combinação de email e onda de avaliação. Para isso, usou-se funções do pacote dplyr.

Código
# Manter apenas as observações não-repetidas
all_waves <- all_waves |>
    # Limpar variável email (letra minúscula e remover espaços)
    dplyr::mutate(email = stringr::str_squish(stringr::str_to_lower(email))) |>
    # Remover entradas de teste
    dplyr::filter(email != "antonelli.thyago@gmail.com") |>
    # Agrupar por tempo e email
    dplyr::group_by(time, email) |>
    # Manter apenas a última observação deste agrupamento
    dplyr::slice_tail(n = 1) |>
    # Desagrupar os dados
    dplyr::ungroup()

Cálculo dos escores

Então, para o cálculo dos escores acontecer corretamente, as variáveis dos itens das escalas tiveram de ser recodificadas. Elas estavam apresentando uma amplitude de 1 a 4, porém para os escores de ambos os instrumentos as respostas devem gradar entre 0 e 3. Para isso, substraímos 1 de cada resposta para que os itens ficassem adequadamente codificados.

Após, os escores da PHQ-9 e GAD-7 foram calculados para todos os sujeitos.

Código
all_waves <- all_waves |>
    # Recodificar a PHQ e GAD pois estavam com range errado (1-4)
    dplyr::mutate(
        dplyr::across(dplyr::matches("^(phq|gad).*$"), \(x) x - 1)
    ) |>
    dplyr::rowwise() |>
    # Calcular escores totais da PHQ-9 e GAD-7
    dplyr::mutate(
        phq_total = phq_01 + phq_02 + phq_03 + phq_04 + phq_05 + phq_06 +
            phq_07 + phq_08 + phq_09,
        gad_total = gad_01 + gad_02 + gad_03 + gad_04 + gad_05 + gad_06 +
            gad_07
    ) |>
    # Desagrupar os dados
    dplyr::ungroup()

Inclusão dos grupos de tratamento

Na sequência das etapas anteriores, usou-se o conjunto de dados carregado como group_table para adicionar através de um dplyr::left_join, uma variável com os grupos de tratamento do estudo (Group e-CBT ou Thrive) associado a cada um dos emails.

Após, removemos todas observações com valores ausentes na nova variável de grupo, chamada group.

Código
all_waves <- all_waves |>
    # Selecionar apenas email, wave, e escores da PHQ e GAD
    dplyr::select(email, time, phq_total, gad_total) |>
    # Inserir coluna com grupos de tratamento (e-CBT ou Thrive)
    dplyr::left_join(group_table, by = dplyr::join_by(email == email)) |>
    # Remover entradas com valor ausente no grupo
    dplyr::filter(!is.na(group))

Filtro de dados per-protocol

Para captar os dados per-protocol, mantivemos apenas as observações que tiveram registros nas quatro ondas de avaliação (ondas 1, 2, 3 e 4). Criamos um vetor de emails (IDs) com estes participantes para filtrar os dados nos chunks seguintes.

Código
email_pp <- all_waves |>
    # Contar quantas entradas tem por email
    dplyr::count(email, sort = TRUE) |>
    # Selecionar apenas as entradas com quatro registros
    dplyr::filter(n == 4) |>
    # Pegar um vetor dos emails
    dplyr::pull(email)

Visualização da trajetória de cada grupo

Na sequência, criamos visualizações para as trajetórias de cada grupo: TCCG e Thrive. Quatro gráficos foram criados: PHQ-9 analisada por intention-to-treat, PHQ-9 per protocol, GAD-7 por intention-to-treat, e GAD-7 per protocol.

No chunk abaixo, o gráfico da PHQ-9 intention-to-treat é criado.

Código
# Calcular média e desvio padrão da PHQ por grupo e tempo (intention-to-treat)
plot_phq_itt <- rcompanion::groupwiseMean(phq_total ~ time + group,
                          data = all_waves,
                          na.rm = TRUE) |>
    dplyr::mutate(group = dplyr::if_else(group == "tccg", "Group e-CBT", "Thrive App")) |>
    ggplot2::ggplot(ggplot2::aes(x = time, y = Mean, group = group,
                                 color = group, label = paste0(
                                                            "n = ",
                                                            as.character(n),
                                                            "\n",
                                                            as.character(Mean),
                                                            " (",
                                                            as.character(Trad.lower),
                                                            ", ",
                                                            as.character(Trad.upper),
                                                            ")"
                                                            )
                                 )
    ) +
    ggplot2::geom_line(linewidth = 0.5,
                       position = ggplot2::position_dodge(0.1)) +
    ggplot2::geom_errorbar(ggplot2::aes(ymin = Trad.lower, ymax = Trad.upper),
                           linewidth = 0.5,
                           position = ggplot2::position_dodge(0.1)) +
    ggrepel::geom_label_repel(size = 1.5, family = "IBM Plex Sans",
                        show.legend = FALSE,
                        position = ggplot2::position_dodge(0.1)) +
    ggsci::scale_color_jama() +
    ggplot2::theme_classic(10, "IBM Plex Sans") +
    ggplot2::labs(x = "Assessment", y = "PHQ-9 score", color = "Treatment") +
    ggplot2::theme(legend.position = "top")

Logo, o PHQ-9 per-protocol também.

Código
# Calcular média e desvio padrão da PHQ por grupo e tempo (per-protocol)
plot_phq_pp <- all_waves |>
    dplyr::filter(email %in% email_pp) |>
    rcompanion::groupwiseMean(phq_total ~ time + group,
                              data = _,
                              na.rm = TRUE) |>
    dplyr::mutate(group = dplyr::if_else(group == "tccg", "Group e-CBT", "Thrive App")) |>
    ggplot2::ggplot(ggplot2::aes(x = time, y = Mean, group = group,
                                 color = group, label = paste0(
                                                            "n = ",
                                                            as.character(n),
                                                            "\n",
                                                            as.character(Mean),
                                                            " (",
                                                            as.character(Trad.lower),
                                                            ", ",
                                                            as.character(Trad.upper),
                                                            ")"
                                                            )
                                 )
    ) +
    ggplot2::geom_line(linewidth = 0.5,
                       position = ggplot2::position_dodge(0.1)) +
    ggplot2::geom_errorbar(ggplot2::aes(ymin = Trad.lower, ymax = Trad.upper),
                           linewidth = 0.5,
                           position = ggplot2::position_dodge(0.1)) +
    ggrepel::geom_label_repel(size = 1.5, family = "IBM Plex Sans",
                        show.legend = FALSE,
                        position = ggplot2::position_dodge(0.1)) +
    ggsci::scale_color_jama() +
    ggplot2::theme_classic(10, "IBM Plex Sans") +
    ggplot2::labs(x = "Assessment", y = "PHQ-9 score", color = "Treatment") +
    ggplot2::theme(legend.position = "top")

Assim como o gráfico da GAD-7 intention-to-treat.

Código
# Calcular média e desvio padrão da PHQ por grupo e tempo (intention-to-treat)
plot_gad_itt <- rcompanion::groupwiseMean(gad_total ~ time + group,
                          data = all_waves,
                          na.rm = TRUE) |>
    dplyr::mutate(group = dplyr::if_else(group == "tccg", "Group e-CBT", "Thrive App")) |>
    ggplot2::ggplot(ggplot2::aes(x = time, y = Mean, group = group,
                                 color = group, label = paste0(
                                                            "n = ",
                                                            as.character(n),
                                                            "\n",
                                                            as.character(Mean),
                                                            " (",
                                                            as.character(Trad.lower),
                                                            ", ",
                                                            as.character(Trad.upper),
                                                            ")"
                                                            )
                                 )
    ) +
    ggplot2::geom_line(linewidth = 0.5,
                       position = ggplot2::position_dodge(0.1)) +
    ggplot2::geom_errorbar(ggplot2::aes(ymin = Trad.lower, ymax = Trad.upper),
                           linewidth = 0.5,
                           position = ggplot2::position_dodge(0.1)) +
    ggrepel::geom_label_repel(size = 1.5, family = "IBM Plex Sans",
                        show.legend = FALSE,
                        position = ggplot2::position_dodge(0.1)) +
    ggsci::scale_color_jama() +
    ggplot2::theme_classic(10, "IBM Plex Sans") +
    ggplot2::labs(x = "Assessment", y = "GAD-7 score", color = "Treatment") +
    ggplot2::theme(legend.position = "top")

E o gráfico da GAD-7 per-protocol.

Código
# Calcular média e desvio padrão da PHQ por grupo e tempo (per-protocol)
plot_gad_pp <- all_waves |>
    dplyr::filter(email %in% email_pp) |>
    rcompanion::groupwiseMean(gad_total ~ time + group,
                              data = _,
                              na.rm = TRUE) |>
    dplyr::mutate(group = dplyr::if_else(group == "tccg", "Group e-CBT", "Thrive App")) |>
    ggplot2::ggplot(ggplot2::aes(x = time, y = Mean, group = group,
                                 color = group, label = paste0(
                                                            "n = ",
                                                            as.character(n),
                                                            "\n",
                                                            as.character(Mean),
                                                            " (",
                                                            as.character(Trad.lower),
                                                            ", ",
                                                            as.character(Trad.upper),
                                                            ")"
                                                            )
                                 )
    ) +
    ggplot2::geom_line(linewidth = 0.5,
                       position = ggplot2::position_dodge(0.1)) +
    ggplot2::geom_errorbar(ggplot2::aes(ymin = Trad.lower, ymax = Trad.upper),
                           linewidth = 0.5,
                           position = ggplot2::position_dodge(0.1)) +
    ggrepel::geom_label_repel(size = 1.5, family = "IBM Plex Sans",
                        show.legend = FALSE,
                        position = ggplot2::position_dodge(0.1)) +
    ggsci::scale_color_jama() +
    ggplot2::theme_classic(10, "IBM Plex Sans") +
    ggplot2::labs(x = "Assessment", y = "GAD-7 score", color = "Treatment") +
    ggplot2::theme(legend.position = "top")

Ao final, os gráficos são unidos da seguinte forma:

  • a1): PHQ-9 ITT
  • a2): PHQ-9 PP
  • b1): GAD-7 ITT
  • b2): GAD-7 PP
Código
library(patchwork)

((plot_phq_itt + plot_phq_pp) + plot_layout(tag_level = "new")) /
    ((plot_gad_itt + plot_gad_pp) + plot_layout(tag_level = "new")) + plot_annotation(tag_levels = c("a", "1"),
                                                   tag_suffix = ")") +
    plot_layout(guides = "collect") & ggplot2::theme(legend.position='bottom')

Modelos lineares mistos para avaliar melhora com o tempo

Código
itt_data <- all_waves |>
    dplyr::mutate(time = as.numeric(time) - 1)

pp_data <- itt_data |>
    dplyr::filter(email %in% email_pp)

itt_phq <- lmerTest::lmer(phq_total ~ time * group + (1 | email),
                          data = itt_data)

pp_phq <- lmerTest::lmer(phq_total ~ time * group + (1 | email),
                          data = pp_data)

itt_phq |> summary()
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: phq_total ~ time * group + (1 | email)
   Data: itt_data

REML criterion at convergence: 1041.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9274 -0.6016 -0.0659  0.5885  3.2064 

Random effects:
 Groups   Name        Variance Std.Dev.
 email    (Intercept)  9.837   3.136   
 Residual             17.601   4.195   
Number of obs: 174, groups:  email, 61

Fixed effects:
                 Estimate Std. Error       df t value   Pr(>|t|)    
(Intercept)       16.3236     1.0393 103.5951  15.706    < 2e-16 ***
time              -2.1090     0.4412 121.7589  -4.780 0.00000496 ***
groupthrive       -2.7633     1.3093 108.9137  -2.110     0.0371 *  
time:groupthrive   1.3235     0.5976 128.7106   2.215     0.0285 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) time   grpthr
time        -0.570              
groupthrive -0.794  0.452       
tim:grpthrv  0.421 -0.738 -0.558
Código
pp_phq |> summary()
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: phq_total ~ time * group + (1 | email)
   Data: pp_data

REML criterion at convergence: 581.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7817 -0.5834 -0.1331  0.6014  2.1646 

Random effects:
 Groups   Name        Variance Std.Dev.
 email    (Intercept) 12.77    3.573   
 Residual             14.81    3.848   
Number of obs: 100, groups:  email, 25

Fixed effects:
                 Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       15.7067     1.2418 43.1496  12.648 4.13e-16 ***
time              -1.9933     0.4443 73.0000  -4.486 2.65e-05 ***
groupthrive       -1.4567     1.9635 43.1496  -0.742    0.462    
time:groupthrive   0.8433     0.7025 73.0000   1.200    0.234    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) time   grpthr
time        -0.537              
groupthrive -0.632  0.339       
tim:grpthrv  0.339 -0.632 -0.537
Código
itt_gad <- lmerTest::lmer(gad_total ~ time * group + (1 | email),
                          data = itt_data)

pp_gad <- lmerTest::lmer(gad_total ~ time * group + (1 | email),
                          data = pp_data)

itt_gad |> summary()
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: gad_total ~ time * group + (1 | email)
   Data: itt_data

REML criterion at convergence: 985

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.11168 -0.55589  0.00863  0.62531  2.88635 

Random effects:
 Groups   Name        Variance Std.Dev.
 email    (Intercept)  6.981   2.642   
 Residual             13.713   3.703   
Number of obs: 172, groups:  email, 61

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)       13.6309     0.8995 102.1111  15.155  < 2e-16 ***
time              -1.5192     0.3912 118.5813  -3.884  0.00017 ***
groupthrive       -1.5596     1.1339 107.6101  -1.375  0.17185    
time:groupthrive   1.0107     0.5286 125.1772   1.912  0.05817 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) time   grpthr
time        -0.578              
groupthrive -0.793  0.458       
tim:grpthrv  0.428 -0.740 -0.566
Código
pp_gad |> summary()
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: gad_total ~ time * group + (1 | email)
   Data: pp_data

REML criterion at convergence: 554.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.15968 -0.46249 -0.01438  0.53032  2.52342 

Random effects:
 Groups   Name        Variance Std.Dev.
 email    (Intercept)  7.028   2.651   
 Residual             12.774   3.574   
Number of obs: 99, groups:  email, 25

Fixed effects:
                 Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       13.8733     1.0318 50.8822  13.445  < 2e-16 ***
time              -1.6600     0.4127 72.1293  -4.022  0.00014 ***
groupthrive       -2.1254     1.6320 50.9208  -1.302  0.19865    
time:groupthrive   1.2079     0.6537 72.1748   1.848  0.06874 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) time   grpthr
time        -0.600              
groupthrive -0.632  0.379       
tim:grpthrv  0.379 -0.631 -0.597