Urban Ipixuna: Descriptive analysis of items

Load required libraries and data

rm(list = ls())
library(day2day)
library(ggplot2)
library(sf)

#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

library(dplyr)

#> 
#> Attaching package: 'dplyr'

#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag

#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

path_main <- git_path()
path_data <- file.path(path_main, "data")
path_processed <- file.path(path_data, "processed")

fidata <- st_read(file.path(path_processed, "fi-items-ipixuna-urban.gpkg"), as_tibble = TRUE)

#> Reading layer `fi-items-ipixuna-urban' from data source 
#>   `/home/rstudio/documents/projects/food-insecurity-mapping/data/processed/fi-items-ipixuna-urban.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 200 features and 36 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -71.70038 ymin: -7.06058 xmax: -71.68109 ymax: -7.03724
#> Geodetic CRS:  WGS 84

Food insecurity score

First, we show the empirical distribution of the food insecurity score for all the households. The score between household with an without children can not be compared given that some items can be answered only for families with children.

# Functions to change plot order of bars, change levels in underlying factor
as_count <- function(x) factor(x, levels = 0:max(x, na.rm = TRUE))
ggplot(fidata, aes(as_count(fi_score))) +
    geom_bar() +
    labs(x = "score", title = "Absolute frequency of food insecurity score") +
    facet_wrap(~ any_children, ncol = 1)

Food insecurity proportion

ggplot(fidata, aes(fi_proportion)) +
    geom_histogram(binwidth = 0.05, color = "white", center = 0.025) +
    labs(x = "proportion", title = "Absolute frequency of food insecurity proportion") +
    facet_wrap(~ any_children, ncol = 1)

fi_proportion_by_children <- fidata |>
    st_set_geometry(NULL) |>
    group_by(any_children) |>
    summarise(mean = mean(fi_proportion))

ggplot(fidata, aes(fi_proportion)) +
    geom_histogram(aes(y = stat(width * density)), binwidth = 0.05,
                   color = "white", center = 0.025) +
    geom_vline(aes(xintercept = mean), fi_proportion_by_children, col = 2) +
    labs(x = "proportion", title = "Relative frequency of food insecurity proportion",
         y = "relative frequency") +
    facet_wrap(~ any_children, ncol = 1)

#> Warning: `stat(width * density)` was deprecated in ggplot2 3.4.0.
#> ℹ Please use `after_stat(width * density)` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

fi_proportion_by_season <- fidata |>
    st_set_geometry(NULL) |>
    group_by(season) |>
    summarise(mean = mean(fi_proportion))

ggplot(fidata, aes(fi_proportion)) +
    geom_histogram(aes(y = stat(width * density)), binwidth = 0.05,
                   color = "white", center = 0.025) +
    geom_vline(aes(xintercept = mean), fi_proportion_by_season, col = 2) +
    labs(x = "proportion", title = "Relative frequency of food insecurity proportion",
         y = "relative frequency") +
    facet_wrap(~ season, ncol = 1)

fi_proportion_by_season_children <- fidata |>
    st_set_geometry(NULL) |>
    group_by(season, any_children) |>
    summarise(mean = mean(fi_proportion)) |>
    ungroup()

#> `summarise()` has grouped output by 'season'. You can override using the `.groups` argument.

ggplot(fidata, aes(fi_proportion)) +
    geom_histogram(aes(y = stat(width * density)), binwidth = 0.05,
                   color = "white", center = 0.025) +
    geom_vline(aes(xintercept = mean), fi_proportion_by_season_children, col = 2) +
    labs(x = "proportion", title = "Relative frequency of food insecurity proportion",
         y = "relative frequency") +
    facet_wrap(~ season + any_children, ncol = 1)

Summary of Items by season

Descriptive analysis of items by season.

fi_items <- fidata |>
    st_set_geometry(NULL) |>
    dplyr::select(matches("^item"), registro, season) %>%
    bind_rows(., .) |>
    within(season[(nrow(fidata)+1):(2*nrow(fidata))] <- "all") |>
    tidyr::pivot_longer(-c(registro, season),
                        names_to = c("item", "section", "question"),
                        names_pattern = "item_([[:digit:]]+)_([[:alpha:]]+)_(.+)",
                        values_to = "varvalue") |>
    mutate(item = as.integer(item), question = gsub("_", " ", question))


fi_items |>
    group_by(season, item, section, question) |>
    summarise(
        n_na = sum(is.na(varvalue)),
        mean = mean(varvalue, na.rm = TRUE),
        median = median(varvalue, na.rm = TRUE),
        variance = var(varvalue, na.rm = TRUE)
    ) |>
    ungroup() |>
    tidyr::pivot_longer(n_na:variance, names_to = "vars", values_to = "value") |>
    mutate(vars = paste0(vars, "-", season)) |>
    dplyr::select(- season) |>
    tidyr::pivot_wider(names_from = vars, values_from = value) |>
    relocate(matches("^n_na"), matches("^mean"), matches("^median"), matches("^variance"),
             .after = question) |>
    knitr::kable()

#> `summarise()` has grouped output by 'season', 'item', 'section'. You can override using the
#> `.groups` argument.

fi_items |>
    group_by(season, item, section, question) |>
    filter(season != "all") |>
    summarise(prop = mean(varvalue, na.rm = TRUE)) |>
    mutate(item = factor(paste(item, section, sep = "-"))) |>
    ggplot() +
        geom_col(aes(item, prop, fill = season), position = position_dodge()) +
        theme_bw()

#> `summarise()` has grouped output by 'season', 'item', 'section'. You can override using the
#> `.groups` argument.

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>   2.725   0.394   2.806