app.Rmd

---
title: "Climate and Conflict Impacts on Child Malnutrition"
runtime: shiny
css: custom.css
output: 
  flexdashboard::flex_dashboard:
    favicon: images/favicon.png
    orientation: rows
    vertical_layout: fill
    # social: menu
    source_code: embed
    theme: 
      bootswatch: united
---

```{r setup}
source("utils.r")
library(shiny)
library(shinyWidgets)
library(flexdashboard)
library(bslib)
library(thematic)
library(leaflet)
library(srvyr)
library(showtext)
library(gtsummary)
library(marginaleffects)

# Match plot themes to dashboard theme 
thematic_rmd(font = "auto")
theme_set(theme_minimal())
knitr::opts_chunk$set(fig.showtext = TRUE)

# Load supplementary georeferenced datasets 
# Write metadata used in the UI 
meta <- here("data/geo") %>% 
  get_shapes() %>%  
  imap(~{
    list(
      CHIRPS = list(
        dat = path(here("data/chirps"), .y, ext = "tif") %>% rast,
        color = "#0098CE", 
        icon = "ion-waterdrop",
        txt = "Total monthly precipitation (mm) averaged",
        source = "https://data.chc.ucsb.edu/products/CHIRPS-2.0/africa_daily",
        info = "CHIRPS v2.0",
        cap = "A 10mm increase in average monthly rainfall"
      ),
      CHIRTS = list(
        dat = path(here("data/chirts"), .y, ext = "tif") %>% rast,
        color = "#D82310", 
        icon = "ion-android-sunny",
        txt = "Maxmimum monthly temperature (C°) averaged",
        source = "https://data.chc.ucsb.edu/experimental",
        info = "CHIRTS-ERA5 (experimental)",
        cap = "A 5 C° increase in the average monthly maximum" 
      ), 
      NDVI =  list(
        dat = path(here("data/ndvi"), .y, ext = "tif") %>% rast, 
        color = "#154D14", 
        icon = "ion-leaf",
        txt = "Maxmimum greenness value (range -1 to 1) averaged",
        source = "https://lpdaac.usgs.gov/products/mod13c1v006",
        info = "MOD13C1 v006",
        cap = "A 0.1 point increase in maximum greenness" 
      ),
      UCDP = list(
        dat = path(here("data/ucdp"), .y, ext = "rds") %>% read_rds,
        color = "#772953", 
        icon = "ion-alert-circled",
        txt = "Count of lethal armed conflict events",
        source = "https://ucdp.uu.se/downloads/olddw.html", 
        info = "GED v19.1",
        cap = "A lethal conflict in the region"
      ),
      WHZ = list(
        color = "#99469D", 
        icon = "ion-woman",
        source = "www.idhsdata.org/idhs-action/variables/HWWHZWHO",
        info = "IPUMS DHS"
      ),
      corr = list(
        pos = list(
          color =  "#2AAD5350",
          icon = "ion-arrow-graph-up-right",
          cap = "<strong>significantly increases</strong> predicted WHZ by"
        ),
        neg = list(
          color = "#EFB73E80",
          icon = "ion-arrow-graph-down-right", 
          cap = "<strong>significantly decreases</strong> predicted WHZ by"
        ),
        insig = list(
          color = "#E9ECEF",
          icon = "ion-android-cancel",
          cap = "has <strong>no significant effect</strong> on WHZ"
        )
      ),
      shapes = .x
    )
  })

# Load survey dataset `dat` 
# Only include samps that have GPS coords and 
# where historical NDVI records are available
dat <- here("data/ipums_dhs.rds") %>% 
  read_rds %>% 
  mutate(sample_lbl = paste(country, sample_yr), iso = iso %>% tolower) %>% 
  filter(
    has_gps, 
    min_date > meta %>% 
      map(~names(.x$NDVI$dat) %>% as_date %>% min) %>% 
      list_c() %>%
      min
  )
```

{.sidebar}
======================================================================

