BYM disease mapping: Germany oral cavity cancer

Introduction

This vignette fits a Besag-York-Mollie (BYM) model to the Germany oral-cavity-cancer data (Besag, York, and Mollié 1991). Counts $y_i$ in 544 districts follow $y_i \sim \mathrm{Poisson}(E_i e^{\eta_i})$ with a structured intrinsic CAR (ICAR) component plus an unstructured iid component.

The example uses the low-level adlaplace pipeline and the random_mult kernel. The generalized log-determinant and rank of the ICAR precision are stored in the bundled germany data; only the precision multipliers enter the AD tape (Sørbye and Rue 2014). The intercept is omitted because the ICAR precision is rank-deficient.

Data

The germany object includes case and expected counts, the district adjacency matrix, the scaled ICAR precision Q_scaled, and the prec list for random_mult. District boundaries are optional: when terra is available, germany$map is loaded from inst/extdata/germany_map.gpkg in the setup chunk.

Latent vector $\gamma = (u_{\mathrm{struct}}, u_{\mathrm{iid}})$ has length $2n$: structured and iid effects on $n$ areas. Observations use poisson_obs with offset $\log E_i$.

n <- length(germany$Y)
n_theta <- 2L
n_gamma <- 2 * n

A_obs <- cbind(
  Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, dims = c(n, n)),
  Matrix::sparseMatrix(i = seq_len(n), j = seq_len(n), x = 1, dims = c(n, n))
)
X_obs <- Matrix::Matrix(0, n, 0)

config <- list(
  transform_theta = TRUE,
  gamma = rep(0, n_gamma),
  theta = log(c(1e-2, 0.1)),
  offset = log(germany$E),
  shards = adlaplace::ad_shards(A_obs, num_shards = 2L),
  verbose = FALSE
)

data_obs <- adlaplace:::ad_data(
  y = germany$Y,
  A = A_obs,
  X = X_obs,
  theta_map = n_theta,
  ad_kind = "observations",
  ad_fun = "poisson_obs"
)

model_struct <- adlaplace:::ad_data(
  gamma_map = Matrix::sparseMatrix(
    i = seq_len(n), j = seq_len(n), x = 1,
    dims = c(n_gamma, n)
  ),
  theta_map = c(1L, n_theta),
  ad_kind = "random",
  ad_fun = "random_mult",
  precision = germany$prec
)


model_iid <- adlaplace:::ad_data(
  gamma_map = Matrix::sparseMatrix(
    i = seq.int(n + 1L, length.out = n), j = seq_len(n), x = 1,
    dims = c(n_gamma, n)
  ),
  theta_map = c(2L, n_theta),
  ad_kind = "random",
  ad_fun = "random_diagonal",
  precision = rep(1, n)
)

Building AD shards

Each model piece is taped separately with ad_fun_ptr(). The random_mult shard for the ICAR term is the slowest step.

ptr_obs <- adlaplace::ad_fun_ptr(data = data_obs, config = config)
ptr_iid <- adlaplace::ad_fun_ptr(data = model_iid, config = config)
ptr_struct <- adlaplace::ad_fun_ptr(data = model_struct, config = config)

ad_fun_plain <- c(ptr_obs, ptr_struct, ptr_iid)
ad_pack <- adlaplace::ad_fun(ad_fun_plain, num_threads = num_threads)

Combining shards

Per-shard log densities sum to the joint log density at a fixed parameter vector.

x_full <- c(config$gamma, config$theta)
shards <- seq.int(from = 0L, length.out = adlaplace:::n_groups(ad_fun_plain))

by_shard <- vapply(
  shards,
  function(s) {
    adlaplace::joint_log_dens(ad_fun_plain, x_full, shards = s, negative = TRUE)
  },
  numeric(1)
)

res_fdf <- adlaplace::fun_obj_fdfh(
  ad_pack,
  config$theta, config$gamma,
  inner = FALSE, verbose = FALSE
)

c(
  sum(by_shard),
  adlaplace::joint_log_dens(ad_fun_plain, x_full, negative = TRUE),
  adlaplace::joint_log_dens(ad_pack, x_full, negative = TRUE),
  res_fdf$f
)

## [1] -1485.755 -1485.755 -1485.755 -1485.755

Joint log density derivatives

Gradients and Hessians on all shards should agree with a subset covering the full model.

