---
title: "Stunting Trends in Mozambique | 2001-2020"
output: 
  flexdashboard::flex_dashboard:
    theme:
      bootswatch: zephyr
    source_code: embed
    orientation: columns
    vertical_layout: fill
---

```{r setup}
knitr::opts_chunk$set(echo = FALSE)
library(flexdashboard)
library(tidyverse)
library(readxl)
library(shiny)
library(shinyWidgets)
library(plotly)
library(ggsci)
library(RColorBrewer)
library(sf)
library(bslib)
library(bsicons)
library(fontawesome)
library(rnaturalearth)
library(rnaturalearthdata)
library(rnaturalearthhires)
library(tmap)
library(tmaptools)
tmap_mode("view")
```

```{r}
my_data <- read_excel("My_Data.xlsx")
stunting1 <- read_excel("UNICEF/UNICEF_WHO_WB_Global_Expanded_Databases_Stunting_May_2023_clean.xlsx")
```

```{r warning=FALSE}
stunting1 <- stunting1 %>%
  mutate(across(.cols = 3:42, .fns = as.numeric))

my_data <- my_data %>%
  mutate(across(c(2:41), round, 1))
```

```{r}
moz_map <- ne_states(geounit = "mozambique", returnclass = "sf")

moz_map$gn_name[moz_map$gn_name == 'Cabo Delgado Province'] <- 'Cabo Delgado'
moz_map$gn_name[moz_map$gn_name == 'Niassa Province'] <- 'Niassa'
moz_map$gn_name[moz_map$gn_name == 'Manica Province'] <- 'Manica'
moz_map$gn_name[moz_map$gn_name == 'Gaza Province'] <- 'Gaza'
moz_map$gn_name[moz_map$gn_name == 'Provincia de Zambezia'] <- 'Zambezia'
moz_map$gn_name[moz_map$gn_name == 'Cabo Delgado Province'] <- 'Cabo Delgado'
moz_map$gn_name[moz_map$gn_name == 'Cidade de Maputo'] <- 'Maputo City'
moz_map$gn_name[moz_map$gn_name == 'Inhambane Province'] <- 'Inhambane'
moz_map$gn_name[moz_map$gn_name == 'Sofala Province'] <- 'Sofala'
moz_map$gn_name[moz_map$gn_name == 'Maputo Province'] <- 'Maputo'

moz_map <- moz_map %>%
  rename(Province = gn_name)
```

```{r}
stunting_map <- stunting1 %>%
  select(2, 32:42) %>%
  pivot_longer(2:12, names_to = "Province", values_to = "Prevalence")

stunting_map$Province[stunting_map$Province == 'Maputo Cidade'] <- 'Maputo City'
stunting_map$Province[stunting_map$Province == 'Maputo Provincia'] <- 'Maputo'

stunting_map <- stunting_map %>%
  left_join(moz_map, stunting_map, by = 'Province')

stunting_map <- stunting_map %>%
  select(c(1:3, 40:41, 124))

stunting_map <- st_as_sf(stunting_map, crs = 4326)
```

```{r}
stunting_range_5m <- stunting1 %>%
  select(-one_of("0 to 23 months", "24 to 59 months")) %>%
  select("Year", c(8:13)) 

stunting_range_5m <- stunting_range_5m %>%        
  pivot_longer(cols = 2:7, names_to = "Range", values_to = "Prevalence") 

desired_order <- c('0 to 5 months', '6 to 11 months', '12 to 23 months', '24 to 35 months', '36 to 47 months', '48 to 59 months')

stunting_range_5m$Range <- factor(stunting_range_5m$Range, levels = desired_order)

stunting_range_5m <- stunting_range_5m
```

```{r}
stunting_range_12m <- stunting1 %>%
  select("Year", c(19, 23, 27, 31))

stunting_range_12m <- stunting_range_12m %>%        
  pivot_longer(cols = 2:5, names_to = "Gender/Range", values_to = "Prevalence")

desired_order12 <- c('Male 0 to 23 months', 'Female 0 to 23 months', 'Male 24 to 59 months', 'Female 24 to 59 months')

stunting_range_12m$`Gender/Range` <- factor(stunting_range_12m$`Gender/Range`, levels = desired_order12)
```

