## ----setup, include = FALSE-----------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, fig.align = "center", out.width = "100%")
knitr::opts_chunk$set(fig.width = 7, fig.height = 4)
knitr::opts_chunk$set(comment = "#>")
old_options <- options(width = 100)

## -------------------------------------------------------------------------------------------------
library(spifa)
library(ggplot2)
library(sf)

## -------------------------------------------------------------------------------------------------
data(ipixuna)
ipixuna

## -------------------------------------------------------------------------------------------------
ipixuna$items |> head()

## -------------------------------------------------------------------------------------------------
nitems <- ncol(ipixuna$items)
nfactors <- 3

## -------------------------------------------------------------------------------------------------
A <- matrix(1, nitems, nfactors)
A[c(4, 8), 1] <- 0
A[c(4, 5, 6, 7, 8, 10), 2] <- 0
A[c(5, 6), 3] <- 0
A

## -------------------------------------------------------------------------------------------------
A_mean <- matrix(0, nitems, nfactors)
A_mean[c(5, 6), 1] <- 1
A_mean[c(4, 8), 3] <- -1
A_mean

A_sd <- matrix(1, nitems, nfactors)
A_sd[A_mean != 0] <- 0.45
A_sd

## -------------------------------------------------------------------------------------------------
phi_mean <- 150
phi_sd <- 0.4

## ----eval = FALSE---------------------------------------------------------------------------------
# samples <- spifa(
#   items ~ wealth, data = ipixuna, nfactors = nfactors,
#   burnin = 10000, niter = 8000, thin = 2,
#   constraints = list(discrimination = A),
#   priors = list(
#     discrimination = list(initial = A_mean, mean = A_mean, sd = A_sd),
#     range = list(initial = phi_mean, mean = phi_mean, sd = phi_sd)
#   )
# )
# samples

## ----echo = FALSE---------------------------------------------------------------------------------
samples <- readRDS(system.file("vigdata", "spifa.rds", package = "spifa"))
samples

## -------------------------------------------------------------------------------------------------
summary(samples, select = "c")

## ----fig.height = 6-------------------------------------------------------------------------------
plot(samples)

## -------------------------------------------------------------------------------------------------
plot(samples, select = "c", nshow = 5)

## -------------------------------------------------------------------------------------------------
plot_trace(samples, select = "Corr")

## -------------------------------------------------------------------------------------------------
plot_trace(samples, select = "phi", facet = FALSE)

## ----fig.height = 5-------------------------------------------------------------------------------
plot_density(samples, select = "c")

## ----fig.height = 5-------------------------------------------------------------------------------
plot_density(samples, select = "B", facet = TRUE, facet_scales = "free_y")

## ----fig.height = 5-------------------------------------------------------------------------------
plot_interval(samples, select = "A")

## ----fig.height = 5-------------------------------------------------------------------------------
plot_interval(samples, select = "phi", horizontal = TRUE)

## ----eval = FALSE---------------------------------------------------------------------------------
# pred_samples <- predict(samples, burnin = 1000, thin = 5)
# posterior::as_draws_df(pred_samples)

## ----echo = FALSE---------------------------------------------------------------------------------
pred_samples <- readRDS(system.file("vigdata", "spifa-pred-observed.rds", package = "spifa"))
attr(pred_samples, "newdata") <- st_geometry(ipixuna)
posterior::as_draws_df(pred_samples)

## -------------------------------------------------------------------------------------------------
plot_predict(pred_samples) +
  scale_colour_distiller(palette = "RdBu")

## -------------------------------------------------------------------------------------------------
bnd <- st_geometry(ipixuna) |>
  st_union() |>
  st_convex_hull() |>
  st_buffer(50)

newdata <- st_sf(wealth = 0, geometry = st_make_grid(bnd, n = c(40, 40))) |>
  st_filter(bnd)
newdata

## ----eval = FALSE---------------------------------------------------------------------------------
# pred_samples <- predict(samples, newdata = newdata, thin = 10)

## ----echo = FALSE---------------------------------------------------------------------------------
pred_grid <- readRDS(system.file("vigdata", "spifa-prediction.rds", package = "spifa"))
newdata <- pred_grid$grid
bnd <- pred_grid$bnd

plot_predict_vig <- function (grid, prefix, boundary = NULL) {
  grid |>
    tidyr::pivot_longer(starts_with(prefix), names_to = "factor", values_to = "value") |>
    dplyr::mutate(factor = sub(paste0(prefix, "(\\d)"), "theta[\\1]", factor)) |>
    ggplot() +
    geom_sf(aes(fill = value), colour = NA) +
    geom_sf(data = boundary, fill = NA, linewidth = 0.6, color = "gray50") +
    facet_wrap(~ factor, labeller = label_parsed) +
    theme_bw(base_size = 9) +
    theme(
      axis.title = element_blank(),
      panel.grid = element_line(linetype = 2, colour = "grey80", linewidth = 0.25),
      panel.background = element_blank(),
      panel.ontop = TRUE,
      legend.position = "bottom",
      legend.title = element_blank(),
      legend.key.height = unit(0.4, "cm"),
      strip.background = element_blank(),
      strip.text = element_text(size = rel(1))
    )
}

## ----eval = FALSE---------------------------------------------------------------------------------
# plot_predict(pred_samples, boundary = bnd) +
#   scale_fill_distiller(palette = "RdBu") +
#   labs(title = "Predictive mean")

## ----echo = FALSE, fig.height = 4-----------------------------------------------------------------
plot_predict_vig(newdata, "theta", bnd) +
  scale_fill_distiller(palette = "RdBu") +
  labs(title = "Predictive mean")

## ----eval = FALSE---------------------------------------------------------------------------------
# plot_predict(pred_samples, boundary = bnd, stat = sd) +
#   scale_fill_viridis_c(option = "magma") +
#   labs(title = "Predictive standard deviation")

## ----echo = FALSE, fig.height = 4-----------------------------------------------------------------
plot_predict_vig(newdata, "sd", bnd) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "Predictive standard deviation")

## ----eval = FALSE---------------------------------------------------------------------------------
# plot_predict(pred_samples, boundary = bnd, stat = function (v) mean(v > 1)) +
#   scale_fill_viridis_c(limits = c(0, 1), direction = -1) +
#   labs(title = expression(paste("Exceedance probability: ", P(theta[j] > 1))))

## ----echo = FALSE, fig.height = 4-----------------------------------------------------------------
plot_predict_vig(newdata, "exceedance", bnd) +
  scale_fill_viridis_c(limits = c(0, 1), direction = -1) +
  labs(title = expression(paste("Exceedance probability: ", P(theta[j] > 1))))

## ----include = FALSE----------------------------------------------------------
options(old_options)