index.Rmd

---
title: "Fatores perinatais e desenvolvimento de transtorno bipolar"
author: "Bruno Montezano"
date: "`r Sys.Date()`"
lang: "pt-br"
output:
    html_document:
        theme: "paper"
        highlight: "kate"
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width=16, fig.height=9)
```

## Introdução

O presente trabalho almeja observar os efeitos de fatores perinatais e
do aleitamento materno no desenvolvimento de transtorno bipolar (TB) aos
18 e 22 anos de idade. Os dados utilizados para tal foram coletados em uma
coorte de nascimento coordenada pelo Centro de Epidemiologia da
Universidade Federal de Pelotas intitulada "Coorte de 1993".

```{r}
# Seed para reprodutibilidade dos resultados
set.seed(1)
# Carrega dados da Coorte de 1993
dados_coorte <- haven::read_sav("data/banco_QI_15_07.sav",
                                user_na = FALSE)
```

## Transtorno bipolar

O diagnóstico de transtorno bipolar será calculado pela MINI (Mini
International Neuropsychiatric Interview) por meio dos critérios
diagnósticos do DSM-IV. Serão calculados diagnósticos de transtorno
bipolar tipo I e transtorno bipolar tipo II.

```{r}
# Número de sujeitos incluídos
dados_coorte |> 
    dplyr::filter((!is.na(MDD_BD1_2_notmood)) |
                  (!is.na(TB1_18anos)) |
                  (!is.na(TBII.p1_18anos))) |> 
    nrow()

# Quem foi avaliado em cada tempo?
#dados_coorte |> 
#    dplyr::filter((!is.na(MDD_BD1_2_notmood)) |
#                  (!is.na(TB1_18anos)) |
#                  (!is.na(TBII.p1_18anos))) |> 
#    dplyr::mutate(
#        av18 = dplyr::if_else(!is.na(MDD_BD1_2_notmood), "Sim", "Não"),
#        av22 = dplyr::if_else(!is.na(kmini1), "Sim", "Não"),
#        av_quais = dplyr::case_when(
#            av18 == "Sim" & av22 == "Sim" ~ "Avaliado em ambos os tempos",
#            av18 == "Sim" & av22 == "Não" ~ "Avaliado aos 18 anos apenas",
#            av18 == "Não" & av22 == "Sim" ~ "Avaliado aos 22 anos apenas"
#        )
#    ) |> 
#    dplyr::count(av_quais)

# Criar variável de desfecho, além de uma variável com os subtipos de TB
dados_com_desfecho <- dados_coorte |> 
    haven::as_factor() |> 
    dplyr::mutate(
        tb1_22 = dplyr::if_else(
            MDD_BD1_2_notmood == "Bipolar Disorder 1",
            "Sim", "Não"
        ),
        tb2_22 = dplyr::if_else(
            MDD_BD1_2_notmood == "Bipolar Disorder 2",
            "Sim", "Não"
        ),
        tb1_18 = dplyr::if_else(
            TB1_18anos == 1,
            "Sim", "Não"
        ),
        tb2_18 = dplyr::if_else(
            TBII.p1_18anos == 1,
            "Sim", "Não"
        ),
        outcome_subtypes = dplyr::case_when(
            tb1_22 == "Sim" ~ "TB I aos 22 anos",
            tb1_18 == "Sim" ~ "TB I aos 18 anos",
            tb2_22 == "Sim" ~ "TB II aos 22 anos",
            tb2_18 == "Sim" ~ "TB II aos 18 anos",
            tb1_22 == "Não" |
            tb2_22 == "Não" |
            tb1_18 == "Não" |
            tb2_18 == "Não" ~ "Sem TB",
            TRUE ~ NA_character_
        ),
        outcome = dplyr::case_when(
            tb1_22 == "Sim" |
            tb2_22 == "Sim" |
            tb1_18 == "Sim" |
            tb2_18 == "Sim" ~ "Sim",
            tb1_22 == "Não" |
            tb2_22 == "Não" |
            tb1_18 == "Não" |
            tb2_18 == "Não" ~ "Não",
            TRUE ~ NA_character_
        )
    ) |> 
    dplyr::filter(!is.na(outcome))
