The method is described in Pichler & Hartig (2021) A new joint species distribution model for faster and more accurate inference of species associations from big community data, https://doi.org/10.1111/2041-210X.13687. The code for producing the results in this paper is available under the subfolder publications in this repo.
The method itself is wrapped into an R package, available under subfolder sjSDM. You can also use it stand-alone under Python (see instructions below). Note: for both the R and the python package, python >= 3.6 and pytorch must be installed (more details below).
Install the package via
install.packages("sjSDM")Depencies for the package can be installed before or after installing the package. Detailed explanations of the dependencies are provided in vignette(“Dependencies”, package = “sjSDM”), source code here. Very briefly, the dependencies can be automatically installed from within R:
sjSDM::install_sjSDM(version = "gpu") # or
sjSDM::install_sjSDM(version = "cpu")To cite sjSDM, please use the following citation:
citation("sjSDM")If you want to install the current (development) version from this repository, run
devtools::install_github("https://github.com/TheoreticalEcology/s-jSDM", subdir = "sjSDM", ref = "devel")Once the dependencies are installed, the following code should run:
Simulate a community and fit model:
library(sjSDM)## ── Attaching sjSDM ──────────────────────────────────────────────────── 1.0.2 ──
## ✔ torch <environment>
## ✔ torch_optimizer
## ✔ pyro
## ✔ madgrad
set.seed(42)
community <- simulate_SDM(sites = 100, species = 10, env = 3, se = TRUE)
Env <- community$env_weights
Occ <- community$response
SP <- matrix(rnorm(200, 0, 0.3), 100, 2) # spatial coordinates (no effect on species occurences)
model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1+X2+X3), spatial = linear(data = SP, formula = ~0+X1:X2), se = TRUE, family=binomial("probit"), sampling = 100L)
summary(model)## LogLik: -511.7972
## Regularization loss: 0
##
## Species-species correlation matrix:
##
## sp1 1.0000
## sp2 -0.3560 1.0000
## sp3 -0.1180 -0.3740 1.0000
## sp4 -0.1220 -0.3350 0.6950 1.0000
## sp5 0.5970 -0.3510 -0.0760 -0.0870 1.0000
## sp6 -0.2140 0.3880 0.1420 0.1350 -0.0730 1.0000
## sp7 0.4640 -0.1540 0.1470 0.1390 0.5120 0.2440 1.0000
## sp8 0.2560 0.1570 -0.4240 -0.4180 0.2150 -0.0320 0.1070 1.0000
## sp9 -0.0730 -0.0470 0.0540 0.0700 -0.3050 -0.2550 -0.2220 -0.1010 1.0000
## sp10 0.1800 0.3920 -0.6130 -0.5760 0.2220 0.1000 0.1090 0.4120 -0.2190 1.0000
##
##
##
## Spatial:
## sp1 sp2 sp3 sp4 sp5 sp6 sp7
## X1:X2 0.2260564 -0.5308802 0.5423007 0.03226206 0.3345731 0.1966482 0.4499481
## sp8 sp9 sp10
## X1:X2 0.4780776 0.2826627 0.2449929
##
##
##
## Estimate Std.Err Z value Pr(>|z|)
## sp1 (Intercept) -0.0290 0.2081 -0.14 0.88925
## sp1 X1 0.8329 0.3936 2.12 0.03433 *
## sp1 X2 -1.6153 0.3787 -4.27 2.0e-05 ***
## sp1 X3 -0.1006 0.3489 -0.29 0.77308
## sp2 (Intercept) -0.0335 0.2164 -0.16 0.87681
## sp2 X1 0.8809 0.4079 2.16 0.03080 *
## sp2 X2 0.1921 0.3994 0.48 0.63055
## sp2 X3 0.4257 0.3576 1.19 0.23383
## sp3 (Intercept) -0.3370 0.2195 -1.54 0.12472
## sp3 X1 1.0309 0.4044 2.55 0.01080 *
## sp3 X2 -0.2759 0.3924 -0.70 0.48190
## sp3 X3 -0.7200 0.3764 -1.91 0.05577 .
## sp4 (Intercept) -0.0535 0.1962 -0.27 0.78510
## sp4 X1 -1.0306 0.3724 -2.77 0.00565 **
## sp4 X2 -1.2974 0.3783 -3.43 0.00061 ***
## sp4 X3 -0.3160 0.3288 -0.96 0.33662
## sp5 (Intercept) -0.1503 0.2044 -0.74 0.46211
## sp5 X1 0.5085 0.3914 1.30 0.19383
## sp5 X2 0.4517 0.3673 1.23 0.21880
## sp5 X3 -0.4063 0.3406 -1.19 0.23291
## sp6 (Intercept) 0.1879 0.2058 0.91 0.36116
## sp6 X1 1.6793 0.4291 3.91 9.1e-05 ***
## sp6 X2 -0.6698 0.3940 -1.70 0.08912 .
## sp6 X3 0.1595 0.3369 0.47 0.63595
## sp7 (Intercept) -0.0188 0.2002 -0.09 0.92511
## sp7 X1 -0.2593 0.3889 -0.67 0.50501
## sp7 X2 0.3116 0.3679 0.85 0.39698
## sp7 X3 -1.0443 0.3462 -3.02 0.00255 **
## sp8 (Intercept) 0.1435 0.1533 0.94 0.34911
## sp8 X1 0.2595 0.2950 0.88 0.37895
## sp8 X2 0.3017 0.2885 1.05 0.29565
## sp8 X3 -1.0177 0.2674 -3.81 0.00014 ***
## sp9 (Intercept) 0.0251 0.1699 0.15 0.88257
## sp9 X1 1.0902 0.3260 3.34 0.00083 ***
## sp9 X2 -0.8162 0.3227 -2.53 0.01143 *
## sp9 X3 0.6126 0.2842 2.16 0.03111 *
## sp10 (Intercept) -0.0960 0.1874 -0.51 0.60858
## sp10 X1 -0.4447 0.3512 -1.27 0.20544
## sp10 X2 -0.9639 0.3397 -2.84 0.00454 **
## sp10 X3 -0.3804 0.3216 -1.18 0.23678
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(model)## LogLik: -511.7972
## Regularization loss: 0
##
## Species-species correlation matrix:
##
## sp1 1.0000
## sp2 -0.3560 1.0000
## sp3 -0.1180 -0.3740 1.0000
## sp4 -0.1220 -0.3350 0.6950 1.0000
## sp5 0.5970 -0.3510 -0.0760 -0.0870 1.0000
## sp6 -0.2140 0.3880 0.1420 0.1350 -0.0730 1.0000
## sp7 0.4640 -0.1540 0.1470 0.1390 0.5120 0.2440 1.0000
## sp8 0.2560 0.1570 -0.4240 -0.4180 0.2150 -0.0320 0.1070 1.0000
## sp9 -0.0730 -0.0470 0.0540 0.0700 -0.3050 -0.2550 -0.2220 -0.1010 1.0000
## sp10 0.1800 0.3920 -0.6130 -0.5760 0.2220 0.1000 0.1090 0.4120 -0.2190 1.0000
##
##
##
## Spatial:
## sp1 sp2 sp3 sp4 sp5 sp6 sp7
## X1:X2 0.2260564 -0.5308802 0.5423007 0.03226206 0.3345731 0.1966482 0.4499481
## sp8 sp9 sp10
## X1:X2 0.4780776 0.2826627 0.2449929
##
##
##
## Estimate Std.Err Z value Pr(>|z|)
## sp1 (Intercept) -0.0290 0.2081 -0.14 0.88925
## sp1 X1 0.8329 0.3936 2.12 0.03433 *
## sp1 X2 -1.6153 0.3787 -4.27 2.0e-05 ***
## sp1 X3 -0.1006 0.3489 -0.29 0.77308
## sp2 (Intercept) -0.0335 0.2164 -0.16 0.87681
## sp2 X1 0.8809 0.4079 2.16 0.03080 *
## sp2 X2 0.1921 0.3994 0.48 0.63055
## sp2 X3 0.4257 0.3576 1.19 0.23383
## sp3 (Intercept) -0.3370 0.2195 -1.54 0.12472
## sp3 X1 1.0309 0.4044 2.55 0.01080 *
## sp3 X2 -0.2759 0.3924 -0.70 0.48190
## sp3 X3 -0.7200 0.3764 -1.91 0.05577 .
## sp4 (Intercept) -0.0535 0.1962 -0.27 0.78510
## sp4 X1 -1.0306 0.3724 -2.77 0.00565 **
## sp4 X2 -1.2974 0.3783 -3.43 0.00061 ***
## sp4 X3 -0.3160 0.3288 -0.96 0.33662
## sp5 (Intercept) -0.1503 0.2044 -0.74 0.46211
## sp5 X1 0.5085 0.3914 1.30 0.19383
## sp5 X2 0.4517 0.3673 1.23 0.21880
## sp5 X3 -0.4063 0.3406 -1.19 0.23291
## sp6 (Intercept) 0.1879 0.2058 0.91 0.36116
## sp6 X1 1.6793 0.4291 3.91 9.1e-05 ***
## sp6 X2 -0.6698 0.3940 -1.70 0.08912 .
## sp6 X3 0.1595 0.3369 0.47 0.63595
## sp7 (Intercept) -0.0188 0.2002 -0.09 0.92511
## sp7 X1 -0.2593 0.3889 -0.67 0.50501
## sp7 X2 0.3116 0.3679 0.85 0.39698
## sp7 X3 -1.0443 0.3462 -3.02 0.00255 **
## sp8 (Intercept) 0.1435 0.1533 0.94 0.34911
## sp8 X1 0.2595 0.2950 0.88 0.37895
## sp8 X2 0.3017 0.2885 1.05 0.29565
## sp8 X3 -1.0177 0.2674 -3.81 0.00014 ***
## sp9 (Intercept) 0.0251 0.1699 0.15 0.88257
## sp9 X1 1.0902 0.3260 3.34 0.00083 ***
## sp9 X2 -0.8162 0.3227 -2.53 0.01143 *
## sp9 X3 0.6126 0.2842 2.16 0.03111 *
## sp10 (Intercept) -0.0960 0.1874 -0.51 0.60858
## sp10 X1 -0.4447 0.3512 -1.27 0.20544
## sp10 X2 -0.9639 0.3397 -2.84 0.00454 **
## sp10 X3 -0.3804 0.3216 -1.18 0.23678
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We also support also other response families: Count data:
model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1+X2+X3), spatial = linear(data = SP, formula = ~0+X1:X2), se = TRUE, family=poisson("log"))Gaussian (normal):
model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1+X2+X3), spatial = linear(data = SP, formula = ~0+X1:X2), se = TRUE, family=gaussian("identity"))Let’s have a look at the importance of the three groups (environment, associations, and space) on the occurences:
imp = importance(model)
print(imp)## sp env spatial biotic
## 1 1 0.8509671 3.718882e-04 0.14866112
## 2 2 0.5491779 3.699944e-03 0.44712225
## 3 3 0.7280060 3.747982e-03 0.26824600
## 4 4 0.9064414 7.932865e-06 0.09355063
## 5 5 0.4953244 2.304945e-03 0.50237066
## 6 6 0.8561825 2.861650e-04 0.14353135
## 7 7 0.5219943 1.883551e-03 0.47612211
## 8 8 0.9194061 4.193856e-03 0.07640003
## 9 9 0.7399929 6.984830e-04 0.25930864
## 10 10 0.9448155 1.054667e-03 0.05412992
plot(imp)
As expected, space has no effect on occurences.
Let’s have a look on community level how the three groups contribute to the overall explained variance
an = anova(model)
print(an)## Analysis of Deviance Table
##
## Terms added sequentially:
##
## Deviance Residual deviance R2 Nagelkerke R2 McFadden
## Biotic 150.78920 521.71684 0.77862 0.1091
## Abiotic 185.05100 336.66584 0.96521 0.2429
## Spatial 2.61524 334.05060 0.96611 0.2448
plot(an)
The anova shows the relative changes in the deviance of the groups and their intersections.
We can also visualize the individual contributions to the species and site
:
plot(an, internal=TRUE)## Registered S3 methods overwritten by 'ggtern':
## method from
## grid.draw.ggplot ggplot2
## plot.ggplot ggplot2
## print.ggplot ggplot2
If it fails, check out the help of ?install_sjSDM, ?installation_help, and vignette(“Dependencies”, package = “sjSDM”).
- Try install_sjSDM()
- New session, if no ‘PyTorch not found’ appears it should work, otherwise see ?installation_help
- If do not get the pkg to run, create an issue issue tracker or write an email to maximilian.pichler at ur.de
pip install sjSDM_pyPython example
import sjSDM_py as fa
import numpy as np
import torch
Env = np.random.randn(100, 5)
Occ = np.random.binomial(1, 0.5, [100, 10])
model = fa.Model_sjSDM(device=torch.device("cpu"), dtype=torch.float32)
model.add_env(5, 10)
model.build(5, optimizer=fa.optimizer_adamax(0.001),scheduler=False)
model.fit(Env, Occ, batch_size = 20, epochs = 10)
# print(model.weights)
# print(model.covariance)## Iter: 0/10 0%| | [00:00, ?it/s]Iter: 0/10 0%| | [00:00, ?it/s, loss=7.328]Iter: 1/10 10%|# | [00:00, 5.28it/s, loss=7.328]Iter: 1/10 10%|# | [00:00, 5.28it/s, loss=7.303]Iter: 1/10 10%|# | [00:00, 5.28it/s, loss=7.305]Iter: 1/10 10%|# | [00:00, 5.28it/s, loss=7.328]Iter: 1/10 10%|# | [00:00, 5.28it/s, loss=7.304]Iter: 1/10 10%|# | [00:00, 5.28it/s, loss=7.292]Iter: 6/10 60%|###### | [00:00, 23.82it/s, loss=7.292]Iter: 6/10 60%|###### | [00:00, 23.82it/s, loss=7.305]Iter: 6/10 60%|###### | [00:00, 23.82it/s, loss=7.285]Iter: 6/10 60%|###### | [00:00, 23.82it/s, loss=7.275]Iter: 6/10 60%|###### | [00:00, 23.82it/s, loss=7.261]Iter: 10/10 100%|##########| [00:00, 26.27it/s, loss=7.261]
Calculate Importance:
Beta = np.transpose(model.env_weights[0])
Sigma = ( model.sigma @ model.sigma.t() + torch.diag(torch.ones([1])) ).data.cpu().numpy()
covX = fa.covariance( torch.tensor(Env).t() ).data.cpu().numpy()
fa.importance(beta=Beta, covX=covX, sigma=Sigma)## {'env': array([[ 3.7089107e-03, 6.4757657e-03, 2.1139428e-03, 7.8506982e-03,
## 5.9689581e-03],
## [ 1.6079791e-02, 3.6896379e-03, 5.2381670e-03, 6.0370858e-03,
## 1.1554083e-03],
## [ 6.1437520e-03, 9.7632296e-03, 7.4961684e-03, -3.9039826e-04,
## 7.9328282e-04],
## [ 9.9422687e-05, 5.3525879e-04, 1.6309794e-03, 1.7054105e-02,
## 1.8254045e-02],
## [ 6.3455396e-04, 9.6351805e-04, 9.7591011e-04, 1.1000523e-04,
## 3.6621885e-03],
## [ 1.5674103e-02, 1.0940240e-03, 1.0806777e-02, 8.0627156e-03,
## 2.0850746e-02],
## [ 1.2024483e-05, 5.0663785e-03, 1.4675152e-02, 1.2784588e-02,
## 5.2863178e-03],
## [ 2.9168164e-03, 2.0089261e-04, 1.3217300e-02, 4.6121329e-04,
## 1.9982974e-05],
## [ 1.0323757e-02, 2.7552131e-03, 8.2153408e-04, 6.8914989e-04,
## 5.9855417e-03],
## [ 7.2189951e-03, 4.6700318e-03, 4.4443756e-03, 2.7086365e-03,
## 2.4578572e-05]], dtype=float32), 'biotic': array([0.9738818 , 0.96779996, 0.9761939 , 0.9624262 , 0.9936538 ,
## 0.94351166, 0.9621755 , 0.9831838 , 0.97942483, 0.9809334 ],
## dtype=float32)}
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent