stLMM formula terms

Formula-term constructors recognized by stLMM.

Details

The symbols documented here are package-native model terms used inside an stLMM() formula. They are not exported as ordinary R functions to call outside a model formula.

IID grouped random effects

iid(group) adds an iid Gaussian random intercept for the observed levels of group. It contributes one column per observed group level to the explicit random-effect design matrix Z.

x:iid(group) adds an iid Gaussian group-varying slope for the numeric covariate x. It contributes one column per observed group level, with row entries scaled by x.

iid(group) + x:iid(group) creates two separate scalar iid random-effect terms with separate variance parameters, typically named iid_1 and iid_2. This is not a correlated intercept/slope block.

IID random effects are sampled through the explicit Z alpha random-effect path. They do not require recover. Prediction from a fit object is available for existing grouping levels. New grouping levels in newdata currently error.

Structured latent process terms

ar1(index) adds a one-dimensional AR1 latent process over the sorted unique values of index. Repeated observations at the same index share a latent process node. The implemented AR1 support uses adjacency in the sorted unique fitted support.

gp(coord1, coord2, ..., cov_model = "exp") adds a dense Gaussian process over unique coordinate rows. Spatial covariance models include cov_model = "exp" with parameter phi and cov_model = "matern" with parameters phi and nu.

nngp(coord1, coord2, ..., m = 15, ordering = "coord", st_scale = 1, cov_model = "exp") adds a nearest-neighbor Gaussian process over unique coordinate rows. Supported ordering values are "coord", "default", "maxmin", "hilbert", "random", or a numeric permutation. The "maxmin" ordering can be very slow for large numbers of unique coordinate rows because it uses an exact iterative max-min construction. The "hilbert" ordering uses only the first two coordinate columns to rank nodes; covariance distances still use all coordinate columns according to the selected cov_model. For space-time covariance models, "hilbert" is currently rejected because the implemented Hilbert ordering is two-dimensional. Space-time covariance models include cov_model = "sep_exp", "multi_res_sep_exp", and "gneiting".

car(area_id, graph = g, car_model = "proper") adds a CAR latent process over areal nodes defined by g, where g is produced by car_graph. Rows with the same area_id share a latent process node. Areas in the graph that do not have observed responses remain in the latent support and can be recovered when they are represented by rows with missing response values. Prediction is implemented for existing graph areas. New CAR areas currently error. For polygon inputs, car_graph() can optionally add auditable nearest-neighbor bridge edges for disconnected island components. car_model = "proper" uses the package's original degree-scaled proper CAR precision. car_model = "leroux" uses the Leroux CAR precision, which bridges iid area effects and intrinsic-CAR smoothing.

dagar(area_id, graph = g, ordering = "coord") adds an ordered directed acyclic graph autoregressive (DAGAR) latent process over the areal nodes in g. The same car_graph object used by car() supplies the undirected adjacency graph; ordering orients each adjacency edge from the earlier area to the later area. Supported orderings are "coord", "default", "maxmin", "hilbert", "random", or a numeric permutation. Coordinate-based orderings require an sf-based graph; for adjacency-matrix graphs, "coord" and "default" use the supplied row order. Prediction is implemented for existing graph areas. New DAGAR areas currently error.

car_time(area_id, time, graph = g, time_model = "ar1", car_model = "proper") adds a separable CAR-by-time latent process over the full Cartesian product of graph areas and sorted unique fitted time values. The latent support is stacked location-major, so all time points for the first area come first. The current precision is sigma_sq^-1 * (Q_space %x% Q_time), with one variance term for the whole space-time process. The default temporal model is time_model = "ar1", which uses adjacency in the sorted fitted time support. time_model = "exp" uses Corr(t_i, t_j) = exp(-lambda * abs(t_i - t_j)); fitted time values must be numeric and finite. New prediction rows must use existing areas and existing fitted time values.

dagar_time(area_id, time, graph = g, ordering = "coord", time_model = "ar1") adds a separable DAGAR-by-time latent process over the full product of ordered DAGAR graph areas and sorted unique fitted time values. The undirected graph supplied by g is first oriented by ordering as in dagar(), and the resulting ordered DAGAR spatial precision is crossed with the temporal precision. The latent support is stacked ordered-area-major, so all time points for the first ordered area come first. The user-facing recovered samples are returned in the original graph-area order crossed with fitted time order. The default temporal model is time_model = "ar1". time_model = "exp" uses Corr(t_i, t_j) = exp(-lambda * abs(t_i - t_j)); fitted time values must be numeric and finite. New prediction rows must use existing graph areas and existing fitted time values.

For space-time covariance models, coordinate order matters and the final coordinate is treated as time. For example, nngp(lon, lat, time, cov_model = "sep_exp") is not equivalent to nngp(time, lon, lat, cov_model = "sep_exp").

For space-time NNGP terms, st_scale controls only the search geometry used to order nodes and choose nearest neighbors. The search distance is equivalent to Euclidean distance after replacing the final coordinate by st_scale * time. Thus st_scale = 1 uses the coordinates as supplied, while st_scale = s treats one unit of the final coordinate as s spatial-coordinate units for graph construction. The original coordinates are still stored on the fitted graph and used for covariance evaluation, so st_scale does not change covariance parameter interpretation or prior calibration. Prediction from recovered NNGP fits inherits the fitted graph's st_scale by default and can override it via predict.stLMM_recovery.