```


### Criação das variáveis de episódio de humor e filtragem para reter variáveis perinatais

```{r}
# Vetor com as variáveis perinatais indicadas pela pesquisadora para
# construção de novas features
vars_vanessa <- c("apartip", "ainduz", "aprobrn", "aerrado",
"aprobrn", "aprobrn1", "aprobrn2", "aprobrn3",
"aprobrn", "aprobrn1", "aprobrn2", "aprobrn3",
"aprobrn", "aprobrn1", "aprobrn2", "aprobrn3",
"aprobrn", "aprobrn1", "aprobrn2", "aprobrn3",
"ahipert", "adiabet", "ainfecur", "aoutinf",
"aanemia", "amotcesa", "apartind", "afumou",
"acompfum", "abebalc", "aalcool", "avivmar",
"asexrn", "apesorn", "acompr", "apcrn", "aescol3",
"aescpai", "aescmae", "aparidad", "aapoio", "aidadmae",
"aidadpai", "apretama", "aapgar1", "aapgar5",
"arenfam")

# Filtrar base para manter subamostra da avaliação aos 22 anos
# Transformar variáveis com rótulos em variáveis tipo fator
# Criar subconjunto com variável diagnóstico e variáveis perinatais
da <- dados_com_desfecho |> 
    dplyr::select(outcome, dplyr::any_of(vars_vanessa))
```

### Criação do diagnóstico de TB

Além da manutenção de variável de TB, a base de dados foi filtrada para
manter somente os sujeitos avaliados aos 18 ou 22 anos (através da filtragem
por meio de uma variável de diagnóstico da doença da avaliação em 2015).

Abaixo, segue um gráfico de barras com a proporção de bipolares e
não-bipolares na amostra:

```{r}
da |> 
    dplyr::count(outcome) |> 
    ggplot2::ggplot(ggplot2::aes(x = outcome, y = n)) +
    ggplot2::geom_col(fill = "royalblue") +
    ggplot2::geom_text(ggplot2::aes(label = n), vjust = -0.5, size = 7) +
    ggplot2::theme_minimal(base_size = 20) +
    ggplot2::labs(x = "\nDiagnóstico",
                  y = "Número de sujeitos\n",
                  title = "Distribuição de sujeitos com diagnóstico e sem diagnóstico de transtorno bipolar",
                  subtitle = "Avaliações aos 18 e 22 anos") +
    ggplot2::theme(
        axis.text.x = ggplot2::element_text(size = 20),
        axis.text.y = ggplot2::element_text(size = 20)
        )
```

E mais um gráfico, agora para demonstrar a distribuição dos
subtipos de transtorno bipolar (tipo I e tipo II) e o momento em que o
diagnóstico foi concretizado:

```{r}
dados_com_desfecho |> 
    dplyr::count(outcome_subtypes) |> 
    ggplot2::ggplot(ggplot2::aes(x = outcome_subtypes, y = n)) +
    ggplot2::geom_col(fill = "darkorange") +
    ggplot2::geom_text(ggplot2::aes(label = n), vjust = -0.5, size = 7) +
    ggplot2::theme_minimal(base_size = 20) +
    ggplot2::labs(x = "\nSubtipo e idade de diagnóstico",
                  y = "Número de sujeitos\n",
                  title = "Distribuição dos subtipos de transtorno bipolar (TB) e idade de diagnóstico",
                  subtitle = "Avaliação do desfecho aos 18 e 22 anos") +
    ggplot2::theme(
        axis.text.x = ggplot2::element_text(size = 20),
        axis.text.y = ggplot2::element_text(size = 20)
        )
```

## Fatores perinatais e transtorno bipolar

### Construção de variáveis e técnicas de imputação

As seguintes variáveis foram construídas:

-   `tb_dic`: Diagnóstico de TB tipo I ou TB tipo II (dicotômica)
- `complicacoes`: Variáveis de prematuridade, traumatismo no parto, disfunção
respiratória no parto, malformação no parto e sofrimento fetal agrupadas
(através do operador OU)
-   `cesarea`: O parto foi cesária?
-   `parto_induzido`: A bolsa rompeu, foi utilizado soro ou ambas as
    coisas?
-   `uti_neo`: O recém-nascido foi para UTI?
-   `hipert_materna`: Hipertensão arterial materna
-   `dm_materna`: Diabetes materna
-   `infec_gest`: Infecção gestacional
-   `anemia_gest`: Anemia gestacional
-   `hemorr_parto`: Apresentou hemorragia no parto
-   `tab_gest`: Contato com tabagismo na gestação (ativo ou passivo)
-   `etil_pesado_gest`: Etilismo pesado na gestação
-   `rup_famil`: Ruptura familiar
-   `baixo_peso`: Baixo peso ao nascer
-   `baixo_comp`: Baixo comprimento ao nascer
-   `per_cef`: Perímetro cefálico (numérico)
-   `escol_paterna`: Escolaridade paterna
-   `escol_materna`: Escolaridade materna
-   `paridade`: Paridade (dicotômica) [4 ou mais e menos que 4]
-   `apoio_familiar`: Teve apoio da famíla/vizinhos/amigos?
-   `idade_mae_num`: Idade da mãe (numérica)
-   `pretende_amamentar`: Mãe pretende amamentar
-   `apgar_1`: Apgar no minuto 1 (escore bruto)
-   `apgar_5`: Apgar no minuto 5 (escore bruto)
-   `renda_familiar`: Renda familiar (em salários mínimos)

As variáveis de idade materna foram utilizadas através das três codificações
pois os autores consideraram a hipótese da relação da idade com o desfecho não
ser puramente linear. Dessa forma, esta foi uma forma na qual se pensou uma
maneira de considerar possível associação.

```{r}
da_proc <- da |>
    dplyr::transmute(
        # Desfecho: transtorno bipolar aos 18 ou 22 anos de idade
        tb_dic = factor(dplyr::if_else(outcome == "Sim",
                                         "Sim",
                                         "Não"
                                         ),
                        levels = c("Não", "Sim")
                           ),
        # Cesárea ou parto normal
        cesarea = as.factor(dplyr::if_else(apartip == "normal",
                                           "Não",
                                           "Sim"
                                           )
                            ),
        # Parto induzido
        parto_induzido = as.factor(dplyr::if_else(
            ainduz %in% c("rompeu bolsa", "soro", "ambos"),
            "Sim",
            "Não"
        )),
        # Ida a UTI
        uti_neo = as.factor(dplyr::if_else(
            aprobrn %in% c("bercario", "aloj.conjunto"),
            "Não",
            dplyr::if_else(aprobrn == "uti", "Sim",
                           NA_character_
                           )
        )),
        # Nascimento prematuro
        prematuro = as.factor(
            dplyr::if_else(aerrado != "prematuro" | is.na(aerrado),
                           "Não",
                           "Sim"
                           )
        ),
        # Traumatismo no parto
        traumatismo = as.factor(
            dplyr::if_else(
                aprobrn1 == "traumatismo" | aprobrn2 == "traumatismo" |
                    aprobrn3 == "traumatismo",
                "Sim",
                "Não",
                missing = "Não"
            )
        ),
        # Estresse respiratório no parto
        disf_resp = as.factor(
            dplyr::if_else(
                aprobrn1 == "stress resp" | aprobrn2 == "stress resp" |
                    aprobrn3 == "stress resp",
                "Sim",
                "Não",
                missing = "Não"
            )
        ),
        # Malformação
        malform = as.factor(
            dplyr::if_else(
                aprobrn1 == "malformacao" | aprobrn2 == "malformacao" |
                    aprobrn3 == "malformacao",
                "Sim",
                "Não",
                missing = "Não"
            )
        ),
        # Sofrimento fetal
        sofrfetal = as.factor(
            dplyr::if_else(
                aprobrn1 == "sofrim fetal" | aprobrn2 == "sofrim fetal" |
                    aprobrn3 == "sofrim fetal",
                "Sim",
                "Não",
                missing = "Não"
            )
        ),
        # Hipertensão materna
        hipert_materna = relevel(
            as.factor(
                dplyr::case_when(
                    as.character(ahipert) == "nao sabe" ~ NA_character_,
                    TRUE ~ as.character(ahipert)
                )
            ), ref = "nao"),
        # Diabetes materna
        dm_materna = as.factor(dplyr::if_else(
            adiabet %in% c("tratado", "nao tratado"),
            "Sim",
            "Não"
        )),
        # Infecção gestacional
        infec_gest = as.factor(
            dplyr::if_else(
                ainfecur %in% c("tratado", "nao tratado") |
                    aoutinf %in% c("tratado", "nao tratado"),
                "Sim",
                "Não"
            )
        ),
        # Anemia gestacional
        anemia_gest = as.factor(dplyr::if_else(
            aanemia %in% c("tratado", "nao tratado"),
            "Sim",
            "Não"
        )),
        # Hemorragia no parto
        hemorr_parto = as.factor(
            dplyr::if_else(
                amotcesa == "hemorragia materna" | apartind == "sangramento",
                "Sim",
                "Não"
            )
        ),
        # Tabagismo gestacional
        tab_gest = as.factor(
            dplyr::if_else(afumou == "sim" | acompfum == "sim",
                           "Sim",
                           "Não"
                           )
        ),
        # Etilismo pesado na gestação
        etil_pesado_gest = as.factor(
            dplyr::if_else(aalcool == "sim",
                           "Sim",
                           "Não"
                           )
        ),
        # Ruptura familiar
        rup_famil = as.factor(
            dplyr::if_else(avivmar == "sim",
                           "Não",
                           "Sim"
                           )
        ),
       # Baixo peso ao nascer
       baixo_peso = as.factor(
           dplyr::if_else(apesorn < 2500,
                          "Sim",
                          "Não"
                          )
       ),
       # Baixo comprimento ao nascer
       baixo_comp = as.factor(
           dplyr::if_else(acompr < 47,
                          "Sim",
                          "Não"
                          )
       ),
       # Perímetro cefálico
       per_cef = apcrn,
       # Escolaridade paterna (em anos)
       escol_paterna = aescpai,
       # Escolaridade materna (em anos)
       escol_materna = aescmae,
       # Paridade (número de filhos)
       paridade = relevel(as.factor(dplyr::if_else(
           aparidad %in% c("0", "1", "2", "3"),
            "Menor que 4",
            "Maior ou igual a 4")
           ), ref = "Menor que 4"), 
       # Apoio familiar
       apoio_familiar = as.factor(
           dplyr::if_else(
               aapoio == "sim",
               "Sim",
               "Não"
           )
       ),
       # Idade da mãe
       idade_mae_num = aidadmae,
       # Pretende amamentar
       pretende_amamentar = as.factor(
           dplyr::if_else(
               apretama == "sim",
               "Sim",
               "Não"
           )
       ),
       # Apgar no primeiro minuto
       apgar_1 = aapgar1,
       # Apgar no minuto 5
       apgar_5 = aapgar5,
       # Renda familiar em salários mínimos (SM)
       renda_familiar = arenfam
    )
```

Nesta etapa abaixo, realizou-se uma checagem em relação aos valores
ausentes em cada uma das variáveis. Posteriormente, elaborou-se um
pequeno fluxo para realizar a imputação das variáveis por meio de mediana
(variáveis numéricas) e moda (variáveis categóricas).

```{r}
# Checar valores ausentes após criação das variáveis propostas
da_proc |> 
    purrr::map_if(is.factor, \(x) table(x, useNA = "always"), .else = summary)

# Agrupar variáveis de complicações perinatais
# O argumento .keep da função dplyr::mutate permite
# que somente as variáveis não utilizadas nas operações
# sejam mantidas na base de dados resultante
da_proc <- da_proc |> 
    dplyr::mutate(
        complicacoes = as.factor(
            dplyr::if_else(
                prematuro == "Sim" | traumatismo == "Sim" |
                disf_resp == "Sim" | malform == "Sim" |
                sofrfetal == "Sim",
                "Sim", "Não"
            )
        ),
        .keep = "unused"
    )

# Criar base com variáveis imputadas pela moda e mediana
da_imp <- da_proc |> 
    recipes::recipe(tb_dic ~ .) |> 
    recipes::step_filter_missing(recipes::all_predictors(), threshold = 0.1) |> 
    recipes::step_impute_mode(recipes::all_nominal_predictors()) |> 
    recipes::step_impute_median(recipes::all_numeric_predictors()) |> 
    recipes::step_nzv(recipes::all_predictors()) |> 
    recipes::prep() |> 
    recipes::bake(new_data = NULL)


# Checar que a imputação funcionou
da_imp |> 
    purrr::map_if(is.factor, \(x) table(x, useNA = "always"), .else = summary)
```

```{r}
# A função abaixo foi criada para verificar a frequência absoluta de
# uma variável independente qualquer estratificada pelo desfecho
ver_freq <- function(var) {
    
    da_imp |> 
        dplyr::group_by(tb_dic) |> 
        dplyr::count({{ var }})
    
}

# Checar frequência agrupada pelo desfecho
da_imp |> 
    purrr::map_if(is.factor, ver_freq, .else = summary)
```

### Quais variáveis foram excluídas e quando?

Na etapa de remoçao de variáveis com mais de 10% de missing,
exclui-se a variável de hemorragia no parto, apgar no primeiro minuto e
apgar no quinto minuto.

Na etapa de retirada de variáveis com pouca ou nenhuma variabilidade,
foram excluídas as seguintes variáveis: UTI neonatal, diabetes materna,
etilismo pesado na gestação, se a mãe pretende amamentar e complicações
no parto.

```{r, eval = FALSE}
da_prep <- da_proc |> 
    recipes::recipe(tb_dic ~ .) |> 
    recipes::step_filter_missing(recipes::all_predictors(), threshold = 0.1) |> 
    recipes::step_impute_mode(recipes::all_nominal_predictors()) |> 
    recipes::step_impute_median(recipes::all_numeric_predictors()) |> 
    recipes::step_nzv(recipes::all_predictors()) |> 
    recipes::prep()

da_prep$steps |> View()
```

### Análise bivariada com variáveis categóricas

```{r}
# Carregar o pacote magrittr para utilizar o placeholder "." na função
# do teste qui-quadrado abaixo
library(magrittr, include.only = "%>%")

cat_biv <- da_imp |>
    purrr::keep(is.factor) |>
    tidyr::pivot_longer(names_to = "variable", values_to = "value", -tb_dic) |>
    dplyr::group_by(variable) %>%
    dplyr::do(chisq.test(.$tb_dic, .$value, simulate.p.value = TRUE) |>
                  broom::tidy()) |> 
    dplyr::select(-c(parameter, method)) |> 
    dplyr::ungroup()

knitr::kable(cat_biv,
             digits = 3,
             col.names = c("Variável",
                           "Estatística",
                           "p-valor"))
```

### Análise bivariada com variáveis numéricas

```{r}
# Realiza análise bivariada com as variáveis numéricas através de teste t
num_biv <- da_imp |>
    dplyr::select(tb_dic, where(is.numeric)) |>
    tidyr::pivot_longer(names_to = "variable", values_to = "value", -tb_dic) |>
    dplyr::group_by(variable) %>%
    dplyr::do(t.test(.$value ~ .$tb_dic) |>
                  broom::tidy()) |> 
    dplyr::ungroup() |> 
    dplyr::select(-method, -alternative)

knitr::kable(num_biv,
             digits = 3)
```

### Análise multivariada

A imputação foi realizada nos passos anteriores para a análise
bivariada, logo, o objeto da base de dados imputada (pré-processada)
será usada para o ajuste do modelo de regressão logística binomial.

Um modelo contendo somente as variáveis
com valor *p* menor que 0,2 na análise bivariada foi construído.

Os valores de odds ratio foram estimados a partir da equação $e^{\beta}$.

#### Modelo de regressão logística binomial

Neste modelo, foram incluídas somente as variáveis independentes que
apresentaram um valor de p menor que 0,2 na análise bivariada (bruta).

```{r}
# Pegar nomes das variáveis que apresentaram p-value menor que 0,2 na análise
# bruta
vars_multivar <- cat_biv |> 
    dplyr::bind_rows(num_biv |> dplyr::select(variable, statistic, p.value)) |> 
    dplyr::ungroup() |> 
    dplyr::filter(p.value < 0.2) |> 
    dplyr::pull(variable)

# Adicionar desfecho na fórmula juntamente aos preditores que forem incluídos
# como argumento da função
add_outcome <- function(predictors) {
    reg_mod <- paste0("tb_dic", " ~ ", predictors)
    return(reg_mod)
}

# Juntar os preditores separados por + (para sintaxe de fórmula) e ao final,
# adicionar o desfecho
mod_filtrado <- vars_multivar |> 
    paste(collapse = " + ") |> 
    add_outcome()

# Especificar modelo a ser ajustado
mod <- parsnip::logistic_reg() |> 
    parsnip::set_mode("classification") |> 
    parsnip::set_engine("glm")

# Ajustar a partir da fórmula construída previamente
mod_filtrado_fit <- mod |> 
    parsnip::fit(as.formula(mod_filtrado), data = da_imp)

# Mostrar resultados tidy
mod_filtrado_fit |> 
    broom::tidy() |> 
    knitr::kable(
        col.names = c("Variável",
                      "Beta",
                      "Erro Padrão",
                      "Estatística",
                      "p-valor"),
        digits = 3,
        caption = "Resumo do modelo final com variáveis que apresentaram p-valor menor que 0,2 na análise bruta."
    )
