Many thanks! using many plot data, I tried to study how similarity of pairwise species (obs.bz) correlated with plot-level species richness (scal.logrich), plot topographic complexity (scal.elevdiff), species abundance (scal.abd and scal.abdpair). ‘scal.’ means standardizing the variable, ‘log’ means log-transforming the relevant variables.

The response ‘obs.bz’ is subject to [0,1) with many values of zero, so I choose zero-inflated beta model to fit the data.

The brm codes I used are:

```
lfm <- bf(obs.bz ~ poly(scal.abd,3)+poly(scal.abdpair,3)+scal.logrich+scal.elevdiff+poly(scal.abd,3):scal.logrich+
poly(scal.abd,2):poly(scal.abdpair,2)+(poly(scal.abd,3)+poly(scal.abdpair,3)||plotcode),
zi ~ poly(scal.abd,1)+poly(scal.abdpair,1)+scal.logrich+scal.elevdiff+poly(scal.abd,1):scal.logrich+poly(scal.abdpair,1):scal.logrich+
poly(scal.abd,1):poly(scal.abdpair,1)+(poly(scal.abd,1)+poly(scal.abdpair,1)||plotcode),
phi ~ poly(scal.abd,1)+poly(scal.abdpair,1)+scal.logrich+scal.elevdiff+poly(scal.abd,3):scal.logrich+poly(scal.abdpair,1):scal.logrich+
poly(scal.abd,2):poly(scal.abdpair,2)+(poly(scal.abd,1)+poly(scal.abdpair,1)||plotcode))
lprior <- c(prior(normal(0,100),class="b"),
prior(normal(0,100),class="Intercept"),
prior(normal(0,100),class="sd"))
bayesfit_zb_pair <- brm(formula=lfm, prior=lprior, data=adat, control = list(adapt_delta = 0.95),
family=zero_inflated_beta(link="probit"),cores=no_cores)
```

the result:

I guess using GAM with brm could help. look forward to commens. Thanks.