log_all <- adlaplace::joint_log_dens(ad_fun_plain, x_full, negative = TRUE)
log_sub <- adlaplace::joint_log_dens(
  ad_fun_plain, x_full,
  shards = 0:2, negative = TRUE
)

grad_all <- adlaplace::grad(ad_fun_plain, x_full, negative = TRUE)
grad_sub <- adlaplace::grad(ad_fun_plain, x_full, shards = 0:2, negative = TRUE)

hess_all <- adlaplace::hessian(ad_fun_plain, x_full, negative = TRUE)
hess_sub <- adlaplace::hessian(ad_fun_plain, x_full, shards = 0:2, negative = TRUE)

rbind(
  log_diff = log_all - res_fdf$f,
  max_abs_grad_diff = max(abs(grad_all - res_fdf$grad)),
  max_abs_hess_diff = max(abs(as.matrix(hess_all) - as.matrix(res_fdf$hessian))),
  log_sub_diff = log_all - log_sub,
  max_abs_grad_sub_diff = max(abs(grad_all - grad_sub)),
  max_abs_hess_sub_diff = max(abs(as.matrix(hess_all) - as.matrix(hess_sub)))
)

##                                [,1]
## log_diff               2.273737e-13
## max_abs_grad_diff      0.000000e+00
## max_abs_hess_diff      0.000000e+00
## log_sub_diff          -7.527037e+02
## max_abs_grad_sub_diff  5.440000e+02
## max_abs_hess_sub_diff  1.000000e+02

Inner optimization

Given hyperparameters $\theta$, inner_opt() maximizes over $\gamma$.

x_outer <- config$theta

inner_res <- adlaplace::inner_opt(
  parameters = log(c(1e-1, 1e-1)),
  gamma = config$gamma,
  ad_fun = ad_pack,
  control = list(
    maxit = 1000L,
    report.level = 4,
    report.freq = 5,
    grad.tol = 1e-8,
    cg.tol = 1e-6
  ),
  deriv = TRUE,
  verbose = FALSE
)

## Beginning optimization
## 
## iter          f         nrm_gr                     status         radCG iter                  CG result
##    5  -353.988972     0.000002   Continuing - TR contract    2.000000     50      Exceeded max CG iters
##   10  -353.988972     0.000002   Continuing - TR contract    0.062500     50      Exceeded max CG iters
##   15  -353.988972     0.000002   Continuing - TR contract    0.001953     50      Exceeded max CG iters
##   20  -353.988972     0.000002   Continuing - TR contract    0.000061     50      Exceeded max CG iters
##   25  -353.988972     0.000002   Continuing - TR contract    0.000002     50      Exceeded max CG iters
##   30  -353.988972     0.000000     Continuing - TR expand    0.000000     16         Intersect TR bound
## 
## Iteration has terminated
##   30  -353.988972     0.000000                    Success

inner_res$opt$fval

## [1] -353.989

summary(inner_res$hessian$trace3)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  -0.37929  -0.10944  -0.06669   0.83530  -0.03906 507.12520

Joint density derivative checks

The same finite-difference idea applies to the full joint log density $-\log p(\gamma, \theta \mid y)$ and its AD gradient and Hessian.

par(mfrow = c(4, 2), mar = c(2, 2, 2, 0), mgp = c(1, 0.5, 0))

x_here <- inner_res$full_parameters
Dpar_dens <- 96 # length(x_here) - 1L
Ngrid <- 11L
shards <- seq.int(from = 0L, length.out = adlaplace:::n_groups(ad_fun_plain))

par_grid <- matrix(x_here, nrow = Ngrid, ncol = length(x_here), byrow = TRUE)
Sx <- x_here[Dpar_dens] + seq(-0.1, 0.1, length.out = Ngrid)
SxD <- Sx[-1] - diff(Sx) / 2
par_grid[, Dpar_dens] <- Sx
x_list <- split(par_grid, row(par_grid))

dens <- vapply(
  x_list,
  adlaplace::joint_log_dens,
  numeric(1),
  ad_fun = ad_fun_plain,
  shards = shards,
  negative = FALSE
)

rbind(
  dens[1:5],
  dens[seq(to = length(dens), length.out = 5)]
)

##             1        2        3        4        5
## [1,] 352.2540 352.8781 353.3638 353.7110 353.9194
## [2,] 353.9194 353.7105 353.3621 352.8739 352.2459