```{r inputs}
selectInput(
  'samp',
  label = 'DHS Sample', 
  choices = dat$sample_lbl, 
  selected = "Kenya 2003"
)

selectInput(
  'exp', 
  label = 'Climate or Conflict Exposure', 
  choices = meta %>% 
    list_transpose %>% 
    names %>% 
    keep(~str_detect(.x, "shapes|corr|WHZ", TRUE)),
  selected = "CHIRTS"
)

sliderInput(
  'space', 
  label = 'Exposure Radius (km)', 
  min = 0, 
  max = 50, 
  value = 10, 
  step = 2.5
)

sliderInput(
  'time', 
  label = 'Exposure Period (months)', 
  min = 0, 
  max = 24, 
  value = c(1, 3),
  step = 1, 
)

materialSwitch(
  "eas", 
  label = "Map Enumeration Areas?",
  value = TRUE,
  status = "primary"
)


```

```{r reactives}
# Adjust figures if user clicks an EA 
ea_clicked <- reactiveVal(NULL)
observeEvent(input$fig1_shape_click, ea_clicked(input$fig1_shape_click$id))
observeEvent(input$fig1_click, ea_clicked(NULL))
observeEvent(input$samp, ea_clicked(NULL))

# Reactive output 
out <- reactive({
  
  # Get user-defined exposure type and period
  samp <- dat %>% filter(sample_lbl == input$samp)
  iso <- samp$iso
  exposure <- meta[[iso]][[input$exp]]$dat
  shapes <- meta[[iso]]$shapes %>% st_as_sf()
  start <- floor_date(samp$first_interview - months(max(input$time)), "month")
  stop <- ceiling_date(samp$first_interview - months(min(input$time)), "month") 
  
  # Get user-defined EA buffers
  km <- input$space
  eas <- list_c(samp$data) %>% 
    filter(!is.na(whz)) %>% 
    count(cluster, geometry, geo) %>% 
    st_as_sf() %>%
    st_transform(crs = samp$epsg) %>% 
    st_buffer(1000 * km) %>%
    st_transform(crs = 4326) 
  
  # National-level Exposure for maps 
  ntl_exposure <- if(input$exp == "UCDP"){
    exposure %>% filter(dates %>% int_overlaps(interval(start, stop))) 
  } else {
    days <- seq(start, stop - days(1), by = "days") %>% as.character()
    months <- exposure[[which(names(exposure) %in% days)]]
    if(input$exp == "NDVI"){max(months)} else {mean(months)}
  }
  
  # Admin1-level summaries (mostly for UCDP maps)
  admin1_exposure <- if(input$exp == "UCDP"){
    idx <- ntl_exposure %>% st_nearest_feature(shapes)
    ntl_exposure$ADMIN_NAME <- shapes$ADMIN_NAME[idx]
    ntl_exposure %>% 
      st_drop_geometry() %>% 
      count(ADMIN_NAME, name = "exposure") %>% 
      left_join(shapes, ., by = "ADMIN_NAME") %>% 
      mutate(exposure = exposure %>% replace_na(0)) %>% 
      select(geo = ADMIN_NAME, exposure)
  } else {
    shapes$exposure <- shapes$geometry %>% 
      map_dbl(~{
        vals <- ntl_exposure %>% crop(vect(.x)) %>% values()
        if(input$exp == "NDVI"){
          max(vals, na.rm = TRUE)
        } else {
          mean(vals, na.rm = TRUE)
        }
      })
    shapes 
  }
  
  # Set `weights = TRUE` for a spatially weighted `fun`, but sacrifice speed
  ea_exposure <- if(input$exp == "UCDP"){
    eas <- eas %>% filter(!st_is_empty(.)) 
    idx <- eas %>% st_nearest_feature(admin1_exposure)
    eas$exposure <- admin1_exposure$exposure[idx]
    eas 
  } else {
    eas$exposure <- ntl_exposure %>%
      zonal(vect(eas), fun = mean, weights = FALSE, na.rm = TRUE) %>%
      unlist()
    eas
  }
  
  # Update sample data with personal exposure 
  personal_exposure <- ea_exposure %>% 
    st_drop_geometry() %>% 
    select(-geo) %>% 
    right_join(list_c(samp$data), by = "cluster", multiple = "all") 
  
  # Calculations for value boxes 1 and 2
  value_boxes <- personal_exposure %>% 
    as_survey_design(weight = weight, id = cluster) %>% 
    mutate(
      acute = whz <= -2,
      severe = whz <= -3
    ) %>% 
    summarise(
      across(
        .names = "{.col}_mean",
        c(exposure, whz, acute, severe),
        ~.x %>% survey_mean(na.rm = TRUE, vartype = "ci") %>% 
          list()
      ), 
      across(
        .names = "{.col}_quant",
        c(exposure, whz, acute, severe),
        ~.x %>% 
          survey_quantile(na.rm = TRUE, vartype = NULL, 
                          quantiles = c(0.25, 0.50, 0.75)) %>% 
          list()
      )
    ) %>% 
    pivot_longer(everything(), names_sep = "_", 
                 names_to = c("measure", "stat")) %>% 
    pivot_wider(names_from = stat, values_from = value) %>% 
    unnest(where(is.list)) %>% 
    rowwise() %>% 
    mutate(across(where(is.double), ~{
      x <- if(measure %in% c("acute", "severe")){
        scales::percent(.x, 0.1)
      } else {
        ifelse(input$exp == "UCDP", scales::number(.x, 1), 
               scales::number(.x, 0.01))
      }
      x %>% list() %>% set_names(measure) 
    }))
  
  # Model for value box 3
  model <- personal_exposure %>% 
    mutate(exposure = if(input$exp == "CHIRPS"){
      exposure / 10
    } else if(input$exp == "CHIRTS"){
      exposure / 5
    } else if(input$exp == "NDVI"){
      exposure * 10
    } else if(input$exp == "UCDP"){
      exposure
    }) %>% 
    as_survey_design(weight = weight, id = cluster) %>% 
    survey::svyglm(
      whz ~ exposure + geo + 
        lz + hh_toilet + hh_electric + hh_drink + hh_natfloor + 
        mom_age + mom_edu +
        kid_agemo + kid_male + kid_diar + kid_fevr + birth_order, 
      design = .
    ) 
  
  output <- list(
    ntl = ntl_exposure, 
    ea = ea_exposure, 
    kid = personal_exposure, 
    model = model,
    value_boxes = value_boxes,
    admin1 = admin1_exposure,
    iso = iso
  )
})
```

