SPIFA Urban Ipixuna MCMC chains: Traceplots and diagnostics

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")
path_fig <- file.path(path_data, "figures")
path_src <- file.path(path_main, "src")

source(file.path(path_src, "ggtheme-publication.R"))
source(file.path(path_src, "mcmc-traceplot-random.R"))
source(file.path(path_src, "mcmc-diagnostics.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 <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-3gp.rds"))
samples_dry <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-dry-3gp.rds"))
samples_wet <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-wet-3gp-restricted.rds"))

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

MCMC traceplots for full model

ggthemr::ggthemr_reset()
pal_aux <- ggthemr:::load_palette.character("fresh")
ggthemr::ggthemr(pal_aux, layout = "clear")

set.seed(2)
plot_mcmc(samples, "c", nshow = 3)

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

ggsave(file.path(path_fig, "mcmc-trace-c.pdf"), width = 7, height = 3)

plot_mcmc(samples, "a", nshow = 5, mat = TRUE)

ggsave(file.path(path_fig, "mcmc-trace-a.pdf"), width = 7, height = 3.5)

set.seed(1)
plot_mcmc(samples, "theta", nshow = 5, mat = TRUE)

ggsave(file.path(path_fig, "mcmc-trace-theta.pdf"), width = 7, height = 3.5)

plot_mcmc(samples, "corr", "R", mat = TRUE)

ggsave(file.path(path_fig, "mcmc-trace-corr.pdf"), width = 7, height = 3)

plot_mcmc(samples, "mgp_sd", mat = TRUE)

ggsave(file.path(path_fig, "mcmc-trace-mgp-sd.pdf"), width = 7, height = 3)

plot_mcmc(samples, "mgp_phi", "phi", mat = TRUE)

ggsave(file.path(path_fig, "mcmc-trace-mgp-phi.pdf"), width = 7, height = 3)

MCMC convergence for full model, dry and wet

diag_labels <- c("rhat" = "hat(R)", "ess_bulk" = "ESS~Bulk", "ess_tail" = "ESS~Tail")

out_full <- spifa_diagnostics(samples, burnin = iter / 5, thin = 1) |>
    pivot_longer(rhat:ess_tail) |>
    mutate(name = factor(name, levels = names(diag_labels), labels = diag_labels))

out_wet <- spifa_diagnostics(samples_wet, burnin = iter / 5, thin = 1) |>
    pivot_longer(rhat:ess_tail) |>
    mutate(name = factor(name, levels = names(diag_labels), labels = diag_labels))

out_dry <- spifa_diagnostics(samples_dry, burnin = iter / 5, thin = 1) |>
    pivot_longer(rhat:ess_tail) |>
    mutate(name = factor(name, levels = names(diag_labels), labels = diag_labels))

bind_rows(out_full, out_wet, out_dry, .id = "model") |>
    mutate(model = factor(model, levels = 1:3, labels = c("Full", "Season: Rainy", "Season: Dry"))) |>
    ggplot() +
    geom_histogram(aes(value), color = 1, linewidth = rel(0.1), bins = 50) +
    facet_grid(model ~ name, scales = "free", labeller = label_parsed) +
    theme_bw() +
    labs(x = "Value", y = "Frequency")

#> Warning: Removed 324 rows containing non-finite outside the scale range (`stat_bin()`).

ggsave(file.path(path_fig, "mcmc-diagnostics-full-wet-dry.pdf"), width = 7, height = 6)

#> Warning: Removed 324 rows containing non-finite outside the scale range (`stat_bin()`).

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  50.518   1.365  51.010