SPIFA Urban Ipixuna Season Missing: summarise results

Load required libraries and data

rm(list = ls())
library(day2day)
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(mirt)

#> Loading required package: stats4

#> Loading required package: lattice

library(spifa)
library(tidyr)
library(sf)

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

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

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

samples1 <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-miss-3gp.rds"))
samples2 <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-miss-dry-3gp.rds"))
samples3 <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-miss-wet-3gp.rds"))

iter <- nrow(samples1[["theta"]])
thin <- 10

MCMC convergence

list(samples1, samples2, samples3) |>
    map(~ gather.spifa(as_tibble(.x, 0, 1, select = "c"))) |>
    bind_rows(.id = "Model") |>
    ggplot(aes(iteration, Value, group = Parameters, col = Parameters)) +
        geom_path(alpha = 0.6, linewidth = rel(0.1)) +
        facet_wrap(~ Model, ncol = 1)

list(samples1, samples2, samples3) |>
    map(~ gather.spifa(as_tibble(.x, 0, 1, select = "a"))) |>
    bind_rows(.id = "Model") |>
    ggplot(aes(iteration, Value, group = Parameters, col = Parameters)) +
        geom_path(alpha = 0.6, linewidth = rel(0.1)) +
        facet_wrap(~ Model, ncol = 1) +
        theme(legend.position = "none")

list(samples1, samples2, samples3) |>
    map(~ gather.spifa(select(as_tibble(.x, 0, 1, select = "theta"), 1:10))) |>
    bind_rows(.id = "Model") |>
    ggplot(aes(iteration, Value, group = Parameters, col = Parameters)) +
        geom_path(alpha = 0.6, linewidth = rel(0.1)) +
        facet_wrap(~ Model, ncol = 1)

list(samples1, samples2, samples3) |>
    map(~ gather.spifa(as_tibble(.x, 0, 1, select = "corr"))) |>
    bind_rows(.id = "Model") |>
    ggplot(aes(iteration, Value, group = Parameters, col = Parameters)) +
        geom_path(alpha = 0.6, linewidth = rel(0.1)) +
        facet_wrap(~ Model, ncol = 1)

list(samples1, samples2, samples3) |>
    map(~ gather.spifa(as_tibble(.x, 0, 1, select = "mgp_sd"))) |>
    bind_rows(.id = "Model") |>
    ggplot(aes(iteration, Value, group = Parameters, col = Parameters)) +
        geom_path(alpha = 0.6, linewidth = rel(0.1)) +
        facet_wrap(~ Model, ncol = 1)

list(samples1, samples2, samples3) |>
    map(~ gather.spifa(as_tibble(.x, 0, 1, select = "mgp_phi"))) |>
    bind_rows(.id = "Model") |>
    ggplot(aes(iteration, Value, group = Parameters, col = Parameters)) +
        geom_path(alpha = 0.6, linewidth = rel(0.1)) +
        facet_wrap(~ Model, ncol = 1)

Credible intervals

list(samples1, samples2, samples3) |>
    map(~ summary(as_tibble(.x, iter/2, select = "a"))) |>
    bind_rows(.id = "Model") |>
    gg_errorbarh(sorted = FALSE) +
    facet_wrap(~ Model)

list(samples1, samples2, samples3) |>
    map(~ summary(as_tibble(.x, iter/2, select = "c"))) |>
    bind_rows(.id = "Model") |>
    gg_errorbarh(sorted = FALSE) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    facet_wrap(~ Model) +
    theme_bw() +
    theme(legend.position = "bottom")

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  80.322   3.316  98.850