# Introduction

According to UNICEF, stunting refers to a child who is too short for his or her age. These children can suffer severe irreversible physical and cognitive damage that accompanies stunted growth. The devastating effects of stunting can last a lifetime and even affect the next generation. Stunting is considered as the main indicator of childhood undernutrition.

Undernutrition in children weakens the immune system and makes children more susceptible to dying from common illnesses such as pneumonia, diarrhea and malaria. Other negative effects include delayed motor development, impaired cognitive function, poor school performance, and lower income earning capacity as an adult.

For this "visualization" exercise, I will focus on Stunting trends of \<5 years old children in Mozambique:

1.  General \<5 Stunting Trends in Mozambique in the first Tab;

2.  Geographical distribution of Stunting in Mozambique over the 2001-2015 period in the second Tab;

3.  Potential correlation between stunting and household's economical background then breastfeeding practices in the third and final Tab.

As this project is entirely focused on visualization, I won't go through extensive data import, cleaning and wrangling processes. I used SQLiteStudio to tidy and transform the [**stunting**]{.underline} data downloaded from UNICEF website: <https://data.unicef.org/wp-content/uploads/2019/04/UNICEF_WHO_WB_Global_Expanded_Databases_Stunting_May_2023.xlsx>

I used SQLiteStudio as well to compile different UNICEF/WHO datasets related to \<5 years old children's information (breastfeeding and maternal health) into one single dataset: [**my_data**]{.underline}

# General Trends

## Column {.tabset .tabset-fade data-width="600"}

Regular prevalence decrease between 2001 and 2020 with Males being more affected than Females (approx. +6-7 % points) and rural areas more affected than urban areas (approx +10-15 % points)

### Countrywide Trend

```{r}
z1 <- stunting1

z1 <- z1 %>%
  plot_ly() %>%
  add_trace(
    x = ~Year,
    y = ~National,
    colors = "steelblue3",
    type = 'scatter',
    mode = 'line',
    line = list()
  ) %>%
  layout(
    title = list(
      text = "MZ Stunting Prevalence of <5 (2001-2020)", font = list(size = 16)
    ),
      xaxis = list(title = "Year"),
      yaxis = list(title = ">5 Stunting Prevalence (%)", range = c(25, 60),
      showlegend = FALSE),
      hovermode = "x unified"
  )

z1
```

### Trend by Gender

```{r}
library(plotly)

stunting_gender <- stunting1 %>%
  pivot_longer(cols = 4:5, names_to = "Gender", values_to = "Prevalence")

acolors <- c("steelblue3", "coral1")

z1.1 <- stunting_gender

z1.1 <- z1.1 %>%
  plot_ly() %>%
  add_trace(
    x = ~Year,
    y = ~Prevalence,
    color = ~Gender,
    colors = acolors,
    type = 'scatter',
    mode = 'line',
    line = list()
  ) %>%
  layout(
    title = list(
      text = "MZ Stunting Prevalence of <5 (2001-2020)", font = list(size = 16)
    ),
    xaxis = list(title = "Year"),
    yaxis = list(title = "<5 Stunting Prevalence (%)", range = c(25, 60)),
    hovermode = "x unified"
  )

z1.1
```

### Rural / Urban

```{r}
library(plotly)

stunting_geo <- stunting1 %>%
  pivot_longer(cols = 6:7, names_to = "Geog", values_to = "Prevalence")

bcolors <- c("coral1", "steelblue3")

z1.2 <- stunting_geo

z1.2 <- z1.2 %>%
  plot_ly %>%
  add_trace(
    x = ~Year,
    y = ~Prevalence,
    color = ~Geog,
    colors = bcolors,
    type = 'scatter',
    mode = 'line',
    line = list()
  ) %>%
  layout(
    title = list(
      text = "MZ Stunting Prevalence of <5 (2001-2020)", font = list(size = 16)
    ),
    xaxis = list(title = "Year"),
    yaxis = list(title = "<5 Stunting Prevalence (%)", range = c(25, 60)),
    hovermode = "x unified"
  )

z1.2
```

