index.qmd

```{r}
#| label: Setup
#| include: false

library(here)

source(here("R", "_setup.R"))
```

<!-- badges: start -->
[![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
[![OSF DOI](https://img.shields.io/badge/OSF-10.17605/OSF.IO/8J94M-1284C5.svg)](https://doi.org/10.17605/OSF.IO/8J94M)
[![License: GPLv3](https://img.shields.io/badge/license-GPLv3-bd0000.svg)](https://www.gnu.org/licenses/gpl-3.0)
[![License: CC BY-NC-SA 4.0](https://img.shields.io/badge/license-CC_BY--NC--SA_4.0-lightgrey.svg)](https://creativecommons.org/licenses/by-nc-sa/4.0/)
<!-- badges: end -->

## Overview

This report provides a reproducible pipeline for processing and analyzing the microdata on nutritional status monitoring in Brazil from the Brazilian Food and Nutrition Surveillance System ([SISVAN](https://sisaps.saude.gov.br/sisvan/)), focusing on the nutritional status of children aged 0–5 years (i.e., younger than 60 months).

If you are working with other age groups, you will need to adapt the code accordingly. We provide some guidance on how to do this along the report.

For instructions on how to run the pipeline, see the repository [README](https://github.com/sustentarea/nutritional-status/blob/main/README.md).

Click [here](https://sustentarea.github.io/nutritional-status-analysis/) to see a report with a longitudinal analysis of the processed data.

## Problem

The Food and Nutrition Surveillance System ([SISVAN](https://sisaps.saude.gov.br/sisvan/)) is a strategic tool for monitoring the nutritional status of the Brazilian population, particularly those served by Brazil's Unified Health System ([SUS](https://www.gov.br/saude/pt-br/sus)). However, despite its broad scope and importance, the anthropometric data recorded in SISVAN often suffer from accessability and quality issues that limit their usefulness for rigorous analyses and evidence-based policymaking [@silva2023a].

Multiple factors contribute to these quality concerns, including the lack of standardized measurement protocols, variability in staff training, inconsistencies in data entry and processing, and incomplete population coverage [@bagni2015; @corsi2017; @perumal2020]. To assess and improve data quality, several indicators have been proposed and applied, such as population coverage [@nascimento2017; @mourao2020], completeness of birth dates and anthropometric measurements [@finaret2018; @nannan2019], digit preference for age, height, and weight [@lyons-amos2017; @bopp2008], the percentage of biologically implausible values [@lawman2015], and the dispersion and distribution of standardized weight and height measurements [@perumal2020; @mei2007].

In light of these challenges, there is a need for an open and reproducible pipeline to process SISVAN microdata. Such a pipeline should facilitate broader access to the data and systematically identify, correct, and remove problematic records, thereby improving the consistency, completeness, and plausibility of the information for research and policymaking.

## Data Availability

::: {style="text-align: left;"}
[![](https://img.shields.io/badge/OSF-10.17605/OSF.IO/8J94M-1284C5.svg)](https://doi.org/10.17605/OSF.IO/8J94M)
:::

The processed data are available in [`csv`](https://en.wikipedia.org/wiki/Comma-separated_values), [`rds`](https://rdrr.io/r/base/readRDS.html), and [`parquet`](https://en.wikipedia.org/wiki/Apache_Parquet) formats via a dedicated repository on the Open Science Framework ([OSF](https://osf.io)), accessible [here](https://doi.org/10.17605/OSF.IO/8J94M). Each dataset is accompanied by a metadata file describing its structure and contents.

You can also retrieve these files directly from [R](https://www.r-project.org/) using the [`osfr`](https://docs.ropensci.org/osfr/) package.

## Methods

### Source of Data

The data used in this report come from the following sources:

- Brazilian Food and Nutrition Surveillance System ([SISVAN](https://sisaps.saude.gov.br/sisvan/)):
  - Microdata on nutritional status monitoring in Brazil [@sisvana], the primary dataset for this pipeline.
- Brazilian Institute of Geography and Statistics ([IBGE](https://www.ibge.gov.br/)):
  - Official codes and metadata for Brazilian municipalities, incorporated via the [`geobr`](https://ipeagit.github.io/geobr/) R package [@pereirab], used to normalize IBGE municipality codes and enrich the analysis with geographic information.
- Department of Informatics of the Brazilian Unified Health System ([DATASUS](https://datasus.saude.gov.br/)):
  - Annual population estimates by municipality, age, and sex for Brazil [@datasusb], used to calculate SISVAN's population coverage.

The [DATASUS](https://datasus.saude.gov.br/) population estimates used in this pipeline are processed through a separate reproducible workflow, available [here](https://sustentarea.github.io/population-estimates/) [@vartanian2025b].

### Data Munging

The data munging follow the data science workflow outlined by @wickham2023e, as illustrated in [@fig-wickham-at-al-2023-figure-1]. All processes were made using
the [Quarto](https://quarto.org/) publishing system, along with the [AWK](https://en.wikipedia.org/wiki/AWK) [@aho2023] and [R](https://www.r-project.org/) [@rcoreteama] programming languages, supported by several R packages.

For data manipulation and workflow, priority was given to packages from the [tidyverse](https://www.tidyverse.org/), [rOpenSci](https://ropensci.org/) and [r-spatial](https://r-spatial.org/) ecosystems, as well as other packages adhering to the tidy tools manifesto [@wickham2023c].

::: {#fig-wickham-at-al-2023-figure-1}
![](images/wickham-at-al-2023-figure-1.png){width=75%}

[Source: Reproduced from @wickham2023e.]{.legend}

Data science workflow created by Wickham, Çetinkaya-Runde, and Grolemund.
:::

### Data Validation

::: {.callout-warning}
The validation steps described below are specifically designed for children aged 0–5 years. If you are working with older children or adolescents (ages 5–19 years), you should adapt the code accordingly. For these age groups, we recommend using the WHO's [`anthroplus`](https://CRAN.R-project.org/package=anthroplus) R package [@dirkschumachera].
:::

Different validation techniques were used to ensure data quality and reliability:

- Duplicate records were removed based on unique combinations of the SISVAN identifier (`id`) and assessment date (`date`). Only the latest record for each individual on a given date was retained.
- Weight and height measurements identified as biologically implausible values (BIVs) according to World Health Organization (WHO) child growth standards [@who2006; @who2008] were set to missing. BIVs were detected by calculating z-scores using the [`anthro_zscores`](https://rdrr.io/cran/anthro/man/anthro_zscores.html) function from the WHO [`anthro`](https://CRAN.R-project.org/package=anthro) R package [@dirkschumacher], based on weight, height, age, and sex. Implausible values were flagged when z-scores exceeded established WHO cutoffs (typically $|z| > 5$). For details, see the [function documentation](https://rdrr.io/cran/anthro/man/anthro_zscores.html).

<!-- @silva2023a quality indicators were also used for validation. Refer to the article for more details. -->

### Data Categorization

Nutritional status categories are ideally determined using z-scores, as recommended by the WHO child growth standards [@who2006, Section C]. However, SISVAN data report age only in years, rather than in days or months as required for accurate z-score calculation. This limitation introduces substantial classification error if z-scores are computed directly. Therefore, we use the nutritional status categories already provided in the SISVAN microdata and set these categories to missing when biologically implausible values (BIVs) were identified.

### Code Style

The Tidyverse [Tidy Tools Manifesto](https://tidyverse.tidyverse.org/articles/manifesto.html) [@wickham2023c], [code style guide](https://style.tidyverse.org/) [@wickhamb] and [design principles](https://design.tidyverse.org/) [@wickhamc] were followed to ensure consistency and enhance readability.

### Reproducibility

The pipeline is fully reproducible and can be run again at any time. To ensure consistent results, the [`renv`](https://rstudio.github.io/renv/) package [@usheya] is used to manage and restore the R environment. See the [README](https://github.com/sustentarea/nutritional-status/blob/main/README.md) file in the code repository to learn how to run it.

## Set Environment

### Load Packages

```{r}
#| label: Set the Environment
#| output: false

library(anthro)
library(brandr)
library(cli)
library(dplyr)
library(forcats)
library(fs)
library(geobr)
library(ggplot2)
library(ggspatial)
library(groomr) # github.com/danielvartan/groomr
library(here)
library(htmltools)
library(httr2)
library(janitor)
library(knitr)
library(labelled)
library(lubridate)
library(nanoparquet)
library(orbis) # github.com/danielvartan/orbis
library(osfr)
library(pal) # gitlab.com/rpkg.dev/pal
library(parallel)
library(quartabs)
library(readr)
library(rutils) # github.com/danielvartan/rutils
library(scales)
library(sf)
library(stringr)
library(tidyr)
library(utils)
library(vroom)
library(zip)
```

### Set Data Directories

```{r}
raw_data_dir <- here("data-raw")
data_dir <- here("data")
```

```{r}
for (i in c(raw_data_dir, data_dir)) {
  if (!dir_exists(i)) dir_create(i, recurse = TRUE)
}
```

### Set Initial Variables

::: {.callout-note}
The `year` variable represent the year of the consolidated [SISVAN](https://sisaps.saude.gov.br/sisvan/) dataset on nutritional status.
:::

```{r}
year <- 2023
```

::: {.callout-note}
The `age_limits` variable define the age range (in years) of individuals to be included in the analysis.
:::

```{r}
age_limits <- c(0, 4) # == Less than 5 years
```

::: {.callout-note}
The `col_selection` variable specifies the columns to be imported from the raw [SISVAN](https://sisaps.saude.gov.br/sisvan/) microdata files.

Click [here](https://s3.sa-east-1.amazonaws.com/ckan.saude.gov.br/SISVAN/estado_nutricional/Dicion%C3%A1rio+de+Dados+-+Estado+Nutricional.pdf) to access the microdata data dictionary (in Portuguese).
:::

```{r}
col_selection <- c(
  "CO_PESSOA_SISVAN",
  "DT_ACOMPANHAMENTO",
  "CO_MUNICIPIO_IBGE",
  "CO_CNES",
  "SG_SEXO",
  "NU_IDADE_ANO",
  "CO_RACA_COR",
  "NU_PESO",
  "NU_ALTURA",
  "PESO X IDADE",
  "PESO X ALTURA",
  "CRI. ALTURA X IDADE",
  "CRI. IMC X IDADE"
)
```

## Download and Import IBGE Municipalities Data

::: {.callout-note}
See the [Source of Data](#source-of-data) section for more information.
:::

```{r}
#| label: Download and Import IBGE Municipalities Data

municipalities_data <- brazil_municipality(year = year)
```

```{r}
municipalities_data |> glimpse()
```

## Download DATASUS Population Estimates

::: {.callout-note}
See the [Source of Data](#source-of-data) section for more information.
:::

### List Files

```{r}
#| label: Download DATASUS Population Estimates

datasus_file_pattern <-
  "datasus-population-estimates-" |>
  paste0(year)
```

```{r}
datasus_file <-
  raw_data_dir |>
  here(paste0(datasus_file_pattern, ".rds"))
```

```{r}
osf_raw_data_id <- "h3pyd"
```

```{r}
osf_raw_data_file <-
  osf_raw_data_id |>
  osf_retrieve_node() |>
  osf_ls_files(
    type = "file",
    pattern = paste0(year, ".rds")
  ) |>
  filter(str_detect(name, paste0("^", year, "\\.rds$")))
```

```{r}
osf_raw_data_file
```

### Download Data

```{r}
osf_raw_data_file |>
  osf_download(
    path = raw_data_dir,
    conflicts = "overwrite"
  ) |>
  pull(local_path)
```

### Rename File

```{r}
if (file_exists(datasus_file)) {
  datasus_file |> file_delete()
}
```

```{r}
raw_data_dir |>
  dir_ls(
    type = "file",
    regexp = paste0(year, "\\.rds$")
  ) |>
  file_move(datasus_file)
```

## Import DATASUS Population Estimates

```{r}
#| label: Import DATASUS Population Estimates

population_estimates_data <- datasus_file |> read_rds()
```

```{r}
population_estimates_data |> glimpse()
```

## Download SISVAN Microdata on Nutritional Status

::: {.callout-note}
See the [Source of Data](#source-of-data) section for more information.
:::

::: {.callout-note}
The microdata files are very large. For practical reasons, some code chunks have `eval: false` set to prevent downloading the data each time the report is rendered. When running the pipeline in a loop or for full automation, remove these lines to enable automatic downloading.
:::

### Download Data

```{r}
#| label: Download SISVAN Microdata on Nutritional Status

file <-
  "sisvan_estado_nutricional_" |>
  paste0(year, ".zip")
```

```{r}
#| eval: false

"https://s3.sa-east-1.amazonaws.com/ckan.saude.gov.br" |>
  path(
    "SISVAN",
    "estado_nutricional",
    file
  ) |>
  request() |>
  req_progress() |>
  req_perform(here(raw_data_dir, file))
```

### Unzip Data

```{r}
#| eval: false

here(raw_data_dir, file) |>
  unzip(exdir = raw_data_dir)
```

### Delete Zip Files

```{r}
#| eval: false

raw_data_dir |>
  dir_ls(type = "file", regexp = "\\.zip$") |>
  file_delete()
```

### Check Data Dimensions

```{r}
file <- file |> str_replace("\\.zip$", "\\.csv")
```

```{r}
raw_data_dir |>
  here(file) |>
  peek_csv_file(
    delim = ";",
    skip = 0,
    has_header = TRUE
  )
```

## Import and Filter Data

::: {.callout-note}
The [`vroom`](https://vroom.r-lib.org/) R package together with the [AWK](https://en.wikipedia.org/wiki/AWK) programming language were use to efficiently handle large datasets and mitigate memory issues. This approach allows the pipeline to run locally on most machines, though we recommend a minimum of 12 GB of RAM for optimal performance. Alternatively, the pipeline can also be executed on cloud platforms such as [Google Colab](https://colab.research.google.com/) and [RStudio Cloud](https://rstudio.cloud/), or using [GitHub Actions](https://docs.github.com/en/actions/concepts/runners/larger-runners) large runners.
:::

### Define Column Names and Schema

```{r}
#| label: Import and Filter Data

col_names <- c(
  "CO_ACOMPANHAMENTO",
  "CO_PESSOA_SISVAN",
  "ST_PARTICIPA_ANDI",
  "CO_MUNICIPIO_IBGE",
  "SG_UF",
  "NO_MUNICIPIO",
  "CO_CNES",
  "NU_IDADE_ANO",
  "NU_FASE_VIDA",
  "DS_FASE_VIDA",
  "SG_SEXO",
  "CO_RACA_COR",
  "DS_RACA_COR",
  "CO_POVO_COMUNIDADE",
  "DS_POVO_COMUNIDADE",
  "CO_ESCOLARIDADE",
  "DS_ESCOLARIDADE",
  "DT_ACOMPANHAMENTO",
  "NU_COMPETENCIA",
  "NU_PESO",
  "NU_ALTURA",
  "DS_IMC",
  "DS_IMC_PRE_GESTACIONAL",
  "PESO X IDADE",
  "PESO X ALTURA",
  "CRI. ALTURA X IDADE",
  "CRI. IMC X IDADE",
  "ADO. ALTURA X IDADE",
  "ADO. IMC X IDADE",
  "CO_ESTADO_NUTRI_ADULTO",
  "CO_ESTADO_NUTRI_IDOSO",
  "CO_ESTADO_NUTRI_IMC_SEMGEST",
  "CO_SISTEMA_ORIGEM_ACOMP",
  "SISTEMA_ORIGEM_ACOMP"
)
```

```{r}
#| code-fold: false

schema <- cols(
  "CO_ACOMPANHAMENTO" = col_character(),
  "CO_PESSOA_SISVAN" = col_character(),
  "ST_PARTICIPA_ANDI" = col_character(),
  "CO_MUNICIPIO_IBGE" = col_integer(),
  "SG_UF" = col_factor(),
  "NO_MUNICIPIO" = col_character(),
  "CO_CNES" = col_integer(),
  "NU_IDADE_ANO" = col_integer(),
  "NU_FASE_VIDA" = col_character(), # decimal mark = "." (double)
  "DS_FASE_VIDA" = col_factor(),
  "SG_SEXO" = col_factor(),
  "CO_RACA_COR" = col_character(),
  "DS_RACA_COR" = col_factor(),
  "CO_POVO_COMUNIDADE" = col_integer(),
  "DS_POVO_COMUNIDADE" = col_factor(),
  "CO_ESCOLARIDADE" = col_character(),
  "DS_ESCOLARIDADE" = col_factor(),
  "DT_ACOMPANHAMENTO" = col_date(),
  "NU_COMPETENCIA" = col_integer(),
  "NU_PESO" = col_double(),
  "NU_ALTURA" = col_integer(),
  "DS_IMC" = col_double(),
  "DS_IMC_PRE_GESTACIONAL" = col_character(), # decimal mark = "." (double)
  "PESO X IDADE" = col_factor(),
  "PESO X ALTURA" = col_factor(),
  "CRI. ALTURA X IDADE" = col_factor(),
  "CRI. IMC X IDADE" = col_factor(),
  "ADO. ALTURA X IDADE" = col_factor(),
  "ADO. IMC X IDADE" = col_factor(),
  "CO_ESTADO_NUTRI_ADULTO" = col_factor(),
  "CO_ESTADO_NUTRI_IDOSO" = col_factor(),
  "CO_ESTADO_NUTRI_IMC_SEMGEST" = col_factor(),
  "CO_SISTEMA_ORIGEM_ACOMP" = col_integer(),
  "SISTEMA_ORIGEM_ACOMP" = col_factor()
)
```

### Import and Filter Data

::: {.callout-important}
You may see warning messages about failed parsing. These warnings are expected due to minor inconsistencies in the SISVAN raw data and do not affect the overall analysis.
:::

```{r}
#| warning: false
#| code-fold: false

data <-
  vroom(
    file = pipe(
      paste0(
        "awk ",
        "-F ", # Field separator
        "';' ",
        "'{", # Program
        "if (",
        "($8 >= ",
        age_limits[1],
        ")",
        " && ",
        "($8 <= ",
        age_limits[2],
        ")",
        ") ",
        "{print}",
        "}' ",
        raw_data_dir |> here(file) # file
      )
    ),
    delim = ";",
    col_names = col_names,
    col_types = schema,
    col_select = all_of(col_selection),
    na = c("", "NA"),
    locale = locale(
      date_names = "pt",
      date_format = "%d/%m/%Y",
      time_format = "%H:%M:%S",
      decimal_mark = ",",
      grouping_mark = ".",
      tz = "America/Sao_Paulo",
      encoding = raw_data_dir |>
        here(file) |>
        guess_encoding() |>
        extract2("encoding") |>
        magrittr::extract(1)
    ),
    guess_max = 100,
    num_threads = detectCores() |>
      multiply_by(0.75) |>
      floor(),
    progress = TRUE
  )
```

```{r}
#| warning: false

data |> glimpse()
```

## Tidy Data

### Rename Columns

```{r}
#| label: Tidy Data

data <-
  data |>
  clean_names() |>
  rename(
    id = co_pessoa_sisvan,
    date = dt_acompanhamento,
    municipality_code = co_municipio_ibge,
    cnes = co_cnes,
    sex = sg_sexo,
    age = nu_idade_ano,
    ethnicity = co_raca_cor,
    weight = nu_peso,
    height = nu_altura,
    weight_for_age = peso_x_idade,
    weight_for_height = peso_x_altura,
    height_for_age = cri_altura_x_idade,
    bmi_for_age = cri_imc_x_idade
  )
```

```{r}
#| code-fold: false

data |> glimpse()
```

### Standardize Columns

```{r}
data <-
  data |>
  mutate(
    sex = sex |>
      as.character() |>
      case_match(
        "F" ~ "Female",
        "M" ~ "Male"
      ) |>
      factor(
        levels = c("Male", "Female"),
        ordered = FALSE
      ),
    ethnicity = ethnicity |>
      as.character() |>
      case_match(
        "01" ~ "White",
        "02" ~ "Black",
        "03" ~ "Yellow",
        "04" ~ "Brown",
        "05" ~ "Indigenous"
      ) |>
      factor(
        levels = c(
          "White",
          "Black",
          "Yellow",
          "Brown",
          "Indigenous"
        ),
        ordered = FALSE
      ),
    weight_for_age = weight_for_age |>
      as.character() |>
      case_match(
        "Muito baixo peso para a idade" ~ "Severely underweight",
        "Baixo peso para a idade" ~ "Underweight",
        "Peso adequado para idade" ~ "Normal",
        "Peso elevado para a idade" ~ "High"
      ) |>
      factor(
        levels = c(
          "Severely underweight",
          "Underweight",
          "Normal",
          "High"
        ),
        ordered = TRUE
      ),
    weight_for_height = weight_for_height |>
      as.character() |>
      case_match(
        "Magreza acentuada" ~ "Severe wasted",
        "Magreza" ~ "Wasted",
        "Peso Adequado ou Eutrofico" ~ "Normal",
        "Risco de sobrepeso" ~ "Possible risk of overweight",
        "Sobrepeso" ~ "Overweight",
        "Obesidade" ~ "Obese"
      ) |>
      factor(
        levels = c(
          "Severe wasted",
          "Wasted",
          "Normal",
          "Possible risk of overweight",
          "Overweight",
          "Obese"
        ),
        ordered = TRUE
      ),
    height_for_age = height_for_age |>
      as.character() |>
      case_match(
        "Muito baixa estatura para idade" ~ "Severely stunted",
        "Baixa estatura para idade" ~ "Stunted",
        "Estatura adequada para a idade" ~ "Normal"
      ) |>
      factor(
        levels = c(
          "Severely stunted",
          "Stunted",
          "Normal"
        ),
        ordered = TRUE
      ),
    bmi_for_age = bmi_for_age |>
      as.character() |>
      case_match(
        "Magreza acentuada" ~ "Severe wasted",
        "Magreza" ~ "Wasted",
        "Eutrofia" ~ "Normal",
        "Risco de sobrepeso" ~ "Possible risk of overweight",
        "Sobrepeso" ~ "Overweight",
        "Obesidade" ~ "Obese"
      ) |>
      factor(
        levels = c(
          "Severe wasted",
          "Wasted",
          "Normal",
          "Possible risk of overweight",
          "Overweight",
          "Obese"
        ),
        ordered = TRUE
      )
  )
```

```{r}
#| code-fold: false

data |> glimpse()
```

## Transform Data

### Remove Duplicates

```{r}
#| label: Transform Data
#| warning: false

data <-
  data |>
  arrange(desc(date)) |>
  distinct(id, date, .keep_all = TRUE)
```

```{r}
#| warning: false

data |> glimpse()
```

### Remove Biological Implausible Values (BVI)

::: {.callout-note}
See the [Data Validation](#data-validation) section for more information.
:::

```{r}
data <-
  data |>
  mutate(
    z_scores = anthro_zscores(
      sex = as.numeric(sex),
      age = age * 12,
      is_age_in_month = TRUE,
      weight = weight,
      lenhei = height,
      measure = "h"
    ),
    weight = if_else(
      (z_scores$fwei == 1) | (z_scores$flen != 1 & z_scores$fwfl == 1),
      NA,
      weight
    ),
    height = if_else(
      z_scores$flen == 1,
      NA,
      height
    ),
    weight_for_age = if_else(
      is.na(weight),
      NA,
      weight_for_age
    ),
    weight_for_height = if_else(
      is.na(weight) | is.na(height),
      NA,
      weight_for_height
    ),
    height_for_age = if_else(
      is.na(height),
      NA,
      height_for_age
    ),
    bmi_for_age = if_else(
      is.na(weight) | is.na(height),
      NA,
      bmi_for_age
    )
  ) |>
  select(-z_scores)
```

```{r}
data |> glimpse()
```

### Fix Municipality Code

```{r}
data <-
  data |>
  rename(municipality_code_6 = municipality_code) |>
  left_join(
    municipalities_data |>
      mutate(
        municipality_code_6 = municipality_code |>
          str_sub(1, 6) |>
          as.integer()
      ) |>
      select(municipality_code, municipality_code_6),
    by = join_by(municipality_code_6)
  ) |>
  select(-municipality_code_6) |>
  relocate(municipality_code, .after = date)
```

```{r}
data |> glimpse()
```

### Arrange Data

```{r}
data <-
  data |>
  arrange(
    date,
    municipality_code,
    cnes,
    sex,
    age
  )
```

```{r}
data |> glimpse()
```

## Create Data Dictionary

### Prepare Metadata

```{r}
#| label: Create Data Dictionary

metadata <-
  data |>
  `var_label<-`(
    list(
      id = "Unique identifier for the individual",
      date = "Date of the individual's nutritional assessment",
      municipality_code = paste0(
        "Institute of Geography and Statistics (IBGE) code of the ",
        "municipality where the assessment was performed"
      ),
      cnes = paste0(
        "National Registry of Health Establishments (CNES) code of the ",
        "health facility where the assessment was performed"
      ),
      sex = "Sex of the individual",
      age = "Age of the individual in years",
      ethnicity = "Self-reported ethnicity/race or color of the individual",
      weight = "Weight of the individual in kilograms",
      height = "Height of the individual in centimeters",
      weight_for_age = paste0(
        "Nutritional status classification (children 0–5) based on ",
        "weight-for-age"
      ),
      weight_for_height = paste0(
        "Nutritional status classification (children 0–5) based on ",
        "weight-for-height"
      ),
      height_for_age = paste0(
        "Nutritional status classification (children 0–10) based on ",
        "height-for-age"
      ),
      bmi_for_age = paste0(
        "Nutritional status classification (children 0–10) based on ",
        "BMI-for-age"
      )
    )
  ) |>
  generate_dictionary(details = "full") |>
  convert_list_columns_to_character()
```

### Visualize Final Data

```{r}
metadata |> glimpse()
```

```{r}
metadata
```

```{r}
data |> glimpse()
```

```{r}
data
```

## Save Data

::: {.callout-note}
The processed data are available in [`csv`](https://en.wikipedia.org/wiki/Comma-separated_values), [`rds`](https://rdrr.io/r/base/readRDS.html) and [`parquet`](https://en.wikipedia.org/wiki/Apache_Parquet) formats through a dedicated repository on the Open Science Framework ([OSF](https://doi.org/10.17605/OSF.IO/8J94M)). See the [Data Availability](#data-availability) section for more information.
:::

### Write Data

```{r}
#| label: Write Data

valid_file_pattern <-
  year |>
  paste0(
    "-age-limits-",
    age_limits[1],
    "-",
    age_limits[2]
  )
```

```{r}
data |>
  write_csv(
    here(data_dir, paste0(valid_file_pattern, ".csv"))
  )
```

```{r}
data |>
  write_rds(
    here(data_dir, paste0(valid_file_pattern, ".rds"))
  )
```

```{r}
data |>
  write_parquet(
    here(data_dir, paste0(valid_file_pattern, ".parquet"))
  )
```

### Write Metadata

```{r}
metadata_file_pattern <-
  "metadata-" |>
  paste0(
    year,
    "-age-limits-",
    age_limits[1],
    "-",
    age_limits[2]
  )
```

```{r}
metadata |>
  write_csv(
    here(data_dir, paste0(metadata_file_pattern, ".csv"))
  )
```

```{r}
metadata |>
  write_rds(
    here(data_dir, paste0(metadata_file_pattern, ".rds"))
  )
```

```{r}
metadata |>
  write_parquet(
    here(data_dir, paste0(metadata_file_pattern, ".parquet"))
  )
```

<!-- ## Upload the Data to OSF (Optional) -->

<!-- Only repository administrators can upload data to OSF. -->
<!-- To enable uploads, set a `OSF_PAT` environment variable with a valid OSF personal access token. -->

```{r}
#| eval: false
#| include: false

Sys.getenv("OSF_PAT") |> osf_auth()
```

```{r}
#| eval: false
#| include: false

osf_data_id <- "q7m9d"
```

```{r}
#| eval: false
#| include: false

osf_data_id |>
  osf_retrieve_node() |>
  osf_upload(
    path = c(
      here(data_dir, paste0(valid_file_pattern, ".csv")),
      here(data_dir, paste0(valid_file_pattern, ".rds")),
      here(data_dir, paste0(valid_file_pattern, ".parquet")),
      here(data_dir, paste0(metadata_file_pattern, ".csv")),
      here(data_dir, paste0(metadata_file_pattern, ".rds")),
      here(data_dir, paste0(metadata_file_pattern, ".parquet"))
    ),
    conflicts = "overwrite"
  )
```

## Explore Data

### Summarize Frequencies

::: {.callout-important}
Some SISVAN data may show discrepancies when compared to official population estimates. For example, the reported number of children under 5 classified as *yellow* in 2023 appears unusually high. If you notice such inconsistencies, check the [SISVAN web reports](https://sisaps.saude.gov.br/sisvan/relatoriopublico/) to confirm before reporting an issue.
:::

```{r}
#| label: Explore Data
#| code-fold: true

vars <- c(
  "sex",
  "age",
  "ethnicity",
  "weight_for_age",
  "weight_for_height",
  "height_for_age",
  "bmi_for_age"
)
```

```{r}
#| code-fold: true

panel_tabset_data <- tibble()

for (i in vars) {
  table <-
    data |>
    arrange(desc(.data[[i]])) |>
    distinct(id, .data[[i]]) |>
    group_by(.data[[i]]) |>
    summarize(n = n(), .groups = "drop") |>
    arrange(desc(n)) |>
    mutate(
      !!i := .data[[i]] |>
        as.factor() |>
        fct_na_value_to_level(level = "(NA)") |>
        fct_reorder(n) |>
        fct_relevel("(NA)", after = 0)
    ) |>
    arrange(desc(.data[[i]])) |>
    mutate(
      n_cum = cumsum(n),
      pct = n |>
        divide_by(sum(n)) |>
        multiply_by(100),
      pct_cum = pct |>
        cumsum() |>
        round(3),
      pct = pct |> round(3),
      across(
        .cols = where(is.numeric),
        .fns = \(x) prettyNum(x, big.mark = ",", decimal.mark = ".")
      )
    ) |>
    rename(
      N = n,
      `N (Cumulative)` = n_cum,
      Percent = pct,
      `Percent (Cumulative)` = pct_cum
    ) |>
    kable()

  panel_tabset_data <-
    panel_tabset_data |>
    bind_rows(
      tibble(
        label = paste0("`", i, "`"),
        table = list(table)
      )
    )
}
```

```{r}
#| code-fold: true
#| output: asis
#| messages: false
#| warning: false

panel_tabset_data |> render_tabset(label, table)
```

### Plot Bar Charts

::: {.callout-important}
Some SISVAN data may show discrepancies when compared to official population estimates. For example, the reported number of children under 5 classified as *yellow* in 2023 appears unusually high. If you notice such inconsistencies, check the [SISVAN web reports](https://sisaps.saude.gov.br/sisvan/relatoriopublico/) to confirm before reporting an issue.
:::

```{r}
#| code-fold: true

panel_tabset_data <- tibble()

for (i in vars) {
  plot <-
    data |>
    mutate(
      !!i := as.character(.data[[i]]),
      !!i := ifelse(
        str_length(.data[[i]]) > 30,
        paste0(str_sub(.data[[i]], 1, 27), "..."),
        .data[[i]]
      ) |>
        fct_na_value_to_level(level = "(NA)")
    ) |>
    group_by(.data[[i]]) |>
    summarize(n = n(), .groups = "drop") |>
    mutate(
      rel = n |>
        divide_by(sum(n)) |>
        round(2),
      !!i := .data[[i]] |>
        fct_reorder(n) |>
        fct_relevel("(NA)", after = 0)
    ) |>
    ggplot(
      aes(
        x = .data[[i]],
        y = n,
        labels = percent(rel)
      )
    ) +
    geom_col(fill = get_brand_color("green")) +
    geom_text(
      hjust = -0.25,
      size = 3
    ) +
    coord_flip() +
    scale_y_continuous(
      expand = expansion(mult = c(0.05, 0.15))
    ) +
    labs(x = NULL, y = NULL) +
    theme(
      axis.text.y = element_text(size = 8)
    )

  panel_tabset_data <-
    panel_tabset_data |>
    bind_rows(
      tibble(
        label = paste0("`", i, "`"),
        plot = list(plot)
      )
    )
}
```

```{r}
#| code-fold: true
#| output: asis
#| messages: false
#| warning: false

panel_tabset_data |> render_tabset(label, plot)
```

### Summarize Descriptive Statistics

```{r}
#| code-fold: true

vars <- c(
  "date",
  "weight",
  "height"
)
```

```{r}
#| code-fold: true

panel_tabset_data <- tibble()

for (i in vars) {
  table <-
    data |>
    arrange(desc(.data[[i]])) |>
    distinct(id, .data[[i]]) |>
    stats_summary(i) |>
    mutate(
      name = c(
        "Class",
        "N",
        "N (Without Missing)",
        "N (Missing)",
        "Mean",
        "Variance",
        "Standard Deviation",
        "Minimum",
        "1st Quartile (Q1)",
        "Median",
        "3rd Quartile (Q3)",
        "Maximum",
        "Interquartile Range (IQR)",
        "Range",
        "Skewness",
        "Kurtosis"
      ),
      value = if_else(
        !str_detect(value, "\\d{4}-\\d{2}-\\d{2}"),
        value |>
          prettyNum(big.mark = ",", decimal.mark = ".") |>
          str_trim(),
        value
      )
    ) |>
    rename(
      Name = name,
      Value = value
    ) |>
    kable()

  panel_tabset_data <-
    panel_tabset_data |>
    bind_rows(
      tibble(
        label = paste0("`", i, "`"),
        table = list(table)
      )
    )
}
```

```{r}
#| code-fold: true
#| output: asis
#| messages: false
#| warning: false

panel_tabset_data |> render_tabset(label, table)
```

### Plot Histograms

```{r}
#| code-fold: true
#| warning: false

panel_tabset_data <- tibble()

for (i in vars) {
  plot <-
    data |>
    tidyr::drop_na(all_of(i)) |>
    ggplot(
      aes(
        x = .data[[i]],
        y = after_stat(count)
      )
    ) +
    geom_histogram(
      bins = 30,
      fill = get_brand_color("gray-d25"),
      color = get_brand_color("white")
    ) +
    # geom_density(
    #   color = "red",
    #   linewidth = 1
    # ) +
    labs(
      title = ifelse(
        i == "date",
        paste("Distribution of", str_to_title(i), "of Nutritional Assessments"),
        paste("Distribution of", str_to_title(i))
      ),
      x = str_to_title(i),
      y = "Frequency"
    ) +
    theme(
      axis.text.x = element_text(margin = margin(0, 0, 10, 0)),
      axis.text.y = element_text(margin = margin(0, 0, 0, 20))
    )

  panel_tabset_data <-
    panel_tabset_data |>
    bind_rows(
      tibble(
        label = paste0("`", i, "`"),
        plot = list(plot)
      )
    )
}
```

```{r}
#| code-fold: true
#| output: asis
#| messages: false
#| warning: false

panel_tabset_data |> render_tabset(label, plot)
```

## Check Relative Coverage

### Transform Data

#### Remove Duplicates by Year

::: {.callout-note}
As described in @silva2023a[p. 4], to calculate SISVAN's total resident population coverage, only the most recent record for each individual within each year is retained for analysis.
:::

```{r}
#| code-fold: true

data <-
  data |>
  mutate(year = year(date)) |>
  arrange(desc(date)) |>
  distinct(id, year, .keep_all = TRUE) |>
  relocate(year, .after = date)
```

```{r}
#| code-fold: true

data |> glimpse()
```

#### Summarize Data by Year

```{r}
#| code-fold: true

data <-
  data |>
  summarize(
    coverage = n(),
    mean_age = age |> mean(na.rm = TRUE),
    mean_weight = weight |> mean(na.rm = TRUE),
    mean_height = height |> mean(na.rm = TRUE),
    .by = c(
      "year",
      "municipality_code"
    )
  )
```

```{r}
#| code-fold: true

data |> glimpse()
```

#### Add Population Estimates

```{r}
#| code-fold: true

data <-
  population_estimates_data |>
  filter(between(age, age_limits[1], age_limits[2])) |>
  summarize(
    n = population |> sum(na.rm = TRUE),
    .by = c(
      "year",
      "municipality_code"
    )
  ) |>
  right_join(
    data,
    by = c(
      "year",
      "municipality_code"
    )
  ) |>
  rename(population = n) |>
  relocate(population, .before = coverage)
```

```{r}
#| code-fold: true

data |> glimpse()
```

#### Add Municipality Data

```{r}
#| code-fold: true

data <-
  data |>
  left_join(
    municipalities_data,
    by = join_by(municipality_code),
    suffix = c(".x", "")
  ) |>
  select(-latitude, -longitude) |>
  relocate(region_code:municipality, .after = year) |>
  relocate(municipality, .after = municipality_code)
```

```{r}
#| code-fold: true

data |> glimpse()
```

### Validate Data

::: {.callout-note}
The population value used here is an estimate. If the SISVAN coverage for a municipality exceeds the estimated population, the population value is adjusted to match the coverage.
:::

```{r}
#| code-fold: true

data <-
  data |>
  mutate(
    population = case_when(
      coverage > population ~ coverage,
      TRUE ~ population
    )
  )
```

```{r}
#| code-fold: true

data |> glimpse()
```

### Calculate Relative Coverage

```{r}
#| code-fold: true

data <-
  data |>
  mutate(
    coverage_pct = coverage |>
      divide_by(population) |>
      multiply_by(100)
  ) |>
  relocate(coverage_pct, .after = coverage)
```

```{r}
#| code-fold: true

data |> glimpse()
```

### Arrange Data

```{r}
#| code-fold: true

data <-
  data |>
  arrange(
    year,
    municipality_code
  )
```

```{r}
#| code-fold: true

data |> glimpse()
```

### Tabulate Relative Coverage by Region

```{r}
#| output: asis
#| code-fold: true

data |>
  mutate(
    region = brazil_region(region_code)
  ) |>
  summarize(
    population = population |> sum(na.rm = TRUE),
    coverage = coverage |> sum(na.rm = TRUE),
    .by = "region"
  ) |>
  slice(c(1, 2, 5, 3, 4)) |>
  mutate(
    coverage_pct = coverage |>
      divide_by(population) |>
      multiply_by(100) |>
      round(3),
    across(
      .cols = where(is.numeric),
      .fns = \(x) prettyNum(x, big.mark = ",", decimal.mark = ".")
    )
  ) |>
  arrange(desc(coverage_pct)) |>
  rename(
    Region = region,
    Population = population,
    `SISVAN Coverage` = coverage,
    `SISVAN Coverage (%)` = coverage_pct
  ) |>
  pipe_table() |>
  cat_lines()
```

### Tabulate Relative Coverage by State

```{r}
#| output: asis
#| code-fold: true

data |>
  mutate(
    state = brazil_state(state_code)
  ) |>
  summarize(
    population = population |> sum(na.rm = TRUE),
    coverage = coverage |> sum(na.rm = TRUE),
    .by = "state"
  ) |>
  arrange(state) |>
  mutate(
    coverage_pct = coverage |>
      divide_by(population) |>
      multiply_by(100) |>
      round(3),
    across(
      .cols = where(is.numeric),
      .fns = \(x) prettyNum(x, big.mark = ",", decimal.mark = ".")
    )
  ) |>
  arrange(desc(coverage_pct)) |>
  rename(
    State = state,
    Population = population,
    `SISVAN Coverage` = coverage,
    `SISVAN Coverage (%)` = coverage_pct
  ) |>
  pipe_table() |>
  cat_lines()
```

### Plot Histogram by Municipality

```{r}
#| code-fold: true

data |>
  tidyr::drop_na(coverage_pct) |>
  ggplot(aes(x = coverage_pct)) +
  geom_histogram(
    aes(y = after_stat(density)),
    bins = 30,
    fill = get_brand_color("gray-d25"),
    color = get_brand_color("white")
  ) +
  geom_density(
    color = "red",
    linewidth = 1
  ) +
  labs(
    title = "SISVAN Coverage by Municipality (%) (Ages 0-5)",
    subtitle = paste0("Year: ", year),
    x = "Coverage (%)",
    y = "Density",
    caption = "Source: SISVAN."
  )
```

### Plot Map by Municipality

#### Set Shape

```{r}
#| code-fold: true

shape <-
  read_municipality(
    year = year |>
      closest_geobr_year(type = "municipality"),
    showProgress = FALSE
  ) |>
  st_transform(st_crs(4326))
```

#### Prepare Plot Data

```{r}
#| code-fold: true

plot_data <-
  data |>
  left_join(
    shape,
    by = join_by(municipality_code == code_muni)
  ) |>
  rename(geometry = geom) |>
  select(
    municipality_code,
    coverage_pct,
    geometry
  ) |>
  tidyr::drop_na(coverage_pct)
```

#### Plot Data

```{r}
#| code-fold: true

brand_div_palette <- function(x) {
  brandr:::make_color_ramp(
    n = x,
    colors = c(
      get_brand_color("dark-red"),
      # get_brand_color("white"),
      get_brand_color_mix(
        position = 950,
        color_1 = "dark-red",
        color_2 = "dark-red-triadic-blue",
        alpha = 0.5
      ),
      get_brand_color("dark-red-triadic-blue")
    )
  )
}
```

```{r}
#| code-fold: true

plot_data |>
  st_as_sf() |>
  ggplot(aes(fill = coverage_pct)) +
  geom_sf(
    color = get_brand_color("gray"),
    linewidth = 0.05
  ) +
  scale_fill_binned(
    breaks = seq(0, 100, 25),
    limits = c(0, 100),
    palette = brand_div_palette,
    na.value = get_brand_color("gray-d25")
  ) +
  annotation_scale(
    aes(),
    location = "br",
    style = "tick",
    height = unit(0.5, "lines")
  ) +
  annotation_north_arrow(
    location = "br",
    height = unit(2, "lines"),
    width = unit(2, "lines"),
    pad_x = unit(0.25, "lines"),
    pad_y = unit(1.25, "lines"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(
    title = "SISVAN Coverage by Municipality (%) (Ages 0-5)",
    subtitle = paste0("Year: ", year),
    fill = NULL,
    caption = "Source: SISVAN."
  )
```

## Citation

::: {.callout-important}
When using this data, you must also cite the original data sources.
:::

To cite this work, please use the following format:

Vartanian, D., Schettino, J. P. J., & Carvalho, A. M. (2025). *A reproducible pipeline for processing and analyzing SISVAN microdata on nutritional status monitoring in Brazil* \[Computer software\]. Sustentarea Research and Extension Group, University of São Paulo. <https://sustentarea.github.io/nutritional-status>

A BibLaTeX entry for LaTeX users is:

```
@software{vartanian2025,
  title = {A reproducible pipeline for processing and analyzing SISVAN microdata on nutritional status monitoring in Brazil},
  author = {{Daniel Vartanian} and {João Pedro Junqueira Schettino} and {Aline Martins de Carvalho}},
  year = {2025},
  address = {São Paulo},
  institution = {Sustentarea Research and Extension Group, University of São Paulo},
  langid = {en},
  url = {https://sustentarea.github.io/nutritional-status}
}
```

## License

[![License: GPLv3](https://img.shields.io/badge/license-GPLv3-bd0000.svg)](https://www.gnu.org/licenses/gpl-3.0)
[![License: CC BY-NC-SA 4.0](https://img.shields.io/badge/license-CC_BY--NC--SA_4.0-lightgrey.svg)](https://creativecommons.org/licenses/by-nc-sa/4.0/)

::: {.callout-important}
The original data sources may be subject to their own licensing terms and conditions.
:::

The code in this repository is licensed under the [GNU General Public License Version 3](https://www.gnu.org/licenses/gpl-3.0), while the report is available under the [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/).

```
Copyright (C) 2025 Sustentarea Research and Extension Group

The code in this report is free software: you can redistribute it and/or
modify it under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your option)
any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program. If not, see <https://www.gnu.org/licenses/>.
```

## Acknowledgments

```{=html}
<br>
```

```{r, results='asis'}
#| echo: false

blocks <- list(
  list(
    logo_link = "https://www.fsp.usp.br/sustentarea/",
    logo_src = "images/sustentarea-icon.svg",
    logo_alt = "Sustentarea Logo",
    logo_width = 100,
    text = 'This work is part of a research project by the <a href="https://www.fsp.usp.br/sustentarea/">Sustentarea</a> Research and Extension Group of the University of São Paulo (<a href="https://usp.br/">USP</a>) titled: <em>Global syndemic: The impact of anthropogenic climate change on the health and nutrition of children under five years old attended by Brazil\'s public health system</em>.'
  ),
  list(
    logo_link = "https://www.gov.br/cnpq/",
    logo_src = "images/cnpq-logo.svg",
    logo_alt = "CNPq Logo",
    logo_width = 160,
    text = 'This work was supported by the Department of Science and Technology of the Secretariat of Science, Technology, and Innovation and of the Health Economic-Industrial Complex (<a href="https://www.gov.br/saude/pt-br/composicao/sectics/">SECTICS</a>) of the <a href="https://www.gov.br/saude/pt-br/composicao/sectics/">Ministry of Health</a> of Brazil, and the National Council for Scientific and Technological Development (<a href="https://www.gov.br/cnpq/">CNPq</a>) (grant no. 444588/2023-0).'
  )
)

blocks |>
  lapply(
    function(x) {
      div(
        style = paste0(
          "display: flex; ",
          "align-items: flex-start; ",
          "margin-bottom: 2em;"
        ),
        div(
          style = paste0(
            "flex: 0 0 30%; ",
            "display: flex; ",
            "justify-content: center; ",
            "margin: auto 0;"
          ),
          tags$a(
            href = x$logo_link,
            tags$img(
              src = x$logo_src,
              alt = x$logo_alt,
              style = paste0(
                "width: ",
                x$logo_width,
                "px; ",
                "height: auto;"
              )
            )
          )
        ),
        div(
          style = paste0(
            "flex: 1; ",
            "padding-left: 1em;"
          ),
          HTML(x$text)
        )
      )
    }
  ) |>
  tagList() |>
  browsable()
```

## References {.unnumbered}

::: {#refs}
:::