Computing approximate Gaussian processes

[katiec]

Hi! I am very new to Bayesian Statistics, so please bear with me…

I am testing the effects of extreme climate events on shorebird productivity (number of chicks). I have 33 years of data (1992-2025, missing 1998 data) across 11 islands. In my baseline models, I’ve included survey date (doy) to control for survey-related effects. I have also included year and terrID as random effects to account for repeated measures.

The model below converges without any warning messages:

pr_ci <- brm(nchicks ~ z_npMHWci + z_npENSOci + z_AHWci + z_PPTci + # climate
               z_doy + I(z_doy^2) + (1|year) + (1|terrID), # base
             data = pr_df,
             family = negbinomial(link = "log"),
             iter = 3000) 

Unsurprisingly, the model residuals have temporal autocorrelation with a one year lag. My colleague suggested I add a Hilbert space Gaussian process approximation with an exponentiated quadratic covariance function following (Riutort-Mayol et al., 2023) to account for this autocorrelation.

I have a few questions about the gp() function in brms…

1. Does the predictor for the GP (in my case, time) need to be centered and/or standardized? My raw time variable is on a continuous scale (1-33). All the continuous predictors in my model are standardized (mean = 0, sd = 1).
2. Is it redundant and/or problematic to include a gp(time) term and (1|year) term in the same model?
3. Adding gp(raw_time, k=20, c=5/4) to the above model results in ~30 divergent iteration warnings… I wasn’t running into these warnings before so I am assuming the issue is my gp() term. Setting k=20 and c=5/4 was recommended here but after attempting to follow the steps outlined in Riutort-Mayol et al., 2023, I don’t think these values are right. But I also don’t trust the values I calculated in Phase A (k=6, c=1.6… doesn’t make sense).

Any help would be greatly appreciated!!!

Thank you :)

[mhollanders]

Hi, I’ll have a crack. Firstly, DOY and survey date are technically different. DOY can be the same between years (e.g., 18 June) but survey date can’t be (18 June 2025 vs 18 June 2026). If you use DOY, that particular date in each year gets the same parameter value. If you use date, each unique date gets its own parameter.

1. Not necessarily, but it does help with setting the length scale prior. The length scale essentially describes the “decay” in correlation on the unit of your input (squared).
2. No, not necessarily. The (1|year)—you could even do gp(year), as a gp() term is just a random effect that includes a marginal scale + length scale for the correlations—captures tendencies of each particular year. gp(time) captures variation within each year.
3. Not sure about this one, are exact GPs on the table, or is your dataset large?

[jsocolar]

@mhollanders

2. No, not necessarily. The (1|year)—you could even do gp(year), as a gp() term is just a random effect that includes a marginal scale + length scale for the correlations—captures tendencies of each particular year. gp(time) captures variation within each year.

I think the OP is already talking about a gp over years, since they say their time variable takes values 1-33 and they have 33 years of data.

@katiec

2. Is it redundant and/or problematic to include a gp(time) term and (1|year) term in the same model?

In general it is not problematic and can often be a good idea. To see that this isn’t problematic, consider a very simple model like y ~ gp(x). In this model, the residual plays the same role that (1 | year) would play in your model.

Two general comments:

1. Your description of an autocorrelative structure with a 1-year lag can be modeled exactly in brms via an autoregressive process, if desired. See brms::ar. If I were you, I would go ahead and try this, given that it matches your description of the generative process and potentially might not display the same computational problems.
2. Are the fluctuations that you observe consistent across all islands, or do you see different islands fluctuating independently of one another? In the latter case, it might be appropriate to fit your gaussian process (or autoregressive process) within islands (i.e. a separate gp for each island). You can achieve this with the by argument to gp or the gr argument to ar.

[katiec]

@mhollanders We used DOY instead of the actual survey date because we wanted to control for when we surveyed the territories within the breeding season as this likely (and does) affects the response variable (e.g., a territory surveyed early in the season (e.g., May 15) is less likely to have chicks than a territory surveyed later in the season (e.g., July 15), regardless of the year these surveys were conducted). Including year as a random effect then accounts for any unexplained year-to-year variation. Exact GPs are possible but I was still running into issues.

@jsocolar is it okay to use ar() even though I am missing all data for 1998 (i.e., time variable jumps from 5 to 7)? I also mistyped and meant my time variable is 0-33, with 0 referring to 1992 and 33 referring to 2025.

I tried ar, but the model is still unhappy (low ESS, crazy high estimates for shape)