## Column {.tabset .tabset-fade data-width="400"}

### Trend by Age Range

```{r}
library(shiny)

ui <- fluidPage(
    
    sidebarLayout(
        sidebarPanel(
            checkboxGroupInput(
                inputId = "rangeInput",
                label = "Age Range in Months",
                choices = unique(stunting_range_5m$`Range`),  
                selected = (stunting_range_5m$`Range`)[c(1:6)]),
                width = 2,
                position = "left"
        ),
    
        mainPanel(
            plotlyOutput("stuntingPlot")
        )
    )
)

server <- function(input, output) {
      output$stuntingPlot <- renderPlotly({
        filtered_data <- stunting_range_5m %>%
            filter(Range %in% input$rangeInput)

        plot_ly(filtered_data, x = ~Year, y = ~Prevalence, color = ~Range, colors = "Set2", type = 'scatter', mode = 'line') %>%
            layout(title = "MZ Stunting Prevalence of <5 (2001-2020)",
                   xaxis = list(title = "Year"),
                   yaxis = list(title = "<5 Stunting Prevalence (%)"),
                   showlegend = FALSE,
                   hovermode = "x unified"
                  )
    })
}

shinyApp(ui = ui, server = server)
```

### Trend by Gender and Age Range

```{r}
library(shiny)

ui <- fluidPage(
    
      sidebarLayout(
        sidebarPanel(
            checkboxGroupInput(
              inputId = "rangeInput",
              label = "Gender & Age Range in Months",
              choices = unique(stunting_range_12m$`Gender/Range`),
              selected = (stunting_range_12m$`Gender/Range`)[c(1:4)]),
              width = 2,
            position = "left"
        ),
        
        mainPanel(
           plotlyOutput("stuntingPlot"),
           width = 10
        )
    )
)

server <- function(input, output) {
    output$stuntingPlot <- renderPlotly({
        filtered_data <- stunting_range_12m %>%
            filter(`Gender/Range` %in% input$rangeInput)

        plot_ly(filtered_data, x = ~Year, y = ~Prevalence, color = ~`Gender/Range`, colors = "Set2", type = 'scatter', mode = 'line') %>%
            layout(title = "MZ Stunting Prevalence of <5 (2001-2020)",
                   xaxis = list(title = "Year"),
                   yaxis = list(title = "<5 Stunting Prevalence (%)"),
                   showlegend = FALSE,
                   hovermode = "x unified"
                   )
    })
}

shinyApp(ui = ui, server = server)
```

# Prevalence Evolution by Province \| Maps

In line with the socio-cultural and economical divide between Southern and Northern Mozambique, the Northern part of the country remains the most affected by stunting.

## Column {data-width="250"}

### Stunting Map \| 2001

```{r message=FALSE}
map1 <- stunting_map %>%
  filter(Year == 2001)

tm_shape(map1) +
  tm_polygons("Prevalence",
              palette = "Blues",
              breaks = seq(0, 70, by = 10),
              id = "Province") +
  tm_text(text = "Prevalence",
          bg.color = "white") +
  tm_view(view.legend.position = c("right", "bottom"))
```

## Column {data-width="250"}

### Stunting Map \| 2008

```{r}
map2 <- stunting_map %>%
  filter(Year == 2008)

tm_shape(map2) +
  tm_polygons("Prevalence",
              palette = "Blues",
              breaks = seq(0, 70, by = 10),
              id = "Province",
              ) +
  tm_text(text = "Prevalence",
          bg.color = "white") +
  tm_view(view.legend.position = c("right", "bottom"))
```

## Column {data-width="250"}

### Stunting Map \| 2015

