SPIFA Credible intervals: full, missing and seasonal

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 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

library(tidyr)
library(purrr)
library(ggplot2)

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"))

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 84

samples_spifa <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-3gp.rds"))
samples_wet <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-wet-3gp-restricted.rds"))
samples_dry <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-dry-3gp.rds"))
samples_miss <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-miss-3gp.rds"))
iter <- nrow(samples_spifa[["theta"]])

Custom function

Create function to plot credible intervals for several models.

ci_intervals <- function(models, param = "c", labels = c("Full", "Rainy", "Dry"),
                         colors = c(rgb(1, 0.5, 0.1), 4, 2), position = position_dodge(0.7)) {
    models |>
        map(~ summary(as_tibble(.x, iter/2, select = param))) |>
        bind_rows(.id = "Model") |>
        mutate(Model = factor(Model, seq_along(labels), labels )) |>
        filter(`50%` != 0) |>
        ggplot(aes(Parameters, `50%`)) +
    geom_hline(yintercept = 0, linetype = 2) +
    geom_errorbar(aes(ymin = `2.5%`, ymax = `97.5%`, group = Model, linetype = "95%"), width = 0,
        position = position) +
    geom_errorbar(aes(ymin = `10%`, ymax = `90%`, color = Model, linetype = "80%"), width = 0,
        linewidth = rel(1.5), position = position) +
    geom_point(aes(group = Model), size = rel(1), position = position) +
    labs(linetype = "Credible interval:", y = "") +
    scale_linetype_manual(values = c(1, 1)) +
    scale_color_manual(values = colors) +
    labs(color = "Model:") +
    theme_bw() +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        legend.position = "bottom")
}

Compare full, dry and wet models

ci_c <- ci_intervals(list(samples_spifa, samples_wet, samples_dry), "c") +
    theme_Publication(base_size = 10)

#> Warning: The `size` argument of `element_line()` 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.

ci_c

ggsave(file.path(path_fig, "results-ci-easiness-full_season.pdf"), ci_c, width = 7, height = 4)

ci_a <- ci_intervals(list(samples_spifa, samples_wet, samples_dry), "a") +
    theme_Publication(base_size = 10) +
    theme(axis.text.x = element_text(size = rel(0.8)))
ci_a

ggsave(file.path(path_fig, "results-ci-discrimination-full_season.pdf"), ci_a, width = 7, height = 4)

Compare full and missing models

ci_c <- ci_intervals(list(samples_spifa, samples_miss), "c", c("Full", "Missing"),
    colors = c(rgb(1, 0.5, 0.1), "gray60"), position = position_dodge(0.4)) +
    theme_Publication(base_size = 10)
ci_c

ggsave(file.path(path_fig, "results-ci-easiness-full_missing.pdf"), ci_c, width = 7, height = 4)

ci_a <- ci_intervals(list(samples_spifa, samples_miss), "a", c("Full", "Missing"),
    colors = c(rgb(1, 0.5, 0.1), "gray60"), position = position_dodge(0.4)) +
    theme_Publication(base_size = 10) +
    theme(axis.text.x = element_text(size = rel(0.8)))
ci_a

ggsave(file.path(path_fig, "results-ci-discrimination-full_missing.pdf"), ci_a, width = 7, height = 4)

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>   7.154   0.545   7.405