inst/doc/spatial-models.R

In tulpa: Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width = 6, fig.height = 4
)
set.seed(20260529)

## ----cran-gate, include = FALSE-----------------------------------------------
# The model fits below are evaluated when the article is built locally, in CI
# and for the pkgdown site (all of which set NOT_CRAN). CRAN's check farm gives
# the whole check a ten-minute budget, which these fits do not fit inside, so
# there the code is shown without being run.
EVAL_FITS <- identical(Sys.getenv("NOT_CRAN"), "true")
knitr::opts_chunk$set(eval = EVAL_FITS)

## ----eval = FALSE-------------------------------------------------------------
# # Current release, built for Windows and macOS
# install.packages("tulpa",
#                  repos = c("https://gcol33.r-universe.dev",
#                            "https://cloud.r-project.org"))
# 
# # Development version from GitHub
# install.packages("pak")
# pak::pak("gcol33/tulpa")
# 
# # Pin a release
# pak::pak("gcol33/tulpa@v0.6.0")

## ----load, message = FALSE----------------------------------------------------
# library(tulpa)

## ----adjacency----------------------------------------------------------------
# K <- 25
# W <- matrix(0, K, K)
# for (i in 1:K) {
#   j <- if (i < K) i + 1 else 1
#   W[i, j] <- W[j, i] <- 1
# }

## ----truth--------------------------------------------------------------------
# field_true <- as.numeric(scale(sin(2 * pi * (1:K) / K)))
# beta0 <- -0.2
# beta1 <-  0.8

## ----simulate-----------------------------------------------------------------
# n      <- 800
# x      <- rnorm(n)
# region <- sample(1:K, n, replace = TRUE)
# eta    <- beta0 + beta1 * x + field_true[region]
# y      <- rbinom(n, size = 1, prob = plogis(eta))
# ds     <- data.frame(y = y, x = x, region = factor(region))

## ----fit----------------------------------------------------------------------
# fit <- tulpa(y ~ x + spatial(region), data = ds, family = "binomial",
#              spatial = list(type = "icar", adjacency = W),
#              mode = "nested_laplace")
# summary(fit)$inference$backend

## ----coef---------------------------------------------------------------------
# coef(fit)

## ----summary------------------------------------------------------------------
# summary(fit)

## ----field-extract------------------------------------------------------------
# g         <- tulpa_outer_grid(fit)
# field_hat <- colSums(g$weights * g$modes)[g$n_fixed + seq_len(K)]
# round(head(field_hat), 3)

## ----field-plot, fig.alt = "Estimated spatial field per region against the simulated truth, near the identity line"----
# plot(field_true - mean(field_true),
#      field_hat  - mean(field_hat),
#      xlab = "true field (centred)", ylab = "estimated field (centred)",
#      pch = 19, col = "steelblue")
# abline(0, 1, lwd = 2, col = "grey40")

## ----field-cor----------------------------------------------------------------
# cor(field_true, field_hat)

## ----tau----------------------------------------------------------------------
# summary(fit)$hyperparameters[, c("parameter", "mean", "ci_lo", "ci_hi")]

## ----confint------------------------------------------------------------------
# confint(fit)

## ----intercept-se-------------------------------------------------------------
# coef(summary(fit))["(Intercept)", "std.error"]

## ----predict------------------------------------------------------------------
# nd <- data.frame(x = seq(-2.5, 2.5, length.out = 50))
# pr <- predict(fit, newdata = nd, type = "link", se.fit = TRUE)
# head(pr, 3)

## ----predict-plot, fig.alt = "Predicted logit response across x with a 95 percent credible band"----
# plot(nd$x, pr$fit, type = "n", xlab = "x", ylab = "predicted logit")
# polygon(c(nd$x, rev(nd$x)), c(pr$lower, rev(pr$upper)),
#         col = adjustcolor("steelblue", 0.25), border = NA)
# lines(nd$x, pr$fit, lwd = 2)

## ----compare------------------------------------------------------------------
# m0 <- tulpa(y ~ x, data = ds, family = "binomial", mode = "laplace")
# cmp <- compare_models(nonspatial = m0, spatial = fit, criterion = "loglik")
# cmp

## ----compare-diff-------------------------------------------------------------
# diff(cmp$logLik)

## ----continuous-spec----------------------------------------------------------
# gp_spec <- spatial_gp(~ lon + lat)
# gp_spec$type

## ----spde-fit-----------------------------------------------------------------
# set.seed(20260531)
# n      <- 300
# coords <- data.frame(lon = runif(n), lat = runif(n))
# field  <- 1.4 * (sin(2.5 * coords$lon) + cos(2.5 * coords$lat))
# field  <- field - mean(field)
# xcov   <- rnorm(n)
# y      <- rpois(n, exp(1.0 + 0.5 * xcov + field))
# 
# spde <- spatial_spde(~ lon + lat, data = coords, max_edge = c(0.2, 0.5))
# fit  <- fit_spde(y = y, X = model.matrix(~ xcov), spatial = spde,
#                  family = "poisson", range = 0.4, sigma = 0.9)
# round(coef(fit), 3)