Spatial Interdependence and Posterior Predictions in a Gaussian Process model

Hi all,

I’m working on a Gaussian Process model, using it to estimate interdependence between observations (military alliances). The challenge is that while most GP applications I can find use the interdependence as a control or examine the impact of various predictors, I want to examine how changing the outcome value in one unit impacts the outcome in all other units.

So far, the best I could come up with is to take some observed data, changing the outcome in one unit, and leaving the rest. Unfortunately, the posterior predictions are exactly the same, despite non-zero correlations between some units in the GP.

Here’s the code- can send data if anyone is interested in tinkering with it directly.

# load packages
library(brms)
library(tidyverse)

# run model 
cred_mod_noar  <- ordbetareg(irt_median ~ bilat + mean_kappa_atop +
                            mean_v2x_polyarchy + us_solsch + log_troops +
                            gp(mean_lat, mean_long,
                               k = 20,
                               cov = "exp_quad",
                               c = 5/4,
                               scale = TRUE,
                               gr = FALSE, iso = TRUE),
                          cores = 4,
                          refresh = 500,
                          backend = "cmdstanr",
                          control = list(adapt_delta = .99,
                                         max_treedepth = 20),
                       data = cred_data)
summary(cred_mod_noar)


# calculate impact of fall in Phillipine credibility (Duterte)
cred_data <- cred_data %>%
              group_by(cred_recip) %>%
              mutate(
                change_cred = irt_median - lag(irt_median)
              ) %>%
              ungroup()
max_fall <- min(cred_data$change_cred, na.rm = TRUE)
max_fall
filter(cred_data, change_cred == max_fall)

cred_orig <- filter(cred_data, year == 2015) 

# predictions with original data
pred_orig <- posterior_epred(cred_mod_noar, newdata = cred_orig)

# new predictions
cred_shift <- cred_orig
cred_shift$irt_median[cred_shift$cred_recip == "Philippines"] <- cred_shift$irt_median[cred_shift$cred_recip == "Philippines"] + max_fall  

pred_change <- posterior_epred(cred_mod_noar, newdata = cred_shift)


# Calculate the change in predictions
effect_change <- pred_change - pred_orig
summary(effect_change)

Should I be working with the covariance matrix for the GP more directly? Clarification if I have misunderstood the nature of how GPs model interdependence is very welcome also.

  • brms Version: 2.21.0

Thanks for your time and help!

I want to examine how changing the outcome value in one unit impacts the outcome in all other units.

That sounds like it might be a good application for the spatial lag model (‘lagsar’ in brms): Spatial simultaneous autoregressive (SAR) structures — sar • brms

If you use it and are unfamiliar with it, you’ll want to see the spatial econometrics literature on interpreting the coefficients (like Introduction to Spatial Econometrics)

Thanks for the quick response. I’ve actually fit some spatial models already- the calculations there are straightforward. I’m therefore wondering if or how to take an analogous approach with GP models.

A quick followup with what is hopefully a quick question for @paul.buerkner, to make sure this is clear.

I can see from this and the brms documentation that posterior predictions (posterior_epred and related) include autocorrelation. Is that also true of GP terms?

Apologies if I’ve missed something here. It seems like the spatial correlations are not showing up in the predictions, but I might also be missing something conceptually.

I don’t fully understand your question. You mean if GPs contain autocorrelation or if posterior predictions include GPs? The answer to the latter question is “yes”. The first question would not quite make sense to me.

Second question is it- thanks!