```

#### Plot com OR do modelo com p \< 0,2 na análise bivariada

Como pode ser observado na tabela acima, traumatismo no parto se mostrou como
um fator de risco para o desenvolvimento de TB aos 18 ou aos 22 anos, porém não
está presente na figura abaixo por ter apresentado um odds ratio que ficaria
muito fora dos limites do gráfico.

Além disso, a variável de traumatismo no parto, como pode ser observado na tabela
abaixo, apresentou um intervalo de confiança superior na razão de chances de 74,8.

```{r, warning = FALSE}
# Seed para reprodutibilidade
set.seed(1)

# Criação dos intervalos de confiança para os parâmetros do modelo ajustado
intervalos <- exp(confint(mod_filtrado_fit$fit)) |>
    tibble::as_tibble(rownames = "term") |> 
    janitor::clean_names()

# Mostrar odds ratio com os respectivos intervalos de confiança
or_table <- mod_filtrado_fit |> 
    broom::tidy() |> 
    dplyr::left_join(intervalos, by = "term") |> 
    dplyr::transmute(
        variavel = term,
        or = exp(estimate),
        lower_int = x2_5_percent,
        upper_int = x97_5_percent,
        pvalue = p.value
        )

or_table |> 
    knitr::kable(
        col.names = c("Variável",
                      "Odds Ratio",
                      "Intervalo Inferior",
                      "Intervalo Superior",
                      "p-valor"),
        caption = "Razão de chances para cada variável incluída no modelo final.",
        digits = 3
    )

# Criar plot com odds ratio listados acima
or_table |> 
    dplyr::filter(variavel != "(Intercept)") |> 
    ggplot2::ggplot() +
    ggplot2::aes(x = or, y = variavel) +
    ggplot2::geom_point(size = 5, color = "royalblue") +
    ggplot2::geom_linerange(
        ggplot2::aes(xmin = lower_int, xmax = upper_int),
        color = "royalblue",
        size = 0.8) +
    ggplot2::geom_vline(xintercept = 1) +
    ggplot2::scale_x_continuous(limits = c(0, 3)) +
    ggplot2::scale_y_discrete(labels = c("Baixo comprimento (sim)",
                                         "Baixo peso ao nascer (sim)",
                                         "Cesária (sim)",
                                         "Escolaridade materna (em anos)",
                                         "Escolaridade paterna (em anos)",
                                         "Idade da mãe",
                                         "Paridade (maior ou igual a 4)",
                                         "Perímetro cefálico",
                                         "Renda familiar (em SM)",
                                         "Tabagismo gestacional (ativo ou passivo)")) +
    ggplot2::labs(x = expression(paste("Razão de chances (", italic("odds ratio"), ")")),
                  y = "Variável") +
    ggplot2::theme_minimal(base_size = 20) +
    ggplot2::theme(
        axis.text.x = ggplot2::element_text(size = 20),
        axis.text.y = ggplot2::element_text(size = 20)
        )
```

```{r, eval = FALSE, warning = FALSE}
# Tabela de odds ratio processada
processed_or_table <- or_table |> 
    dplyr::filter(variavel %in% c("paridadeMaior ou igual a 4", "tab_gestSim",
                                  "escol_materna", "escol_paterna",
                                  "idade_mae_num", "per_cef",
                                  "renda_familiar")) |> 
    dplyr::mutate(
        dplyr::across(dplyr::ends_with("int"), \(x) as.character(format(x, digits = 2))),
        intervals = stringr::str_glue("({lower_int}-{upper_int})"),
        variavel = dplyr::case_when(
            variavel == "paridadeMaior ou igual a 4" ~ "Parity",
            variavel == "tab_gestSim" ~ "Tobbaco in pregnancy",
            variavel == "escol_materna" ~ "Maternal years of education",
            variavel == "escol_paterna" ~ "Paternal years of education",
            variavel == "idade_mae_num" ~ "Maternal age",
            variavel == "per_cef" ~ "Head circumference",
            variavel == "renda_familiar" ~ "Family income"
        ),
        pvalue = format(pvalue, digits = 3)
    ) |> 
    dplyr::select(-dplyr::ends_with("int")) |> 
    dplyr::select(variable = variavel, adjusted_or = or, intervals, pvalue)

# Exportar tabela para Word
processed_or_table |> 
    flextable::flextable() |> 
    flextable::save_as_docx(path = "output/tabela_or.docx")