Sample Analysis  
======================================================================

Row 
-----------------------------------------------------------------------

### Value 1

```{r vb1}
renderValueBox({
  vals <- out()$value_boxes
  valueBox(
    value = make_vb_icon("vb_exp", vals$coef$exposure), 
    icon =  meta[[out()$iso]][[input$exp]]$icon,
    color =  meta[[out()$iso]][[input$exp]]$color %>% paste0("96"),
    caption = paste(
      "<strong>Estimated Mean Exposure</strong><br>",
      "Interquartile range:", vals$`_q25`$exposure, "to",
      vals$`_q75`$exposure, "<br>",
      "Median value:",  vals$`_q50`$exposure
    )
  )
})
```

### Value 2 

```{r vb2}
renderValueBox({
  vals <- out()$value_boxes
  valueBox(
    value = make_vb_icon("vb_whz",  vals$coef$whz),
    icon = meta[[out()$iso]]$WHZ$icon,
    color = meta[[out()$iso]]$WHZ$color %>% paste0("50"),
    caption = paste(
      "<strong>Estimated Mean WHZ</strong><br>", 
      "Acutely malnourished:", vals$coef$acute, "<br>",
      "Severely malnourished:", vals$coef$severe
    )
  )
})
```

### Value 3 

```{r vb3}
renderValueBox({
  val <- out()$model %>% 
    broom::tidy(conf.int = TRUE) %>%
    filter(term == "exposure") %>% 
    mutate(sig = gtools::stars.pval(p.value))
  est <- val$estimate 
  sig <- val$`p.value` <= 0.05
  pos <- est > 0
  est <- est %>% scales::number(0.01) 
  corr <- meta[[out()$iso]]$corr
  
  if(sig){
    if(pos){
      corr <- corr$pos
    } else {
      corr <- corr$neg
    }
  } else {
    corr <- corr$insig 
  }

  valueBox(
    value =  make_vb_icon("vb_cor", est),
    icon = corr$icon ,
    color = corr$color,
    caption = meta[[out()$iso]][[input$exp]]$cap %>% 
      paste(corr$cap) %>% 
      paste(est)
  )
})

```

