The default SPI model uses a negative binomial distribution for
counts and no additional explanatory variables (covariates). The
negative binomial distribution allows counts to vary more than a Poisson
distribution would. This article shows how to assess the effects of
covariates, observation models, spatial structure, and the training
period. These options apply to the model-based SPI
(spi_index()), not to the simple SPI.
The examples use the bundled synthetic data. Model fitting requires INLA and can take several minutes.
Adjusting for covariates
Covariates are optional and are not part of the default SPI used in
the accompanying study. Adding covariates changes the expected counts
and therefore what SPI compares. For example, including travel time
makes the expected count depend partly on access to care. Pass
district-year covariates to spi_index() or
spi_expected() through the covariates
argument. The example data include DTP3 coverage, urban population
share, and travel time to care:
synth$covariates |>
dplyr::slice_head(n = 6)
#> adm2_guid year dtp3 urban_prop travel_time_min
#> 1 {01325AA0-BEA1-66FE-9B5C-88AA603382F3} 2015 84.8 0.566 15.9
#> 2 {01325AA0-BEA1-66FE-9B5C-88AA603382F3} 2016 89.9 0.566 16.5
#> 3 {01325AA0-BEA1-66FE-9B5C-88AA603382F3} 2017 89.9 0.566 17.1
#> 4 {01325AA0-BEA1-66FE-9B5C-88AA603382F3} 2018 87.3 0.566 23.2
#> 5 {01325AA0-BEA1-66FE-9B5C-88AA603382F3} 2019 89.9 0.566 19.3
#> 6 {01325AA0-BEA1-66FE-9B5C-88AA603382F3} 2020 87.6 0.566 22.1To inspect the covariate effects, fit the model once with
spi_expected(), log-transforming the skewed travel-time
column:
fit_adj <- spi_expected(
cases = synth$cases,
population = synth$population,
adjacency = synth$boundaries,
covariates = synth$covariates,
log_transform = "travel_time_min",
id_col = "adm2_guid",
n_draws = 500,
seed = 42,
verbose = FALSE
)Effects are reported as rate ratios for a one-standard-deviation increase in each covariate. The model first rescales each covariate to have mean 0 and standard deviation 1. Here are the three covariate effects, excluding the seasonal terms:
summary(fit_adj)$effects |>
dplyr::filter(covariate %in% c("dtp3", "urban_prop", "travel_time_min")) |>
dplyr::select(covariate, rr_median, rr_q025, rr_q975, signif)
#> # A tibble: 3 x 5
#> covariate rr_median rr_q025 rr_q975 signif
#> <chr> <dbl> <dbl> <dbl> <lgl>
#> 1 dtp3 1.00 0.953 1.05 FALSE
#> 2 urban_prop 1.03 0.938 1.13 FALSE
#> 3 travel_time_min 0.991 0.968 1.01 FALSEThese coefficients describe associations with reporting after accounting for the other model terms. They do not show that changing a covariate would cause reporting to change.
Compare overdispersion options
Compare Poisson, Poisson with an observation-level random effect, and negative binomial models:
od <- spi_compare_overdispersion(
cases = synth$cases,
population = synth$population,
adjacency = synth$boundaries,
id_col = "adm2_guid",
specs = c("none", "iid", "nb"),
n_draws = 100,
verbose = FALSE
)
od$summary
#> # A tibble: 3 x 11
#> spec n_obs dic waic p_eff sd_spatial phi_spatial phi_pegged sd_extra
#> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl> <dbl>
#> 1 none 28320 81917. 81996. 241. 0.865 0.928 FALSE NA
#> 2 iid 28320 79563. 79528. 5667. 0.864 0.926 FALSE 0.423
#> 3 nb 28320 80299. 80296. 240. 0.867 0.934 FALSE NA
#> # i 2 more variables: cpo_valid <dbl>, pit_ks <dbl>od$summary has one row per model. DIC (deviance
information criterion) and WAIC (widely applicable information
criterion) compare fit while allowing for model complexity; lower values
are preferred. cpo_valid is the share of observations with
a valid conditional predictive ordinate, a check of how well each
observation is predicted when left out of the fit. pit_ks
measures how far the probability integral transform (PIT) values are
from a uniform distribution; smaller values indicate closer agreement
with that reference.
od$recommendation gives the model selected by the
package’s diagnostic rule, with its reasoning:
od$recommendation$choice
#> [1] "none"
od$recommendation$reasoning
#> [1] "Excluded: effective parameter count exceeds 20% of observations (iid). Smallest PIT distance among remaining models (none, nb): none."How the diagnostic rule selects a model
WAIC can favour a model with one random effect per observation even
when that model overfits. spi_compare_overdispersion()
therefore first excludes any model with valid CPO values for fewer than
half of the observations, an effective number of parameters above 20% of
the observations, or a BYM2 phi below 0.02 or above 0.98.
Among the remaining models, it selects the one with the smallest PIT
distance.
To skip the manual step,
spi_expected(overdispersion = "auto") runs this same
comparison internally and refits with the selected model. Here the
diagnostic rule selects overdispersion = "none". Negative
binomial remains the package default because it is the pre-specified
model used in the accompanying study;
overdispersion = "auto" is an optional diagnostic workflow
for other applications.
Spatial structure
spi_expected() builds the neighbour graph itself when
adjacency is an sf object. To inspect or reuse
the graph, build it with spi_adjacency():
adj <- spi_adjacency(synth$boundaries, id_col = "adm2_guid")
adj
#> Neighbour list object:
#> Number of regions: 236
#> Number of nonzero links: 1,320
#> Percentage nonzero weights: 2.3700
#> Average number of links: 5.5932Districts are neighbours when their boundaries touch at any point
(contiguity = "queen"); contiguity = "rook"
requires a shared edge. With handle_islands = TRUE, the
default, a district with no neighbours is linked to its nearest district
by centroid distance. A group of connected districts that has no link to
the rest remains separate. The printed graph calls these groups
“disjoint connected subgraphs”.
The spatial term is a BYM2 random effect (Riebler et al. 2016). It
splits district-level variation into a spatially structured part, shared
with neighbours, and an unstructured part specific to each district;
phi is the share of that variation that is spatially
structured. When a district has little data, the model uses information
from other districts, including its neighbours, to help estimate its
expected count. This makes the estimate more stable, but it can pull an
unusual district towards its neighbours.
Sensitivity checks
Results can depend on the training period, observation model, covariates, centring, and reporting year. Each check below calculates SPI for 2022 to 2024 with one setting changed, using 100 posterior draws to keep the run time down. Each check counts district-years that cover all 12 months and have an SPI below 1. The count is a simple illustration of how classification changes between settings, not a measure of model performance:
sensitivity_index <- function(...) {
spi_index(
cases = synth$cases,
population = synth$population,
adjacency = synth$boundaries,
first_assessment = 2022,
last_assessment = 2024,
id_col = "adm2_guid",
n_draws = 100,
seed = 42,
verbose = FALSE,
...
)
}
count_below_1 <- function(index) {
index |>
tibble::as_tibble() |>
dplyr::filter(n_months == 12) |>
dplyr::summarise(
district_years = dplyr::n(),
below_1 = sum(spi_median < 1)
)
}
dplyr::bind_rows(
default = count_below_1(sensitivity_index()),
centre_none = count_below_1(sensitivity_index(centre = "none")),
may_to_april = count_below_1(sensitivity_index(year_end_month = 4)),
covariates = count_below_1(sensitivity_index(
covariates = synth$covariates,
log_transform = "travel_time_min"
)),
poisson_iid = count_below_1(sensitivity_index(overdispersion = "iid")),
.id = "setting"
)
#> # A tibble: 5 x 3
#> setting district_years below_1
#> <chr> <int> <int>
#> 1 default 708 410
#> 2 centre_none 708 337
#> 3 may_to_april 708 402
#> 4 covariates 708 409
#> 5 poisson_iid 708 404Starting assessment in 2021 skips the earliest assessment years. It does not change the training history for years already included, because each fit still uses all preceding years:
spi_later <- spi_index(
cases = synth$cases,
population = synth$population,
adjacency = synth$boundaries,
first_assessment = 2021,
id_col = "adm2_guid"
)Keep the chosen settings the same when comparing years, and report them with the results.