Package {SDALGCP2}

SDALGCP2: Fast Spatially Discrete Approximation to Log-Gaussian Cox Processes

Description

A faster, modernised successor to the SDALGCP package for fitting a spatially discrete approximation to a log-Gaussian Cox process (SDA-LGCP) to spatially aggregated disease counts. Performance-critical steps are implemented in C++ via RcppArmadillo. See DESIGN.md in the package source for the full design and statistical roadmap.

Author(s)

Maintainer: Olatunji Johnson olatunjijohnson21111@gmail.com

Authors:

• Emanuele Giorgi

• Peter Diggle

See Also

Useful links:

• https://github.com/olatunjijohnson/SDALGCP2

• https://olatunjijohnson.github.io/SDALGCP2/

• Report bugs at https://github.com/olatunjijohnson/SDALGCP2/issues

Fit a spatial SDA-LGCP model

Description

End-to-end user entry point: generates candidate points inside each region, assembles the aggregated region-level correlation array (C++), and estimates parameters by Monte Carlo maximum likelihood.

Usage

SDALGCP2(
  formula,
  data,
  my_shp,
  delta,
  phi = NULL,
  method = 1L,
  weighted = FALSE,
  pop_shp = NULL,
  kappa = 0.5,
  par0 = NULL,
  control.mcmc = NULL,
  phi_method = c("grid", "direct"),
  nugget = FALSE,
  confounding = c("none", "restricted"),
  reanchor = 0L,
  rho = 0.55,
  giveup = 1000L,
  nthreads = 0L,
  messages = FALSE
)

Arguments

Value

an object of class "SDALGCP2".

See Also

mcml_fit, precompute_corr, sda_points

Examples

library(sf)
## ---- simulate a lattice of regions and aggregated counts ----
set.seed(1)
bound <- st_as_sfc(st_bbox(c(xmin = 0, ymin = 0, xmax = 20, ymax = 20)))
shp   <- st_sf(geometry = st_make_grid(bound, n = c(8, 8)))
N     <- nrow(shp)

pts   <- sda_points(shp, delta = 1.2, method = 3)      # regular grid points
phi_grid <- seq(1, 5, length.out = 8)
corr  <- precompute_corr(pts, phi_grid)
Sig   <- 0.5 * corr$R[, , which.min(abs(phi_grid - 2.5))]
x1    <- as.numeric(scale(st_coordinates(st_centroid(shp))[, 1]))
pop   <- round(runif(N, 500, 3000))
y     <- rpois(N, pop * exp(cbind(1, x1) %*% c(-6, 0.5) +
                            as.numeric(t(chol(Sig)) %*% rnorm(N))))
dat   <- data.frame(y = y, x1 = x1, pop = pop)

## ---- fit ----
ctrl <- control_mcmc(n.sim = 6000, burnin = 1500, thin = 6, h = 1.65 / N^(1/6))
fit  <- SDALGCP2(y ~ x1 + offset(log(pop)), dat, shp, delta = 1.2,
                 phi = phi_grid, method = 3, control.mcmc = ctrl)
summary(fit)

## ---- predict ----
pred_d <- predict(fit, type = "discrete",   sampler = "mcmc",    control.mcmc = ctrl)
pred_c <- predict(fit, type = "continuous", sampler = "laplace", cellsize = 1,
                  control.mcmc = ctrl)

Fit a spatio-temporal SDA-LGCP model (Kronecker-free)

Description

Separable space-time SDA-LGCP for aggregated counts observed over the same N regions at T times. The spatial scale phi is profiled on a grid; the temporal Matern range nu is estimated continuously. The likelihood never forms the (NT)\times(NT) covariance.

Usage

SDALGCP2_ST(
  formula,
  data,
  my_shp,
  times,
  delta,
  phi = NULL,
  kappa = 0.5,
  kappa_t = 0.5,
  method = 3L,
  weighted = FALSE,
  pop_shp = NULL,
  control.mcmc = NULL,
  reanchor = 0L,
  rasters = NULL,
  covariates = NULL,
  confounding = c("none", "restricted"),
  berkson = TRUE,
  max_iter = 10L,
  tol = 0.001,
  messages = FALSE
)