```{r modal_windows}
# Modal dialogue for value boxes

observeEvent(input$vb_exp, {
  showModal(
    modalDialog(
      title = input$exp,
      easyClose = TRUE,
      size = "l",
      footer = a(
        href = meta[[out()$iso]][[input$exp]]$source, 
        "Source:", meta[[out()$iso]][[input$exp]]$info
      ),
      HTML(
        meta[[out()$iso]][[input$exp]]$txt,
        "<strong>", paste0(min(input$time),"-", max(input$time)), 
        "months</strong> before the first interview for the <strong>",
        input$samp, "</strong>sample. <br><br>This value represents", 
        ifelse(
          input$exp == "UCDP", 
          paste(
            "all reported lethal conflicts in the admin 1 area where the",
            "child resides (including events that may have crossed borders)."
          ),
          paste(
            "a spatial average for all grid cells at 0.05 arc-degree",
            "resolution within a <strong>", input$space, "km</strong> radius",
            "from the center of the enumeration area where the child resides."
          )
        ),
        "<br><br>"
      ),
      out()$ea %>% 
        st_drop_geometry() %>% 
        mutate(exposure = exposure %>% scales::number(0.01)) %>% 
        rename(
          "DHS Cluster ID" = cluster,
          "Admin 1" = geo, 
          "Sampled children age 0-5" = n,
          "EA-level exposure" = exposure
        ) %>% 
        gt::gt() %>% 
        make_gt()
    )
  )
})


observeEvent(input$vb_whz, {
  showModal(modalDialog(
    title = "WHZ for Children Aged 0-5 Years",
    easyClose = TRUE,
    size = "l",
    footer = a(
      href = meta[[out()$iso]]$WHZ$source, 
      "Source:", meta[[out()$iso]]$WHZ$info
    ),
    HTML(
      "Difference between the child's weight and the median weight of a", 
      "healthy child of the same height and sex (defined by the WHO",
      "Global Database on Child Growth and Malnutrition), expressed in units",
      "equal to <strong>one standard deviation</strong>",
      "of the reference population's distribution (known as a Z-score).",
      "<br><br>The WHO database define <strong>acute malnutrition</strong>",
      "as a WHZ score <strong>below -2</strong>, and",
      "<strong>severe malnutrition</strong>",
      "as a WHZ score <strong>below -3</strong>.<br><br>",
      "Weighted mean WHZ scores for each Admin 1 derived from the <strong>", 
      input$samp, "</strong> sample are estimated below",
      "(cluster-robust 95% confidence intervals are shown).<br><br>"
    ),
    out()$kid %>% 
      as_survey_design(weight = weight, id = cluster) %>%
      group_by(geo) %>%
      summarise(Mean = survey_mean(whz, na.rm = TRUE, vartype = "ci")) %>%
      mutate(geo = case_when(
        str_detect(geo, "Southern Nations") ~ "SNNP",
        .default = geo
      )) %>%
      mutate(across(where(is.double), ~.x %>% scales::number(0.01))) %>% 
      rename(`Admin 1` = geo) %>% 
      rename_with(~.x %>% str_replace("Mean_", "95% CI "), matches("_")) %>% 
      gt::gt() %>% 
      make_gt()
  ))
})


observeEvent(input$vb_cor, {
  showModal(modalDialog(
    title = "Correlation",
    easyClose = TRUE,
    size = "l",
    footer = NULL,
    HTML(
      "The relationship between <strong>", input$exp, "</strong> and WHZ",
      "is modeled via linear regression with cluster-robust standard",
      "errors accounting for similarities between households in the same",
      "enumeration area.<br><br>Model controls are shown in the table",
      "below (Admin 1 geographies are included as fixed-effects, but",
      "not shown)." 
    ),
    out()$model %>% 
      tbl_regression(
        include = -c("geo"),
        show_single_row = where(is.logical),
        intercept = FALSE,
        label = c(
          exposure = input$exp, 
          lz = "Livelihood Zone (FEWSNET)",
          hh_electric = "Household has electricity access",
          hh_natfloor = "Household has natural flooring",
          hh_toilet = "Household toilet type",
          hh_drink = "Household drinking water type",
          mom_age = "Mother's age (years)",
          mom_edu = "Mother's education",
          kid_agemo = "Child's age (months)",
          kid_male = "Child is male",
          kid_diar = "Child had recent diarrhea",
          kid_fevr = "Child had recent fever",
          birth_order = "Child's place in birth order"
        )
      ) %>% 
      bold_labels() %>%
      add_significance_stars(
        hide_ci = FALSE,
        hide_se = TRUE
      ) %>%
      modify_header(
        update = list(
          label ~ " ",
          estimate = "**Estimate**"
        )
      ) %>%
      gtsummary::as_gt() %>% 
      gt::opt_interactive(
        use_search = TRUE, 
        use_highlight = TRUE,
        use_compact_mode = TRUE, 
        use_page_size_select = TRUE,
        page_size_default = 10,
        page_size_values =  c(5, 10, 15, 20)
      ) %>% 
      gt::render_gt()
  ))
})

```

