How do I run a Nearest Neighbors Gaussian Process spatial model in tinyVAST and compare the output to SPDE models

[MorganJohnston4]

Troubleshooting report: pomega1_g returns constant value with NNGP domain despite non-trivial spatial parameters

Summary

When fitting a tinyVAST model with spatial_domain set to an NNGP domain (via make_nngp_domain()), the spatial random effect component (predict(..., what = "pomega1_g")) returns a constant value across all locations — both at the original fitted data locations and on new prediction grids. This happens despite the model reporting a non-trivial log_kappa estimate (implying a finite, identified spatial range) and despite mu_g varying appropriately across locations. The exact same model structure (single-variable SEM spatial term, no time dimension) works correctly when run using the package's own vignette example with a continuous SPDE mesh (fm_mesh_2d()) instead of an NNGP domain.

Environment

• tinyVAST installed from vast-lib/tinyVAST@dev (per the packageVersion("tinyVAST") < "1.6.0" guard in my script)
• Please let me know if you want exact sessionInfo() output

Steps to reproduce / debugging history

I'm including the full debugging path since each step ruled out a plausible alternative explanation.

1. Initial model

nngp_domain <- make_nngp_domain(sf_grid, n = 4)

dat_nngp$var <- "density_m2"

sem_nngp <- "
  density_m2 <-> density_m2, spatial_sd
"

fit_nngp <- tinyVAST(
  formula = density_m2 ~ 0 + Substrate2 + O2 + s(Depth),
  data = dat_nngp,
  spatial_domain = nngp_domain,
  family = tweedie("log"),
  space_columns = c("X","Y"),
  space_term = sem_nngp,
  control = tinyVASTcontrol(gmrf_parameterization = "projection")
)

Note: I initially also passed time_column = "Year", but my data has no Year/time column at all — removing this had no effect on the issue described below.
Also, I have no idea if specifying the space_term is needed for the NNGP model. It seemed like a potential spot that might be restricting the spatial term output.

2. Model converges with sensible diagnostics

fit_nngp$sdrep

              Estimate Std. Error
alpha_j    -10.2100205 0.93357297
alpha_j     -9.3866305 0.92725634
alpha_j    -15.3674782 0.92328090
alpha_j     -0.6970693 0.31208113
theta_z      6.4097543 0.55930298
log_lambda  -3.7071005 0.84432905
log_sigma    1.5732367 0.02639787
log_sigma    1.3221383 0.01273130
log_kappa    1.4458246 0.19988266
Maximum gradient component: 0.00740102

Max gradient is small, SEs look reasonable, and log_kappa has a tight, well-identified estimate. Note spatial_sd (the name I assigned in the SEM string) does not appear in the output at all — only the generic theta_z etc. labels — so it's unclear whether theta_z corresponds to spatial_sd or whether the custom name is simply not being used in sdreport()'s printed output.

3. pomega1_g is constant at original data locations

pred_pomega <- predict(fit_nngp, what = "pomega1_g")  # no newdata
range(pred_pomega)

[1] 1.687459 1.687459

This is constant across all 37,816 observations.

4. mu_g (full predicted mean) varies normally

pred_mu <- predict(fit_nngp, what = "mu_g")
range(pred_mu)

[1] 5.284952e-08 8.935022e-04

5. palpha1_g (fixed-effects-only linear predictor) also varies, and on its own does not span a wide enough range to explain mu_g's full variation alone

pred_palpha <- predict(fit_nngp, what = "palpha1_g")
range(pred_palpha)

[1] -16.064548  -9.386631

So mu_g varies, palpha1_g varies, but pomega1_g — the spatial-only component — is flat. This seems inconsistent: if there were truly no spatial signal, I'd expect pomega1_g to hover near zero, not sit at 1.69.

6. Raw TMB report shows no per-observation spatial random effect array

rep_nngp <- fit_nngp$obj$report()
names(rep_nngp)

[1] "deviance"    "mu_i"        "nll"         "log_tau"     "p2_i"        "H"           "negloglik_i"
[8] "Q_ss"        "p_i"         "devresid_i"

No omega-named or per-location field of length matching the data appears here.

7. Sanity check: the package's own vignette example (continuous SPDE mesh, not NNGP) works correctly

Using the exact code from the "Spatial modeling" vignette (spatial.Rmd), with fm_mesh_2d() instead of an NNGP domain:

sem = "
  density <-> density, spatial_sd
"
out = tinyVAST(data = Data, formula = n ~ s(w), spatial_domain = mesh,
               control = tinyVASTcontrol(getsd = FALSE), space_term = sem)

pred_toy <- predict(out, what = "pomega1_g")
range(pred_toy)

[1] -2.073545  1.822922

This varies correctly. So the bug (if it is one) appears specific to the NNGP domain code path, not the SEM string syntax or the pomega1_g extraction mechanism in general.

In addition, the sdmTMB model I ran as a comparison did report spatial vairaion, the fixed effects were similar to the fixed efects for NNGP, but the predicted density between the two models was off by a factor of 100x.

Other notes from debugging

• I confirmed my prediction grid and the NNGP domain's sf_areal are on the same coordinate scale (UTM km) before suspecting a units issue, and confirmed mesh/grid resolution between my sdmTMB SPDE mesh (~4.2 km median edge length) and the NNGP sf_areal grid (2 km cells) are broadly comparable, so I don't think this is a resolution mismatch.
• I believe I figured out how to at least measure cell size and mesh triangle edge length to try and make them match closely.
• Increasing NNGP neighbor count (n = 4, 8, 16, 32) did not change log_kappa meaningfully (1.486 to 1.493), so the constant pomega1_g doesn't appear to be a too-few-neighbors artifact either.
• Attached is my full debugging code

Question for you

Is pomega1_g expected to work for NNGP domains in the current dev branch, or is this a known gap/in-progress feature?

testing_spde_vs_nngp.docx

[James-Thorson]

I caught a similar bug for other non-SPDE domains recently, and fixed it in the tinyVAST@dev branch. However, I can't tell from my NEWS text whether it also applied to NNGP. Do you mind testing again using tinyVAST@dev branch rather than CRAN verison?

[MorganJohnston4]

I am pretty sure I was running the dev branch. But I updated it again and made sure I was on the version 1.6.0.9000. I am still running into errors when predicting:
pred_nngp <- predict(fit_nngp, what = "pomega1_g")
Error in add_predictions(object = object, newdata = newdata, remove_origdata = remove_origdata) :
object 'newdatadata' not found

Some how before I was able to get past that error, but it still would not give me more than one value for Omega

Thanks

[James-Thorson]

Thanks for reporting this bug! I confirm that the tinyVAST@dev branch had a typo, which I've now fixed. The predict function now appears to be working for NNGP, based on exploration here.

I'm closing for now, but could you please report back with any further concerns?

[MorganJohnston4]

That worked! The scale of NNGP and SPDE omega's is still off...but at least we are now getting a different value for each predicted tile. Thanks!