Areal kernel demonstration

Abstract

Task We briefly demonstrate the fitting of areal kernel models in Stan.

Code to fit centroid kernel

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

sf <- mw

centroid_distance <- function(sf) {
  cent <- sf::st_centroid(sf)
  D <- sf::st_distance(cent, cent)
  return(D)
}

D <- centroid_distance(sf)

## Warning in st_centroid.sf(sf): st_centroid assumes attributes are constant over geometries of x

dat <- list(n = nrow(sf),
            y = sf$y,
            m = sf$n_obs,
            mu = rep(0, nrow(sf)),
            D = D)

nsim_warm <- 100
nsim_iter <- 1000

fit_centroid <- rstan::stan(
  "centroid.stan",
  data = dat,
  warmup = nsim_warm,
  iter = nsim_iter
)

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

Code to fit integrated kernel

sf <- mw

The number of integration points within each area:

L <- 10

The method for selecting integration points:

type <- "hexagonal"

The number of areas:

n <- nrow(sf)

Draw \(L\) samples from each area according to method type:

samples <- sf::st_sample(sf, type = type, size = rep(L, n))

plot_samples <- function(samples){
  ggplot(sf) +
    geom_sf(fill = "lightgrey") +
    geom_sf(data = samples, alpha = 0.5, shape = 4) +
    labs(x = "Longitude", y = "Latitude") +
    theme_minimal() +
    labs(fill = "") +
    theme(axis.title = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank())
}

plot_samples(samples)

Construct an \(nL \times nL\) matrix containing the Euclidean distance between each sample. Note that this the \(L\) actually obtained will not be exact for some settings of type (hexagonal, regular):

S <- sf::st_distance(samples, samples)
dim(S)

## [1] 280 280

We use a start index data structure for unequal number of points in each area. Note that other possible ways to do this include a padded data structure (matrix(unlist(padded_sample_index), nrow = n, ncol = max_length, byrow = TRUE)) and a database structure group_db <- data.frame(i = unlist(sample_index), group_id = rep(seq(length(sample_index)), lengths(sample_index)))).

sample_index <- sf::st_intersects(sf, samples)
sample_lengths <- lengths(sample_index)
start_index <- sapply(sample_index, function(x) x[1])

dat <- list(n = nrow(sf),
            y = sf$y,
            m = sf$n_obs,
            mu = rep(0, nrow(sf)),
            sample_lengths = sample_lengths,
            total_samples = sum(sample_lengths),
            start_index = start_index,
            S = S)

fit_integrated <- rstan::stan(
  "integrated.stan",
  data = dat,
  warmup = nsim_warm,
  iter = nsim_iter
)

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

Comparison

Sampling speed

samples_per_second <- function(fit) {
  # Outside of warm-up phase
  times <- rstan::get_elapsed_time(fit)
  # Number of samples in first chain
  nsim_iter <- fit@stan_args[[1]]$iter
  nsim_iter / mean(times[, 2])
}
samples_per_second(fit_centroid)

## [1] 196.4958

samples_per_second(fit_integrated)

## [1] 3.540992

Convergence analysis

Trace-plots:

bayesplot::mcmc_trace(fit_centroid, pars = "l")

bayesplot::mcmc_trace(fit_centroid, pars = "beta_0")

bayesplot::mcmc_trace(fit_centroid, pars = "rho[28]") # For example

bayesplot::mcmc_trace(fit_integrated, pars = "l")

bayesplot::mcmc_trace(fit_integrated, pars = "beta_0")

bayesplot::mcmc_trace(fit_integrated, pars = "rho[28]")

Rank histograms:

bayesplot::mcmc_rank_overlay(fit_centroid, pars = "l")

bayesplot::mcmc_rank_overlay(fit_centroid, pars = "beta_0")

bayesplot::mcmc_rank_overlay(fit_centroid, pars = "rho[28]")

bayesplot::mcmc_rank_overlay(fit_integrated, pars = "l")

bayesplot::mcmc_rank_overlay(fit_integrated, pars = "beta_0")

bayesplot::mcmc_rank_overlay(fit_integrated, pars = "rho[28]")

All Rhat < 1.1 necessary but not sufficient:

which(rstan::summary(fit_centroid)[["summary"]][, "Rhat"] > 1.1)

## named integer(0)

which(rstan::summary(fit_integrated)[["summary"]][, "Rhat"] > 1.1)

## named integer(0)

Fitted values comparison:

fitted_centroid <- summary(fit_centroid)$summary[paste0("rho[", 1:28, "]"), "mean"]
fitted_integrated <- summary(fit_integrated)$summary[paste0("rho[", 1:28, "]"), "mean"]
plot(fitted_centroid, fitted_integrated)