grad <- do.call(
  cbind,
  lapply(
    x_list,
    adlaplace::grad,
    ad_fun = ad_fun_plain,
    inner = FALSE,
    shards = shards,
    negative = FALSE
  )
)

plot(Sx, grad[Dpar_dens, ], type = "l", ylab = "AD gradient")
points(SxD, diff(dens) / diff(Sx), pch = 16)
legend("topright", legend = c("AD", "finite diff"), lty = c(1, NA), pch = c(NA, 1), bty = "n")

hes <- array(
  dim = c(length(x_here), length(x_here), Ngrid),
  dimnames = list(NULL, NULL, NULL)
)
for (i in seq_len(Ngrid)) {
  hes[, , i] <- as.matrix(adlaplace::hessian(
    ad_fun_plain,
    x_list[[i]],
    inner = FALSE,
    shards = shards,
    negative = FALSE
  ))
}

grad_slope <- apply(grad, 1, diff) / mean(diff(Sx))
s_par2 <- sort(unique(c(Dpar_dens + 0:5, Dpar_dens + length(germany$Y), seq(to = length(x_here), length.out = 5))))
s_par2 <- s_par2[seq(1, min(c(prod(par("mfrow")) - 1L, length(s_par2))))]
for (Dpar2 in s_par2) {
  plot(
    Sx, hes[Dpar_dens, Dpar2, ],
    type = "l",
    ylab = paste0("H[", Dpar_dens, ",", Dpar2, "]"),
    ylim = range(c(hes[Dpar_dens, Dpar2, ], grad_slope[, Dpar2]), na.rm = TRUE)
  )
  points(SxD, grad_slope[, Dpar2], pch = 16)
}

hes_sub <- hes[1:length(config$gamma), 1:length(config$gamma), ]
hes_eval <- apply(hes_sub, 3, eigen)

Profile likelihood

log_lik_laplace() runs the inner optimizer and returns the Laplace approximation to the marginal log likelihood.

res_noderiv <- adlaplace::log_lik_laplace(
  x = x_outer,
  ad_fun = ad_pack,
  config = config,
  deriv = FALSE
)

