I have a dataset with populations and several binary features for each population. I also have latitude and longitude data for the populations.

Imagine the data looks as follows:

```
A B C lat long
1 0 0 12 0
1 1 0 11 0.5
0 0 1 15 5
```

I know that the presence or absence of some of these features are correlated. Most importantly, I know that there is some borrowing between populations of the presence or absence of features . I also have some additional controls. I am not interested in the correlations between A, B, C, etc. What I am interested in is estimating the amount of contact between the populations based on the presence or absence of all the features.

If I was only looking at one of the features, in BRMS I could do something like :

```
brm(bf(A ~ 1 + B + C + gp(lat, long)),
family = bernoulli,
...)
```

To estimate all the model combinations, I could do:

```
brm( bf(A ~ 1 + B + C + gp(latitude, longitude)) +
bf(B ~ 1 + A + C + gp(latitude, longitude)) +
bf(C ~ 1 + B + A + gp(latitude, longitude)),
family = bernoulli,
...)
```

But as I understand it, this estimates independent GPs for each model. From what Iâ€™ve read, it is not possible using BRMS to share the same predictor for multiple models.

Looking at the STAN code, I think that changing:

```
...
vector[N_A] gp_pred_A_1 = gp(Xgp_A_1, sdgp_A_1[1], lscale_hasb_1[1], zgp_A_1);
vector[N_B] gp_pred_B_1 = gp(Xgp_B_1, sdgp_B_1[1], lscale_hasb_1[1], zgp_B_1);
vector[N_C] gp_pred_C_1 = gp(Xgp_C_1, sdgp_C_1[1], lscale_hasb_1[1], zgp_C_1);
vector[N_A] mu_hasb = Intercept_A + rep_vector(0.0, N_A) + gp_pred_A_1;
vector[N_B] mu_hasd = Intercept_B + rep_vector(0.0, N_B) + gp_pred_B_1;
vector[N_C] mu_hasf = Intercept_C + rep_vector(0.0, N_C) + gp_pred_C_1;
...
```

To:

```
...
vector[N_A] gp_pred_A_1 = gp(Xgp_A_1, sdgp_A_1[1], lscale_hasb_1[1], zgp_A_1);
vector[N_A] mu_hasb = Intercept_A + rep_vector(0.0, N_A) + gp_pred_A_1;
vector[N_B] mu_hasd = Intercept_B + rep_vector(0.0, N_B) + gp_pred_A_1;
vector[N_C] mu_hasf = Intercept_C + rep_vector(0.0, N_C) + gp_pred_A_1;
...
```

(Xgp_A_1, Xgp_B_1, etc. are the same for all three models)

will share the same GP across all three models. So, my first question would be: **is this a valid model?**

Would it maybe also make sense to fix the scale to be the same for all three models, but allow the *sdgp*â€™s to vary by model? does that make any sense? Intuitively I would like to also specify that a model which allows different GPs by model, but where they inform each other and the *lscale* and *sdgp* cannot vary freely.

The intuition I have is that while between population x1 and x2 there is only a fixed amount of contact, feature A may be more easily borrowed than feature B, for example, so a single GP might be too strong of a constraint.

Final bonus question: is it possible to specify the GP to only apply to the presence of the response variable? So to state that the presence of a feature is spatially correlated but its absence is not.