Advanced features

Model

We assume that the response variable \(Y_{it}\) for region \(i\) at time \(t\) follows a normal distribution given the partition \(M\): \[ Y_{it} \mid \mu_{it}, M \stackrel{ind}{\sim} \text{Normal}(\mu_{it}, \sigma^2), \] where the mean \(\mu_{it}\) is modeled based on the cluster assignment \(c_i\): \[ \mu_{it} = \alpha_{c_i} + f_{c_i,t}, \] where \(\alpha_{c_i}\) is the intercept for cluster \(c_i\), and \(f_{c_i,t}\) is a latent random effect modeled as a random walk process. Specifically, we impose the condition: \[ f_{c_i,t} - f_{c_i,t-1} \sim N(0, \nu^{-1}), \quad \text{for } t = 2, \dots, n. \]

The prior for the hyperparameter \(\nu_c\) is defined as: \[ \log(\nu_c) \sim \text{Normal}(-2, 1). \]

Initial clustering

sfclust uses an undirected graph to represent connections between regions and proposes spatial clusters using this graph through minimum spanning trees (MST). By default, sfclust accepts the argument graphdata, which should include:

• An undirected igraph object representing spatial connections,
• An MST of that graph, and
• A membership vector indicating the initial cluster assignments for each region.

For simplicity, you can use the genclust function to generate an initial random partitioning with a specified number of clusters. In this example, we create a partition with 20 clusters:

set.seed(123)
initial_cluster <- genclust(stgaus, nclust = 20)
names(initial_cluster)

#> [1] "graph"      "mst"        "membership" "valid_ids"

Now, let’s visualize how the regions were randomly clustered. Panel (A) shows the full adjacency graph (all neighbor connections), panel (B) shows the MST derived from it (the backbone used for cluster proposals), and panel (C) shows the resulting initial partition:

gg1 <- ggraph(initial_cluster$graph, layout = st_coordinates(st_centroid(st_geometry(stgaus)))) +
    geom_edge_fan(linetype = 1, color = 2) +
    geom_node_point(size = 1.5, color = 1) +
    geom_sf(data = st_geometry(stgaus), fill = NA, color = 1, linewidth = 0.5) +
    labs(subtitle = "(A)") +
    theme_void()
gg2 <- ggraph(initial_cluster$mst, layout = st_coordinates(st_centroid(st_geometry(stgaus)))) +
    geom_edge_fan(linetype = 1, color = 2) +
    geom_node_point(size = 1.5, color = 1) +
    geom_sf(data = st_geometry(stgaus), fill = NA, color = 1, linewidth = 0.5) +
    labs(subtitle = "(B)") +
    theme_void()
gg3 <- st_sf(st_geometry(stgaus), cluster = factor(initial_cluster$membership)) |>
  ggplot() +
    geom_sf(aes(fill = cluster), color = 1) +
    labs(subtitle = "(C)") +
    theme_void() +
    theme(legend.position = "none")
gg1 + gg2 + gg3 & theme(plot.margin = margin(0, 0, 0, 0))

Sampling with sfclust

We will initiate the Bayesian clustering algorithm using our generated partition (initial_cluster). Additionally, we use the capabilites of INLA::inla to define a model that includes a random walk process in the linear predictor, and custom priors for the scale hyperparameter. To begin, we will run the algorithm for only 50 iterations.

The logpen argument sets the log-penalty on the number of clusters, specifically log(1-q) where q is the prior probability of keeping a cluster. A large negative value such as -50 favors fewer clusters, acting as a strong parsimony prior. The burnin argument discards the first iterations before saving samples, and thin retains only every thin-th iteration, reducing autocorrelation in the chain.

result0 <- sfclust(stgaus, graphdata = initial_cluster, logpen = -50,
  formula = y ~ f(id_time, model = "rw1",
                  hyper = list(prec = list(prior = "normal", param = c(-2, 1)))),
  niter = 50, burnin = 10, thin = 2, nmessage = 10
)
result0

#> Within-cluster formula:
#> y ~ f(id_time, model = "rw1", hyper = list(prec = list(prior = "normal", 
#>     param = c(-2, 1))))
#> 
#> Clustering hyperparameters:
#> log(1-q)    birth    death   change    hyper 
#>  -50.000    0.425    0.425    0.100    0.050 
#> 
#> Clustering movement counts:
#>  births  deaths changes  hypers 
#>      11       8       1       2 
#> 
#> Log marginal likelihood (sample 25 out of 25): 2570.063

The output indicates that after starting with 20 clusters, the algorithm created 11 new clusters (births), removed 8 clusters (deaths), changed the membership of 1 cluster, and modified the minimum spanning tree (MST) 2 times.

The plot method allows us to select which graph to produce. For example, we can visualize only the log marginal likelihood to diagnose convergence. With just 50 iterations, we can see that the log marginal likelihood has not yet achieved convergence.

plot(result0, which = 3)

Continue sampling

To achieve convergence, we can continue the sampling using the update function and specify the number of new iterations with the niter argument:

result <- update(result0, niter = 1000)
result

#> Within-cluster formula:
#> y ~ f(id_time, model = "rw1", hyper = list(prec = list(prior = "normal", 
#>     param = c(-2, 1))))
#> 
#> Clustering hyperparameters:
#> log(1-q)    birth    death   change    hyper 
#>  -50.000    0.425    0.425    0.100    0.050 
#> 
#> Clustering movement counts:
#>  births  deaths changes  hypers 
#>      33      44       7      43 
#> 
#> Log marginal likelihood (sample 1000 out of 1000): 12602.37

With 1000 additional iterations, there have been many clustering movements. Furthermore, when visualizing the results, we can observe that the log marginal likelihood has achieved convergence.

plot(result, which = 3)

The side-by-side comparison below highlights the contrast: panel (A) shows the unstable log marginal likelihood from the short initial run, while panel (B) shows the stabilized trajectory after the full run, indicating convergence.

gg1 <- plot(result0, which = 3) + labs(subtitle = "(A)")
gg2 <- plot(result, which = 3) + labs(subtitle = "(B)")
gg1 + gg2

if (save_figures) {
  ggsave(file.path(path_figures, "stgaus-marginal-likelihood.pdf"), width = 10, height = 4.5,
    device = cairo_pdf)
}

The update function can also be used to select a different past iteration as the active clustering, without running additional MCMC. This is useful when you want to inspect or use a specific sample rather than the last one:

result_other <- update(result, sample = 750)
result_other

#> Within-cluster formula:
#> y ~ f(id_time, model = "rw1", hyper = list(prec = list(prior = "normal", 
#>     param = c(-2, 1))))
#> 
#> Clustering hyperparameters:
#> log(1-q)    birth    death   change    hyper 
#>  -50.000    0.425    0.425    0.100    0.050 
#> 
#> Clustering movement counts:
#>  births  deaths changes  hypers 
#>      33      44       7      43 
#> 
#> Log marginal likelihood (sample 750 out of 1000): 12319.21