## Beginning optimization
## 
## iter           f         nrm_gr                     status         radCG iter                  CG result
##   1  -1544.148145    10.761614     Continuing - TR expand    2.000000     50      Exceeded max CG iters
##   2  -1544.281934     9.892267                 Continuing    2.000000     50      Exceeded max CG iters
##   3  -1544.300928     1.776112                 Continuing    2.000000     50      Exceeded max CG iters
##   4  -1544.304080     1.751663                 Continuing    2.000000     50      Exceeded max CG iters
##   5  -1544.304701     0.386926                 Continuing    2.000000     50      Exceeded max CG iters
##   6  -1544.304829     0.372286                 Continuing    2.000000     50      Exceeded max CG iters
##   7  -1544.304856     0.086241                 Continuing    2.000000     50      Exceeded max CG iters
##   8  -1544.304862     0.081523                 Continuing    2.000000     50      Exceeded max CG iters
##   9  -1544.304863     0.018936                 Continuing    2.000000     50      Exceeded max CG iters
##  10  -1544.304864     0.017943                 Continuing    2.000000     50      Exceeded max CG iters
## 
## iter           f         nrm_gr                     status         radCG iter                  CG result
##  11  -1544.304864     0.004168                 Continuing    2.000000     50      Exceeded max CG iters
##  12  -1544.304864     0.003951                 Continuing    2.000000     50      Exceeded max CG iters
##  13  -1544.304864     0.000916                 Continuing    2.000000     50      Exceeded max CG iters
##  14  -1544.304864     0.000868                 Continuing    2.000000     50      Exceeded max CG iters
##  15  -1544.304864     0.000201                 Continuing    2.000000     50      Exceeded max CG iters
##  16  -1544.304864     0.000190                 Continuing    2.000000     50      Exceeded max CG iters
##  17  -1544.304864     0.000190   Continuing - TR contract    1.000000     50      Exceeded max CG iters
##  18  -1544.304864     0.000190   Continuing - TR contract    0.500000     50      Exceeded max CG iters
##  19  -1544.304864     0.000190   Continuing - TR contract    0.250000     50      Exceeded max CG iters
##  20  -1544.304864     0.000190   Continuing - TR contract    0.125000     50      Exceeded max CG iters
## 
## iter           f         nrm_gr                     status         radCG iter                  CG result
##  21  -1544.304864     0.000190   Continuing - TR contract    0.062500     50      Exceeded max CG iters
##  22  -1544.304864     0.000190   Continuing - TR contract    0.031250     50      Exceeded max CG iters
##  23  -1544.304864     0.000190   Continuing - TR contract    0.015625     50      Exceeded max CG iters
##  24  -1544.304864     0.000190   Continuing - TR contract    0.007812     50      Exceeded max CG iters
##  25  -1544.304864     0.000190   Continuing - TR contract    0.003906     50      Exceeded max CG iters
##  26  -1544.304864     0.000190   Continuing - TR contract    0.001953     50      Exceeded max CG iters
##  27  -1544.304864     0.000190   Continuing - TR contract    0.000977     50      Exceeded max CG iters
##  28  -1544.304864     0.000190   Continuing - TR contract    0.000488     50      Exceeded max CG iters
##  29  -1544.304864     0.000190   Continuing - TR contract    0.000244     50      Exceeded max CG iters
##  30  -1544.304864     0.000190   Continuing - TR contract    0.000122     50      Exceeded max CG iters
## 
## iter           f         nrm_gr                     status         radCG iter                  CG result
##  31  -1544.304864     0.000190   Continuing - TR contract    0.000061     50      Exceeded max CG iters
##  32  -1544.304864     0.000190   Continuing - TR contract    0.000031     50      Exceeded max CG iters
##  33  -1544.304864     0.000190   Continuing - TR contract    0.000015     50      Exceeded max CG iters
##  34  -1544.304864     0.000190   Continuing - TR contract    0.000008     50      Exceeded max CG iters
##  35  -1544.304864     0.000190   Continuing - TR contract    0.000004     50      Exceeded max CG iters
##  36  -1544.304864     0.000190   Continuing - TR contract    0.000002     50      Exceeded max CG iters
##  37  -1544.304864     0.000190   Continuing - TR contract    0.000001     50      Exceeded max CG iters
##  38  -1544.304864     0.000190   Continuing - TR contract    0.000000     50      Exceeded max CG iters
##  39  -1544.304864     0.000190   Continuing - TR contract    0.000000     50         Intersect TR bound
##  40  -1544.304864     0.000190   Continuing - TR contract    0.000000     22         Intersect TR bound
## 
## iter           f         nrm_gr                     status         radCG iter                  CG result
##  41  -1544.304864     0.000190   Continuing - TR contract    0.000000     12         Intersect TR bound
##  42  -1544.304864     0.000190   Continuing - TR contract    0.000000      7         Intersect TR bound
##  43  -1544.304864     0.000190   Continuing - TR contract    0.000000      4         Intersect TR bound
##  44  -1544.304864     0.000190   Continuing - TR contract    0.000000      3         Intersect TR bound
## 
## Iteration has terminated
##  44  -1544.304864     0.000190Radius of trust region is less than stop.trust.radius

res_deriv <- adlaplace::log_lik_laplace(
  x = x_outer,
  ad_fun = ad_pack,
  config = config,
  deriv = TRUE
)