pr_gp <- brm(nchicks ~ z_npMHWci + z_npENSOci + z_AHWci + z_PPTci + # climate
               ar(p = 1) + 
               z_doy + I(z_doy^2) + (1|year) + (1|terrID), # base
             data = pr_df,
             family = negbinomial(link = "log"),
             iter = 3000)
# Warning messages:
#   1: There were 24 divergent transitions after warmup. See
# https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
# to find out why this is a problem and how to eliminate them. 
# 2: Examine the pairs() plot to diagnose sampling problems
# 
# 3: 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 

> summary(pr_gp)
 Family: negbinomial 
  Links: mu = log 
Formula: nchicks ~ z_npMHWci + z_npENSOci + z_AHWci + z_PPTci + ar(p = 1) + z_doy + I(z_doy^2) + (1 | year) + (1 | terrID) 
   Data: pr_df (Number of observations: 1185) 
  Draws: 4 chains, each with iter = 3000; warmup = 1500; thin = 1;
         total post-warmup draws = 6000

Correlation Structures:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
ar[1]     0.46      0.42    -0.71     0.98 1.00      494      713
sderr     0.38      0.20     0.02     0.71 1.01      180      793

Multilevel Hyperparameters:
~terrID (Number of levels: 106) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.21      0.11     0.02     0.41 1.00      866     1196

~year (Number of levels: 33) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.33      0.10     0.14     0.54 1.00     1465     1636

Regression Coefficients:
           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     -0.85      0.13    -1.11    -0.61 1.00      316     1371
z_npMHWci     -0.04      0.09    -0.22     0.15 1.00     3279     3488
z_npENSOci     0.15      0.08    -0.00     0.32 1.00     3174     3664
z_AHWci        0.07      0.09    -0.12     0.24 1.00     2742     2893
z_PPTci        0.00      0.09    -0.18     0.18 1.00     2976     3729
z_doy         -0.09      0.08    -0.25     0.07 1.00     3332     4271
Iz_doyE2      -0.14      0.05    -0.23    -0.06 1.00     4575     4473

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape  2752.36  86433.46     1.25  3177.69 1.01      211     1229

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Warning message:
There were 24 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup 

I am not sure how to check if the temporal autocorrelation is consistent across all islands. Some islands only had one territory and/or were not checked every year. Below is how I tested for temporal autocorrelation by time.

pr_ci <- brm(nchicks ~ z_npMHWci + z_npENSOci + z_AHWci + z_PPTci + # climate
               z_doy + I(z_doy^2) + (1|year) + (1|terrID), # base
             data = pr_df,
             family = negbinomial(link = "log"),
             iter = 3000)

res <- createDHARMa(
  simulatedResponse = t(posterior_predict(pr_ci)),
  observedResponse = pr_df$nchicks,
  fittedPredictedResponse = apply(t(posterior_epred(pr_ci)), 1, mean),
  integerResponse = TRUE)
plot(res) 

res_agg <- recalculateResiduals(res, group = pr_df$t)
testTemporalAutocorrelation(res_agg, time = unique(pr_df$t))

Thank you both so much for replying and providing guidance!

[avehtari]

@katiec

1. Does the predictor for the GP (in my case, time) need to be centered and/or standardized? My raw time variable is on a continuous scale (1-33). All the continuous predictors in my model are standardized (mean = 0, sd = 1).

brms centers automatically, but doesn’t standardize. Now I can’t remember how the box is set when using brms gp(..., k=...), but you can easily test whether the result changes if you do standardize time given to gp(..., k=...)

@katiec

2. Is it redundant and/or problematic to include a gp(time) term and (1|year) term in the same model?

This will make the posterior shape challenging. You could improve the posterior shape by adding sum-to-zero-constraint, but that is not currently possible in brms

@katiec

3. the above model results in ~30 divergent iteration warnings

This is quite a small number and may be solved by increasing adapt_delta a little bit to something like 0.9 or 0.95. If that doesn’t solve it, you can also test with posterior SBC whether that small number of divergences affect inference for the quantity of interest

[katiec]

I run into fewer issues when I use the raw time variables, but I am still not sure what the best values of k and c are despite reading Riutort-Mayol et al., 2023. Step A1 confuses me.

Increasing adapt_delta produces more warnings, so I will look into the posterior SBC. Thanks!

[avehtari]

For example, if your centered time would be from -10 to 10, then S would be 10. If you don’t have any initial guess for the lengthscale, you can start with lengthscale guess which makes l/S \in [0.5,1], which would be in this happen with l \in [5,10]. Using that guess, you can get first guess for k and c. After first MCMC run, you can then update k and c if needed