```{r}
map3 <- stunting_map %>%
  filter(Year == 2015)

tm_shape(map3) +
  tm_polygons("Prevalence",
              palette = "Blues",
              breaks = seq(0, 70, by = 10),
              id = "Province") +
  tm_text(text = "Prevalence",
          bg.color = "white") +
  tm_view(view.legend.position = c("right", "bottom"))
```

## Column {data-width="250"}

### Stunting Map \| 2020

```{r}
map4 <- stunting_map %>%
  filter(Year == 2020)

tm_shape(map4) +
  tm_polygons("Prevalence",
              palette = "Blues",
              breaks = seq(0, 70, by = 10),
              id = "Province") +
  tm_text(text = "Prevalence",
          bg.color = "white") +
  tm_view(view.legend.position = c("right", "bottom"))
```

# Potential Correlations

## Column {data-width="500"}

To few data available to confirm any clear correlation. It is nevertheless quite clear that the poorer a household is the higher the probability to have stunted children is.

### Stunting & Household Wealth Quintile (WQ)

```{r}
stunting_wq <- my_data %>%
  select(1:2, 7:11) %>%
  rename(c("Stunting_Prevalence" = "stun_total", "Poorest" = "stun_q1", "Wealth_Quintile_2" = "stun_q2", "Wealth_Quintile_3" = "stun_q3", "Wealth_Quintile_4" = "stun_q4", "Richest" = "stun_q5"))
```

```{r}
z3.1 <- stunting_wq %>%
  plot_ly() %>%
  add_trace(x = ~Year,
            y = ~Stunting_Prevalence,
            type = 'scatter',
            mode = 'lines+markers',
            legendgroup = 'group1',
            name = 'MZ <5 Stunting %')
z3.1 <- z3.1 %>% 
  add_trace(x = ~Year,
            y = ~Poorest,
            type = "scatter",
            mode = "markers",
            size = 1.2,
            legendgroup = 'group2',
            name = 'WQ 1 - Poorest')
z3.1 <- z3.1 %>% 
  add_trace(x = ~Year,
            y = ~Wealth_Quintile_2,
            type = "scatter",
            mode = "markers",
            size = 1.2,
            legendgroup = 'group2',
            name = 'WQ 2')
z3.1 <- z3.1 %>% 
  add_trace(x = ~Year,
            y = ~Wealth_Quintile_3,
            type = "scatter",
            mode = "markers",
            size = 1.2,
            legendgroup = 'group2',
            name = 'WQ 3')
z3.1 <- z3.1 %>% 
  add_trace(x = ~Year,
            y = ~Wealth_Quintile_4,
            type = "scatter",
            mode = "markers",
            size = 1.2,
            legendgroup = 'group2',
            name = 'WQ 4')
z3.1 <- z3.1 %>% 
  add_trace(x = ~Year,
            y = ~Richest,
            type = "scatter",
            mode = "markers",
            size = 1.2,
            legendgroup = 'group2',
            name = 'WQ 5 - Richest')

z3.1 <- z3.1 %>% 
  layout(
    hovermode = "x unified")
z3.1
```

## Column {.tabset .tabset-fade data-width="500"}

Exclusive breastfeeding (for the first - months) has regularly increased over the 2004-2015 period despite important differences between gender, rural/urban origines and wealth. Too little data to establish clear correlation.

```{r}
stunting_bf <- my_data %>% 
  select(1:2, 32:41) %>% 
  rename(c("Stunting_Prevalence" = "stun_total", "BF_Total" = "bf_total", "BF_Male" = "bf_m", "BF_Female" = "bf_f",  "BF_Rural" = "bf_rur", "BF_Urban" = "bf_urb", "BF_WQ1" = "bf_q1", "BF_WQ2" = "bf_q2", "BF_WQ3" = "bf_q3", "BF_WQ4" = "bf_q4", "BF_WQ5" = "bf_q5")) %>% 
  filter(!row_number() %in% c(3))

```

### Exclusive Breastfeeding Prevalence