Arguments

Details

With rasters or covariates the covariate surface is taken to be constant over time (time-varying covariates can still be supplied as ordinary columns of data). confounding = "restricted" constrains the space-time random effect to the orthogonal complement of the fixed-effect design and is fitted by an analytic Laplace-marginal likelihood; it reduces to the spatial restricted fit when T = 1 and is not currently combined with rasters/covariates.

Value

an object of class c("SDALGCP2_ST","SDALGCP2") with phi_opt, nu_opt, coefficient table and covariance.

See Also

sdalgcp (friendly wrapper), predict.SDALGCP2_ST

Examples

data(sdalgcp_data)
shp <- sdalgcp_data
## build a 3-time panel (data frame, N*T rows ordered by time then region)
times <- 1:3
dat <- do.call(rbind, lapply(times, function(t) {
  d <- sf::st_drop_geometry(shp); d$time <- t
  d$cases <- rpois(nrow(d), d$pop * exp(-6 + 0.6 * d$x1 + 0.1 * (t - 2)))
  d
}))
fit <- SDALGCP2_ST(cases ~ x1 + offset(log(pop)), dat, shp, times = times,
                   delta = 1.5,
                   control.mcmc = control_mcmc(n.sim = 2000, burnin = 500, thin = 5))
fit$phi_opt; fit$nu_opt

## restricted spatial regression against space-time confounding
fit_c <- SDALGCP2_ST(cases ~ x1 + offset(log(pop)), dat, shp, times = times,
                     delta = 1.5, phi = c(2, 4, 6), confounding = "restricted")
fit_c$beta_opt

## a spatially continuous (raster) covariate, aggregated on the intensity scale
r <- terra::rast(terra::ext(0, 20, 0, 20), resolution = 0.5)
terra::values(r) <- as.numeric(scale(terra::crds(r)[, 1])); names(r) <- "z"
fit_r <- SDALGCP2_ST(cases ~ z + offset(log(pop)), dat, shp, times = times,
                     delta = 1.5, phi = c(2, 4, 6), rasters = r, max_iter = 4,
                     control.mcmc = control_mcmc(n.sim = 2000, burnin = 500, thin = 5))
fit_r$beta_opt

Fit an SDA-LGCP with covariates measured on a different support

Description

Covariates observed on a different support from the outcome (e.g. air- quality monitors at point locations) are kriged to the candidate points and enter the model on the intensity scale with a Berkson correction that propagates the prediction uncertainty (see math/confounding-and-misalignment.pdf).

Usage

SDALGCP2_misaligned(
  formula,
  data,
  delta,
  covariates,
  phi = NULL,
  method = 3L,
  weighted = FALSE,
  pop_shp = NULL,
  berkson = TRUE,
  control.mcmc = NULL,
  max_iter = 10L,
  tol = 0.001,
  messages = FALSE
)

Arguments

Value

an object of class "SDALGCP2" with misaligned = TRUE.

See Also

SDALGCP2_raster, sdalgcp

Examples

data(sdalgcp_data)
set.seed(1)
## a covariate z observed only at 40 scattered monitor points (a different support)
mon <- sf::st_as_sf(data.frame(x = runif(40, 0, 20), y = runif(40, 0, 20)),
                    coords = c("x", "y"))
mon$z <- scale(sf::st_coordinates(mon)[, 1])[, 1]
fit <- SDALGCP2_misaligned(cases ~ z + offset(log(pop)), sdalgcp_data, delta = 1.5,
                           covariates = list(z = mon),
                           control.mcmc = control_mcmc(n.sim = 2000, burnin = 500, thin = 5))
summary(fit)

Fit an SDA-LGCP with spatially continuous (raster) covariates

Description

Covariates supplied as rasters enter the model at the candidate-point level and are aggregated on the intensity (exp) scale via a log-sum-exp offset b_i(\beta)=\log\sum_k w_{ik}\exp(z(x_{ik})^\top\beta) – the statistically correct alternative to averaging the predictor over each polygon. Estimation is a Gauss-Newton fixed point that reuses mcml_fit with the intensity-tilted effective design.

Usage

