Semi-variogram of the latent factors with CIFA

Compute summaries for the latent abilities.

theta_summary <- as_tibble(samples, iter/thin/2, select = "theta") |>
  summary() |>
  tidyr::extract(Parameters, c("subject", "factor"),
                 "\\[([[:digit:]]+),([[:digit:]]+)\\]", convert = TRUE)

theta_median <- theta_summary |>
    dplyr::select(subject, factor, `50%`) |>
    pivot_wider(names_from = factor, values_from = `50%`, names_prefix = "F") |>
    dplyr::select(-subject)

factorsdata <- fidata |>
    bind_cols(theta_median) |>
    dplyr::select(registro, any_children, matches("^F[1-3]$")) |>
    tidyr::pivot_longer(matches("^F[1-3]$"), names_to = "factor_label",
                        values_to = "factor_value")

Variogram of latent factors

Computing extent of coordinates.

extent <- factorsdata |>
  st_distance() |>
  max() |>
  as.numeric()

Computing variograms for each factor.

factorsdata[c("longitude_3857", "latitude_3857")] <-
    st_coordinates(st_transform(factorsdata, 3857))

# empirical variogram
gg_variogs <- factorsdata |>
  st_set_geometry(NULL) |>
  group_by(factor_label) |>
  tidyr::nest() |>
  mutate(variogs = map(data,
    ~ gstat::variogram(factor_value ~ 1, ~ longitude_3857 + latitude_3857, .,
                       # cutoff = extent * 0.7, width = extent / 200)
                       cutoff = extent * 0.5, width = extent / 40)
    )
  ) |>
  ungroup()

# prepare empirical variogram
gg_variog_emp <- dplyr::select(gg_variogs, factor_label, variogs) |>
  tidyr::unnest(col = variogs)

# variogram fitting
gg_variogs <- gg_variogs |>
  mutate(psill = 1, nugget = 0.2, range = c(50, 15, 60)) |>
  mutate(model = pmap(list(psill, range, nugget), function(...) vgm(model = "Exp", ...))) |>
  mutate(vgfit = map2(variogs, model, ~ fit.variogram(.x, fit.method = 1, model = .y)))

# prepare theoretical variogram
gg_variog_theo <- gg_variogs |>
    mutate(vgline = map(vgfit, ~ variogramLine(.x, maxdist = extent * 0.5, n = 1000))) |>
    dplyr::select(factor_label, vgline) |>
    tidyr::unnest(col = vgline)

Visualizing variograms per factor.

# Names for the strips
# factor_labels <- c("F1" = "(A)~First~factor",
#         "F2" = "(B)~Second~factor",
#         "F3" = "(C)~Third~factor")
factor_labels <- c(
    "F1" = "atop('(F1) Severe food insecurity', 'including children')",
    "F2" = "atop('(F2) Anxiety and relying on', 'social relations')",
    "F3" = "atop('(F3) Adults eating less and', 'going hungry')"
)
gg_variog_emp$factor_label <- factor_labels[gg_variog_emp$factor_label]
gg_variog_theo$factor_label <- factor_labels[gg_variog_theo$factor_label]

ggthemr::ggthemr_reset()
pal_aux <- ggthemr:::load_palette.character("fresh")
pal_aux$swatch[c(2, 3)] <- pal_aux$swatch[c(3, 2)]
ggthemr::ggthemr(pal_aux, layout = "clear")

gg <- gg_variog_emp |>
  ggplot(aes(dist, gamma)) +
    geom_point(aes(size = np, alpha = np)) +
    geom_line(aes(dist, gamma, color = "Exponential model"), data = gg_variog_theo) +
    expand_limits(y = 0, x = 0) +
    scale_x_continuous(limits = c(0, extent * 0.5)) +
    scale_y_continuous(limits = c(0.5, 1.2)) +
    scale_size_continuous(range = c(1, 2.5), breaks = seq(250, 1250, 250)) +
    facet_wrap(~factor_label, ncol = 3, labeller = label_parsed) +
    labs(y = expression(widetilde(gamma)(u):Semi-variogram),
    x = expression(u:Distance~"in"~meters), size = "Number of pairs:", alpha = "Number of pairs:", color = "") +
     guides(color= guide_legend(), size=guide_legend()) +
    # labs(x = "distance", y = "semi-variogram") +
    theme(legend.position = "bottom") +
    theme_Publication()

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

gg

#> Warning: The `scale_name` argument of `discrete_scale()` is deprecated as of ggplot2 3.5.0.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

Save graph of variogram.

ggsave(file.path(path_data, "figures", "cifa_variogs.pdf"), gg, width = 7, height = 3.7)