```{r}
z3.2 <- stunting_bf

z3.2 <- z3.2 %>%
  plot_ly() %>%
  add_trace(
    x = ~Year,
    y = ~BF_Total,
    colors = "steelblue3",
    type = 'scatter',
    mode = 'lines+markers',
    name = 'BF Total'
  ) %>%
  layout(
    title = list(
      text = "MZ Exclusive Brestfeeding Prevalence (2003-2015)", font = list(size = 16)
                ),
      xaxis = list(title = "Year"),
      yaxis = list(title = "Excl. Breastfeeding Prevalence (%)", range = c(15, 60),
      showlegend = FALSE),
      hovermode = "x unified"
        )

z3.2 <- z3.2 %>% 
  add_trace(
    x = ~Year,
    y = ~BF_Male,
    type = "scatter",
    mode = 'lines+markers', 
    line = list(color = 'tomato3', width = .7, dash = 'dash'),
    legendgroup = 'group2',
    name = 'Male'
    )
z3.2 <- z3.2 %>% 
  add_trace(
    x = ~Year,
    y = ~BF_Female,
    type = "scatter",
    mode = 'lines+markers', 
    line = list(color = 'dodgerblue', width = .7, dash = 'dash'),
    legendgroup = 'group2',
    name = 'Female'
    )
z3.2 <- z3.2 %>% 
  add_trace(
    x = ~Year,
    y = ~BF_Rural,
    type = "scatter",
    mode = 'lines+markers', 
    line = list(color = 'palegreen4', width = .7, dash = 'dot'),
    legendgroup = 'group3',
    name = 'Rural'
    )
z3.2 <- z3.2 %>% 
  add_trace(
    x = ~Year,
    y = ~BF_Urban,
    type = "scatter",
    mode = 'lines+markers', 
    line = list(color = 'tan4', width = .7, dash = 'dot'),
    legendgroup = 'group3',
    name = 'Urban'
    )
z3.2 <- z3.2 %>% 
  add_trace(x = ~Year,
            y = ~BF_WQ1,
            type = "scatter",
            mode = "markers",
            size = 1.2,
            legendgroup = 'group4',
            name = 'WQ 1 - Poorest')
z3.2 <- z3.2 %>% 
  add_trace(x = ~Year,
            y = ~BF_WQ2,
            type = "scatter",
            mode = "markers",
            size = 1.2,
            legendgroup = 'group4',
            name = 'WQ 2')
z3.2 <- z3.2 %>% 
  add_trace(x = ~Year,
            y = ~BF_WQ3,
            type = "scatter",
            mode = "markers",
            size = 1.2,
            legendgroup = 'group4',
            name = 'WQ 3')
z3.2 <- z3.2 %>% 
  add_trace(x = ~Year,
            y = ~BF_WQ4,
            type = "scatter",
            mode = "markers",
            size = 1.2,
            legendgroup = 'group4',
            name = 'WQ 4')
z3.2 <- z3.2 %>% 
  add_trace(x = ~Year,
            y = ~BF_WQ5,
            type = "scatter",
            mode = "markers",
            size = 1.2,
            legendgroup = 'group4',
            name = 'WQ 5 - Richest')
z3.2
```

### Stunting & Exclusive Breastfeeding

```{r}
z3.3 <- stunting_bf

z3.3 <- z3.3 %>%
  plot_ly() %>%
  add_trace(
    x = ~Year,
    y = ~Stunting_Prevalence,
    type = 'scatter',
    mode = 'lines+markers',
    name = 'Stunting'
  ) %>%
  layout(
    title = list(
      text = "Stunting & Exclusive Breastfeeding", font = list(size = 16)
                ),
      xaxis = list(title = "Year"),
      yaxis = list(title = "Prevalence (%)", range = c(25, 50),
      showlegend = FALSE),
      hovermode = "x unified"
        )

z3.3 <- z3.3 %>% 
  add_trace(x = ~Year,
            y = ~BF_Total,
            type = 'scatter',
    mode = 'lines+markers',
            name = 'Exclusive BF')

z3.3
```