SDALGCP2_raster(
  formula,
  data,
  my_shp,
  delta,
  rasters,
  phi = NULL,
  method = 3L,
  weighted = FALSE,
  pop_shp = NULL,
  kappa = 0.5,
  tilt_spatial = FALSE,
  control.mcmc = NULL,
  max_iter = 10L,
  tol = 0.001,
  messages = FALSE
)

Arguments

Value

an object of class "SDALGCP2" (as mcml_fit) with extra fields raster = TRUE and n_iter.

See Also

SDALGCP2, mcml_fit

Examples

data(sdalgcp_data)
## a spatially continuous covariate supplied as a raster layer named "z"
r <- terra::rast(terra::ext(0, 20, 0, 20), resolution = 0.5)
terra::values(r) <- as.numeric(scale(terra::crds(r)[, 1]))   # west-east gradient
names(r) <- "z"
df <- sf::st_drop_geometry(sdalgcp_data)
fit <- SDALGCP2_raster(cases ~ z + offset(log(pop)), df, sdalgcp_data,
                       delta = 1.5, rasters = r,
                       control.mcmc = control_mcmc(n.sim = 2000, burnin = 500, thin = 5))
summary(fit)

Coefficient plot of fixed effects (and sigma^2) with confidence intervals

Description

Coefficient plot of fixed effects (and sigma^2) with confidence intervals

Usage

coef_plot(object, level = 0.95, intercept = FALSE)

Arguments

Value

a ggplot object.

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
coef_plot(fit)

Wald confidence intervals for an SDALGCP2 fit

Description

Wald confidence intervals for an SDALGCP2 fit

Usage

## S3 method for class 'SDALGCP2'
confint(object, parm, level = 0.95, ...)

Arguments

Value

a matrix of lower/upper confidence limits.

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
confint(fit)

MCMC control settings for the MALA sampler

Description

MCMC control settings for the MALA sampler

Usage

control_mcmc(
  n.sim = 10000,
  burnin = 2000,
  thin = 8,
  h = NULL,
  c1.h = 0.01,
  c2.h = 1e-04
)

Arguments

Value

a named list consumed by laplace_sampling / the fit.

Examples

## 1000 retained draws (5000 iterations, 2000 burn-in, thin every 3)
ctrl <- control_mcmc(n.sim = 5000, burnin = 2000, thin = 3)
str(ctrl)

Aggregated correlation array (C++)

Description

Aggregated correlation array (C++)

Usage

corr_aggregate_cpp(coords, weights, phi, kappa, weighted, nthreads = 0L)

Arguments

Value

a numeric array of dimension N x N x length(phi).

Aggregated correlation and its phi-derivatives at one phi (C++, Matern)

Description

Aggregated correlation and its phi-derivatives at one phi (C++, Matern)

Usage

corr_and_grad_cpp(coords, weights, phi, kappa, weighted, nthreads = 0L)

Arguments

Value

list with N x N matrices R, dR (dR/dphi) and d2R (d2R/dphi2).

Aggregated cross-covariance between prediction points and regions (C++)

Description

Mirrors compute_cross_cov() for continuous prediction (bottleneck B5): returns an n_pred x N matrix with entries sum_l w_jl * matern(||x_pred - x_jl||, phi, kappa) (weighted) or the mean.

Usage

cross_cov_cpp(pred, coords, weights, phi, kappa, weighted, nthreads = 0L)

Arguments

Exceedance probabilities P(risk > threshold)

Description

Exceedance probabilities P(risk > threshold)

Usage

exceedance(object, thresholds, which = c("adjusted_rr", "relative_risk"))

Arguments

Value

a matrix of exceedance probabilities (locations x thresholds).

See Also

map_exceedance to map them.

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
pr <- predict(fit, type = "discrete")
## P(adjusted relative risk > 1) and > 1.5 for every region
ex <- exceedance(pr, thresholds = c(1, 1.5), which = "adjusted_rr")
head(ex)

Conditional mode and Laplace covariance for [S | Y], Poisson, non-nested (C++)

Description

Conditional mode and Laplace covariance for [S | Y], Poisson, non-nested (C++)

Usage

laplace_mode_poisson_cpp(y, m, mu, Sigma, tol = 1e-08, maxit = 100L)