Row
-----------------------------------------------------------------------

### 

```{r}
output$fig1 <- renderLeaflet({
  ea_col <- meta[[out()$iso]]$WHZ$color
  map <- leaflet(options = leafletOptions(minZoom = 5, maxZoom = 8)) %>% 
    addProviderTiles(provider = providers$CartoDB.PositronNoLabels)
  
  map <- if(input$exp == "UCDP"){
    pal <- colorNumeric(
      palette = colorRamp(c("transparent", meta[[out()$iso]]$UCDP$color), 
                          interpolate = "linear"),
      domain = out()$admin1$exposure
    )
    map %>% 
      addMarkers(
        clusterOptions = markerClusterOptions(
          iconCreateFunction=JS("function (cluster) {
            var childCount = cluster.getChildCount();
            var c = ' marker-cluster-';
            if (childCount < 5) {  c += 'small';
            } else if (childCount < 10) {
              c += 'medium';
            } else {
              c += 'large';
            }
            return new L.DivIcon({
              html: '<div><span>' +
              childCount + '</span></div>',
              className: 'marker-cluster' + c,
              iconSize: new L.Point(40, 40)
            });
          }"),
          showCoverageOnHover = FALSE,
          singleMarkerMode = TRUE
        ),
        data = out()$ntl %>%
          rowwise() %>%
          mutate(
            kbl = pick(
              -c(country, index, geometry, ADMIN_NAME)) %>%
              relocate(admin2, .before = source_article) %>%
              relocate(admin1 = admin, .before = admin2) %>%
              mutate(across(everything(), as.character)) %>%
              pivot_longer(everything()) %>%
              mutate(name = name %>% str_replace_all("_", " ") %>%
                       str_to_title) %>%
              knitr::kable(format = "html", col.names = NULL) 
          ),
        popup = ~kbl
      ) %>%
      addPolygons(
        data = out()$admin1,
        color = "black",
        weight = 1,
        fillColor = ~pal(out()$admin1$exposure),
        fillOpacity = 1
      )
  } else {
    map %>% 
      addRasterImage(
        colors = c("transparent", meta[[out()$iso]][[input$exp]]$color),
        out()$ntl
      ) 
  }
  map <- map %>%
    addPolygons(
      layerId = ~cluster,
      group = "eas",
      data = out()$ea,
      weight = 0,
      fillColor = "black",
      fillOpacity = 0.1,
      popup = ~paste(
        "<strong>Admin 1:</strong>", 
        geo, 
        "<br><strong>Exposure:</strong>",
        exposure %>% scales::number(0.01),
        "<br><strong style='color:",  
        ea_col,
        "'>Sampled Children:</strong>", 
        n
      )
    )
  if(!input$eas){map <- map %>% hideGroup("eas")}
  return(map)
})

leafletOutput("fig1")
```

### 