## Beginning optimization
## 
## iter           f         nrm_gr                     status         radCG iter                  CG result
##   1  -1544.148145    10.761614     Continuing - TR expand    2.000000     50      Exceeded max CG iters
##   2  -1544.281934     9.892267                 Continuing    2.000000     50      Exceeded max CG iters
##   3  -1544.300928     1.776112                 Continuing    2.000000     50      Exceeded max CG iters
##   4  -1544.304080     1.751663                 Continuing    2.000000     50      Exceeded max CG iters
##   5  -1544.304701     0.386926                 Continuing    2.000000     50      Exceeded max CG iters
##   6  -1544.304829     0.372286                 Continuing    2.000000     50      Exceeded max CG iters
##   7  -1544.304856     0.086241                 Continuing    2.000000     50      Exceeded max CG iters
##   8  -1544.304862     0.081523                 Continuing    2.000000     50      Exceeded max CG iters
##   9  -1544.304863     0.018936                 Continuing    2.000000     50      Exceeded max CG iters
##  10  -1544.304864     0.017943                 Continuing    2.000000     50      Exceeded max CG iters
## 
## iter           f         nrm_gr                     status         radCG iter                  CG result
##  11  -1544.304864     0.004168                 Continuing    2.000000     50      Exceeded max CG iters
##  12  -1544.304864     0.003951                 Continuing    2.000000     50      Exceeded max CG iters
##  13  -1544.304864     0.000916                 Continuing    2.000000     50      Exceeded max CG iters
##  14  -1544.304864     0.000868                 Continuing    2.000000     50      Exceeded max CG iters
##  15  -1544.304864     0.000201                 Continuing    2.000000     50      Exceeded max CG iters
##  16  -1544.304864     0.000190                 Continuing    2.000000     50      Exceeded max CG iters
##  17  -1544.304864     0.000190   Continuing - TR contract    1.000000     50      Exceeded max CG iters
##  18  -1544.304864     0.000190   Continuing - TR contract    0.500000     50      Exceeded max CG iters
##  19  -1544.304864     0.000190   Continuing - TR contract    0.250000     50      Exceeded max CG iters
##  20  -1544.304864     0.000190   Continuing - TR contract    0.125000     50      Exceeded max CG iters
## 
## iter           f         nrm_gr                     status         radCG iter                  CG result
##  21  -1544.304864     0.000190   Continuing - TR contract    0.062500     50      Exceeded max CG iters
##  22  -1544.304864     0.000190   Continuing - TR contract    0.031250     50      Exceeded max CG iters
##  23  -1544.304864     0.000190   Continuing - TR contract    0.015625     50      Exceeded max CG iters
##  24  -1544.304864     0.000190   Continuing - TR contract    0.007812     50      Exceeded max CG iters
##  25  -1544.304864     0.000190   Continuing - TR contract    0.003906     50      Exceeded max CG iters
##  26  -1544.304864     0.000190   Continuing - TR contract    0.001953     50      Exceeded max CG iters
##  27  -1544.304864     0.000190   Continuing - TR contract    0.000977     50      Exceeded max CG iters
##  28  -1544.304864     0.000190   Continuing - TR contract    0.000488     50      Exceeded max CG iters
##  29  -1544.304864     0.000190   Continuing - TR contract    0.000244     50      Exceeded max CG iters
##  30  -1544.304864     0.000190   Continuing - TR contract    0.000122     50      Exceeded max CG iters
## 
## iter           f         nrm_gr                     status         radCG iter                  CG result
##  31  -1544.304864     0.000190   Continuing - TR contract    0.000061     50      Exceeded max CG iters
##  32  -1544.304864     0.000190   Continuing - TR contract    0.000031     50      Exceeded max CG iters
##  33  -1544.304864     0.000190   Continuing - TR contract    0.000015     50      Exceeded max CG iters
##  34  -1544.304864     0.000190   Continuing - TR contract    0.000008     50      Exceeded max CG iters
##  35  -1544.304864     0.000190   Continuing - TR contract    0.000004     50      Exceeded max CG iters
##  36  -1544.304864     0.000190   Continuing - TR contract    0.000002     50      Exceeded max CG iters
##  37  -1544.304864     0.000190   Continuing - TR contract    0.000001     50      Exceeded max CG iters
##  38  -1544.304864     0.000190   Continuing - TR contract    0.000000     50      Exceeded max CG iters
##  39  -1544.304864     0.000190   Continuing - TR contract    0.000000     50         Intersect TR bound
##  40  -1544.304864     0.000190   Continuing - TR contract    0.000000     22         Intersect TR bound
## 
## iter           f         nrm_gr                     status         radCG iter                  CG result
##  41  -1544.304864     0.000190   Continuing - TR contract    0.000000     12         Intersect TR bound
##  42  -1544.304864     0.000190   Continuing - TR contract    0.000000      7         Intersect TR bound
##  43  -1544.304864     0.000190   Continuing - TR contract    0.000000      4         Intersect TR bound
##  44  -1544.304864     0.000190   Continuing - TR contract    0.000000      3         Intersect TR bound
## 
## Iteration has terminated
##  44  -1544.304864     0.000190Radius of trust region is less than stop.trust.radius