Arguments

Value

list with mode and Sigma_tilde.

Sample the latent field [S | Y] (Poisson, non-nested) via C++ MALA

Description

Draws posterior samples of the latent Gaussian field for the Poisson, non-nested case. The Laplace mode (Newton step) and the adaptive Metropolis- adjusted Langevin (MALA) loop both run in C++ for speed, with a fixed-seed path for reproducibility.

Usage

laplace_sampling(mu, Sigma, y, units.m, control.mcmc)

Arguments

Value

list with samples (kept x n matrix) and h (step sizes).

Examples

## sample [S | Y] for a tiny 10-unit Poisson example
set.seed(1)
n <- 10
D <- as.matrix(dist(cbind(runif(n), runif(n))))
Sigma <- 0.4 * exp(-D / 0.3)
mu <- rep(log(2), n); m <- rep(100, n)
y <- rpois(n, m * exp(mu + as.numeric(t(chol(Sigma)) %*% rnorm(n))))
out <- laplace_sampling(mu, Sigma, y, m, control_mcmc(n.sim = 2000, burnin = 500, thin = 3))
dim(out$samples)     # (retained draws) x n

Primary biliary cirrhosis incidence in North East England

Description

A real aggregated disease-count dataset: incident primary biliary cirrhosis (a chronic liver disease) cases by Lower-layer Super Output Area (LSOA) in the Newcastle and Gateshead area of North East England, with population and area deprivation covariates. This is the case study of Johnson et al. (2019) and a realistic test bed for the spatial model: cases ~ deprivation + offset(log(pop)).

Usage

liver

Format

An sf object of 545 LSOA polygons (British National Grid, EPSG:27700) with columns:

lsoa

LSOA 2004 census code.

cases

observed incident case count in the LSOA.

pop

population at risk (the offset; use offset(log(pop))).

IMD

Index of Multiple Deprivation score (higher = more deprived).

Income

income-deprivation score.

Employment

employment-deprivation score.

geometry

the LSOA polygon.

Source

Johnson, O., Diggle, P. and Giorgi, E. (2019), "A spatially discrete approximation to log-Gaussian Cox processes for modelling aggregated disease count data", Statistics in Medicine, 38(24), 4871-4884. doi:10.1002/sim.8339. Population and area-deprivation covariates are from the 2004 English indices of deprivation (Lower-layer Super Output Area level). See data-raw/liver.R in the package sources.

See Also

sdalgcp_data for a small simulated example.

Examples

data(liver)
summary(liver$cases)
plot(liver["IMD"])

Adaptive MALA sampler for [S | Y], Poisson, non-nested (C++)

Description

Draw order per iteration is d normals then one uniform, giving reproducible results under a common seed and the same mode/Sigma.tilde.

Usage

mala_poisson_cpp(
  y,
  m,
  mu,
  Sigma,
  mode,
  Sigma_tilde,
  n_sim,
  burnin,
  thin,
  h_init,
  c1,
  c2
)

Arguments

Value

list with samples (kept x d matrix of S draws) and h.

Map exceedance probabilities P(risk > threshold)

Description

Map exceedance probabilities P(risk > threshold)

Usage

map_exceedance(
  x,
  threshold = 1,
  which = c("adjusted_rr", "relative_risk"),
  bound = NULL,
  ...
)

Arguments

Value

a ggplot object.

See Also

exceedance for the underlying probabilities.

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
pr <- predict(fit, type = "discrete")
map_exceedance(pr, threshold = 1.5)           # P(adjusted RR > 1.5)

Importance-sampling diagnostics for an MCML fit

Description

The MCML estimate reweights latent samples drawn at the anchor towards the optimum. When the optimum is far from the anchor the weights become uneven and the estimate unreliable. This reports the effective sample size of the importance weights at the maximiser and a Monte Carlo standard error for the maximised log-likelihood, \mathrm{SE}\approx\sqrt{1/\mathrm{ESS}-1/B}.

Usage

mc_diagnostics(object, warn_frac = 0.1)

Arguments

Value

invisibly, a list with B, ESS, ESS_frac and se_loglik.

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
d <- mc_diagnostics(fit)
d$ESS_frac           # importance-sampling ESS as a fraction of the draws

