Food insecurity maps: full model
Maps of the prediction of the latent abilities and the spatial component.
Load required libraries and data
rm(list = ls())
library(day2day)
library(spifa)
library(xtable)
library(tibble)
library(sf)#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUElibrary(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, unionlibrary(tidyr)
library(purrr)
library(ggplot2)
library(patchwork)
path_main <- git_path()
path_data <- file.path(path_main, "data")
path_raw <- file.path(path_data, "raw")
path_processed <- file.path(path_data, "processed")
path_modelled <- file.path(path_data, "modelled")
path_fig <- file.path(path_data, "figures")
path_src <- file.path(path_main, "src")
source(file.path(path_src, "itemnames.R"))
source(file.path(path_src, "ggtheme-publication.R"))
source(file.path(path_src, "maps.R"))
fidata <- file.path(path_processed, "fi-items-ipixuna-urban.gpkg") |>
st_read(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 84ipi_osm <- readRDS(file.path(path_processed, "ipixuna-osm.rds"))
pred_spifa <- readRDS(file.path(path_modelled, "spifa-prediction-ipixuna-urban-3gp.rds"))
pred_spifa_sp <- readRDS(file.path(path_modelled, "spifa-predictionsp-ipixuna-urban-3gp.rds"))Get summary
# custom function
fun_sum <- function(data, regex, fun){
dplyr::select(data, matches(regex)) |>
st_drop_geometry() |>
apply(1, fun)
}
# compute summaries
pred_summary <- pred_spifa["id"]
pred_summary[paste0("F", 1:3, "-exceedance")] <- lapply(c("^F1-", "^F2-", "^F3-"),
function(z) fun_sum(pred_spifa, z, function(x) sum(x > 0) / length(x)))
pred_summary[paste0("F", 1:3, "-median")] <- lapply(c("^F1-", "^F2-", "^F3-"),
function(z) fun_sum(pred_spifa, z, median))
# compute summaries
pred_summary_sp <- pred_spifa_sp["id"]
pred_summary_sp[paste0("F", 1:3, "-exceedance")] <- lapply(c("^F1-", "^F2-", "^F3-"),
function(z) fun_sum(pred_spifa_sp, z, function(x) sum(x > 0) / length(x)))
pred_summary_sp[paste0("F", 1:3, "-median")] <- lapply(c("^F1-", "^F2-", "^F3-"),
function(z) fun_sum(pred_spifa_sp, z, median))Detect high insecure areas
mylabs <- c(
# "F1-exceedance" = "(F1)~Severe~\n~food~insecurity~including~children",
"F1-exceedance" = "atop('(F1) Severe food insecurity', 'including children')",
# "F2-exceedance" = "(F2)~Anxiety~and~relying~on~social~relations",
"F2-exceedance" = "atop('(F2) Anxiety and relying on', 'social relations')",
"F3-exceedance" = "atop('(F3) Adults eating less and', 'going hungry')"
)
pred_summary_sp |>
tidyr::pivot_longer(matches("exceedance$"), names_to = "factor_exceed") |>
mutate(factor_exceed = mylabs[factor_exceed]) |>
mutate(value = ifelse(value >= 0.75, 1, NA)) |>
fi_ipi_exceed(osm = ipi_osm, title = "title", license = TRUE) +
theme_map_Publication2(base_size = 9) +
theme(legend.position = "none")#> Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
#> ℹ Please use the `linewidth` argument instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
Text to add
text_df <- tibble(
# text = c(toupper(letters[1:10])),
text = c("C2", "D2", "B", "A", "C", "F2", "G2", "D", "I2", "J2"),
x1 = c(-71.695, -71.69, -71.6975, -71.6837, -71.6915, -71.686, -71.698, -71.694,
-71.6875, -71.6875),
y1 = c(-7.045, -7.045, -7.051, -7.0585, -7.0488, -7.048, -7.057, -7.045, -7.0452,
-7.0558),
x2 = c(-71.6975, -71.6845, -71.6993, -71.6817, -71.699, -71.683, -71.699, -71.6975,
-71.684, -71.682),
y2 = c(-7.045, -7.044, -7.048, -7.055, -7.049, -7.049, -7.06, -7.045, -7.045,
-7.055),
factor = c(rep("F1-median", 4), rep("F2-median", 3), rep("F3-median", 3)),
factor_exceed = c(rep(mylabs[[1]], 4), rep(mylabs[[2]], 3),
rep(mylabs[[3]], 3)),
) |>
dplyr::filter(text %in% toupper(letters[1:4])) |>
st_as_sf(coords = c("x1", "y1"), crs = st_crs(4326)) |>
st_transform(crs = st_crs(3857))
text_df[c("x1", "y1")] <- st_coordinates(text_df)
text_df <- text_df |>
st_set_geometry(NULL) |>
st_as_sf(coords = c("x2", "y2"), crs = st_crs(4326)) |>
st_transform(crs = st_crs(3857))
text_df[c("x2", "y2")] <- st_coordinates(text_df)
text_df <- text_df |> st_set_geometry(NULL)
text_df#> # A tibble: 4 × 7
#> text factor factor_exceed x1 y1 x2 y2
#> * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 B F1-median atop('(F1) Severe food insecurity', 'including ch… -7.98e6 -7.87e5 -7.98e6 -7.87e5
#> 2 A F1-median atop('(F1) Severe food insecurity', 'including ch… -7.98e6 -7.88e5 -7.98e6 -7.87e5
#> 3 C F2-median atop('(F2) Anxiety and relying on', 'social relat… -7.98e6 -7.87e5 -7.98e6 -7.87e5
#> 4 D F3-median atop('(F3) Adults eating less and', 'going hungry… -7.98e6 -7.86e5 -7.98e6 -7.86e5Visualize prediction of latent abilities
fi_ipi(pred_summary, `F1-median`, ipi_osm, title = mylabs[[1]], text_df = text_df) +
fi_ipi(pred_summary, `F2-median`, ipi_osm, title = mylabs[[2]], text_df = text_df, axis.y = FALSE) +
fi_ipi(pred_summary, `F3-median`, ipi_osm, title = mylabs[[3]], text_df = text_df, axis.y = FALSE,
license = TRUE) &
theme_map_Publication(base_size = 9) &
theme(strip.text.x = element_text(margin=margin(l=0, r=0, t=2, b=2)))#> # A tibble: 2 × 7
#> text factor factor_exceed x1 y1 x2 y2
#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 B F1-median atop('(F1) Severe food insecurity', 'including ch… -7.98e6 -7.87e5 -7.98e6 -7.87e5
#> 2 A F1-median atop('(F1) Severe food insecurity', 'including ch… -7.98e6 -7.88e5 -7.98e6 -7.87e5#> Warning: `aes_()` was deprecated in ggplot2 3.0.0.
#> ℹ Please use tidy evaluation idioms with `aes()`
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.#> # A tibble: 1 × 7
#> text factor factor_exceed x1 y1 x2 y2
#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 C F2-median atop('(F2) Anxiety and relying on', 'social relat… -7.98e6 -7.87e5 -7.98e6 -7.87e5
#> # A tibble: 1 × 7
#> text factor factor_exceed x1 y1 x2 y2
#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 D F3-median atop('(F3) Adults eating less and', 'going hungry… -7.98e6 -7.86e5 -7.98e6 -7.86e5
ggsave(file.path(path_fig, "osmmap-theta-median.pdf"), width = 7, height = 3.5)
pred_summary |>
tidyr::pivot_longer(matches("exceedance$"), names_to = "factor_exceed") |>
mutate(factor_exceed = mylabs[factor_exceed]) |>
within({value[value < 0.5] <- NA}) |>
fi_ipi_exceed(osm = ipi_osm, title = "title", text_df = text_df, license = TRUE) +
theme_map_Publication2(base_size = 9)
ggsave(file.path(path_fig, "osmmap-theta-exceedance.pdf"), width = 7, height = 3.6)Visualize prediction of the spatial component of latent abilities
fi_ipi(pred_summary_sp, `F1-median`, ipi_osm, title = mylabs[[1]], text_df = text_df) +
fi_ipi(pred_summary_sp, `F2-median`, ipi_osm, title = mylabs[[2]], text_df = text_df, axis.y = FALSE) +
fi_ipi(pred_summary_sp, `F3-median`, ipi_osm, title = mylabs[[3]], text_df = text_df, axis.y = FALSE,
license = TRUE) &
theme_map_Publication(base_size = 9) &
theme(strip.text.x = element_text(margin=margin(l=0, r=0, t=2, b=2)))#> # A tibble: 2 × 7
#> text factor factor_exceed x1 y1 x2 y2
#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 B F1-median atop('(F1) Severe food insecurity', 'including ch… -7.98e6 -7.87e5 -7.98e6 -7.87e5
#> 2 A F1-median atop('(F1) Severe food insecurity', 'including ch… -7.98e6 -7.88e5 -7.98e6 -7.87e5
#> # A tibble: 1 × 7
#> text factor factor_exceed x1 y1 x2 y2
#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 C F2-median atop('(F2) Anxiety and relying on', 'social relat… -7.98e6 -7.87e5 -7.98e6 -7.87e5
#> # A tibble: 1 × 7
#> text factor factor_exceed x1 y1 x2 y2
#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 D F3-median atop('(F3) Adults eating less and', 'going hungry… -7.98e6 -7.86e5 -7.98e6 -7.86e5
ggsave(file.path(path_fig, "osmmap-thetaspatial-median.pdf"), width = 7, height = 3.6)
pred_summary_sp |>
tidyr::pivot_longer(matches("exceedance$"), names_to = "factor_exceed") |>
mutate(factor_exceed = mylabs[factor_exceed]) |>
within({value[value < 0.5] <- NA}) |>
fi_ipi_exceed(osm = ipi_osm, title = "title", text_df = text_df, license = TRUE) +
theme_map_Publication2(base_size = 9)
ggsave(file.path(path_fig, "osmmap-thetaspatial-exceedance.pdf"), width = 7, height = 3.6)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 17.383 1.039 17.478