Process terms can be covariate-scaled with the same colon convention used for spatially varying coefficients, for example x:nngp(lon, lat) or x:ar1(day). x:car(area_id, graph = g) is the CAR spatially varying coefficient form. x:car_time(area_id, time, graph = g) is the CAR-time spatial-temporal varying coefficient form. x:dagar_time(area_id, time, graph = g) is the DAGAR-time spatial-temporal varying coefficient form.

Structured process terms are collapsed during parameter updates. Gaussian process fits use recover after fitting to draw latent process samples before using process contributions in fitted.stLMM or predict.stLMM_recovery. Polya-Gamma process fits save in-chain process draws by default as save_process = list(start = 1, thin = 1) and recover selects from those draws.

Covariance functions

Dense GP and NNGP terms use the process variance sigma_sq and the theta parameters named by get_cor_models.

For cov_model = "exp", with spatial distance \(h\), $$ C(h) = \sigma^2 \exp(-\phi h). $$

For cov_model = "matern", $$ C(h) = \sigma^2 \frac{(\phi h)^\nu}{2^{\nu - 1}\Gamma(\nu)} K_\nu(\phi h), $$ where \(K_\nu(\cdot)\) is the modified Bessel function of the second kind.

For cov_model = "sep_exp", with spatial distance \(h\) and temporal distance \(u\), $$ C(h, u) = \sigma^2 \exp(-\phi h)\exp(-\lambda u). $$

For cov_model = "multi_res_sep_exp", $$ C(h, u) = \sigma^2\left[ \alpha \exp(-\phi_1 h)\exp(-\lambda_1 u) + (1 - \alpha)\exp(-\phi_2 h)\exp(-\lambda_2 u) \right]. $$ This two-scale separable covariance follows the spatial-temporal forest inventory model used by May and Finley (2025).

For cov_model = "gneiting", $$ C(h, u) = \sigma^2 (1 + a |u|^{2\alpha})^{-(\delta + d/2)} \exp\left\{ -\frac{c h^{2\gamma}}{(1 + a |u|^{2\alpha})^{\beta\gamma}} \right\}. $$ Here \(h\) is Euclidean spatial distance, \(u\) is temporal distance, and \(d\) is the number of spatial coordinate columns before the final time coordinate. This is Gneiting (2002, equation 12). Parameter supports are \(a > 0\), \(c > 0\), \(0 < \alpha \le 1\), \(0 \le \beta \le 1\), \(0 < \gamma \le 1\), and \(\delta \ge 0\). Free parameters use the strict interior of these supports where needed by the Metropolis transformation; boundary values such as beta = fixed(0), alpha = fixed(1), gamma = fixed(1), or delta = fixed(0) can be fixed explicitly.

The covariance used by Datta et al. (2016, equation 6.1) is recovered in two spatial dimensions by setting alpha = fixed(1), gamma = fixed(0.5), delta = fixed(0), and using beta as Datta's interaction parameter. A conservative starting point is to estimate a, c, and optionally beta, while fixing alpha, gamma, and delta.

For ar1(index), the fitted support is the sorted unique values of index. The correlation between support positions is $$ \mathrm{Cor}(w_i, w_j) = \phi^{|i-j|}. $$ This is an ordered-support AR1, not a numeric time-gap model.

For car(area, graph = g), the proper CAR precision is $$ Q = \sigma^{-2}(D - \rho G), $$ where \(G\) is the binary adjacency matrix stored by car_graph() and \(D\) is the diagonal degree matrix. The graph is not stored as a row-standardized matrix; rather, this degree-scaled precision is equivalent to a CAR conditional specification whose mean uses the neighbor average, $$ E(w_i \mid w_{-i}) = \frac{\rho}{d_i}\sum_{j \sim i} w_j, $$ with conditional variance \(\sigma^2/d_i\). The current implementation restricts \(\rho\) to \((0,1)\).

For car(area, graph = g, car_model = "leroux"), the Leroux CAR precision is $$ Q = \sigma^{-2}\{(1-\rho)I + \rho(D - G)\}. $$ This model has the same graph input and parameter names as the proper CAR model. As \(\rho\) approaches zero, the prior approaches iid area effects with variance \(\sigma^2\); as \(\rho\) approaches one, it approaches an intrinsic-CAR smoothing structure. The implementation keeps \(\rho\) strictly inside \((0,1)\).