Monte Carlo maximum likelihood estimation for the spatial SDA-LGCP

Description

Vectorised, Cholesky-based MCML estimation. Simulates the latent field at an anchor, then profiles the importance-sampling MCML objective over the supplied phi grid.

Usage

mcml_fit(
  formula,
  data,
  corr,
  par0 = NULL,
  control.mcmc = NULL,
  phi_method = c("grid", "direct"),
  nugget = FALSE,
  reanchor = 0L,
  reanchor_tol = 0.01,
  messages = FALSE
)

Arguments

Value

an object of class "SDALGCP2" (estimates, covariance, profile, latent samples and metadata).

See Also

SDALGCP2 (the end-to-end wrapper), precompute_corr

Examples

data(sdalgcp_data)
df  <- sf::st_drop_geometry(sdalgcp_data)
pts <- sda_points(sdalgcp_data, delta = 1.2, method = 3)
cc  <- precompute_corr(pts, phi = seq(2, 8, length.out = 6))
fit <- mcml_fit(cases ~ x1 + offset(log(pop)), df, cc,
                control.mcmc = control_mcmc(n.sim = 2000, burnin = 500, thin = 5))
summary(fit)

Posterior-predictive model checking for an SDALGCP2 fit

Description

Compares observed counts with fitted Poisson means, returns Pearson residuals, and tests them for residual spatial autocorrelation with Moran's I. A non-significant Moran's I indicates the spatial random effect has absorbed the spatial structure.

Usage

model_check(object, pred = NULL, nsim = 999, plot = TRUE)

Arguments

Value

invisibly, a list with fitted, residuals and moran.

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
chk <- model_check(fit, plot = FALSE)
chk$moran            # residual Moran's I and its permutation p-value

Profile likelihood and confidence interval for the spatial scale phi

Description

Spline-smoothed profile deviance for phi, with the coverage-level confidence interval where the deviance crosses the chi-squared cutoff.

Usage

phi_profile(object, coverage = 0.95, plot = TRUE)

Arguments

Value

invisibly, a list with the interval and the smoothed profile; a ggplot is drawn when plot = TRUE.

Examples

data(sdalgcp_data)
## profile phi on a grid (scale = "grid") so there is a deviance curve to draw
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(scale = "grid", n_sim = 2000,
                                         burnin = 500, thin = 5, reanchor = 0))
phi_profile(fit)

Plot an SDALGCP2 fit (the phi profile deviance)

Description

Plot an SDALGCP2 fit (the phi profile deviance)

Usage

## S3 method for class 'SDALGCP2'
plot(x, ...)

Arguments

Value

invisibly, the profile (see phi_profile).

Examples

data(sdalgcp_data)
fit <- SDALGCP2(cases ~ x1 + offset(log(pop)),
                sf::st_drop_geometry(sdalgcp_data), sdalgcp_data, delta = 1.2,
                control.mcmc = control_mcmc(n.sim = 2000, burnin = 500, thin = 5))
plot(fit)   # profile deviance for the spatial scale phi

Map a spatio-temporal prediction for one time

Description

Maps a chosen quantity ("relative_risk", "adjusted_rr", "relative_risk_se", "adjusted_rr_se" or "exceedance") for a selected time slice of a spatio-temporal prediction.

Usage

## S3 method for class 'SDALGCP2_ST_pred'
plot(
  x,
  time = attr(x, "times")[1],
  what = c("relative_risk", "adjusted_rr", "relative_risk_se", "adjusted_rr_se",
    "exceedance"),
  threshold = 1,
  which = c("adjusted_rr", "relative_risk"),
  ...
)

Arguments

Value

a ggplot object.

Examples

data(sdalgcp_data)
times <- 1:3
panel <- do.call(rbind, lapply(times, function(t) {
  d <- sdalgcp_data; d$time <- t
  d$cases <- rpois(nrow(d), d$pop * exp(-6 + 0.6 * d$x1 + 0.1 * (t - 2)))
  d
}))
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = panel, time = "time",
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
pr <- predict(fit)
plot(pr, time = 2)                 # one time slice
plot(pr, time = NULL)              # facet all times
plot(pr, what = "exceedance", threshold = 1.2, time = 3)