Stan code

centroid.stan

// centroid.stan

functions {
  real xbinomial_logit_lpdf(real y, real m, real eta) {
    real dens = lchoose(m, y) + y * log(inv_logit(eta)) + (m - y) * log(1 - inv_logit(eta));
    real jacobian = log(inv_logit(eta)) + log(inv_logit(-eta)); // Alternative: eta - 2 * log(1 + exp(eta))
    return dens + jacobian;
  }

  matrix cov_matern32(matrix D, real l) {
    int n = rows(D);
    matrix[n, n] K;
    real norm_K;
    real sqrt3;
    sqrt3 = sqrt(3.0);

    // Diagonal entries
    for(i in 1:n){
      K[i, i] = 1;
    }

    for(i in 1:(n - 1)){
      for(j in (i + 1):n){
        norm_K = D[i, j] / l;
        K[i, j] = (1 + sqrt3 * norm_K) * exp(-sqrt3 * norm_K); // Fill lower triangle
        K[j, i] = K[i, j]; // Fill upper triangle
      }
    }
    return(K);
  }
}

data {
  int<lower=0> n; // Number of regions
  vector[n] y; // Vector of observed responses
  vector[n] m; // Vector of sample sizes
  vector[n] mu; // Prior mean vector
  matrix[n, n] D; // Distances between centroids
}

parameters {
  real beta_0; // Intercept
  vector[n] phi; // Spatial effects
  real<lower=0> sigma_phi; // Standard deviation of spatial effects
  real<lower=0> l; // Kernel lengthscale
}

transformed parameters {
  vector[n] eta = beta_0 + sigma_phi * phi;
}

model {
  matrix[n, n] K = cov_matern32(D, l);
  // I could do this?
    // matrix[n, n] L = cholesky_decompose(K);
    // y ~ multi_normal_cholesky(mu, L);
    l ~ gamma(1, 1);
    sigma_phi ~ normal(0, 2.5); // Weakly informative prior
    beta_0 ~ normal(-2, 1);
    phi ~ multi_normal(mu, K);
    for(i in 1:n) {
      y[i] ~ xbinomial_logit(m[i], eta[i]);
    }
}

generated quantities {
  real tau_phi = 1 / sigma_phi^2; // Precision of spatial effects
  vector[n] rho = inv_logit(beta_0 + sigma_phi * phi);
  vector[n] log_lik;
  for (i in 1:n) {
    log_lik[i] = xbinomial_logit_lpdf(y[i] | m[i], eta[i]);
  }
}

integrated.stan

// integrated.stan

functions {
  real xbinomial_logit_lpdf(real y, real m, real eta) {
    real dens = lchoose(m, y) + y * log(inv_logit(eta)) + (m - y) * log(1 - inv_logit(eta));
    real jacobian = log(inv_logit(eta)) + log(inv_logit(-eta)); // Alternative: eta - 2 * log(1 + exp(eta))
    return dens + jacobian;
  }

  matrix cov_matern32(matrix D, real l) {
    int n = rows(D);
    matrix[n, n] K;
    real norm_K;
    real sqrt3;
    sqrt3 = sqrt(3.0);

    for(i in 1:(n - 1)){
      // Diagonal entries (apart from the last)
      K[i, i] = 1;
      for(j in (i + 1):n){
        // Off-diagonal entries
        norm_K = D[i, j] / l;
        K[i, j] = (1 + sqrt3 * norm_K) * exp(-sqrt3 * norm_K); // Fill lower triangle
        K[j, i] = K[i, j]; // Fill upper triangle
      }
    }
    K[n, n] = 1;

    return(K);
  }

  matrix cov_sample_average(matrix S, real l, int n, int[] start_index, int[] sample_lengths, int total_samples) {
    matrix[n, n] K;
    matrix[total_samples, total_samples] kS = cov_matern32(S, l);

    for(i in 1:(n - 1)) {
      // Diagonal entries (apart from the last)
      K[i, i] = mean(block(kS, start_index[i], start_index[i], sample_lengths[i], sample_lengths[i]));
      for(j in (i + 1):n) {
        // Off-diagonal entries
        K[i, j] = mean(block(kS, start_index[i], start_index[j], sample_lengths[i], sample_lengths[j]));
        K[j, i] = K[i, j];
      }
    }
    K[n, n] = mean(block(kS, start_index[n], start_index[n], sample_lengths[n], sample_lengths[n]));

    return(K);
  }
}

