---
title: "Bayesian Spatial Item Factor Analysis"
output: rmarkdown::html_vignette
bibliography: references.bib
vignette: >
  %\VignetteIndexEntry{Bayesian Spatial Item Factor Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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)
```

In this vignette, we show how to use the **spifa** package to fit a spatial item
factor model, summarise and visualize the posterior samples, and predict and map
the latent factors at new locations.

## Introduction

We will focus on the spatial item factor model proposed in @chaconmontalvan2025. Let
$Y_{ij}$ be the binary response for item $j$ in individual $i$. The model can be defined
using an auxiliary variable $Z_{ij}$ such that:

$$
\begin{aligned}
  {Y}_{ij}  & =
  \begin{cases}
    1, & \text{if} ~ {Z}_{ij} > 0\\
    0, & \text{otherwise}
  \end{cases}\\
  {Z}_{ij} & = c_j + \boldsymbol{a}_j^\intercal\boldsymbol{\theta}_i + \epsilon_{ij},
  ~~ \epsilon_{ij} \sim {N}(0, 1),
\end{aligned}
$$

where the easiness parameters $c_j$ define how common is to endorse the item $j$, and the
discrimination parameters $\boldsymbol{a}_j$ define how important is item $j$ to
discriminate the latent factors/abilities $\boldsymbol{\theta}_i$.

The vector of latent abilities $\boldsymbol{\theta}_i$ is modelled in terms of
predictors $\boldsymbol{x}_i$, a vector of Gaussian processes $\boldsymbol{w}(s_i)$, and a
multivariate non-spatial term $\boldsymbol{v}_i$:
$$
\boldsymbol{\theta}_i = \boldsymbol{B} \boldsymbol{x}_i + \boldsymbol{T} \boldsymbol{w}(s_i) + \boldsymbol{v}_i,
$$
where $\boldsymbol{B}$ is the matrix of multivariate effects of the predictors,
$\boldsymbol{T}$ defines the relationship between the spatial processes and the latent
abilities, and $\boldsymbol{v}_i$ is allowed to have a correlations structure
$\boldsymbol{R}$.

This model is unidentifiable, so it is important to restrict the discrimination
parameters $\boldsymbol{a}_j$ and set informative priors for the spatial range of the
Gaussian processes $\boldsymbol{w}(s_i)$ before fitting it, to facilitate convergence
of the posterior sampling algorithm -- in practice, both are best informed by
exploratory analysis (e.g. a preliminary non-spatial item factor analysis for the
former, an empirical variogram for the latter), though we use them directly below
without walking through that derivation. The linear transformation $\boldsymbol{T}$,
by contrast, is not determined from exploratory analysis -- it is a structural choice
made by the user (usually diagonal, i.e. one independent Gaussian process per
factor), though alternative structures can be specified and compared.


## Ipixuna data

```{r}
library(spifa)
library(ggplot2)
library(sf)
```

The `ipixuna` data under analysis is included in the **spifa** package: a simulated
dataset of 10 `items` responses on 100 household locations, on the same area of study
as in @chaconmontalvan2025. It is a `sf` object containing the `items` responses as a
matrix, locations under the `geometry` column, and a predictor `wealth`.

```{r}
data(ipixuna)
ipixuna
```

The `items` can be easily obtained as an element of the data `ipixuna`:

```{r}
ipixuna$items |> head()
```

## Sample from the posterior

As mentioned before, an identifiable spifa model requires restricting the
discrimination parameters, along with adequate initial values and priors for the
spatial parameters (mainly the range). In practice, these are best obtained from a
preliminary analysis, as shown in `vignette("vg04-integrated-workflow")`. Here, we
instead assume the discrimination structure and some initial knowledge of the spatial
parameters are already available, and focus on defining, fitting, and diagnosing the
model.

First, we obtain the number of items (10) and define the number of factors (3) we will
use:

```{r}
nitems <- ncol(ipixuna$items)
nfactors <- 3
```

### Restrictions

We define the discrimination matrix (10x3): items 4 and 8 do not inform factor 1,
items 4, 5, 6, 7, 8, and 10 do not inform factor 2, and items 5 and 6 do not inform
factor 3.

```{r}
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
```

### Priors

In factor analysis, the sign of the factors is not identifiable, so we need to fix the
direction of each factor (which determines how we interpret it). A common way to do
this is to set a strong positive or negative prior on a few discrimination parameters.
The discrimination parameters are assumed to have a Normal prior, so `mean`/`sd` are
directly its mean and standard deviation. Here, we set the mean prior to 1 for items 5
and 6 on factor 1, and to -1 for items 4 and 8 on factor 3, each with a standard
deviation prior of 0.45:

```{r}
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
```

Similarly, we will define the prior for the range parameter of the Gaussian processes. By
default the `spifa()` function introduces 1 GP per factor, meaning that 3 range parameters
will be used in the model. Since the range must be positive, it is instead assumed to
have a Log-Normal prior: `mean`/`sd` are the mean and standard deviation of its
logarithm, not of the range itself. The range is in the same units as the coordinates
(metres here); `150` is a plausible order of magnitude for the spacing between
households in this study area. If we want to set the same prior for all three, we can
define a single mean and standard deviation that will be recycled across them:

```{r}
phi_mean <- 150
phi_sd <- 0.4
```

### Sampling

The function `spifa()` is used to sample from the posterior distribution of the model.
The main arguments are:

- `formula`: expression to define the response items and predictors to include in the
  latent abilities,
- `data`: a `sf` object containing the terms defined in `formula`,
- `nfactors`: the number of factors for the model,
- `niter`, `thin`, `burnin`: usual sampling arguments,
- `constraints`: a named list of constraints for certain parameters,
- `priors`: a named list of prior information including initial values.

In the following, we fit a spifa model with `3` latent factors, explained by `wealth`,
constraints for the *discrimination* parameters and priors for the *discrimination* and
*range* parameters. By default, it will introduce a GP for each latent factor.

```{r, 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
```

```{r, echo = FALSE}
samples <- readRDS(system.file("vigdata", "spifa.rds", package = "spifa"))
samples
```

The `print()` method shows metadata of the model such as *model type*, *formula*,
*dimensions*, MCMC dimensions including the number of saved *samples*. Additionally, it
shows the summary and diagnostics of the item model parameters (easiness and
discrimination), and the factor model parameters (effects, loading, range, and
correlation). Among the diagnostics, `ess_bulk` is the bulk effective sample size and
`rhat` is the potential scale reduction factor on split chains, which should be close
to 1 at convergence.

## Summarise and visualize samples

### Summary

While `print()` gives a quick overview, `summary()` computes the full set of posterior
statistics (mean, median, sd, mad, quantiles, `ess_bulk`, `ess_tail`, and `rhat`) for a
chosen parameter group, returned as a tidy tibble. Its main arguments are `select` (the
parameter group, e.g. `"c"`, `"A"`, `"phi"`), and `burnin`/`thin`, in case we want to
discard some initial samples or thin them further before summarising:

```{r}
summary(samples, select = "c")
```

### Plot

`plot()` gives a quick overview of a fitted model: traceplots and density plots, side
by side, for the `easiness` and `discrimination` parameters by default.

```{r, fig.height = 6}
plot(samples)
```

Use the `select` argument for a specific parameter group and `nshow` to define the
number of random parameters to show:

```{r}
plot(samples, select = "c", nshow = 5)
```

### Plotting functions

For more flexibility, **spifa** provides three individual plotting functions:

- `plot_trace()`,
- `plot_density()`,
- `plot_interval()`

They take as input the spifa `samples` with the `select` argument to define the parameter
group.

#### Trace plots

`plot_trace()` draws the values against iteration per parameter. By default, it plots
one panel per parameter (`facet = TRUE`), or all of them in a single panel when
`facet = FALSE`.

Traceplots of the residual correlation:

```{r}
plot_trace(samples, select = "Corr")
```

Traceplots of the range parameter of the GPs:

```{r}
plot_trace(samples, select = "phi", facet = FALSE)
```

#### Density plots

`plot_density()` draws the posterior density of each parameter in the selected group.
Similar to `plot_trace()`, you can control the `facet` argument.

Comparing the easiness parameters:

```{r, fig.height = 5}
plot_density(samples, select = "c")
```

Comparing the predictor effects:

```{r, fig.height = 5}
plot_density(samples, select = "B", facet = TRUE, facet_scales = "free_y")
```

#### Interval plots

`plot_interval()` draws the posterior mean/median and credible interval of
each parameter in the selected group.

Compare the discrimination parameters:

```{r, fig.height = 5}
plot_interval(samples, select = "A")
```

Compare the range parameters:

```{r, fig.height = 5}
plot_interval(samples, select = "phi", horizontal = TRUE)
```

## Predict the latent factors

**spifa** provides Bayesian prediction of the latent factors at (i) the observed locations
and (ii) new locations through the `predict()` function.

### Observed locations

By default `predict()` will simply return the samples of the latent factors at observed
locations as a `posterior::draws_array()` object:

```{r, eval = FALSE}
pred_samples <- predict(samples, burnin = 1000, thin = 5)
posterior::as_draws_df(pred_samples)
```

```{r, 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)
```

We can visualize a summary (`mean` by default) of these predictive samples using
`plot_predict()` and integrate it with `ggplot2` functionality:

```{r}
plot_predict(pred_samples) +
  scale_colour_distiller(palette = "RdBu")
```

### New locations

If prediction is desired for new location/profiles, then the argument `newdata` can be
provided, which is expected to be a `sf`/`sfc` object containing the new
location/profiles.

First, we create a grid where prediction will be performed. Then, we create a `newdata`
`sf` object using a predictor profile representing an individual with average
(standardized) wealth:

```{r}
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
```

We perform prediction for this `newdata` and thin the posterior samples to reduce
computational cost:

```{r, eval = FALSE}
pred_samples <- predict(samples, newdata = newdata, thin = 10)
```

```{r, 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))
    )
}
```

Now, we use the `plot_predict()` function to visualize the predictive `mean` by default.

```{r, eval = FALSE}
plot_predict(pred_samples, boundary = bnd) +
  scale_fill_distiller(palette = "RdBu") +
  labs(title = "Predictive mean")
```

```{r, echo = FALSE, fig.height = 4}
plot_predict_vig(newdata, "theta", bnd) +
  scale_fill_distiller(palette = "RdBu") +
  labs(title = "Predictive mean")
```

Other summaries can be mapped by simply providing a function in the `stat` argument:

```{r, eval = FALSE}
plot_predict(pred_samples, boundary = bnd, stat = sd) +
  scale_fill_viridis_c(option = "magma") +
  labs(title = "Predictive standard deviation")
```

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

Custom functions can easily be provided, for example, we can plot the probability of
exceeding the value of 1 to identify hotspots:

```{r, 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))))
```

```{r, 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))))
```

## References

```{r, include = FALSE}
options(old_options)
```