Map a fitted SDALGCP2 prediction

Description

Maps any of the four predicted quantities from predict.SDALGCP2 – the relative risk "relative_risk", the covariate-adjusted relative risk "adjusted_rr", or their standard errors "relative_risk_se"/"adjusted_rr_se" – for either discrete (choropleth) or continuous (raster) predictions.

Usage

## S3 method for class 'SDALGCP2_pred'
plot(
  x,
  variable = c("relative_risk", "adjusted_rr", "relative_risk_se", "adjusted_rr_se"),
  bound = NULL,
  midpoint = NULL,
  title = NULL,
  ...
)

Arguments

Value

a ggplot object.

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
pr <- predict(fit, type = "discrete")
plot(pr, variable = "relative_risk")          # choropleth of relative risk
plot(pr, variable = "adjusted_rr_se")         # its uncertainty

Map an sdalgcp fit

Description

Predicts and maps a chosen quantity. Works for spatial fits (discrete or continuous) and spatio-temporal fits (select a time).

Usage

## S3 method for class 'sdalgcp'
plot(
  x,
  what = c("relative_risk", "adjusted_rr", "relative_risk_se", "adjusted_rr_se",
    "exceedance"),
  type = c("discrete", "continuous"),
  time = NULL,
  threshold = 1,
  which = c("adjusted_rr", "relative_risk"),
  cellsize = NULL,
  sampler = c("mcmc", "laplace"),
  ...
)

Arguments

Value

a ggplot object.

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
plot(fit)                              # relative-risk map (predicts internally)
plot(fit, what = "exceedance", threshold = 1.5)

Precompute aggregated region-level correlation matrices

Description

Builds the N \times N \times length(phi) array of region-level correlations used by the SDA-LGCP model, where

R(\phi)_{ij} = \sum_{k,l} w_{ik} w_{jl}\, C(\lVert x_{ik}-x_{jl}\rVert; \phi, \kappa)

(population-weighted) or the unweighted mean over candidate-point pairs. The heavy reduction runs in C++ (OpenMP-parallel over region pairs).

Usage

precompute_corr(points, phi, kappa = 0.5, weighted = NULL, nthreads = 0L)

Arguments

Value

a list with R (the correlation array) and phi, carrying weighted, my_shp and S_coord attributes on R.

See Also

sda_points, mcml_fit

Examples

data(sdalgcp_data)
pts <- sda_points(sdalgcp_data, delta = 1.2, method = 3)
cc  <- precompute_corr(pts, phi = c(2, 4, 6))
dim(cc$R)            # N x N x length(phi)

Reference (pure-R) aggregated correlation builder

Description

Slow but dependency-light implementation kept for correctness testing and benchmarking against precompute_corr. Computes the exponential (kappa = 0.5) aggregated correlation only.

Usage

precompute_corr_ref(points, phi, weighted = NULL)

Arguments

Value

the N \times N \times length(phi) correlation array.

Predict relative risk from a fitted SDALGCP2 model

Description

Predict relative risk from a fitted SDALGCP2 model

Usage

## S3 method for class 'SDALGCP2'
predict(
  object,
  type = c("discrete", "continuous"),
  sampler = c("mcmc", "laplace"),
  cellsize = NULL,
  pred.loc = NULL,
  control.mcmc = NULL,
  ...
)

Arguments

Value

an sf (class c("SDALGCP2_pred", "sf", "data.frame")) with one row per location – polygons for type = "discrete", grid-cell points for type = "continuous" – carrying the posterior mean and standard error of two relative-risk quantities:

relative_risk, relative_risk_se

the relative risk \exp(d'\beta + S) – the fitted risk relative to the offset baseline, combining the covariate effect and the residual spatial variation. This is the headline disease-mapping quantity.

adjusted_rr, adjusted_rr_se

the covariate-adjusted relative risk \exp(S) – the purely spatial relative risk that remains after holding the covariates fixed (the spatial signal the covariates do not explain).

The full posterior draws are retained as object attributes so that exceedance and map_exceedance can be computed for either quantity. Map a column with plot.SDALGCP2_pred.

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))