For dagar(area, graph = g), let \(\pi\) be the ordering of the graph areas and let $$ N_\pi(i) = \{j : i \sim j,\ \pi^{-1}(j) < \pi^{-1}(i)\} $$ be the directed parent set for area \(i\). Let \(n_\pi(i)=|N_\pi(i)|\). The ordered DAGAR conditional specification is $$ w_i \mid w_{<i,\pi} \sim N\left( \frac{\rho}{1 + \{n_\pi(i)-1\}\rho^2} \sum_{j \in N_\pi(i)} w_j,\, \frac{1 + \{n_\pi(i)-1\}\rho^2}{1-\rho^2} \right), $$ where the second normal argument is precision. Nodes with no directed parents have conditional mean zero. In stLMM scale convention, $$ Q = \sigma^{-2} L_\pi' F_\pi L_\pi, $$ where \(L_\pi\) is unit lower triangular and \(F_\pi\) contains the conditional precisions. The ordering is part of the model specification, not a CHOLMOD or display-only ordering. This is the ordered spatial DAGAR model of Datta, Banerjee, Hodges, and Gao (2019).

For car_time(area, time, graph = g), the separable precision is $$ Q = \sigma^{-2}\{Q_{\mathrm{space}}(\rho) \otimes Q_{\mathrm{time}}\}. $$ car_model selects whether Q_space is the proper CAR precision \(D - \rho G\) or the Leroux precision \((1-\rho)I + \rho(D-G)\). With time_model = "ar1", Q_time is the ordered-support AR1 precision with parameter phi. With time_model = "exp", the temporal correlation is $$ \mathrm{Cor}(w_t, w_{t'}) = \exp(-\lambda |t - t'|), $$ represented internally through its Markov precision on sorted numeric fitted time values.

For dagar_time(area, time, graph = g), the separable precision is $$ Q = \sigma^{-2}\{Q_{\mathrm{dagar}}(\rho) \otimes Q_{\mathrm{time}}\}. $$ Q_dagar is the scale-free ordered DAGAR precision defined above, after orienting graph edges by ordering. Q_time is the same AR1 or exponential-time Markov precision used by car_time(). The log determinant is evaluated by the Kronecker identity, $$ \log|Q| = -ST\log(\sigma^2) + T\log|Q_{\mathrm{dagar}}(\rho)| + S\log|Q_{\mathrm{time}}|, $$ where \(S\) is the number of graph areas and \(T\) is the number of fitted time values.

Controls and priors

Term-specific covariance controls use generated names such as iid_1, ar1_1, gp_1, nngp_1, car_1, dagar_1, car_time_1, and dagar_time_1. The names are assigned by term type in formula order.

IID variance priors are inverse-gamma prior objects, for example:

priors = list(iid_1 = list(sigma_sq = ig(shape, scale)))

Process terms require a process variance prior plus finite-support priors for their theta parameters, for example:

priors = list(
  nngp_1 = list(
    sigma_sq = ig(shape, scale),
    phi = uniform(lower, upper)
  )
)

CAR terms use a process variance prior and CAR rho bounds:

priors = list(
  car_1 = list(
    sigma_sq = ig(shape, scale),
    rho = uniform(lower, upper)
  )
)

DAGAR terms use the same process variance and rho prior structure:

priors = list(
  dagar_1 = list(
    sigma_sq = ig(shape, scale),
    rho = uniform(lower, upper)
  )
)

CAR-time terms use one process variance, spatial CAR rho, and a temporal parameter. For time_model = "ar1", the temporal parameter is AR1 phi:

priors = list(
  car_time_1 = list(
    sigma_sq = ig(shape, scale),
    rho = uniform(lower, upper),
    phi = uniform(lower, upper)
  )
)

For time_model = "exp", replace phi with positive temporal decay lambda.

DAGAR-time terms use the same process variance, ordered-DAGAR rho, and temporal parameter structure:

priors = list(
  dagar_time_1 = list(
    sigma_sq = ig(shape, scale),
    rho = uniform(lower, upper),
    phi = uniform(lower, upper)
  )
)

For time_model = "exp", replace phi with positive temporal decay lambda.

Use get_cor_models to inspect supported cov_model values and their correlation parameter names.

References

Gneiting, T. (2002). Nonseparable, stationary covariance functions for space-time data. Journal of the American Statistical Association, 97(458), 590–600. doi:10.1198/016214502760047113.

May, P. B. and Finley, A. O. (2025). Spatial-temporal prediction of forest attributes using latent Gaussian models and inventory data. Spatial Statistics, 69, 100917. doi:10.1016/j.spasta.2025.100917.

See also

stLMM, recover, fitted.stLMM, predict.stLMM_recovery, get_cor_models, stLMM_priors

Examples

adj <- matrix(
  c(
    0, 1, 0,
    1, 0, 1,
    0, 1, 0
  ),
  nrow = 3,
  byrow = TRUE,
  dimnames = list(c("a", "b", "c"), c("a", "b", "c"))
)
g <- car_graph(adj)

set.seed(1)
dat <- data.frame(
  y = rnorm(3),
  x = rnorm(3),
  area = c("a", "b", "c")
)

fit <- stLMM(
  y ~ x + car(area, graph = g),
  data = dat,
  priors = list(
    resid = list(tau_sq = ig(2, 1)),
    car_1 = list(sigma_sq = ig(2, 1), rho = uniform(0.05, 0.95))
  ),
  n_samples = 8,
  warmup = FALSE,
  verbose = FALSE
)
rec <- recover(fit)
fitted(rec)
#> [1] -0.17575214  0.05087614 -0.69700934