```{r}
renderPlot({
  if(input$exp == "NDVI"){
    jitter = 0.01
  } else if(input$exp == "UCDP"){
    jitter = 0.1
  } else if(input$exp == "CHIRTS"){
    jitter = 0.2
  } else {
    jitter = 0.5
  }
  fig2 <- out()$model %>% 
    predictions(
      type = "response", 
      variables = "exposure",
      newdata = datagrid()
    ) %>% 
    ggplot(aes(x = exposure, y = estimate)) + 
    geom_line(color = meta[[out()$iso]][[input$exp]]$color) +
    geom_ribbon(alpha = 0.2, aes(ymin = conf.low, ymax = conf.high)) +
    labs(y = "WHZ", x = paste("Exposure \n Jittered +/-", jitter))
  if(is.null(ea_clicked())){
    fig2 <- fig2 +
      geom_point(
        data = out()$model$data,
        mapping = aes(x = exposure, y = whz),
        color = "#00000005", size = 3, stroke = 0,
        position = position_jitter(width = jitter, height = 0, seed = 1)
      )
  } else {
    fig2 <- fig2 +
      geom_point(
        data = out()$model$data,
        aes(x = exposure, y = whz, color = cluster == ea_clicked()),
        size = 3, stroke = 0, show.legend = FALSE,
        position = position_jitter(width = jitter, height = 0, seed = 1),
      ) +
      scale_color_manual(values = c(
        "TRUE" = meta[[out()$iso]]$WHZ$color,
        "FALSE" = "#00000002"
      ))
  }
  
  if(input$exp == "CHIRPS"){
    fig2 <- fig2 + scale_x_continuous(labels = ~.x * 10)
  } else if(input$exp == "CHIRTS"){
    fig2 <- fig2 + scale_x_continuous(labels = ~.x * 5)
  } else if(input$exp == "NDVI"){
    fig2 <- fig2 + scale_x_continuous(labels = ~.x / 10)
  }
  
  return(fig2)
})
```




Differential Exposure 
============================================================================



### {data-width=300}

Exposure to climate and conflict factors may have a greater or lesser impact on malnutrition for children in different situations. Put another way: exposure may "interact" with other factors. 

In the figures below, **parallel lines** suggest that exposure tends to act the same way for children in each group. 

However, **intersecting lines** suggest differential exposure experienced by children in different groups. 

Select an "interaction" to test whether children in different groups are equally susceptible to exposure: 

<div style="margin-top: 2rem;">
```{r}
radioButtons(
  "interact", 
  label = NULL,
  choices = c(
    "Admin 1" = 'geo', 
    "Livelihood Zone" = 'lz', 
    "Household Toilet Type" = 'hh_toilet', 
    "Household has Electricity" = 'hh_electric', 
    "Household Drinking Water" = 'hh_drink', 
    "Household has Natural Flooring" = 'hh_natfloor',
    "Mother's age (years)" = 'mom_age', 
    "Mother's Education" = 'mom_edu', 
    "Child's age (months)" = 'kid_agemo', 
    "Child's Sex" = 'kid_male', 
    "Child had Recent Diarrhea" = 'kid_diar', 
    "Child had Recent Fever" = 'kid_fevr',
    "Child's place in birth order" = 'birth_order'
  ),
  selected = "kid_male"
)
```
</div>


### 

```{r}
renderPlot({
  interact <- paste("whz ~ exposure *", input$interact) 
  controls <- c( 'geo',  'lz', 'hh_toilet', 'hh_electric','hh_drink',
                 'hh_natfloor', 'mom_age', 'mom_edu', 'kid_agemo', 'kid_male',
                 'kid_diar', 'kid_fevr','birth_order') %>%
    keep(~.x != input$interact) %>%
    paste(collapse = " + ")
  full <- paste(interact, " + ", controls) %>% as.formula()
  out()$kid %>%
    mutate(kid_male = ifelse(kid_male, "Male", "Female")) %>% 
    as_survey_design(weight = weight, id = cluster) %>%
    survey::svyglm(full, design = .) %>%
    plot_predictions(condition = c("exposure", input$interact)) + 
    labs(color = NULL, fill = NULL, y = "WHZ", x = "Exposure")
})
```