## region-level (discrete) prediction: an sf you can map or st_write()
pr <- predict(fit, type = "discrete")
head(pr)                       # relative_risk / adjusted_rr (+ standard errors)
plot(pr, variable = "relative_risk")

## continuous surface on a grid
pr_c <- predict(fit, type = "continuous", cellsize = 1)

Discrete (region x time) prediction for a spatio-temporal fit

Description

Draws the latent field at the fitted optimum and returns posterior mean and SD of the incidence relative risk \exp(\mu+S) and covariate-adjusted relative risk \exp(S) for every region and time.

Usage

## S3 method for class 'SDALGCP2_ST'
predict(object, control.mcmc = NULL, ...)

Arguments

Value

a long sf of class c("SDALGCP2_ST_pred", "sf", "data.frame") with one row per region and time (ordered region-fastest within each time block) and columns region, time, relative_risk, relative_risk_se (\exp(\mu+S)), adjusted_rr and adjusted_rr_se (\exp(S)) – the same column names as the spatial predict.SDALGCP2. The posterior draws are kept in object attributes (for exceedance); map a time slice with plot.SDALGCP2_ST_pred.

Examples

data(sdalgcp_data)
## stack the spatial example into a 3-time panel with a mild temporal trend
times <- 1:3
panel <- do.call(rbind, lapply(times, function(t) {
  d <- sdalgcp_data; d$time <- t
  d$cases <- rpois(nrow(d), d$pop * exp(-6 + 0.6 * d$x1 + 0.1 * (t - 2)))
  d
}))
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = panel, time = "time",
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
pr <- predict(fit)        # a long sf: region x time
head(pr)
plot(pr, time = 2)        # map the relative risk at time 2

Predict relative risk from an sdalgcp fit

Description

Returns a prediction object carrying, for every location, the posterior mean and standard error of the relative risk relative_risk (\exp(\eta)=\exp(d'\beta+S)) and the covariate-adjusted relative risk adjusted_rr (\exp(S)). Map it with plot() and get hotspot probabilities with exceedance.

Usage

## S3 method for class 'sdalgcp'
predict(
  object,
  type = c("discrete", "continuous"),
  sampler = c("mcmc", "laplace"),
  cellsize = NULL,
  ...
)

Arguments

Value

for a spatial fit, an sf of class "SDALGCP2_pred" with relative_risk, relative_risk_se, adjusted_rr and adjusted_rr_se columns (polygons for type = "discrete", grid points for "continuous"); for a spatio-temporal fit, an "SDALGCP2_ST_pred" object (see predict.SDALGCP2_ST).

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
pr <- predict(fit)            # discrete by default; an sf of relative risks
head(pr)

Print an SDALGCP2 fit

Description

Print an SDALGCP2 fit

Usage

## S3 method for class 'SDALGCP2'
print(x, ...)

Arguments

Value

x, invisibly.

Print a summary of an SDALGCP2 fit

Description

Print a summary of an SDALGCP2 fit

Usage

## S3 method for class 'summary.SDALGCP2'
print(x, ...)

Arguments

Value

x, invisibly.

One-call panel of post-fit graphics

Description

Returns the maps and summaries an analyst usually wants after fitting: relative-risk and uncertainty maps, an exceedance map, the coefficient plot and the phi profile. The pieces are returned as a named list of ggplot objects so they can be arranged or printed individually.

Usage

report(object, pred = NULL, threshold = 1.5, ...)

Arguments

Value

a named list of ggplot objects.

Examples

data(sdalgcp_data)
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = sdalgcp_data,
               control = sdalgcp_control(n_sim = 2000, burnin = 500, thin = 5,
                                         reanchor = 0))
figs <- report(fit, threshold = 1.5)
names(figs)          # relative_risk, uncertainty, exceedance, coefficients, ...
figs$relative_risk   # print one of the maps

Generate candidate sampling points inside each region

Description

For every polygon feature in my_shp it produces candidate points and aggregation weights, in the list format consumed by precompute_corr.

Usage

sda_points(
  my_shp,
  delta,
  method = 1L,
  weighted = FALSE,
  pop_shp = NULL,
  rho = 0.55,
  giveup = 1000L
)

Arguments

Value

a list of length nrow(my_shp); each element has xy and weight. Carries "weighted" and "my_shp" attributes.

See Also

precompute_corr, which consumes this output.

Examples

data(sdalgcp_data)
pts <- sda_points(sdalgcp_data, delta = 1.2, method = 3)   # regular grid points
length(pts)          # one entry per region
str(pts[[1]])        # $xy candidate coordinates and $weight

Fit a spatially discrete LGCP model for aggregated counts

Description

The main user interface, designed to feel like glm: give a formula and an sf data object and it does the rest. The same call covers three settings, chosen from the arguments you supply:

• spatial (default): sdalgcp(y ~ x + offset(log(pop)), data);

• raster covariates: add rasters = a SpatRaster whose layers are named in the formula – these enter on the intensity scale (see SDALGCP2_raster);

• spatio-temporal: add time = the name of a time column.

Usage

sdalgcp(
  formula,
  data,
  time = NULL,
  rasters = NULL,
  covariates = NULL,
  popden = NULL,
  control = sdalgcp_control(),
  verbose = FALSE
)

Arguments

Value

a fitted model object of class c("sdalgcp", ...) with print, summary, confint, predict and plot methods.

See Also

predict.sdalgcp, sdalgcp_control, SDALGCP2, SDALGCP2_raster, SDALGCP2_ST

Examples

library(sf)
set.seed(1)
grid <- st_make_grid(st_as_sfc(st_bbox(c(xmin = 0, ymin = 0, xmax = 20, ymax = 20))),
                     n = c(8, 8))
regions <- st_sf(geometry = grid)
regions$x1  <- as.numeric(scale(st_coordinates(st_centroid(regions))[, 1]))
regions$pop <- round(runif(nrow(regions), 500, 3000))
regions$cases <- rpois(nrow(regions), regions$pop * exp(-6 + 0.5 * regions$x1))

fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = regions)  # that's it
summary(fit)
rr <- predict(fit)          # an sf you can plot() directly
plot(fit)                   # default relative-risk map