data.frame(
  log_lik = c(res_noderiv$log_lik, res_deriv$log_lik),
  neg_log_lik = c(res_noderiv$neg_log_lik, res_deriv$neg_log_lik),
  inner_fval = c(res_deriv$opt$fval, res_deriv$opt$fval),
  row.names = c("deriv = FALSE", "deriv = TRUE")
)

##                 log_lik neg_log_lik inner_fval
## deriv = FALSE -1446.652    1446.652  -1544.305
## deriv = TRUE  -1446.652    1446.652  -1544.305

res_deriv$grad

## [1]  -3.231589 -27.001930

Outer optimization

outer_fn() and outer_gr() wrap log_lik_laplace() for use with stats::optim().

cache <- new.env(parent = emptyenv())
cache$gamma <- rep(0, n_gamma)

control_inner <- list(
  maxit = 100L,
  report.level = 0,
  report.freq = 0
)

opt <- optim(
  par = config$theta,
  fn = adlaplace::outer_fn,
  gr = adlaplace::outer_gr,
  method = "L-BFGS-B",
  lower = rep(-15, n_theta),
  control = list(maxit = 200L, factr = 1e3, trace = 1, REPORT = 1),
  config = config,
  ad_fun = ad_pack,
  cache = cache,
  control_inner = control_inner
)

## iter    1 value 1438.546746
## iter    2 value 1435.053755
## iter    3 value 1434.646017
## iter    4 value 1434.383745
## iter    5 value 1423.213290
## iter    6 value 1423.182386
## iter    7 value 1416.787819
## iter    8 value 1403.801739
## iter    9 value 1400.771909
## iter   10 value 1398.424226
## iter   11 value 1397.683740
## iter   12 value 1397.123150
## iter   13 value 1396.812895
## iter   14 value 1396.642743
## iter   15 value 1396.618185
## iter   16 value 1396.507839
## iter   17 value 1396.480000
## iter   18 value 1396.463077
## iter   19 value 1396.456788
## iter   20 value 1396.447679
## iter   21 value 1396.442234
## iter   22 value 1396.439843
## iter   23 value 1396.438560
## iter   24 value 1396.437934
## iter   25 value 1396.437930
## iter   26 value 1396.437616
## iter   27 value 1396.437458
## iter   28 value 1396.437380
## iter   29 value 1396.437341
## iter   30 value 1396.437321
## iter   31 value 1396.437320
## iter   32 value 1396.437310
## iter   33 value 1396.437306
## iter   34 value 1396.437304
## iter   35 value 1396.437303
## iter   36 value 1396.437302
## iter   37 value 1396.437302
## iter   38 value 1396.437302
## iter   39 value 1396.437302
## iter   40 value 1396.437302
## iter   41 value 1396.437302
## iter   42 value 1396.437302
## iter   43 value 1396.437302
## iter   44 value 1396.437302
## iter   45 value 1396.437302
## iter   46 value 1396.437302
## iter   47 value 1396.437302
## final  value 1396.437302 
## converged

res <- adlaplace::log_lik_laplace(
  x = opt$par,
  gamma = cache$gamma,
  ad_fun = ad_pack,
  config = config,
  control = control_inner,
  deriv = FALSE
)

u_hat <- res$opt$solution
eta_hat <- u_hat[seq_len(n)] + u_hat[seq.int(n + 1L, length.out = n)]
sd_hat <- exp(opt$par)
data.frame(
  component = c("SD structured (ICAR)", "SD unstructured (iid)"),
  estimate = sd_hat
)

##               component     estimate
## 1  SD structured (ICAR) 2.003913e-01
## 2 SD unstructured (iid) 2.725064e-06

Spatial effects

if (!is.null(germany$map)) {
  germany$map$eta_hat <- eta_hat
  terra::plot(
    germany$map, "eta_hat",
    main = expression(hat(eta) ~ "(adlaplace mode)")
  )
} else {
  par(mfrow = c(1, 1), mar = c(4, 4, 2, 1), mgp = c(2.2, 0.7, 0))
  plot(
    seq_len(n), eta_hat,
    type = "h", xlab = "district", ylab = expression(hat(eta)),
    main = "Total spatial effect (adlaplace mode; install terra for map)"
  )
}

## Warning: [is.lonlat] coordinates are out of range for lon/lat

The structured (ICAR) component typically dominates; the unstructured component is often near zero. Results are Laplace-approximate MLEs and conditional modes, not fully Bayesian posterior summaries.