```

```{r, eval = FALSE}
# Criar plot com odds ratio, porém apenas variáveis significativas da bivariada
plot_apenas_sig <- or_table |> 
    dplyr::filter(variavel %in% c("paridadeMaior ou igual a 4", "tab_gestSim",
                                  "escol_materna", "escol_paterna",
                                  "idade_mae_num", "per_cef",
                                  "renda_familiar")) |> 
    ggplot2::ggplot() +
    ggplot2::aes(x = or, y = variavel) +
    ggplot2::geom_point(size = 5, color = "royalblue") +
    ggplot2::geom_linerange(
        ggplot2::aes(xmin = lower_int, xmax = upper_int),
        color = "royalblue",
        size = 0.8) +
    ggplot2::geom_vline(xintercept = 1) +
    ggplot2::scale_x_continuous(limits = c(0, 3)) +
    ggplot2::scale_y_discrete(labels = c("Maternal education level",
                                         "Paternal education level",
                                         "Maternal age",
                                         "Multiparity",
                                         "Head circumference",
                                         "Family income",
                                         "Tobacco in pregnancy")) +
    ggplot2::labs(x = "Adjusted OR",
                  y = "Exposure") +
    ggplot2::theme_minimal(base_size = 20) +
    ggplot2::theme(
        axis.text.x = ggplot2::element_text(size = 20),
        axis.text.y = ggplot2::element_text(size = 20)
        )

ggplot2::ggsave(plot = plot_apenas_sig,
                filename = "output/or_plot.png",
                dpi = 300)

ggplot2::ggsave(plot = plot_apenas_sig,
                filename = "output/or_plot.pdf",
                dpi = 300)
```



## Tabela descritiva com dados perinatais

```{r, warning = FALSE}
pvalue <- function(x, ...) {
    # Construir vetores com y (desfecho), e grupos (estratificação) chamado g
    y <- unlist(x)
    g <- factor(rep(1:length(x), times=sapply(x, length)))
    if (is.numeric(y)) {
        # Para variáveis numéricas, performar um teste t de duas amostras
        p <- t.test(y ~ g)$p.value
    } else {
        # For categorical variables, perform a chi-squared test of independence
        # Para variáveis categóricas, performar um teste qui-quadrado de independência
        p <- chisq.test(table(y, g))$p.value
    }
    # Formatar o p-valor usando sintaxe em HTML para o símbolo de "menor que"
    # A string vazia coloca o output do p-valor na linha abaixo do rótulo da variável
    c("", sub("<", "&lt;", format.pval(p, digits=3, eps=0.001)))
}

# Gerar objeto da tabela com dados perinatais
# O underscore (_) indica o placeholder do pipe nativo do R
tabela_peri <- da_proc |> 
    table1::table1(~ cesarea + parto_induzido + uti_neo + hipert_materna +
                       dm_materna + infec_gest + anemia_gest + hemorr_parto +
                       tab_gest + etil_pesado_gest + rup_famil + baixo_peso +
                       baixo_comp + per_cef + escol_paterna + escol_materna +
                       paridade + apoio_familiar + idade_mae_num +
                       pretende_amamentar + apgar_1 + apgar_5 + renda_familiar +
                       complicacoes | tb_dic, data = _,
                   extra.col = list("p-value" = pvalue),
                   overall = FALSE)

# Transformar para formato kable para melhor visualização
tabela_peri |> 
    table1::t1kable()
```

```{r, eval = FALSE}
# Exportar tabela para Word
tabela_peri |> 
    table1::t1flex() |> 
    flextable::save_as_docx(path = "output/tabela_peri.docx")
```

## Tabela descritiva com dados do follow-up

```{r, warning = FALSE}
# ld023a = estado civil
#lescoljovemano = anos de estudo

# Gerar tabela descritiva com dados do follow-up
tabela_follow <- dados_com_desfecho |> 
    dplyr::rename(
        cocaina_lifetime = ldc007, # Necessita recode
        maconha_lifetime = ldc006, # Necessita recode
        status_socio = labep3,
        trabalha_atualmente = ld011,
        sexo = sexonovo,
        estado_civil = ld023a,
        anos_de_estudo = lescoljovemano) |> 
    dplyr::mutate(
        dplyr::across(dplyr::ends_with("lifetime"),
                      \(x) as.factor(dplyr::if_else(
                          as.character(x) == "Nunca usei",
                          "Não", "Sim"
                      )))
    ) |> 
    table1::table1(~ sexo + status_socio + estado_civil + anos_de_estudo + trabalha_atualmente + maconha_lifetime + cocaina_lifetime | outcome,
                   data = _, # Placeholder do pipe nativo
                   extra.col = list("p-value" = pvalue),
                   overall = FALSE)