Control settings for sdalgcp

Description

Bundles the technical knobs so that a default fit needs none of them.

Usage

sdalgcp_control(
  delta = NULL,
  points_per_region = 16,
  point_method = c("regular", "uniform", "ssi"),
  scale = c("continuous", "grid"),
  phi = NULL,
  kappa = 0.5,
  kappa_t = 0.5,
  nugget = FALSE,
  confounding = c("none", "restricted"),
  reanchor = 2L,
  n_sim = 10000L,
  burnin = 2000L,
  thin = 8L,
  tilt_spatial = FALSE,
  nthreads = 0L
)

Arguments

Value

a list of control settings.

See Also

sdalgcp

Examples

## defaults, then a faster grid-based fit with a nugget term
str(sdalgcp_control())
ctrl <- sdalgcp_control(scale = "grid", nugget = FALSE, n_sim = 4000,
                        burnin = 1000, thin = 6)

Simulated aggregated disease-count data

Description

A small, self-contained example dataset used throughout the help pages and vignettes. It is simulated from the model the package fits: an 8x8 lattice of regions, a spatially structured covariate, a latent Gaussian spatial field with exponential covariance, and Poisson counts with a population offset. The true fixed effects are (Intercept) = -6 and x1 = 0.6; the latent field has variance \sigma^2 = 0.3 and exponential scale \phi = 4.

Usage

sdalgcp_data

Format

An sf object of 64 POLYGON regions with columns:

region

integer region identifier (1-64).

cases

observed disease count in the region.

x1

a standardised, spatially structured covariate.

pop

population at risk (the offset; use offset(log(pop))).

geometry

the region polygon.

Source

Simulated; see data-raw/sdalgcp_data.R in the package sources.

See Also

liver for a real disease-count example.

Examples

data(sdalgcp_data)
summary(sdalgcp_data$cases)
plot(sdalgcp_data["cases"])

Summary of an SDALGCP2 fit

Description

Summary of an SDALGCP2 fit

Usage

## S3 method for class 'SDALGCP2'
summary(object, ...)

Arguments

Value

an object of class "summary.SDALGCP2" with a coefficient table.