Point estimates: CIFA, SPIFA
Compare observed data with point estimates from CIFA and SPIFA models.
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_cifa <- readRDS(file.path(path_modelled, "cifa-ipixuna-urban.rds"))
samples_spifa <- readRDS(file.path(path_modelled, "spifa-ipixuna-urban-3gp.rds"))
iter <- nrow(samples_spifa[["theta"]])Compute summaries from data
items_summary <- st_drop_geometry(fidata) |>
dplyr::select(matches("^item_[0-9]+_"), registro) |>
pivot_longer(-registro, names_to = "Item", values_to = "varvalue", names_ptypes = factor()) |>
group_by(Item) |>
summarise(`NA` = sum(is.na(varvalue)), `$\\pi$` = mean(varvalue, na.rm = TRUE)) |>
mutate(
Question = gsub("[A-D]-[0-9]+: ", "", clean_item_names(Item)),
Item = as.character(1:n()),
Section = map2(LETTERS[1:4], c(3, 4, 6, 5),
~ c(sprintf("\\multirow{%i}[3]{*}{%s}", .y, .x), rep("", .y - 1))) %>%
Reduce(c, .)
) |>
dplyr::select(Item, Section, Question, `NA`, `$\\pi$`) |>
print(n = Inf)#> # A tibble: 18 × 5
#> Item Section Question `NA` `$\\pi$`
#> <chr> <chr> <chr> <int> <dbl>
#> 1 1 "\\multirow{3}[3]{*}{A}" worried that food ends 0 0.565
#> 2 2 "" run out of food 0 0.52
#> 3 3 "" ate few food types 0 0.64
#> 4 4 "\\multirow{4}[3]{*}{B}" skipped a meal 0 0.305
#> 5 5 "" ate less than required 0 0.41
#> 6 6 "" hungry but did not eat 0 0.24
#> 7 7 "" at most one meal per day 0 0.26
#> 8 8 "\\multirow{6}[3]{*}{C}" ate few food types 25 0.486
#> 9 9 "" ate less than required 25 0.314
#> 10 10 "" decreased food quantity 25 0.36
#> 11 11 "" skipped a meal 25 0.229
#> 12 12 "" hungry but did not eat 25 0.2
#> 13 13 "" at most one meal per day 25 0.177
#> 14 14 "\\multirow{5}[3]{*}{D}" food just with farinha 0 0.165
#> 15 15 "" credit for eating 0 0.675
#> 16 16 "" borrowed food 0 0.14
#> 17 17 "" had meals at neighbors 0 0.175
#> 18 18 "" reduced meat or fish 0 0.535Model summaries
discrimination_summary <- summary(as_tibble(samples_cifa, iter/5, select = c("a"))) |>
bind_rows(summary(as_tibble(samples_spifa, iter/5, select = c("a"))), .id = "model") |>
select(model, Parameters, `50%`) |>
within({`50%`[`50%` == 0] <- NA}) |>
separate_wider_regex(Parameters,
patterns = c(".+\\[", Item = "[0-9]+", ",", factor = "[1-3]+", "\\]")) |>
pivot_wider(names_from = c("model", "factor"), values_from = `50%`,
names_glue = "M{model}-F{factor}")
difficulty_summary <-
# summary(as_tibble(samples_cifa, iter/5, select = c("c"))) |>
bind_rows(summary(as_tibble(samples_spifa, iter/5, select = c("c"))), .id = "model") |>
select(model, Parameters, `50%`) |>
separate_wider_regex(Parameters,
patterns = c(".+\\[", Item = "[0-9]+", "\\]")) |>
pivot_wider(names_from = "model", values_from = `50%`, names_glue = "c")
full_summary <- left_join(items_summary, discrimination_summary) |>
left_join(difficulty_summary)#> Joining with `by = join_by(Item)`
#> Joining with `by = join_by(Item)`full_summary#> # A tibble: 18 × 12
#> Item Section Question `NA` `$\\pi$` `M1-F1` `M1-F2` `M1-F3` `M2-F1` `M2-F2` `M2-F3` c
#> <chr> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 "\\multiro… worried… 0 0.565 NA 1.63 NA NA 1.77 NA 0.439
#> 2 2 "" run out… 0 0.52 NA NA 1.51 NA NA 1.88 0.216
#> 3 3 "" ate few… 0 0.64 1.68 NA NA 1.81 NA NA 1.46
#> 4 4 "\\multiro… skipped… 0 0.305 1.47 NA 1.02 1.92 NA 0.998 -0.628
#> 5 5 "" ate les… 0 0.41 0.880 NA 1.75 1.39 NA 1.50 -0.0743
#> 6 6 "" hungry … 0 0.24 1.26 NA 1.49 1.81 NA 1.54 -1.44
#> 7 7 "" at most… 0 0.26 1.57 NA NA 1.80 NA NA -0.551
#> 8 8 "\\multiro… ate few… 25 0.486 1.71 NA NA 1.88 NA NA 0.588
#> 9 9 "" ate les… 25 0.314 1.90 NA NA 2.23 NA NA -0.360
#> 10 10 "" decreas… 25 0.36 2.17 NA NA 2.53 NA NA -0.0331
#> 11 11 "" skipped… 25 0.229 1.99 NA NA 2.52 NA NA -1.05
#> 12 12 "" hungry … 25 0.2 2.11 NA NA 2.59 NA NA -1.33
#> 13 13 "" at most… 25 0.177 1.93 NA NA 2.43 NA NA -1.51
#> 14 14 "\\multiro… food ju… 0 0.165 0.343 NA 1.25 0.626 NA 1.28 -1.62
#> 15 15 "" credit … 0 0.675 NA 0.716 NA NA 0.796 NA 0.629
#> 16 16 "" borrowe… 0 0.14 NA 1.44 NA NA 1.61 NA -1.90
#> 17 17 "" had mea… 0 0.175 NA 0.982 NA NA 1.04 NA -1.26
#> 18 18 "" reduced… 0 0.535 1.29 NA NA 1.43 NA NA 0.751Create and save latex table
# convert to latex
full_table <- xtable(
full_summary,
align = c("c", "c", "c", "l", "c", "c", "c", "c", "c", "c", "c", "c", "r"),
caption = paste0(
"Summary of the food insecurity items: ",
"i) the number of missing values ($\\#$NA) and the proportion of endorsement ",
"($\\pi$) are shown for the descriptive analysis, ",
"ii) the posterior median of the discrimination parameters $\\{\\m{\\hat{A}}_{\\cdot 1}, \\m{\\hat{A}}_{\\cdot 2}, \\m{\\hat{A}}_{\\cdot 3}\\}$ are shown for the ",
"confirmatory factor analysis (CIFA), and ",
"iii) the posterior median of the discrimination and easiness parameters ",
" $\\{\\m{\\hat{A}}_{\\cdot 1}, \\m{\\hat{A}}_{\\cdot 2}, \\m{\\hat{A}}_{\\cdot 3}, \\hat{\\ve{c}}\\}$ are shown for the ",
"spatial item factor analysis (SPIFA)."),
label = "tab:items"
)
# improve row names and add separators
addtorow <- list()
addtorow$pos <- list(0, 0, nrow(full_table))
addtorow$command <- c(
"\\toprule[.15em] \\multirow{2}[3]{*}{\\textbf{Item}} & \\multirow{2}[3]{*}{\\textbf{Section}} & \\multirow{2}[3]{*}{\\textbf{Question}} & \\multicolumn{2}{c}{\\textbf{Descriptive}} & \\multicolumn{3}{c}{\\textbf{CIFA}} & \\multicolumn{4}{c}{\\textbf{SPIFA}} \\\\\n \\cmidrule(rl){4-5} \\cmidrule(rl){6-8} \\cmidrule(rl){9-12}",
" & & & $\\#$NA & $\\pi$ & ${\\m{\\hat{A}}_{\\cdot 1}}$ & ${\\m{\\hat{A}}_{\\cdot 2}}$ & ${\\m{\\hat{A}}_{\\cdot 3}}$ & ${\\m{\\hat{A}}_{\\cdot 1}}$ & ${\\m{\\hat{A}}_{\\cdot 2}}$ & ${\\m{\\hat{A}}_{\\cdot 3}}$ & \\multicolumn{1}{c}{$\\hat{\\ve{c}}$} \\\\\n \\midrule[.1em]",
"\\bottomrule[.15em]"
)
# print latex with additional formating
print(full_table,
file = file.path(path_fig, "results-summary-cifa-spifa.tex"),
caption.placement = "top",
include.rownames = FALSE,
size = "\\scriptsize",
sanitize.text.function = identity,
add.to.row = addtorow,
include.colnames = FALSE,
hline.after = c(3, 7, 13),
NA.string = "$\\cdot$"
)Time to execute the task
Only useful when executed with Rscript.
proc.time()#> user system elapsed
#> 3.951 0.411 4.046