data {
  int<lower=0> n; // Number of region
  vector[n] y; // Vector of observed responses
  vector[n] m; // Vector of sample sizes
  vector[n] mu; // Prior mean vector

  int sample_lengths[n]; // Number of Monte Carlo samples in each area
  int<lower=0> total_samples; // sum(sample_lengths)
  int start_index[n]; // Start indicies for each group of samples
  matrix[total_samples, total_samples] S; // Distances between all points (could be sparser!)
}

parameters {
  real beta_0; // Intercept
  vector[n] phi; // Spatial effects
  real<lower=0> sigma_phi; // Standard deviation of spatial effects
  real<lower=0> l; // Kernel lengthscale
}

transformed parameters {
  vector[n] eta = beta_0 + sigma_phi * phi;
}

model {
  matrix[n, n] K = cov_sample_average(S, l, n, start_index, sample_lengths, total_samples);
  // I could do this?
  // matrix[n, n] L = cholesky_decompose(K);
  // y ~ multi_normal_cholesky(mu, L);
  l ~ gamma(1, 1);
  sigma_phi ~ normal(0, 2.5); // Weakly informative prior
  beta_0 ~ normal(-2, 1);
  phi ~ multi_normal(mu, K);
  for(i in 1:n) {
   y[i] ~ xbinomial_logit(m[i], eta[i]);
  }
}

generated quantities {
  real tau_phi = 1 / sigma_phi^2; // Precision of spatial effects
  vector[n] rho = inv_logit(beta_0 + sigma_phi * phi);
  vector[n] log_lik;
  for (i in 1:n) {
    log_lik[i] = xbinomial_logit_lpdf(y[i] | m[i], eta[i]);
  }
}

Original computing environment

sessionInfo()

## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Monterey 12.6
## 
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tidyr_1.2.0          magrittr_2.0.3       dplyr_1.0.9          rmarkdown_2.14       bayesplot_1.9.0      rstan_2.21.5         ggplot2_3.3.6        StanHeaders_2.21.0-7
##  [9] sf_1.0-7             bsae_0.2.7          
## 
## loaded via a namespace (and not attached):
##  [1] viridis_0.6.2        sass_0.4.1           bit64_4.0.5          jsonlite_1.8.0       viridisLite_0.4.0    ids_1.0.1            bslib_0.3.1          RcppParallel_5.1.5  
##  [9] assertthat_0.2.1     distributional_0.3.0 posterior_1.2.2      askpass_1.1          highr_0.9            stats4_4.2.0         tensorA_0.36.2       blob_1.2.3          
## [17] yaml_2.3.5           backports_1.4.1      pillar_1.7.0         RSQLite_2.2.14       lattice_0.20-45      glue_1.6.2           RcppEigen_0.3.3.9.2  uuid_1.1-0          
## [25] digest_0.6.29        checkmate_2.1.0      colorspace_2.0-3     cowplot_1.1.1        htmltools_0.5.2      Matrix_1.5-1         plyr_1.8.7           pkgconfig_2.0.3     
## [33] purrr_0.3.4          scales_1.2.0         processx_3.5.3       tibble_3.1.7         openssl_2.0.1        proxy_0.4-26         generics_0.1.2       farver_2.1.0        
## [41] ellipsis_0.3.2       cachem_1.0.6         withr_2.5.0          cli_3.3.0            crayon_1.5.1         memoise_2.0.1        evaluate_0.15        ps_1.7.0            
## [49] fansi_1.0.3          class_7.3-20         pkgbuild_1.3.1       tools_4.2.0          loo_2.5.1            prettyunits_1.1.1    lifecycle_1.0.1      matrixStats_0.62.0  
## [57] stringr_1.4.0        munsell_0.5.0        callr_3.7.0          compiler_4.2.0       jquerylib_0.1.4      e1071_1.7-9          rlang_1.0.2          classInt_0.4-3      
## [65] units_0.8-0          grid_4.2.0           ggridges_0.5.3       rstudioapi_0.13      labeling_0.4.2       orderly_1.4.3        gtable_0.3.0         codetools_0.2-18    
## [73] abind_1.4-5          inline_0.3.19        DBI_1.1.2            reshape2_1.4.4       R6_2.5.1             gridExtra_2.3        knitr_1.39           bit_4.0.4           
## [81] fastmap_1.1.0        utf8_1.2.2           rprojroot_2.0.3      BH_1.78.0-0          KernSmooth_2.23-20   stringi_1.7.6        parallel_4.2.0       Rcpp_1.0.8.3        
## [89] vctrs_0.4.1          tidyselect_1.1.2     xfun_0.31