index.qmd

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

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

## Overview

This report explores potential associations between childhood undernutrition and the Standardized Precipitation Evapotranspiration Index ([SPEI](https://en.wikipedia.org/wiki/Standardised_Precipitation_Evapotranspiration_Index)) in Brazilian municipalities (2008–2019). It is part of the [Sustentarea](https://www.fsp.usp.br/sustentarea) Research and Extension Group's project titled *Global syndemic: The impact of anthropogenic climate change on the health and nutrition of children under five years old served by Brazil's public health system (SUS)*.

::: {.callout-warning}
Please note that this report is designed to support decision-making and may not include all the details of the analysis.

Any additional analyses or updates will be incorporated into the report in future revisions, if necessary.
:::

## Question

This analysis seeks to address the following question:

::: {style="text-align: center; font-size: 1.1em; font-weight: bold; color: #ED6B4D; padding-top: 1em; padding-bottom: 1em;"}
Does the Standardized Precipitation Evapotranspiration Index (SPEI) significantly improve the prediction of childhood undernutrition in Brazilian municipalities?
:::

## Methods

### Approach and Procedure Method

This study employed the hypothetical-deductive method, also known as the method of conjecture and refutation [@popper1979, p. 164], as its problem-solving approach. Procedurally, it applied an enhanced version of Null Hypothesis Significance Testing (NHST), grounded on the original ideas of Neyman-Pearson framework for data testing [@neyman1928; @neyman1928a; @perezgonzalez2015].

The analysis was conducted using Generalized Additive Models ([GAMs](https://en.wikipedia.org/wiki/Generalized_additive_model)) to model and control for potential effects, alongside visual inspections of the data. Results are also presented for each cluster of the Revised Multidimensional Index for Sustainable Food Systems (MISFS-R) [@carvalho2021a; @norde2023].

### Source of Data/Information

The data used in this analysis have as sources:

- The Brazilian Institute of Geography and Statistics ([IBGE](https://www.ibge.gov.br/)) Automatic Retrieval System ([SIDRA](https://sidra.ibge.gov.br/)), for data on GDP per capita and the Gini Index [@ibgen; @ibgel].
- Brazil's Food and Nutrition Surveillance System ([SISVAN](https://sisaps.saude.gov.br/sisvan/)), for data on malnutrition [@sisvana].
- [WorldClim](https://worldclim.org/data/index.html), for data on bioclimatic variables, which allowed us to calculate the Standardised Precipitation Evapotranspiration Index (SPEI) for municipalities in Brazil [@harris2020; @fick2017].

Some data are imported directly from the source, while others rely on external data files available in the [data](https://github.com/sustentarea/gs-data-analysis-report-3/tree/main/data) directory of the code repository.

### Data Wrangling

Data wrangling and analysis followed the data science framework outlined by @wickham2023e, as illustrated in [@fig-wickham-at-al-2024-figure-1]. All processes were made using the R programming language [@rcoreteama], RStudio IDE [@positteama], and several R packages.

The [tidyverse](https://www.tidyverse.org/) and [rOpenSci](https://ropensci.org/) peer-reviewed package ecosystem and other R packages adherents of the tidy tools manifesto [@wickham2023e] were prioritized. All processes were made in order to provide result reproducibility and to be in accordance with the FAIR principles [@wilkinson2016].

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

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

Model of the data science process created by Wickham, Çetinkaya-Runde, and Grolemund.
:::

The Tidyverse [code style guide](https://style.tidyverse.org/) and [design principles](https://design.tidyverse.org/) were followed to ensure consistency and enhance readability.

All the analyses are 100% reproducible and can be run again at any time. See the [README](https://github.com/sustentarea/gs-data-analysis-report-3/blob/main/README.md) file in the code repository to learn how to run them.

### Model Parameters

The models were built using the `mgcv` R package [@wood]. Since the dependent variables are relative frequencies (continuous), we used the beta distribution family (@fig-beta-distributions) with a logit link function (@eq-logit-link-function)(@fig-logit-funciton)(see @casella2002[p. 591). The [REML](https://en.wikipedia.org/wiki/Restricted_maximum_likelihood) (Restricted Maximum Likelihood) method was used to estimate the smoothing parameters.

$$
\text{logit}(P) = \ln\left(\frac{P}{1 - P}\right) = \beta_{0} + \beta_{1} X_{1} + \cdots + \beta_{k} X_{k}
$$ {#eq-logit-link-function}

::: {#fig-beta-distributions}
```{r}
list <- 
  dplyr::tibble(
    alpha = c(0.5, 5, 1, 2, 2),
    beta = c(0.5, 1, 3, 2, 5),
    color = c(
      "a = b == 0.5",
      "a = 5, b = 1",
      "a = 1, b = 3",
      "a = 2, b = 2",
      "a = 2, b = 5"
    )
  ) %>%
  split(., seq(nrow(.)))

plot <- 
  ggplot2::ggplot(NULL, ggplot2::aes(x = x, color = color)) +
  ggplot2::labs(
    x = "x", 
    y = "Probability Density Function (PDF)",
    color = "Parameters"
  ) +
  ggplot2::scale_y_continuous(limits = c(0, 2.5))+
  scale_color_brand_d()

for (i in list) {
  plot <-
    plot +
    ggplot2::stat_function(
      data = dplyr::tibble(x = 0:1, color = factor(i$color)),
      fun = stats::dbeta,
      args = list(shape1 = i$alpha, shape2 = i$beta),
      n = 1000,
      linewidth = 1.5
    )
}

plot |> print() |> rutils::shush()
```

[Source: Created by the authors.]{.legend}

Examples of beta distributions with different parameters.
:::

::: {#fig-logit-funciton}
```{r}
ggplot2::ggplot() +
  ggplot2::stat_function(
    data = dplyr::tibble(x = 0:1, color = factor(i$color)),
    fun = stats::qlogis,
    args = list(location = 0, scale = 1, log = FALSE),
    n = 1000,
    linewidth = 1.5,
    color = get_brand_color("red")
  ) +
  ggplot2::geom_hline(
    yintercept = 0, 
    linewidth = 0.25,
    linetype = "dashed",
    color = get_brand_color("grey")
  ) +
  ggplot2::geom_vline(
    xintercept = 0.5, 
    linewidth = 0.25,
    linetype = "dashed",
    color = get_brand_color("grey")
  ) +
  ggplot2::lims(
    x = c(0, 1),
    y = c(-6, 6)
  ) +
  ggplot2::labs(
    x = "Probability",
    y = "Logit(0,1)"
  )
```

[Source: Created by the authors.]{.legend}

Example of a logit function (i.e., the quantile function of the logistic distribution).
:::

### Hypothesis Testing

We tested whether SPEI significantly improves model fit when predicting MBEPR & BEIPR (stunting) and MAPER & MPEPR (wasting) with nested models. We compared a restricted model (excluding SPEI) with a full model (including SPEI). To ensure practical significance, we applied a Minimum Effect Size (MES) criterion, following the original Neyman-Pearson framework for hypothesis testing [@neyman1928; @neyman1928a; @perezgonzalez2015].

The MES was set at Cohen's threshold for small effects ($f^2 = 0.02$, equivalent to $\text{R}^2 = 0.01960784$). Thus, SPEI was considered significant only if its inclusion accounted for at least $1.960784\%$ of the variance in the dependent variable.

The test was structured as follows:

- **Null hypothesis** ($\text{H}_{0}$): Adding *SPEI* does not meaningfully improve the model’s fit, indicated by a negligible change in the adjusted $\text{R}^{2}$ or a non-significant F-test (with a Type I error probability ($\alpha$) of $0.05$).

- **Alternative Hypothesis** ($\text{H}_{a}$): Adding *SPEI* meaningfully improves the model’s fit, indicated by an increase in the adjusted $\text{R}^{2}$ exceeding the MES and a significant F-test (with $\alpha < 0.05$).

Formally:

$$
\begin{cases}
\text{H}_{0}: \Delta \ \text{Adjusted} \ \text{R}^{2} \leq \text{MES} \quad \text{or} \quad \text{F-test is not significant} \ (\alpha \geq 0.05) \\
\text{H}_{a}: \Delta \ \text{Adjusted} \ \text{R}^{2} > \text{MES} \quad \text{and} \quad \text{F-test is significant} \ (\alpha < 0.05)
\end{cases}
$$

Where:

$$
\Delta \ \text{Adjusted} \ \text{R}^{2} = \text{Adjusted} \ \text{R}^{2}_{\text{full}} - \text{Adjusted} \ \text{R}^{2}_{\text{restricted}}
$$

The restricted model is the same as the first model presented in this document, minus the SPEI variable: ~ `te(gini_index, gdp_per_capita)` + `s(year)` (Continuous `year`).

### Interpretation of Results

#### Standardized Precipitation Evapotranspiration Index (SPEI)

Since the original SPEI authors [@vicente-serrano2010] did not establish definitive thresholds for SPEI values and their corresponding drought conditions, we adopt the benchmark values provided by @mehr2020 (@tbl-mehr-et-al-2020-table-2).

::: {#tbl-mehr-et-al-2020-table-2}
| Classification         | SPI Threshold      | SPEI Threshold     |
|-------------------------|--------------------|--------------------|
| Extremely wet          | 2.0 ≤ SPI          | 1.83 ≤ SPEI        |
| Severely wet           | 2.0 > SPI ≥ 1.5    | 1.82 > SPEI ≥ 1.43 |
| Moderately wet         | 1.49 > SPI ≥ 1.0   | 1.42 > SPEI ≥ 1.0  |
| Near normal            | -1.0 ≤ SPI ≤ 1.0   | -1.0 ≤ SPEI ≤ 1.0  |
| Moderate drought (MoD) | -1.49 ≤ SPI < -1.0 | -1.42 ≤ SPEI < -1.0|
| Severe drought (SD)    | -2.0 ≤ SPI < -1.5  | -1.82 ≤ SPEI < -1.43|
| Extreme drought (ED)   | SPI < -2.0         | SPEI < -1.83       |

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

Classifications of wet and dry conditions using SPI and SPEI indices.
:::

#### Revised Multidimensional Index for Sustainable Food Systems (MISFS-R)

Results are also presented for each cluster of the Revised Multidimensional Index for Sustainable Food Systems (MISFS-R) ([@fig-norde-2023-figure-6]).

The MISFS is a tool designed to assess the sustainability of food systems at a subnational level in Brazil, incorporating local behaviors and practices. The MISFS-R is a revised version that introduces new indicators and a refined methodology for calculating the index ([@fig-norde-2023-figure-6]). For more details, see @carvalho2021a and @norde2023.

::: {#fig-norde-2023-figure-6}
![](images/norde-2023-figure-6.jpg){width=90%}

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

Dendrogram for cluster analysis between Brazilian states considering all the Revised Multidimensional Index for Sustainable Food Systems (MISFS-R) indicators and geographical location of each cluster.
:::

#### Pratical Significance

To ensure practical significance, the adjusted $\text{R}^2$ of the models are analysed for their effect sizes considering a confidence interval of $95\%$. We use @cohen1988a benchmark for interpretation.

A $\text{R}^2$ less than $\approx 0.0196$ is considered negligeble.

> SMALL EFFECT SIZE: $f^2 = .02$. Translated into $\text{R}^{2}$ or partial $\text{R}^{2}$ for Case 1, this gives $.02 / (1 + .02) = .0196$. We thus define a small effect as one that accounts for 2% of the $\text{Y}$ variance (in contrast with 1% for $r$), and translate to an $\text{R} = \sqrt{0196} = .14$ (compared to .10 for $r$). This is a modest enough amount, **just barely escaping triviality** and (alas!) all too frequently in practice represents the true order of magnitude of the effect being tested [@cohen1988a, p. 413].

> [...] in many circumstances, all that is intended by "proving" the null hypothesis is that the ES [Effect Size] is not necessarily zero but **small enough to be negligible** [...]. [@cohen1988a, p. 461].

## Setting the Enviroment

```{r}
#| eval: false

library(beepr)
library(broom)
library(cli)
library(clipr)
library(colorspace)
library(dplyr)
library(effectsize)
library(GGally)
library(geobr)
library(ggplot2)
library(ggspatial)
library(glmmTMB)
library(here)
library(janitor)
library(lme4)
library(lubridate)
library(lubritime) # github.com/danielvartan/lubritime
library(magrittr)
library(mgcv)
library(mgcViz)
library(MuMIn)
library(pal) # gitlab.com/rpkg.dev/pal
library(patchwork)
library(performance)
library(polyglotr)
library(prettycheck) # github.com/danielvartan/prettycheck
library(psychometric)
library(purrr)
library(r2glmm)
library(ragg)
library(RColorBrewer)
library(readr)
library(readxl)
library(rutils) # github.com/danielvartan/rutils
library(sidrar)
library(stats)
library(stringr)
library(summarytools)
library(tidyr)
```

```{r}
#| include: false

library(magrittr)
library(mgcv)
library(GGally)
library(ggplot2)
library(summarytools)
```

```{r}
source(here::here("R", "cohens_f_squared.R"))
source(here::here("R", "get_and_aggregate_sidra_by_year.R"))
source(here::here("R", "gam_misfs.R"))
source(here::here("R", "plot_brazil.R"))
source(here::here("R", "plot_dist.R"))
source(here::here("R", "plot_gam.R"))
source(here::here("R", "plot_ggally.R"))
source(here::here("R", "summarise_coefs.R"))
source(here::here("R", "summarise_r2.R"))
source(here::here("R", "tabset_panel_brazil_municipality.R"))
source(here::here("R", "tabset_panel_gam.R"))
source(here::here("R", "tabset_panel_gam_by_misfs.R"))
source(here::here("R", "tabset_panel_var_distribution.R"))
source(here::here("R", "tabset_panel_var_distribution_by_misfs.R"))
source(here::here("R", "utils.R"))
source(here::here("R", "utils-plots.R"))
```

## Importing and Tidying the Data

### Nutrition Data

Based on SISVAN [Nutritional Status](https://opendatasus.saude.gov.br/dataset/sisvan-estado-nutricional) dataset [@sisvana].

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

test_nutrition_data <-
  here::here("data", "Banco_dados_malnutritio_clima - Adaptado.csv") |>
  readr::read_csv(col_types = readr::cols(.default = "c")) |>
  janitor::clean_names() |>
  dplyr::rename(
    year = ano,
    municipality_code = code_muni,
  ) |>
  dplyr::select(year, municipality_code, mbepr, beipr, maper, mpepr) |>
  dplyr::mutate(
    dplyr::across(
      .cols = -dplyr::all_of(c("year", "municipality_code")),
      .fns = as.numeric
    ),
    dplyr::across(
      .cols = dplyr::all_of(c("year", "municipality_code")),
      .fns = as.integer
    )
  )
```

```{r}
#| output: false

nutrition_data <-
  here::here("data", "Banco_dados_malnutritio_clima - Adaptado.csv") |>
  readr::read_csv(col_types = readr::cols(.default = "c")) |>
  janitor::clean_names() |>
  dplyr::rename(
    year = ano,
    municipality_code = code_muni,
    sisvan_cover = cobrs,
    number_of_children = n_ao_de_criana_as_municipio_x,
    n_mbepr = muito_baixa_e_i_n_x,
    n_beipr = baixa_e_i_n_x,
    n_maper = magreza_acentuada_p_e_n_x,
    n_mpepr = magreza_p_e_n_x
  ) |>
  dplyr::select(
    year, municipality_code, misf, number_of_children, sisvan_cover,
    n_mbepr, n_beipr, n_maper, n_mpepr
  ) |>
  dplyr::mutate(
    dplyr::across(
      .cols = dplyr::all_of(
        c("number_of_children", "n_mbepr", "n_beipr", "n_maper", "n_mpepr")
      ),
      .fns = ~ 
        dplyr::case_when(
          !(as.numeric(.x) %% 1 == 0) ~ stringr::str_remove(.x, "\\."),
          TRUE ~ .x
        ) |>
        as.integer()
    )
  ) |>
  dplyr::mutate(
    # year = as.integer(year),
    year =
      year |>
      factor(
        levels = year |> unique() |> sort(),
        ordered = TRUE
      ),
    municipality_code = as.integer(municipality_code),
    misf = factor(misf, levels = c("A", "B", "C", "D"), ordered = FALSE),
    sisvan_cover = as.numeric(sisvan_cover),
    mbepr = n_mbepr,
    beipr = n_beipr,
    n_mbepr_beipr = n_mbepr + n_beipr,
    mbepr_beipr = n_mbepr_beipr,
    maper = n_maper,
    mpepr = n_mpepr,
    n_maper_mpepr = n_maper + n_mpepr,
    maper_mpepr = n_maper_mpepr
  ) |>
  dplyr::filter(
    dplyr::between(sisvan_cover, 0, 1),
    number_of_children >= 0,
    number_of_children > n_mbepr,
    number_of_children > n_beipr,
    number_of_children > mbepr_beipr,
    number_of_children > n_maper,
    number_of_children > n_mpepr,
    number_of_children > maper_mpepr
  ) |>
  dplyr::mutate(
    dplyr::across(
      .cols = dplyr::all_of(
        c("mbepr", "beipr", "mbepr_beipr", "maper", "mpepr", "maper_mpepr")
      ),
      .fns = ~
        .x  %>%
      `/`(number_of_children * sisvan_cover)
    ),
    dplyr::across(
      .cols = dplyr::all_of(
        c("mbepr", "beipr", "mbepr_beipr", "maper", "mpepr", "maper_mpepr")
      ),
      .fns = ~ pmax(0.00001, pmin(.x, 0.99999))
    )
  ) |>
  dplyr::select(
    year, municipality_code, misf, number_of_children, sisvan_cover, 
    n_mbepr, mbepr, n_beipr, beipr, n_mbepr_beipr, mbepr_beipr,
    n_maper, maper, n_mpepr, mpepr, n_maper_mpepr, maper_mpepr
  )
```

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

var <- "maper"
diff <- (nutrition_data[[var]] - test_nutrition_data[[var]])

summary(diff)
diff |> hist(breaks = seq(-1, 1, 0.01))
```

### GDP Per Capita Data

Source: IBGE-SIDRA [Table 5938](https://sidra.ibge.gov.br/tabela/5938) – Gross domestic product at current prices, taxes net of subsidies on products at current prices, and gross value added at current prices, total and by economic activity, and their respective shares - **Reference year 2010** [@ibgen].

```{r}
#| include: false

file_path <- here::here("data", "ibge-table-5938.rds")

if (prettycheck:::test_file_exists(file_path)) {
  ibge_table_5938 <- readRDS(file_path)
} else {
  cli::cli_abort("Data file not found.")
}
```

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

tictoc::tic()

ibge_table_5938 <- 
  get_and_aggregate_sidra_by_year(
    years = nutrition_data$year |> unique(),
    api_start = "/t/5938/n6/all/v/37/p/",
    api_end = "/d/v37%200"
  ) |>
  regularize_col_names()

tictoc::toc()

beepr::beep(1)
Sys.sleep(3)
```

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

# Last download: 2025-01-22 (~6m1s)

ibge_table_5938 |> 
  saveRDS(here::here("data", "ibge-table-5938.rds"))

ibge_table_5938 |> 
  readr::write_csv(here::here("data", "ibge-table-5938.csv"))
```

Source: IBGE-SIDRA [Table 6579](https://sidra.ibge.gov.br/tabela/6579) – Estimated resident population [@ibgel].

```{r}
#| include: false

file_path <- here::here("data", "ibge-table-6579.rds")

if (prettycheck:::test_file_exists(file_path)) {
  ibge_table_6579 <- readRDS(file_path)
} else {
  cli::cli_abort("Data file not found.")
}
```

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

tictoc::tic()

ibge_table_6579 <- 
  get_and_aggregate_sidra_by_year(
    years = nutrition_data$year |> unique(),
    api_start = "/t/6579/n6/all/v/all/p/",
    api_end = ""
  ) |>
  regularize_col_names()

tictoc::toc()

beepr::beep(1)
Sys.sleep(3)
```

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

# Last download: 2025-01-22 (~4m42s)

ibge_table_6579 |> 
  saveRDS(here::here("data", "ibge-table-6579.rds"))

ibge_table_6579 |> 
  readr::write_csv(here::here("data", "ibge-table-6579.csv"))
```

```{r}
gdp_data <- 
  ibge_table_5938 |>
  dplyr::left_join(
    ibge_table_6579,
    by = c("year", "municipality_code")
  ) |>
  dplyr::select(year, municipality_code, municipality.x, value.x, value.y) |>
  dplyr::rename(
    gdp = value.x,
    population = value.y,
    municipality = municipality.x
  ) |>
  dplyr::filter(gdp >= 0, population >= 0) |>
  dplyr::mutate(
    # year = as.integer(year),
    year =
      year |>
      factor(
        levels = year |> unique() |> sort(),
        ordered = TRUE
      ),
    gdp = gdp * 1000,
    gdp_per_capita = gdp / population,
    municipality_code = municipality_code |> as.integer()
  )
```

### Gini Index Data

::: {.callout-note}
The Brazilian Institute of Geography and Statistics ([IBGE](https://www.ibge.gov.br/)) Automatic Retrieval System ([SIDRA](https://sidra.ibge.gov.br/)) provides GINI data by municipality only for the year 1991 ([Table 115](https://sidra.ibge.gov.br/tabela/155)) [@ibgeo]. However, we found data processed by the Institute for Applied Economic Research ([IPEA](https://www.ipea.gov.br)), which had access to census data for the years 2000 and 2010 [@ibgep].
:::

Source: IBGE/IPEA – Gini Index of per capita household income by municipality: Period: 1991, 2000, and 2010 [@ibgel].

```{r}
gini_data <-
  "http://tabnet.datasus.gov.br/cgi/ibge/censo/bases/ginibr.csv" |>
  readr::read_delim(
    delim = ";",
    col_names = FALSE,
    col_types = readr::cols(.default = "c"),
    # locale = readr::locale(decimal_mark = ","),
    trim_ws = TRUE,
    skip = 3
  ) |>
  dplyr::rename_with(~ c("municipality", "x1991", "x2000", "x2010")) |>
  dplyr::slice(1:5565) |>
  dplyr::mutate(
    dplyr::across(
      .cols = dplyr::everything(),
      .fns = ~ iconv(.x, from = "latin1", to = "UTF-8")
    ),
    dplyr::across(
      .cols = dplyr::starts_with("x"),
      .fns = ~ dplyr::case_when(
        .x == "..." ~ NA,
        TRUE ~ .x |> stringr::str_replace_all(",", ".")
      )
    ),
    dplyr::across(
      .cols = dplyr::starts_with("x"),
      .fns = as.numeric
    )
  ) |>
  tidyr::pivot_longer(
    cols = starts_with("x"),
    names_to = "year",
    values_to = "gini_index"
  ) |>
  dplyr::mutate(
    year = year |> stringr::str_remove("x") |> as.integer(),
    municipality_code = 
      municipality |> 
      stringr::str_extract("\\d*") |> 
      as.integer(),
    municipality = 
      municipality |> 
      stringr::str_remove("\\d*") |> 
      stringr::str_trim()
  )|>
  dplyr::relocate(year, municipality_code, .before = municipality)
```

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

# Last download: 2025-01-29 (~5s)

gini_data |> 
  saveRDS(here::here("data", "ibge_gini_index.rds"))

gini_data |> 
  readr::write_csv(here::here("data", "ibge_gini_index.csv"))
```

### SPEI Data

Based on WorldClim 2.1 [Historical Monthly Weather]((https://worldclim.org/data/monthlywth.html)) dataset [@harris2020; @fick2017].

```{r}
#| output: false

spei_data <- 
  here::here("data", "spei_Extreme_drought_event_municipality_year2.csv") |>
  readr::read_csv(col_types = readr::cols(.default = "c")) |>
  janitor::clean_names() |>
  dplyr::rename(municipality_code = code_muni) |>
  dplyr::select(
    municipality_code, 
    dplyr::all_of(paste0("spei_12m_", nutrition_data$year |> unique()))
  ) |>
  dplyr::mutate(
    dplyr::across(
      .cols = dplyr::starts_with("spei_12m"),
      .fns = as.numeric
    )
  ) |>
  tidyr::pivot_longer(
    cols = starts_with("spei_12m"),
    names_to = "year",
    values_to = "spei_12m"
  ) |>
  dplyr::mutate(
    year =
      year |>
      stringr::str_remove("spei_12m_") |>
      as.integer(),
    year = 
      year |>
      factor(
        levels = year |> unique() |> sort(),
        ordered = TRUE
      ),
    municipality_code = municipality_code |> as.integer()
  )
```

### Analysis Data

<!-- https://geoftp.ibge.gov.br/organizacao_do_territorio/estrutura_territorial/evolucao_da_divisao_territorial_do_brasil/evolucao_da_divisao_territorial_do_brasil_1872_2010/notas_tecnicas.pdf -->

```{r}
#| output: false

data <- 
  nutrition_data |>
  dplyr::mutate(year = year |> as.character()) |>
  # In order to adapt the data to the municipalities in the nutrition and
  # gini data, since they do not have the check digit.
  dplyr::rename(temp_municipality_code = municipality_code) |>
  dplyr::left_join(
    y = gdp_data |> 
      dplyr::mutate(
        year = year |> as.character(),
        temp_municipality_code =
          municipality_code |>
          stringr::str_sub(end = -2) |>
          as.integer()
      ),
    by = c("year", "temp_municipality_code")
  ) |>
  dplyr::left_join(
    y = 
      gini_data |> 
      dplyr::filter(year == 2010) |>
      dplyr::select(municipality_code, gini_index) |>
      dplyr::rename(temp_municipality_code = municipality_code),
    by = c("temp_municipality_code"),
    relationship = "many-to-one"
  ) |>
  dplyr::left_join(
    y = spei_data |>
      dplyr::mutate(
        year = year |> as.character(),
        temp_municipality_code =
          municipality_code |>
          stringr::str_sub(end = -2) |>
          as.integer()
      ),
    by = c("year", "temp_municipality_code")
  ) |>
  dplyr::rename(
    municipality_code = municipality_code.x
  ) |>
  dplyr::select(
    year, municipality_code, municipality, misf, 
    number_of_children, sisvan_cover, 
    mbepr, beipr, mbepr_beipr, maper, mpepr, maper_mpepr,
    gini_index, gdp_per_capita, spei_12m
  ) |>
  dplyr::mutate(
    # year = as.integer(year),
    year =
      year |>
      factor(
        levels = year |> unique() |> sort(),
        ordered = TRUE
      )
  ) |>
  dplyr::filter(sisvan_cover >= 0.05) |>
  tidyr::drop_na(number_of_children, sisvan_cover)

data
```

```{r}
#| output: false

dplyr::tibble(
  n = data |> nrow(),
  n_sisvan_cover_less_than_0_05 = 
    data |> 
    dplyr::filter(sisvan_cover < 0.05) |> 
    nrow(),
  n_sisvan_cover_more_than_1 =
    data |> 
    dplyr::filter(sisvan_cover > 1) |> 
    nrow(),
  n_brazil_municipalities =
    geobr::read_municipality(year = 2022, showProgress = FALSE) |> 
    dplyr::pull(code_muni) |> 
    length() |>
    rutils::shush(),
  missing_municipalities = 
    geobr::read_municipality(year = 2022, showProgress = FALSE) |> 
    dplyr::pull(code_muni) |> 
    as.integer() |>
    setdiff(data |> dplyr::pull(municipality_code)) |>
    length() |>
    rutils::shush()
) |>
  tidyr::pivot_longer(cols = dplyr::everything())
```

#### Dictionary

- `year`: Year of data collection.
- `municipality_code`: Brazilian Institute of Geography and Statistics (IBGE) municipality code.
- `misf`: Cluster of the Revised Multidimensional Index for Sustainable Food Systems (MISFS-R) (A, B, C, or D).
- `number_of_children`: Number of children under five years old in the municipality.
- `sisvan_cover`: Proportion of children under five covered by the Brazilian Food and Nutrition Surveillance System (SISVAN) in the municipality.
- `mbepr`: Relative frequency of children under five years old with Very Short Stature for Age (*Muito Baixa Estatura Para a Idade*) in the municipality.
- `beipr`: Relative frequency of children under five years old with Short Stature for Age (*Baixa Estatura Para Idade*) in the municipality.
- `mbepr_beipr`: Relative frequency of children under five years old with Very Short/Short Stature for Age (*Muito Baixa/Baixa Estatura Para Idade*) in the municipality.
- `maper`: Relative frequency of children under five years old with Severe Thinness for Height (*Magreza Acentuada Para a Estatura*) in the municipality.
- `mpepr`: Relative frequency of children under five years old with Thinness for Height (*Magreza Por Estatura*) in the municipality.
- `maper_mpepr`: Relative frequency of children under five years old with Severe/Moderate Thinness for Height (*Magreza Acentuada/Moderada Para a Estatura*) in the municipality.
- `gini_index`: Gini index of per capita household income in the municipality at 2010. Due to limited data availability at the municipal level, all values are based on the 2010 Brazilian Census.
- `gdp_per_capita`: Gross Domestic Product (GDP) per capita in the municipality in Brazilian Reais (BRL) based at current prices on 2010.
- `spei_12m`: Standardised Precipitation Evapotranspiration Index (SPEI) in a 12-month times cale for the municipality.

## Checking SISVAN Cover

```{r}
brand_div_palette <- function(x) {
  make_color_vector(
    n_prop = 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}
#| eval: false
#| include: false

plot <-
  data |>
  dplyr::mutate(sisvan_cover = sisvan_cover * 100) |>
  dplyr::filter(year == 2019) |>
  plot_brazil_municipality(
    col_fill = "sisvan_cover",
    comparable_areas = TRUE,
    transform = "identity",
    breaks = seq(0, 100, 25),
    reverse = FALSE,
    limits = c(0, 100),
    print = FALSE,
    palette = brand_div_palette
  ) +
  ggplot2::labs(
    title = "Year: 2019"
  )

plot |> print() |> rutils::shush()
```

```{r}
#| eval: false

# Run this chunk to produce the plots and the animation.

data |> 
  dplyr::mutate(sisvan_cover = sisvan_cover * 100) |>
  animate_plot_brazil_municipality(
    col_fill = "sisvan_cover",
    col_group = "year",
    group_label = "Year",
    comparable_areas = TRUE
    suffix = NULL,
    width = 1344,
    height = 960,
    dpi = 150,
    transform = "identity",
    breaks = seq(0, 100, 25),
    reverse = FALSE,
    limits = c(0, 100),
    palette = brand_div_palette
  )
```

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

# Run this chunk to produce the panel-tabset for the plots.

tabset_panel_brazil_municipality(
  col = "sisvan_cover",
  caption = long_string(
      "
      Percentage of children under five monitored by Brazil's Food and 
      Nutrition Surveillance System (SISVAN) across Brazil's 
      historically comparable municipalities
      (minimum comparable areas [@ehrl2017]).
      "
  ),
  source = "Created by the authors based on data from @sisvana and @ibgel.",
  years = 2008:2019,
  heading = "###"
)
```

{{< include ./qmd/_panel-tabset-brazil-municipality-sisvan-cover.qmd >}}

## Checking SPEI Variations

### Year-to-Year comparison

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

plot <-
  data |>
  dplyr::filter(year == 2019) |>
  plot_brazil_municipality(
    col_fill = "spei_12m",
    comparable_areas = TRUE,
    transform = "identity",
    # breaks = c(-2, -1.83, -1.43, -1, 1, 1.42, 1.82, 2),
    breaks = seq(-2 , 2, 0.5),
    reverse = FALSE,
    limits = c(-2, 2),
    print = FALSE,
    palette = brand_div_palette
  ) +
  ggplot2::labs(
    title = "Year: 2019"
  )

plot |> print() |> rutils::shush()
```

```{r}
#| eval: false

# Run this chunk to produce the plots and the animation.

data |> 
  animate_plot_brazil_municipality(
    col_fill = "spei_12m",
    col_group = "year",
    group_label = "Year",
    comparable_areas = TRUE
    suffix = NULL,
    width = 1344,
    height = 960,
    dpi = 150,
    transform = "identity",
    # breaks = c(-2, -1.83, -1.43, -1, 1, 1.42, 1.82, 2),
    breaks = seq(-2 , 2, 0.5),
    reverse = FALSE,
    limits = c(-2, 2),
    palette = brand_div_palette
  )
```

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

# Run this chunk to produce the panel-tabset for the plots.

tabset_panel_brazil_municipality(
  col = "spei_12m",
  caption = long_string(
      "
      Standardized Precipitation Evapotranspiration Index (SPEI)
      across Brazil's historically comparable municipalities 
      (minimum comparable areas [@ehrl2017]).
      "
  ),
  source = long_string(
      "
      Created by the authors using data from the WorldClim 2.1
      [Historical Monthly Weather](https://worldclim.org/data/monthlywth.html) 
      dataset [@harris2020; @fick2017].
      "
  ),
  years = 2008:2019,
  heading = "####"
)
```

{{< include ./qmd/_panel-tabset-brazil-municipality-spei-12m.qmd >}}

### Comparison with External Data

@fig-fao-2025-figure-annex-6-1 presents an independent analysis for comparison, while our results are shown in @fig-spei-comparison.

::: {#fig-fao-2025-figure-annex-6-1}
![](images/fao-2025-figure-annex-6-1.png){width=90%}

[Source: Reproduced from @fao2025[Figure VI.1, p. 185].]{.legend}

Average Standardized Precipitation Evapotranspiration Index (SPEI) across Brazil's historically comparable municipalities (minimum comparable areas [@ehrl2017]).
:::

::: {#fig-spei-comparison}
```{r}
plot <- 
  data |>
  dplyr::filter(year >= 2010, year <= 2019) |>
  plot_brazil_municipality(
    col_fill = "spei_12m",
    comparable_areas = TRUE,
    transform = "identity",
    binned = FALSE,
    breaks = seq(-2 , 2, 0.5),
    reverse = FALSE,
    limits = c(-2, 2),
    print = FALSE,
    quiet = TRUE,
    palette = function(x) {
      make_color_vector(
        n_prop = x, 
        colors = c(
          "#6F0322",
          "#BC2D32",
          "#E28668",
          "#F9D0BC",
          "#F1F1F3",
          "#BBDAE9",
          "#65A9D1",
          "#266DB1",
          "#0D2648"
        )
      )
    }
  ) +
  ggplot2::labs(title = "2010-2019") +
  ggplot2::theme(plot.title = element_text(hjust = 0.5))

plot |> print() |> rutils::shush()
```

[Source: Created by the authors using data from the WorldClim 2.1 [Historical Monthly Weather](https://worldclim.org/data/monthlywth.html) dataset [@harris2020; @fick2017].]{.legend}

Average Standardized Precipitation Evapotranspiration Index (SPEI) across Brazil's historically comparable municipalities (minimum comparable areas [@ehrl2017]).
:::

## Checking Distributions

### General Distributions

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

# Run this chunk to produce the panel-tabset for the distributions.

data |>
  tabset_panel_var_distribution(
    cols = 
      names(data) |> 
      stringr::str_subset("^municipality", negate = TRUE),
    jitter = FALSE,
    heading = "####"
  )
```

{{< include ./qmd/_panel-tabset-var-distribution.qmd >}}

### By MISFS-R Clusters

```{r}
categories <- c("A", "B", "C", "D")

for (i in categories) {
  i_data_name  <- paste0("data_misfs_", tolower(i))
  
  assign(
    i_data_name,
    data |> dplyr::filter(misf == i)
  )
}
```

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

# Run this chunk to produce the panel-tabset for the distributions.

data |>
  tabset_panel_var_distribution_by_misfs(
    cols = 
      names(data) |> 
      stringr::str_subset("^municipality", negate = TRUE),
    jitter = FALSE,
    heading = "####"
  )
```

{{< include ./qmd/_panel-tabset-var-distribution-by-misfs.qmd >}}

## Checking Correlations

::: {#fig-correlation-matrix-1}
```{r}
#| fig-width: 10
#| fig-height: 7

data |>
  plot_ggally(
    cols = c(
      "year", "misf", "sisvan_cover", "gini_index", "gdp_per_capita", "spei_12m"
    ),
    mapping = ggplot2::aes(colour = misf)
  )
```

[Source: Created by the authors.]{.legend}

Correlation matrix of the secondary variables colored by the Revised Multidimensional Index for Sustainable Food Systems (MISFS-R) clusters.
:::

::: {#fig-correlation-matrix-2}
```{r}
#| fig-width: 10
#| fig-height: 10

data |>
  plot_ggally(
    cols = c(
      "year", "misf", "mbepr", "beipr", "mbepr_beipr", 
      "gini_index", "gdp_per_capita", "spei_12m"
    ),
    mapping = ggplot2::aes(colour = misf)
  )
```

[Source: Created by the authors.]{.legend}

Correlation matrix of the stature variables colored by the Revised Multidimensional Index for Sustainable Food Systems (MISFS-R) clusters.
:::

::: {#fig-correlation-matrix-3}
```{r}
#| fig-width: 10
#| fig-height: 10

data |>
  plot_ggally(
    cols = c(
      "year", "misf", "maper",  "mpepr", "maper_mpepr", 
      "gini_index", "gdp_per_capita", "spei_12m"
    ),
    mapping = ggplot2::aes(colour = misf)
  )
```

[Source: Created by the authors.]{.legend}

Correlation matrix of the thinness variables colored by the Revised Multidimensional Index for Sustainable Food Systems (MISFS-R) clusters.
:::

## Modeling the Data

::: {.callout-warning}
Please note that in some graphics, the x-axis is inverted to improve data comprehension.
:::

### By `s(spei_12m)` + `te(gini_index, gdp_per_capita)` + `s(year)` (Continuous `year`)

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 1

data |>
  dplyr::mutate(year = year |> as.character() |> as.integer()) |>
  tabset_panel_gam(
    rhs_formula = "s(spei_12m) + te(gini_index, gdp_per_capita) + s(year)",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-1.qmd >}}

### By `s(spei_12m)` + `s(gini_index)` + `s(gdp_per_capita)` + `s(year)` (Ordered `year`)

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 2

data |>
  dplyr::mutate(year = year |> as.character() |> as.integer()) |>
  tabset_panel_gam(
    rhs_formula = "s(spei_12m) + s(gini_index) + s(gdp_per_capita) + s(year)",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-2.qmd >}}

### By `s(spei_12m)` + `te(gini_index, gdp_per_capita)` + `year` (Unordered `year`)

::: {.callout-note}
In this model, the `year` variable is treated as a **unordered** categorical variable.
:::

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 3

data |>
  dplyr::mutate(year = factor(year, ordered = FALSE)) |>
  tabset_panel_gam(
    rhs_formula = "s(spei_12m) + te(gini_index, gdp_per_capita) + year",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-3.qmd >}}

### By `s(year)` (Continuous)

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 4

data |>
  dplyr::mutate(year = year |> as.character() |> as.integer()) |>
  tabset_panel_gam(
    rhs_formula = "s(year)",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-4.qmd >}}

### By `year` (Ordered)

::: {.callout-note}
In this model, the `year` variable is treated as a **ordered** categorical variable.

`.L`, `.Q`, and `.C` are, respectively, the coefficients for the ordered factor coded with linear, quadratic, and cubic contrasts.
:::

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 5

data |>
  tabset_panel_gam(
    rhs_formula = "year",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-5.qmd >}}

### By `year` (Unordered)

::: {.callout-note}
In this model, the `year` variable is treated as a **unordered** categorical variable.
:::

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 6

data |>
  dplyr::mutate(year = factor(year, ordered = FALSE)) |>
  tabset_panel_gam(
    rhs_formula = "year",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-6.qmd >}}

## Modeling the Data by MISFS-R Clusters

### By `s(spei_12m)` + `te(gini_index, gdp_per_capita)` + `s(year)` (Continuous `year`)

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 1

data |>
  dplyr::mutate(year = year |> as.character() |> as.integer()) |>
  tabset_panel_gam_by_misfs(
    rhs_formula = "s(spei_12m) + te(gini_index, gdp_per_capita) + s(year)",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-1-by-misfs.qmd >}}

### By `s(spei_12m)` + `te(gini_index, gdp_per_capita)` + `year` (Unordered `year`)

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 2

data |>
  dplyr::mutate(year = year |> as.character() |> as.integer()) |>
  tabset_panel_gam_by_misfs(
    rhs_formula = "s(spei_12m) + s(gini_index) + s(gdp_per_capita) + s(year)",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-2-by-misfs.qmd >}}

### By `s(spei_12m)` + `te(gini_index, gdp_per_capita)` + `year` (Unordered `year`)

::: {.callout-note}
In this model, the `year` variable is treated as a **unordered** categorical variable.
:::

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 3

data |>
  dplyr::mutate(year = factor(year, ordered = FALSE)) |>
  tabset_panel_gam_by_misfs(
    rhs_formula = "s(spei_12m) + te(gini_index, gdp_per_capita) + year",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-3-by-misfs.qmd >}}

### By `s(year)` (Continuous)

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 4

data |>
  dplyr::mutate(year = year |> as.character() |> as.integer()) |>
  tabset_panel_gam_by_misfs(
    rhs_formula = "s(year)",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-4-by-misfs.qmd >}}

### By `year` (Ordered)

::: {.callout-note}
In this model, the `year` variable is treated as a **ordered** categorical variable.

`.L`, `.Q`, and `.C` are, respectively, the coefficients for the ordered factor coded with linear, quadratic, and cubic contrasts.
:::

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 5

data |>
  tabset_panel_gam_by_misfs(
    rhs_formula = "year",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-5-by-misfs.qmd >}}

### By `year` (Unordered)

::: {.callout-note}
In this model, the `year` variable is treated as a **unordered** categorical variable.
:::

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

# Run this chunk to produce the panel-tabset for the models.

model_id <- 6

data |>
  dplyr::mutate(year = factor(year, ordered = FALSE)) |>
  tabset_panel_gam_by_misfs(
    rhs_formula = "year",
    type = model_id
  )
```

{{< include ./qmd/_panel-tabset-gam-6-by-misfs.qmd >}}

## Testing the Hypothesis

::: {.callout-note}
The rationale behind the test can be found in the Methods section.
:::

The hypothesis test was designed as follows:

$$
\begin{cases}
\text{H}_{0}: \Delta \ \text{Adjusted} \ \text{R}^{2} \leq \text{MES} \quad \text{or} \quad \text{F-test is not significant} \ (\alpha \geq 0.05) \\
\text{H}_{a}: \Delta \ \text{Adjusted} \ \text{R}^{2} > \text{MES} \quad \text{and} \quad \text{F-test is significant} \ (\alpha < 0.05)
\end{cases}
$$

The restricted model is identical to the first model presented in this document, except for the exclusion of the SPEI variable: ~ `te(gini_index, gdp_per_capita)` + `s(year)` (Continuous `year`).

### MBEPR & BEIPR – Very Short/Short Stature for Age (*Muito Baixa/Baixa Estatura Para Idade*)

```{r}
mbepr_beipr_restricted <- mgcv::gam(
  mbepr_beipr ~ te(gini_index, gdp_per_capita) + s(year),
  data = dplyr::mutate(data, year = as.integer(as.character(year))),
  family = mgcv::betar(link = "logit"),
  method = "REML"
)
```

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

mbepr_beipr_restricted |> summary()
```

```{r}
#| output: asis

adj_r_squared_restricted <-
  mbepr_beipr_restricted |>
  summarise_r2(
    dplyr::mutate(data, year = as.integer(as.character(year))) |>
      nrow()
  )

adj_r_squared_restricted|> md_named_tibble()
```

```{r}
mbepr_beipr_full <- mgcv::gam(
  mbepr_beipr ~ s(spei_12m) + te(gini_index, gdp_per_capita) + s(year),
  data = dplyr::mutate(data, year = as.integer(as.character(year))),
  family = mgcv::betar(link = "logit"),
  method = "REML"
)
```

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

mbepr_beipr_full |> summary()
```

```{r}
#| output: asis

adj_r_squared_full <-
  mbepr_beipr_full |>
  summarise_r2(
    dplyr::mutate(data, year = as.integer(as.character(year))) |>
      nrow()
  )

adj_r_squared_full|> md_named_tibble()
```

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

stats::anova(mbepr_beipr_restricted, mbepr_beipr_full, test = "F")
```

```{r}
effect_size <- cohens_f_squared_summary(
  adj_r_squared_restricted, 
  adj_r_squared_full
  )

effect_size |> list_as_tibble()
```

### MAPER & MPEPR – Severe/Moderate Thinness for Height (*Magreza Acentuada/Moderada Para a Estatura*)

```{r}
maper_mpepr_restricted <- mgcv::gam(
  maper_mpepr ~ te(gini_index, gdp_per_capita) + s(year),
  data = dplyr::mutate(data, year = as.integer(as.character(year))),
  family = mgcv::betar(link = "logit"),
  method = "REML"
)
```

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

maper_mpepr_restricted |> summary()
```

```{r}
#| output: asis

adj_r_squared_restricted <-
  maper_mpepr_restricted |>
  summarise_r2(
    dplyr::mutate(data, year = as.integer(as.character(year))) |>
      nrow()
  )

adj_r_squared_restricted|> md_named_tibble()
```

```{r}
maper_mpepr_full <- mgcv::gam(
  maper_mpepr ~ s(spei_12m) + te(gini_index, gdp_per_capita) + s(year),
  data = dplyr::mutate(data, year = as.integer(as.character(year))),
  family = mgcv::betar(link = "logit"),
  method = "REML"
)
```

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

maper_mpepr_full |> summary()
```

```{r}
#| output: asis

adj_r_squared_full <-
  maper_mpepr_full |>
  summarise_r2(
    dplyr::mutate(data, year = as.integer(as.character(year))) |>
      nrow()
  )

adj_r_squared_full|> md_named_tibble()
```

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

stats::anova(maper_mpepr_restricted, maper_mpepr_full, test = "F")
```

```{r}
effect_size <- cohens_f_squared_summary(
  adj_r_squared_restricted, 
  adj_r_squared_full
  )

effect_size |> list_as_tibble()
```

## Conclusion

Does the Standardized Precipitation Evapotranspiration Index (SPEI) significantly improve the prediction of childhood undernutrition in Brazilian municipalities?

Our analysis show that the SPEI **[does not]{.brand-red}** significantly improve the prediction of childhood undernutrition in Brazilian municipalities. The adjusted R-squared values indicate that the inclusion of the SPEI in the models does not significantly improve the prediction of the nutritional status of children under five years old in the municipalities.

## Acknowledgments

![](images/sustentarea-icon-rgb-150-dpi.png){style="width: 17.5%;"}

This analysis is part of the [Sustentarea](https://www.fsp.usp.br/sustentarea) Research and Extension Group's project: *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 (SUS)*.

![](images/cnpq-logo.png){style="width: 30%;"}

This research was supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico - Brazil ([CNPq](https://www.gov.br/cnpq/)).

## How to Cite

To cite this work, please use the following format:

Magalhães, A. R., Vartanian, D, & Carvalho, A. M. (2025). *Global syndemic project data analysis: Report 3: Exploring potential associations between childhood undernutrition and the Standardized Precipitation Evapotranspiration Index (SPEI) in Brazilian municipalities (2008–2019)*. Sustentarea Research and Extension Group at the University of São Paulo. https://sustentarea.github.io/gs-data-analysis-report-3

A BibTeX entry for LaTeX users is

```         
@techreport{magalhaes2025,
  title = {Global syndemic project data analysis: Report 3: Exploring potential associations between childhood undernutrition and the Standardized Precipitation Evapotranspiration Index (SPEI) in Brazilian municipalities (2008–2019)},
  author = {{Arthur Ramalho Magalhães} and {Daniel Vartanian} and {Aline Martins de Carvalho}},
  year = {2025},
  address = {São Paulo},
  institution = {Sustentarea Research and Extension Group at the University of São Paulo},
  langid = {en},
  url = {https://sustentarea.github.io/gs-data-analysis-report-3}
}
```

## References {.unnumbered}

::: {#refs}
:::