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

can_run <- requireNamespace("mirt", quietly = TRUE) &&
  requireNamespace("gstat", quietly = TRUE)
knitr::opts_chunk$set(eval = can_run)

## ----eval = !can_run, echo = FALSE, results = "asis"----------------------------------------------
# cat("This vignette needs the mirt and gstat packages, which are not installed.")

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

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

## -------------------------------------------------------------------------------------------------
items <- unclass(ipixuna$items)
colnames(items) <- paste0("item", seq_len(nitems))
head(items)

## -------------------------------------------------------------------------------------------------
set.seed(11)
efa <- mirt(items, model = nfactors, itemtype = "2PL", verbose = FALSE)
loadings <- summary(efa, rotate = "varimax", verbose = FALSE)$rotF
loadings

## -------------------------------------------------------------------------------------------------
A <- (abs(loadings) > 0.24) * 1
A

## -------------------------------------------------------------------------------------------------
fscore <- fscores(efa, rotate = "varimax")
colnames(fscore) <- paste0("F", seq_len(nfactors))
head(fscore)

## -------------------------------------------------------------------------------------------------
fscore_sf <- st_sf(
  id = ipixuna$id, wealth = ipixuna$wealth, as.data.frame(fscore),
  geometry = st_transform(st_geometry(ipixuna), 32719)
)
fscore_sf

## ----fig.height = 3-------------------------------------------------------------------------------
fscore_sf |>
  tidyr::pivot_longer(F1:F3) |>
  ggplot(aes(wealth, value)) +
    geom_point() +
    geom_smooth(method = 'lm', formula = y ~ x) +
    facet_wrap(~ name) +
    theme_bw()

## ----fig.height = 5-------------------------------------------------------------------------------
vg1 <- variogram(F1 ~ wealth, fscore_sf, width = 20, cutoff = 1000)
vg1_fit <- fit.variogram(vg1, model = vgm(psill = 0.6, "Exp", range = 100, nugget = 0.1))
vg1_fit
plot(vg1, vg1_fit, pch = 19, main = "F1")

## ----fig.height = 5-------------------------------------------------------------------------------
vg2 <- variogram(F2 ~ wealth, fscore_sf, width = 50, cutoff = 1000)
vg2_fit <- fit.variogram(vg2, model = vgm(psill = 0.6, "Exp", range = 100, nugget = 0.1))
vg2_fit
plot(vg2, vg2_fit, pch = 19, main = "F2")

## ----fig.height = 5-------------------------------------------------------------------------------
vg3 <- variogram(F3 ~ wealth, fscore_sf, width = 50, cutoff = 800)
vg3_fit <- fit.variogram(vg3, model = vgm(psill = 0.6, "Exp", range = 100, nugget = 0.1))
vg3_fit
plot(vg3, vg3_fit, pch = 19, main = "F3")

## -------------------------------------------------------------------------------------------------
set.seed(1)
samples <- spifa(
  items ~ wealth, data = ipixuna, nfactors = nfactors,
  burnin = 1000, niter = 4000,
  constraints = list(discrimination = A, loading = diag(nfactors)),
  priors = list(range = list(initial = 150, mean = 150, sd = 0.4)))
samples

## ----fig.height = 6-------------------------------------------------------------------------------
plot(samples, c("c", "A"))
plot(samples, select = c("B", "T", "Corr"))
plot(samples, select = c("phi"))

## -------------------------------------------------------------------------------------------------
pred_samples <- predict(samples)
plot_predict(pred_samples) +
  scale_colour_distiller(palette = "RdYlBu")

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