SPIFA estimates: all, rainy and dry data
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)
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"))
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 84samples_spifa <- 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_spifa[["theta"]])Compute summaries from data
summ_all <- summary(as_tibble(samples_spifa, iter/5, select = c("mgp_sd", "mgp_phi", "corr")))
summ_dry <- summary(as_tibble(samples_dry, iter/5, select = c("mgp_sd", "mgp_phi", "corr")))
summ_wet <- summary(as_tibble(samples_wet, iter/5, select = c("mgp_sd", "mgp_phi", "corr")))
map(list(full = samples_spifa, rainy = samples_wet, dry = samples_dry),
~ summary(as_tibble(.x, iter/5, select = c("mgp_sd")))) |>
bind_rows(.id = "model")#> # A tibble: 9 × 7
#> model Parameters `2.5%` `10%` `50%` `90%` `97.5%`
#> <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 full T[1,1] 0.292 0.347 0.468 0.611 0.698
#> 2 full T[2,2] 0.0885 0.117 0.199 0.324 0.413
#> 3 full T[3,3] 0.152 0.190 0.289 0.411 0.493
#> 4 rainy T[1,1] 0.270 0.338 0.487 0.687 0.823
#> 5 rainy T[2,2] 0.0767 0.101 0.165 0.269 0.336
#> 6 rainy T[3,3] 0.142 0.186 0.295 0.444 0.546
#> 7 dry T[1,1] 0.207 0.269 0.435 0.645 0.771
#> 8 dry T[2,2] 0.0825 0.109 0.184 0.307 0.391
#> 9 dry T[3,3] 0.136 0.177 0.276 0.414 0.507map(list(full = samples_spifa, rainy = samples_wet, dry = samples_dry),
~ summary(as_tibble(.x, iter/5, select = c("mgp_phi")))) |>
bind_rows(.id = "model")#> # A tibble: 9 × 7
#> model Parameters `2.5%` `10%` `50%` `90%` `97.5%`
#> <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 full mgp_phi[1] 127. 153. 216. 303. 365.
#> 2 full mgp_phi[2] 45.1 55.9 83.0 121. 148.
#> 3 full mgp_phi[3] 43.9 52.8 76.8 112. 136.
#> 4 rainy mgp_phi[1] 115. 139. 197. 278. 331.
#> 5 rainy mgp_phi[2] 44.1 55.0 80.7 119. 147.
#> 6 rainy mgp_phi[3] 46.6 57.0 83.8 122. 148.
#> 7 dry mgp_phi[1] 98.2 122. 178. 259. 314.
#> 8 dry mgp_phi[2] 43.9 54.3 79.6 117. 144.
#> 9 dry mgp_phi[3] 42.8 52.4 76.2 113. 139.map(list(full = samples_spifa, rainy = samples_wet, dry = samples_dry),
~ summary(as_tibble(.x, iter/5, select = c("corr")))) |>
bind_rows(.id = "model")#> # A tibble: 9 × 7
#> model Parameters `2.5%` `10%` `50%` `90%` `97.5%`
#> <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 full Corr[2,1] 0.615 0.667 0.756 0.838 0.876
#> 2 full Corr[3,1] 0.712 0.766 0.853 0.921 0.949
#> 3 full Corr[3,2] 0.637 0.697 0.796 0.877 0.912
#> 4 rainy Corr[2,1] 0.598 0.677 0.799 0.894 0.932
#> 5 rainy Corr[3,1] 0.587 0.671 0.797 0.896 0.936
#> 6 rainy Corr[3,2] 0.450 0.548 0.713 0.843 0.896
#> 7 dry Corr[2,1] 0.415 0.491 0.631 0.753 0.815
#> 8 dry Corr[3,1] 0.531 0.617 0.764 0.876 0.919
#> 9 dry Corr[3,2] 0.472 0.568 0.717 0.837 0.885Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 4.006 0.398 4.088