Modelling species responses with HMSC
We applied HMSC from the Hmsc R package (Tikhonov et al. 2020) to
fit a joint species distribution model to plant community data
simultaneously including information on traits (genome size),
environmental covariates, and phylogenetic relationships in a single
model (Pollock et al. 2012; Warton et al. 2015). We
included MAP, MAT, altitude, soil N, and soil P as fixed effects (i.e.,
predictors) and used the sampling site as the sampling unit and the
random effect. HMSC estimates species responses (slope parameters) to
environmental covariates and uses these responses as the species’
functional niche to estimate the strength and direction in which these
niches moderate multiple species responses to environmental filtering.
To account for the phylogenetic correlations among all the species, we
included a phylogenetic correlation matrix (construction details of the
phylogenetic trees see above) in the model’s covariance structure. HMSC
model can indicate the strength of the phylogenetic signal based on
parameter ρ (varying from 0 to 1, a higher value indicating stronger
phylogenetic signal) (Ovaskainen et al. 2017).
We fitted the model to the plant community (with a probit-link for the
presence/absence data) with Bayesian inference, using the posterior
sampling scheme (Ovaskainen et al. 2016). We used four Markov
Chain Monte Carlo (MCMC) chains, each of which consisted of 15,000
iterations, out of which we discarded the first 5,000 as the burn-in and
thinned the remaining by 10 to yield in total 1000 posterior samples per
chain. We assessed the convergence of the MCMC chains by examining the
distribution of the potential scale reduction factor over the parameters
that measure the responses of the species to the fixed effects included
in the model. The MCMC convergence of the HMSC models was satisfactory
as the potential scale reduction factors (psrf) are close to one (the
psrf of the 𝛽-parameters that measure species response to environmental
covariates (Ovaskainen et al. 2017) were on average 1.16). We
examined the explanatory and predictive powers of the HMSC models
through species-specific Tjur’s R2 values (Tjur 2009).
Twofold cross-validation was performed to assess the predictive power of
the model.
Based on the output of HMSC, we first characterized species’ responses
to each of the environmental variables as a mean parameter and its 95%
posterior probability. Second, we predicted the community weighted mean
genome size based on 161 species considered in the model as responses to
the included environmental covariates. Note that genome size was
log-transformed. Third, we used HMSC to capture the species-to-species
associations (positive or negative co-occurrences) with latent
variables. We considered statistical support for these associations
based on 90% posterior probabilities. Association patterns were
classified as positive or negative, with this result used to estimate
the frequency of two taxa co-existing (or not) compared to the frequency
expected based on shared habitats or random occurrences.