# Transformar tabela em formato kable para melhor visualização no Rmd compilado
tabela_follow |> 
    table1::t1kable()
```

```{r, eval = FALSE}
# Exportar tabela para Word
tabela_follow |> 
    table1::t1flex() |> 
    flextable::save_as_docx(path = "output/tabela_follow.docx")
```

## Ocitocina e TB

```{r}
# ainduz = parto foi induzido (rompeu bolsa, soro, ambos, nao)

dados_oci <- dados_com_desfecho |> 
    dplyr::mutate(ocitocina = dplyr::if_else(
        ainduz == "soro", "Sim", "Não"
    ))

janitor::tabyl(
    dat = dados_oci,
    var1 = ocitocina,
    var2 = outcome,
    show_na = TRUE
)

chi_oci <- chisq.test(x = dados_oci$ocitocina,
           y = dados_oci$outcome,
           correct = TRUE)
```

O p-valor do teste qui-quadrado entre ocitocina e o desfecho de TB no início da
vida adulta foi de `r chi_oci$p.value`

# Possível mediação da renda familiar no efeito da paridade

Abaixo, cria-se um modelo de modelagem de equações estruturais para verificar
a possível mediação da renda familiar no efeito da paridade (número de filhos)
no desfecho de desenvolvimento de transtorno bipolar na idade adulta.

```{r}
set.seed(1)

da_med <- da_imp |> 
    dplyr::select(tb_dic, paridade, renda_familiar) |> 
    dplyr::mutate(
        paridade = dplyr::if_else(
            paridade == "Menor que 4", 0, 1
        )
    )

med_mod <- "# Mediação da renda para paridade como preditor de TB
renda_familiar ~ a*paridade
tb_dic ~ b*renda_familiar + c*paridade
# Efeito indireto
indireto := a * b
# Efeito direto
direto := c
# Efeito total
total := c + (a * b)"

med_fit <- lavaan::sem(
    med_mod,
    data = da_med,
    ordered = "tb_dic"
)

mc_mediation_ci <- semTools::monteCarloCI(med_fit, nRep = 20000)

mc_mediation_ci |>
    tibble::rownames_to_column("Efeito") |>
    dplyr::rename("Beta" = est,
                  "IC 2.5%" = ci.lower,
                  "IC 97.5%" = ci.upper) |>
    knitr::kable(caption = "Beta e intervalos de confiança de cada efeito do modelo de mediação baseado no método Monte Carlo com 20.000 amostras.",
                 digits = 3)

broom::tidy(med_fit, conf.int = TRUE) |> 
    dplyr::select("Termo" = term,
                  "Rótulo" = label,
                  "Beta" = estimate,
                  "Erro" = std.error,
                  "p-valor" = p.value,
                  "IC 2.5%" = conf.low,
                  "IC 97.5%" = conf.high,
                  "Beta_std" = std.all) |> 
    knitr::kable(caption = "Tabela com parâmetros do modelo de modelagem de equações estruturais e respectivos rótulos.", digits = 3, format = "html")
```

```{r, eval = FALSE}
resumo_mediacao <- broom::tidy(med_fit, conf.int = TRUE) |> 
    dplyr::select("Termo" = term,
                  "Rótulo" = label,
                  "Beta" = estimate,
                  "Erro" = std.error,
                  "p-valor" = p.value,
                  "IC 2.5%" = conf.low,
                  "IC 97.5%" = conf.high,
                  "Beta_std" = std.all)

resumo_mediacao |> 
    dplyr::filter(Rótulo %in% c("indireto", "direto", "total")) |> 
    dplyr::select(-Termo) |> 
    dplyr::rename(B = Beta,
                  "SE B" = Erro,
                  pvalue = `p-valor`,
                  beta = Beta_std,
                  ic_lower = `IC 2.5%`,
                  ic_upper = `IC 97.5%`
                  ) |> 
    dplyr::mutate(
        dplyr::across(c(dplyr::starts_with("ic"), pvalue, B, beta, `SE B`),
                      \(x) as.character(format(x, digits = 2))),
        intervals = stringr::str_glue("({ic_lower}-{ic_upper})"),
    ) |> 
    dplyr::select(
        Effect = Rótulo,
        B,
        intervals,
        `SE B`,
        beta,
        pvalue
    ) |> 
    flextable::flextable() |> 
    flextable::save_as_docx(path = "output/tabela_mediacao.docx")


```