Repository: TheoreticalEcology/s-jSDM (language: python)

Issue:
Readme Files
Hi Max,

at the moment, there are 2 Readme.md files, one in the top level, one in the package, and they are copies. Do we really need both? I'm not sure if the README.md in the package folder is really used for anything?

F
importance
in the vignette it says importance is deprecated, then I would also erase it from the vignette.

Nevertheless, let's maybe have a quick chat about this, as it seems it has different code than anova, so it's not clear to me if it is the same and just deprecated, or if it is not, then why not keep it and document it properly?


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 97cae9c..642a7e5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,3 +33,5 @@ sjSDM.Rproj
 sjSDM/vignettes/sjSDM_Introduction.R
 .DS_Store
 .DS_Store
+.DS_Store
+.DS_Store
diff --git a/Code/Theory/VP_Anova.R b/Code/Theory/VP_Anova.R
new file mode 100644
index 0000000..0328cd5
--- /dev/null
+++ b/Code/Theory/VP_Anova.R
@@ -0,0 +1,196 @@
+n = 50
+X1 = runif(n)
+X2 = X1 + rnorm(n, sd = 0.2)
+X3 = runif(n)
+cor(X1, X2)
+
+Y = X1 + X3+ rnorm(n, sd = 1.0)
+df = data.frame(Y = Y, X1= X1, X2 = X2, X3 = X3)
+
+m = lm(Y~X1+X2+X3, data = df)
+summary(m)
+car::Anova(m)
+anova(m)
+df = data.frame(Y = Y, X1= X1, X2 = X2, X3 = X3)
+
+imp = function(df) {
+  model =  lm(Y~X1+X2+X3, data= df) 
+  res = rep(0, 3)
+  for(i in 1:3){ 
+    model_tmp = model
+    model_tmp$coefficients[2:4] = model_tmp$coefficients[2:4] * diag(1.0, 3)[i,]
+    pred = predict(model_tmp, newdata = df)
+    res[i] = 1-sum((df$Y - pred  )**2)/sum((df$Y)**2)
+  }
+  return(res)
+}
+
+imp(df)
+
+
+# Otso
+beta = coef(m)
+beta
+cM = cov(cbind(0, df[,-1]))
+
+Rs = 
+sapply(2:4, function(i) {
+t(beta[c(1, i)] ) %*% cM[c(1, i), c(1,i)] %*% beta[c(1, i)] 
+})
+Rs / sum(Rs)
+
+
+
+Rs2 = Rs / summary(m)$coef[-1,2]
+Rs2 / sum(Rs2)
+
+
+relaimpo::calc.relimp(m, type = c("lmg", "last", "first", "betasq", "pratt"), rela = TRUE )
+
+car::Anova(m)
+
+
+
+X1 = runif(1000, -1, 1)
+X2 = X1 + rnorm(1000, sd = 0.3)
+X3 = runif(1000, -1 ,1 )
+cor(X1, X2)
+X1 = scale(X1)
+X2 = scale(X2)
+X3 = scale(X3)
+prob = plogis(5*X1 + 5*X3)
+Y = rbinom(1000, 1, prob)
+
+m = glm(Y~X1+X2+X3, family = binomial(), data = df)
+rsq::rsq.partial(m, type = "n")
+car::Anova(m)
+anova(m)
+summary(m)
+df = data.frame(Y = Y, X1= X1, X2 = X2, X3 = X3)
+
+imp_glm = function(df) {
+  model =  glm(Y~X1+X2+X3, family = binomial(), data = df)
+  null = glm(Y~1, family = binomial(), data = df)
+  pred_null = predict(null, newdata = df, type = "response")
+  f_L = sum(dbinom(df$Y, 1, pred_null, TRUE))
+  
+  res = rep(0, 3)
+  for(i in 1:3){ 
+    model_tmp = model
+    model_tmp$coefficients[2:4] = model$coefficients[2:4] * diag(1.0, 3)[i,]
+    model_tmp$coefficients[1] = 0
+    pred = predict(model_tmp, newdata = df, type = "response")
+    f_m = sum(dbinom(df$Y, 1, pred, TRUE))
+    res[i] = 1 -  f_m / f_L # pseudo Rsquared (McFadden)
+  }
+  res = c(res, 1- sum(dbinom(df$Y, 1, predict(model, newdata = df, type = "response"), TRUE))/f_L)
+  
+  return(res)
+}
+(imp_glm(df))
+
+
+model =  glm(Y~X3, family = binomial(), data = df)
+model
+
+library(lme4)
+df$site = as.factor(sample.int(10, size = 1000, replace = TRUE))
+mm = glmer(Y~X1+X2+X3+(1|site), data = df, family = binomial())
+rr = partR2::partR2(mm, partvars = c("X1", "X2", "X3"), R2_type="conditional")
+rr$BW
+
+library(dominanceanalysis)
+dominanceAnalysis(m)
+
+
+YY = binomial()$linkfun(scales::rescale(as.numeric(df$Y), to = c(0.001, 1-0.001) ))
+coef(m)
+
+imp2 = function(df) {
+  model =  glm(Y~X1+X2+X3, data= df, family = binomial()) 
+  res = rep(0, 3)
+  for(i in 1:3){ 
+    model_tmp = model
+    model_tmp$coefficients[2:4] = model_tmp$coefficients[2:4] * diag(1.0, 3)[i,]
+    pred = predict(model_tmp, newdata = df)
+    res[i] = 1-sum((YY - pred  )**2)/sum((YY)**2)
+  }
+  return(res)
+}
+imp2(df)
+
+library(iml)
+Pred = Predictor$new(model = m, data = df[,-c(1, 5)])
+Pred$task = "classif"
+res = sapply(1:1000, function(i) iml::Shapley$new(Pred, x.interest = df[i,2:4])$results[,2])
+
+
+imp_glm2 = function(df) {
+  model =  glm(Y~X1+X2+X3, family = binomial(), data = df)
+  null = glm(Y~1, family = binomial(), data = df)
+  pred_null = predict(null, newdata = df, type = "response")
+  f_L = sum(dbinom(df$Y, 1, pred_null, TRUE))
+  
+  res = rep(0, 3)
+  for(i in 1:3){ 
+    df_tmp = df
+    df_tmp[,i+1] = 0
+    pred = predict(model, newdata = df_tmp, type = "response")
+    f_m = sum(dbinom(df$Y, 1, pred, TRUE))
+    res[i] = 1 -  f_m / f_L # pseudo Rsquared (McFadden)
+  }
+  res = c(res, 1- sum(dbinom(df$Y, 1, predict(model, newdata = df, type = "response"), TRUE))/f_L)
+  res[1:3] = res[4] - res[1:3]
+  
+  return(res)
+}
+(imp_glm2(df))
+
+
+library(Hmsc)
+
+n = 200
+X1 = runif(n, -1, 1)
+X2 = 0*X1 + rnorm(n, sd = 0.3)
+X3 = runif(n, -1 ,1 )
+cor(X1, X2)
+X1 = scale(X1)
+X2 = scale(X2)
+X3 = scale(X3)
+prob = plogis(-2*X1 + 2*X3)
+Y = rbinom(n, 1, prob)
+df = data.frame(Y = Y, X1= X1, X2 = X2, X3 = X3)
+m = glm(Y~X1+X2+X3, family = binomial(), data = df)
+rsq::rsq.partial(m, type = "n")
+car::Anova(glm(Y~X1+X2+X3, family = binomial(), data = df), error.estimate=c("deviance"), test.statistic = "F")
+anova(m)
+summary(m)
+df = data.frame(Y = Y, X1= X1, X2 = X2, X3 = X3)
+
+
+Y = as.matrix(df$Y)
+XData = df[,-1]
+m = Hmsc(Y = Y, XData = XData, XFormula = ~X1+X2+X3, distr = "normal")
+nChains = 2
+test.run = FALSE
+if (test.run){
+  #with this option, the vignette runs fast but results are not reliable
+  thin = 1
+  samples = 10
+  transient = 5
+  verbose = 5
+} else {
+  #with this option, the vignette evaluates slow but it reproduces the results of the
+  #.pdf version
+  thin = 5
+  samples = 1000
+  transient = 500*thin
+  verbose = 500*thin
+}
+
+m = sampleMcmc(m, thin = thin, samples = samples, transient = transient,
+               nChains = nChains, verbose = verbose)
+Hmsc::computeVariancePartitioning(m)$vals
+
+
+
diff --git a/README.Rmd b/README.Rmd
index f1d3b0a..cdb84d6 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -11,17 +11,25 @@ knitr::opts_chunk$set(echo = TRUE)
 
 # s-jSDM - Fast and accurate Joint Species Distribution Modeling
 
-## About the method
+## About sjSDM
 
-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://www.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 sjSDM package is an R package for estimating so-called joint species distribution models. A jSDM is a GLMM that models a multivariate (i.e. a many-species) response to the environment, space and a covariance term that models conditional (on the other terms) correlations between the outputs (i.e. species). 
 
-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.7 and pytorch must be installed (more details below).
+![image](sjSDM/vignettes/jSDM-structure.png)
 
-## Installing the R / Python package
+A big challenge in jSDM implementation is computational speed. The goal of the sjSDM (which stands for "scalable joint species distribution models") is to make jSDM computations fast and scalable. Unlike many other packages, which use a latent-variable approximation to make estimating jSDMs faster, sjSDM fits a full covariance matrix in the likelihood, which is, however, numerically approximated via simulations. 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://www.doi.org/10.1111/2041-210X.13687. 
 
-### R-package
+The core code of sjSDM is implemented in Python / PyTorch, which is then wrapped into an R package. In principle, you can also use it stand-alone under Python (see instructions below). Note: for both the R and the python package, python \>= 3.7 and pytorch must be installed (more details below). However, for most users, it will be more convenient to use sjSDM via the sjSDM R package, which also provides a large number of downstream functionalities. 
 
-Install the package via
+To get citation info for sjSDM when you use it for your reseach, type  
+
+```{r,eval=FALSE}
+citation("sjSDM")
+```
+
+## Installing the R package
+
+sjSDM is distributed via [CRAN](https://cran.rstudio.com/web/packages/sjSDM/index.html). For most users, it will be best to install the package from CRAN
 
 ```{r,eval=FALSE}
 install.packages("sjSDM")
@@ -34,62 +42,64 @@ sjSDM::install_sjSDM(version = "gpu") # or
 sjSDM::install_sjSDM(version = "cpu")
 ```
 
-To cite sjSDM, please use the following citation:
+For advanced users: if you want to install the current (development) version from this repository, run
 
 ```{r,eval=FALSE}
-citation("sjSDM")
+devtools::install_github("https://github.com/TheoreticalEcology/s-jSDM", subdir = "sjSDM", ref = "master")
 ```
 
-### Development
+dependencies should be installed as above. If the installation fails, check out the help of ?install_sjSDM, ?installation_help, and vignette("Dependencies", package = "sjSDM").
 
-If you want to install the current (development) version from this repository, run
+1.  Try install_sjSDM()
+2.  New session, if no 'PyTorch not found' appears it should work, otherwise see ?installation_help
+3.  If do not get the pkg to run, create an issue [issue tracker](https://github.com/TheoreticalEcology/s-jSDM/issues) or write an email to maximilian.pichler at ur.de
 
-```{r,eval=FALSE}
-devtools::install_github("https://github.com/TheoreticalEcology/s-jSDM", subdir = "sjSDM", ref = "master")
-```
+## Basic Workflow
 
-Once the dependencies are installed, the following code should run:
+Load the package
 
-## Workflow
+```{r, message = F}
+library(sjSDM)
+```
 
-Simulate a community and fit a sjSDM model:
+Simulate some community data 
 
 ```{r}
-library(sjSDM)
 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)
+```
+
+This fits the standard SDM with environmental, spatial and covariance terms 
 
+```{r}
 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)
-plot(model)
 ```
 
-We support other distributions:
-
--   Count data with Poisson:
+We can also plot the niche estimates, i.e the estimates in the environmental component
 
-```{r,eval=FALSE}
-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"))
+```{r, echo=F}
+plot(model)
 ```
 
--   Count data with negative Binomial (which is still experimental, if you run into errors/problems, please let us know):
+We can also visualize the species-species association matrix
 
-```{r,eval=FALSE}
-model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1+X2+X3), spatial = linear(data = SP, formula = ~0+X1:X2), se = TR, family="nbinom")
+```{r}
+image(getCor(model))
 ```
 
--   Gaussian (normal):
 
-```{r,eval=FALSE}
-model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1+X2+X3), spatial = linear(data = SP, formula = ~0+X1:X2), se = TR, family=gaussian("identity"))
-```
+## Anova / Variation partitioning
 
-### Anova
+### Global ANOVA
 
-ANOVA can be used to partition the three components (abiotic, biotic, and spatial):
+As in other models, it can be interesting to analyze how much variation is explained by which parts of hte model. 
+
+![image](sjSDM/vignettes/jSDM-ANOVA.png)
+For the Env, Spatial, Covariance terms, this is implemented in 
 
 ```{r,fig.height=7, fig.width=6.3}
 an = anova(model)
@@ -113,15 +123,9 @@ The plot function returns the results for the internal metacommunity structure:
 print(results$data$Species)
 ```
 
-## Installation trouble shooting
-
-If the installation fails, check out the help of ?install_sjSDM, ?installation_help, and vignette("Dependencies", package = "sjSDM").
-
-1.  Try install_sjSDM()
-2.  New session, if no 'PyTorch not found' appears it should work, otherwise see ?installation_help
-3.  If do not get the pkg to run, create an issue [issue tracker](https://github.com/TheoreticalEcology/s-jSDM/issues) or write an email to maximilian.pichler at ur.de
+## Python Package
 
-### Python Package
+If you want to use sjSDM from python (as said, not encouraged because all help and downstream functions are in R), install via 
 
 ```{bash,eval=FALSE}
 pip install sjSDM_py
diff --git a/README.md b/README.md
index dc68ca5..aa6e022 100644
--- a/README.md
+++ b/README.md
@@ -3,34 +3,58 @@
 
 [![Project Status: Active – The project has reached a stable, usable
 state and is being actively
-developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
+developed.](http://www.repostatus.org/badges/latest/active.svg)](http://www.repostatus.org/#active)
 [![License: GPL
 v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
-[![CRAN_Status_Badge](https://www.r-pkg.org/badges/version/sjSDM)](https://cran.r-project.org/package=sjSDM)
+[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/sjSDM)](https://cran.r-project.org/package=sjSDM)
 ![R-CMD-check](https://github.com/TheoreticalEcology/s-jSDM/workflows/R-CMD-check/badge.svg?branch=master)
 [![Publication](https://img.shields.io/badge/Publication-10.1111/2041-green.svg)](https://www.doi.org/10.1111/2041-210X.13687)
 
 # s-jSDM - Fast and accurate Joint Species Distribution Modeling
 
-## About the method
-
+## About sjSDM
+
+The sjSDM package is an R package for estimating so-called joint species
+distribution models. A jSDM is a GLMM that models a multivariate (i.e. a
+many-species) response to the environment, space and a covariance term
+that models conditional (on the other terms) correlations between the
+outputs (i.e. species).
+
+<figure>
+<img src="sjSDM/vignettes/jSDM-structure.png" alt="image" />
+<figcaption aria-hidden="true">image</figcaption>
+</figure>
+
+A big challenge in jSDM implementation is computational speed. The goal
+of the sjSDM (which stands for “scalable joint species distribution
+models”) is to make jSDM computations fast and scalable. Unlike many
+other packages, which use a latent-variable approximation to make
+estimating jSDMs faster, sjSDM fits a full covariance matrix in the
+likelihood, which is, however, numerically approximated via simulations.
 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.
+<https://www.doi.org/10.1111/2041-210X.13687>.
+
+The core code of sjSDM is implemented in Python / PyTorch, which is then
+wrapped into an R package. In principle, you can also use it stand-alone
+under Python (see instructions below). Note: for both the R and the
+python package, python \>= 3.7 and pytorch must be installed (more
+details below). However, for most users, it will be more convenient to
+use sjSDM via the sjSDM R package, which also provides a large number of
+downstream functionalities.
 
-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.7 and pytorch must be installed (more details below).
+To get citation info for sjSDM when you use it for your reseach, type
 
-## Installing the R / Python package
+``` r
+citation("sjSDM")
+```
 
-### R-package
+## Installing the R package
 
-Install the package via
+sjSDM is distributed via
+[CRAN](https://cran.rstudio.com/web/packages/sjSDM/index.html). For most
+users, it will be best to install the package from CRAN
 
 ``` r
 install.packages("sjSDM")
@@ -48,37 +72,33 @@ sjSDM::install_sjSDM(version = "gpu") # or
 sjSDM::install_sjSDM(version = "cpu")
 ```
 
-To cite sjSDM, please use the following citation:
-
-``` r
-citation("sjSDM")
-```
-
-### Development
-
-If you want to install the current (development) version from this
-repository, run
+For advanced users: if you want to install the current (development)
+version from this repository, run
 
 ``` r
 devtools::install_github("https://github.com/TheoreticalEcology/s-jSDM", subdir = "sjSDM", ref = "master")
 ```
 
-Once the dependencies are installed, the following code should run:
+dependencies should be installed as above. If the installation fails,
+check out the help of ?install_sjSDM, ?installation_help, and
+vignette(“Dependencies”, package = “sjSDM”).
 
-## Workflow
+1.  Try install_sjSDM()
+2.  New session, if no ‘PyTorch not found’ appears it should work,
+    otherwise see ?installation_help
+3.  If do not get the pkg to run, create an issue [issue
+    tracker](https://github.com/TheoreticalEcology/s-jSDM/issues) or
+    write an email to maximilian.pichler at ur.de
+
+## Basic Workflow
 
-Simulate a community and fit a sjSDM model:
+Load the package
 
 ``` r
 library(sjSDM)
 ```
 
-    ## ── Attaching sjSDM ──────────────────────────────────────────────────── 1.0.4 ──
-
-    ## ✔ torch <environment> 
-    ## ✔ torch_optimizer  
-    ## ✔ pyro  
-    ## ✔ madgrad
+Simulate some community data
 
 ``` r
 set.seed(42)
@@ -86,186 +106,180 @@ 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)
+```
+
+This fits the standard SDM with environmental, spatial and covariance
+terms
 
+``` r
 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)
 ```
 
     ## Family:  binomial 
     ## 
-    ## LogLik:  -510.9816 
+    ## LogLik:  -505.3381 
     ## Regularization loss:  0 
     ## 
     ## Species-species correlation matrix: 
     ## 
     ##  sp1  1.0000                                 
-    ##  sp2 -0.3780  1.0000                             
-    ##  sp3 -0.2050 -0.4070  1.0000                         
-    ##  sp4 -0.1850 -0.3860  0.8220  1.0000                     
-    ##  sp5  0.6820 -0.4090 -0.1240 -0.0730  1.0000                 
-    ##  sp6 -0.3050  0.4870  0.1630  0.1510 -0.1220  1.0000             
-    ##  sp7  0.5830 -0.1190  0.0960  0.1200  0.5520  0.2450  1.0000         
-    ##  sp8  0.3140  0.1690 -0.5280 -0.5460  0.2330 -0.0480  0.1300  1.0000     
-    ##  sp9 -0.0620 -0.0250  0.0840  0.0640 -0.4010 -0.3430 -0.2060 -0.1380  1.0000 
-    ##  sp10     0.2080  0.4750 -0.7140 -0.6490  0.2540  0.1410  0.1480  0.4560 -0.2850  1.0000
+    ##  sp2 -0.3700  1.0000                             
+    ##  sp3 -0.1980 -0.4260  1.0000                         
+    ##  sp4 -0.1670 -0.3830  0.8330  1.0000                     
+    ##  sp5  0.6670 -0.3620 -0.1330 -0.1040  1.0000                 
+    ##  sp6 -0.2910  0.4780  0.1730  0.1860 -0.1020  1.0000             
+    ##  sp7  0.5740 -0.1110  0.1270  0.1790  0.5370  0.2740  1.0000         
+    ##  sp8  0.2870  0.2150 -0.5160 -0.5250  0.2000 -0.0060  0.1100  1.0000     
+    ##  sp9 -0.0610 -0.0560  0.0510  0.0580 -0.3950 -0.3590 -0.2320 -0.1310  1.0000 
+    ##  sp10     0.2050  0.5050 -0.7150 -0.6640  0.2550  0.1480  0.1380  0.4670 -0.2610  1.0000
     ## 
     ## 
     ## 
     ## Spatial: 
-    ##            sp1       sp2      sp3       sp4      sp5      sp6      sp7      sp8
-    ## X1:X2 2.103188 -4.041381 3.452883 0.2332844 2.681165 1.325118 3.126471 1.928931
-    ##             sp9     sp10
-    ## X1:X2 0.9001696 1.262238
+    ##           sp1       sp2      sp3       sp4      sp5       sp6      sp7      sp8
+    ## X1:X2 2.10835 -4.061844 3.446407 0.4750176 2.757261 0.9577486 3.384753 2.053963
+    ##            sp9     sp10
+    ## X1:X2 1.003981 1.293781
     ## 
     ## 
     ## 
     ##                  Estimate Std.Err Z value Pr(>|z|)    
-    ## sp1 (Intercept)   -0.0847  0.2671   -0.32  0.75124    
-    ## sp1 X1             1.3854  0.5241    2.64  0.00820 ** 
-    ## sp1 X2            -2.4736  0.4839   -5.11  3.2e-07 ***
-    ## sp1 X3            -0.2583  0.4362   -0.59  0.55385    
-    ## sp2 (Intercept)   -0.0145  0.2601   -0.06  0.95560    
-    ## sp2 X1             1.2578  0.5233    2.40  0.01625 *  
-    ## sp2 X2             0.2357  0.4909    0.48  0.63112    
-    ## sp2 X3             0.6825  0.4302    1.59  0.11266    
-    ## sp3 (Intercept)   -0.5653  0.2861   -1.98  0.04819 *  
-    ## sp3 X1             1.4285  0.5099    2.80  0.00509 ** 
-    ## sp3 X2            -0.4155  0.5096   -0.82  0.41489    
-    ## sp3 X3            -1.1364  0.4898   -2.32  0.02034 *  
-    ## sp4 (Intercept)   -0.1156  0.2580   -0.45  0.65406    
-    ## sp4 X1            -1.5792  0.4921   -3.21  0.00133 ** 
-    ## sp4 X2            -1.9313  0.5088   -3.80  0.00015 ***
-    ## sp4 X3            -0.4306  0.4314   -1.00  0.31822    
-    ## sp5 (Intercept)   -0.2109  0.2526   -0.83  0.40378    
-    ## sp5 X1             0.7425  0.4843    1.53  0.12525    
-    ## sp5 X2             0.5624  0.4582    1.23  0.21969    
-    ## sp5 X3            -0.7171  0.4154   -1.73  0.08433 .  
-    ## sp6 (Intercept)    0.2184  0.2707    0.81  0.41973    
-    ## sp6 X1             2.6087  0.5552    4.70  2.6e-06 ***
-    ## sp6 X2            -1.1176  0.5271   -2.12  0.03400 *  
-    ## sp6 X3             0.2021  0.4461    0.45  0.65049    
-    ## sp7 (Intercept)   -0.0719  0.2448   -0.29  0.76903    
-    ## sp7 X1            -0.3372  0.4899   -0.69  0.49132    
-    ## sp7 X2             0.3403  0.4328    0.79  0.43175    
-    ## sp7 X3            -1.4822  0.4269   -3.47  0.00052 ***
-    ## sp8 (Intercept)    0.1574  0.1625    0.97  0.33270    
-    ## sp8 X1             0.3657  0.3158    1.16  0.24688    
-    ## sp8 X2             0.3236  0.3102    1.04  0.29688    
-    ## sp8 X3            -1.2363  0.2850   -4.34  1.4e-05 ***
-    ## sp9 (Intercept)    0.0235  0.2003    0.12  0.90667    
-    ## sp9 X1             1.4160  0.3943    3.59  0.00033 ***
-    ## sp9 X2            -1.0606  0.3755   -2.82  0.00473 ** 
-    ## sp9 X3             0.7943  0.3444    2.31  0.02111 *  
-    ## sp10 (Intercept)  -0.0825  0.2076   -0.40  0.69104    
-    ## sp10 X1           -0.5510  0.3781   -1.46  0.14505    
-    ## sp10 X2           -1.3145  0.3777   -3.48  0.00050 ***
-    ## sp10 X3           -0.5257  0.3590   -1.46  0.14310    
+    ## sp1 (Intercept)   -0.1038  0.2614   -0.40  0.69129    
+    ## sp1 X1             1.3685  0.4894    2.80  0.00517 ** 
+    ## sp1 X2            -2.5386  0.4626   -5.49  4.1e-08 ***
+    ## sp1 X3            -0.2941  0.4331   -0.68  0.49713    
+    ## sp2 (Intercept)   -0.0106  0.2760   -0.04  0.96949    
+    ## sp2 X1             1.2541  0.5173    2.42  0.01534 *  
+    ## sp2 X2             0.2723  0.5312    0.51  0.60824    
+    ## sp2 X3             0.7237  0.4605    1.57  0.11603    
+    ## sp3 (Intercept)   -0.5153  0.2854   -1.81  0.07100 .  
+    ## sp3 X1             1.5114  0.5174    2.92  0.00349 ** 
+    ## sp3 X2            -0.4924  0.5080   -0.97  0.33235    
+    ## sp3 X3            -1.0819  0.4862   -2.23  0.02606 *  
+    ## sp4 (Intercept)   -0.0771  0.2559   -0.30  0.76318    
+    ## sp4 X1            -1.5116  0.4940   -3.06  0.00222 ** 
+    ## sp4 X2            -1.9738  0.4985   -3.96  7.5e-05 ***
+    ## sp4 X3            -0.3837  0.4295   -0.89  0.37164    
+    ## sp5 (Intercept)   -0.2368  0.2424   -0.98  0.32864    
+    ## sp5 X1             0.7438  0.4670    1.59  0.11121    
+    ## sp5 X2             0.5777  0.4317    1.34  0.18081    
+    ## sp5 X3            -0.7728  0.3984   -1.94  0.05244 .  
+    ## sp6 (Intercept)    0.3047  0.2753    1.11  0.26847    
+    ## sp6 X1             2.5735  0.6142    4.19  2.8e-05 ***
+    ## sp6 X2            -1.0934  0.5190   -2.11  0.03513 *  
+    ## sp6 X3             0.1742  0.4538    0.38  0.70103    
+    ## sp7 (Intercept)   -0.0224  0.2574   -0.09  0.93054    
+    ## sp7 X1            -0.3184  0.5029   -0.63  0.52667    
+    ## sp7 X2             0.3480  0.4616    0.75  0.45090    
+    ## sp7 X3            -1.5991  0.4387   -3.64  0.00027 ***
+    ## sp8 (Intercept)    0.1415  0.1673    0.85  0.39759    
+    ## sp8 X1             0.3254  0.3263    1.00  0.31864    
+    ## sp8 X2             0.3401  0.3154    1.08  0.28092    
+    ## sp8 X3            -1.2411  0.2907   -4.27  2.0e-05 ***
+    ## sp9 (Intercept)    0.0200  0.1955    0.10  0.91852    
+    ## sp9 X1             1.3218  0.3775    3.50  0.00046 ***
+    ## sp9 X2            -1.0367  0.3700   -2.80  0.00508 ** 
+    ## sp9 X3             0.7911  0.3337    2.37  0.01775 *  
+    ## sp10 (Intercept)  -0.0983  0.2091   -0.47  0.63821    
+    ## sp10 X1           -0.5550  0.3861   -1.44  0.15064    
+    ## sp10 X2           -1.2310  0.3852   -3.20  0.00139 ** 
+    ## sp10 X3           -0.5567  0.3514   -1.58  0.11315    
     ## ---
     ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 
-``` r
-plot(model)
-```
+We can also plot the nich estimates, i.e the estimates in the
+environmental component
 
     ## Family:  binomial 
     ## 
-    ## LogLik:  -510.9816 
+    ## LogLik:  -505.3381 
     ## Regularization loss:  0 
     ## 
     ## Species-species correlation matrix: 
     ## 
     ##  sp1  1.0000                                 
-    ##  sp2 -0.3780  1.0000                             
-    ##  sp3 -0.2050 -0.4070  1.0000                         
-    ##  sp4 -0.1850 -0.3860  0.8220  1.0000                     
-    ##  sp5  0.6820 -0.4090 -0.1240 -0.0730  1.0000                 
-    ##  sp6 -0.3050  0.4870  0.1630  0.1510 -0.1220  1.0000             
-    ##  sp7  0.5830 -0.1190  0.0960  0.1200  0.5520  0.2450  1.0000         
-    ##  sp8  0.3140  0.1690 -0.5280 -0.5460  0.2330 -0.0480  0.1300  1.0000     
-    ##  sp9 -0.0620 -0.0250  0.0840  0.0640 -0.4010 -0.3430 -0.2060 -0.1380  1.0000 
-    ##  sp10     0.2080  0.4750 -0.7140 -0.6490  0.2540  0.1410  0.1480  0.4560 -0.2850  1.0000
+    ##  sp2 -0.3700  1.0000                             
+    ##  sp3 -0.1980 -0.4260  1.0000                         
+    ##  sp4 -0.1670 -0.3830  0.8330  1.0000                     
+    ##  sp5  0.6670 -0.3620 -0.1330 -0.1040  1.0000                 
+    ##  sp6 -0.2910  0.4780  0.1730  0.1860 -0.1020  1.0000             
+    ##  sp7  0.5740 -0.1110  0.1270  0.1790  0.5370  0.2740  1.0000         
+    ##  sp8  0.2870  0.2150 -0.5160 -0.5250  0.2000 -0.0060  0.1100  1.0000     
+    ##  sp9 -0.0610 -0.0560  0.0510  0.0580 -0.3950 -0.3590 -0.2320 -0.1310  1.0000 
+    ##  sp10     0.2050  0.5050 -0.7150 -0.6640  0.2550  0.1480  0.1380  0.4670 -0.2610  1.0000
     ## 
     ## 
     ## 
     ## Spatial: 
-    ##            sp1       sp2      sp3       sp4      sp5      sp6      sp7      sp8
-    ## X1:X2 2.103188 -4.041381 3.452883 0.2332844 2.681165 1.325118 3.126471 1.928931
-    ##             sp9     sp10
-    ## X1:X2 0.9001696 1.262238
+    ##           sp1       sp2      sp3       sp4      sp5       sp6      sp7      sp8
+    ## X1:X2 2.10835 -4.061844 3.446407 0.4750176 2.757261 0.9577486 3.384753 2.053963
+    ##            sp9     sp10
+    ## X1:X2 1.003981 1.293781
     ## 
     ## 
     ## 
     ##                  Estimate Std.Err Z value Pr(>|z|)    
-    ## sp1 (Intercept)   -0.0847  0.2671   -0.32  0.75124    
-    ## sp1 X1             1.3854  0.5241    2.64  0.00820 ** 
-    ## sp1 X2            -2.4736  0.4839   -5.11  3.2e-07 ***
-    ## sp1 X3            -0.2583  0.4362   -0.59  0.55385    
-    ## sp2 (Intercept)   -0.0145  0.2601   -0.06  0.95560    
-    ## sp2 X1             1.2578  0.5233    2.40  0.01625 *  
-    ## sp2 X2             0.2357  0.4909    0.48  0.63112    
-    ## sp2 X3             0.6825  0.4302    1.59  0.11266    
-    ## sp3 (Intercept)   -0.5653  0.2861   -1.98  0.04819 *  
-    ## sp3 X1             1.4285  0.5099    2.80  0.00509 ** 
-    ## sp3 X2            -0.4155  0.5096   -0.82  0.41489    
-    ## sp3 X3            -1.1364  0.4898   -2.32  0.02034 *  
-    ## sp4 (Intercept)   -0.1156  0.2580   -0.45  0.65406    
-    ## sp4 X1            -1.5792  0.4921   -3.21  0.00133 ** 
-    ## sp4 X2            -1.9313  0.5088   -3.80  0.00015 ***
-    ## sp4 X3            -0.4306  0.4314   -1.00  0.31822    
-    ## sp5 (Intercept)   -0.2109  0.2526   -0.83  0.40378    
-    ## sp5 X1             0.7425  0.4843    1.53  0.12525    
-    ## sp5 X2             0.5624  0.4582    1.23  0.21969    
-    ## sp5 X3            -0.7171  0.4154   -1.73  0.08433 .  
-    ## sp6 (Intercept)    0.2184  0.2707    0.81  0.41973    
-    ## sp6 X1             2.6087  0.5552    4.70  2.6e-06 ***
-    ## sp6 X2            -1.1176  0.5271   -2.12  0.03400 *  
-    ## sp6 X3             0.2021  0.4461    0.45  0.65049    
-    ## sp7 (Intercept)   -0.0719  0.2448   -0.29  0.76903    
-    ## sp7 X1            -0.3372  0.4899   -0.69  0.49132    
-    ## sp7 X2             0.3403  0.4328    0.79  0.43175    
-    ## sp7 X3            -1.4822  0.4269   -3.47  0.00052 ***
-    ## sp8 (Intercept)    0.1574  0.1625    0.97  0.33270    
-    ## sp8 X1             0.3657  0.3158    1.16  0.24688    
-    ## sp8 X2             0.3236  0.3102    1.04  0.29688    
-    ## sp8 X3            -1.2363  0.2850   -4.34  1.4e-05 ***
-    ## sp9 (Intercept)    0.0235  0.2003    0.12  0.90667    
-    ## sp9 X1             1.4160  0.3943    3.59  0.00033 ***
-    ## sp9 X2            -1.0606  0.3755   -2.82  0.00473 ** 
-    ## sp9 X3             0.7943  0.3444    2.31  0.02111 *  
-    ## sp10 (Intercept)  -0.0825  0.2076   -0.40  0.69104    
-    ## sp10 X1           -0.5510  0.3781   -1.46  0.14505    
-    ## sp10 X2           -1.3145  0.3777   -3.48  0.00050 ***
-    ## sp10 X3           -0.5257  0.3590   -1.46  0.14310    
+    ## sp1 (Intercept)   -0.1038  0.2614   -0.40  0.69129    
+    ## sp1 X1             1.3685  0.4894    2.80  0.00517 ** 
+    ## sp1 X2            -2.5386  0.4626   -5.49  4.1e-08 ***
+    ## sp1 X3            -0.2941  0.4331   -0.68  0.49713    
+    ## sp2 (Intercept)   -0.0106  0.2760   -0.04  0.96949    
+    ## sp2 X1             1.2541  0.5173    2.42  0.01534 *  
+    ## sp2 X2             0.2723  0.5312    0.51  0.60824    
+    ## sp2 X3             0.7237  0.4605    1.57  0.11603    
+    ## sp3 (Intercept)   -0.5153  0.2854   -1.81  0.07100 .  
+    ## sp3 X1             1.5114  0.5174    2.92  0.00349 ** 
+    ## sp3 X2            -0.4924  0.5080   -0.97  0.33235    
+    ## sp3 X3            -1.0819  0.4862   -2.23  0.02606 *  
+    ## sp4 (Intercept)   -0.0771  0.2559   -0.30  0.76318    
+    ## sp4 X1            -1.5116  0.4940   -3.06  0.00222 ** 
+    ## sp4 X2            -1.9738  0.4985   -3.96  7.5e-05 ***
+    ## sp4 X3            -0.3837  0.4295   -0.89  0.37164    
+    ## sp5 (Intercept)   -0.2368  0.2424   -0.98  0.32864    
+    ## sp5 X1             0.7438  0.4670    1.59  0.11121    
+    ## sp5 X2             0.5777  0.4317    1.34  0.18081    
+    ## sp5 X3            -0.7728  0.3984   -1.94  0.05244 .  
+    ## sp6 (Intercept)    0.3047  0.2753    1.11  0.26847    
+    ## sp6 X1             2.5735  0.6142    4.19  2.8e-05 ***
+    ## sp6 X2            -1.0934  0.5190   -2.11  0.03513 *  
+    ## sp6 X3             0.1742  0.4538    0.38  0.70103    
+    ## sp7 (Intercept)   -0.0224  0.2574   -0.09  0.93054    
+    ## sp7 X1            -0.3184  0.5029   -0.63  0.52667    
+    ## sp7 X2             0.3480  0.4616    0.75  0.45090    
+    ## sp7 X3            -1.5991  0.4387   -3.64  0.00027 ***
+    ## sp8 (Intercept)    0.1415  0.1673    0.85  0.39759    
+    ## sp8 X1             0.3254  0.3263    1.00  0.31864    
+    ## sp8 X2             0.3401  0.3154    1.08  0.28092    
+    ## sp8 X3            -1.2411  0.2907   -4.27  2.0e-05 ***
+    ## sp9 (Intercept)    0.0200  0.1955    0.10  0.91852    
+    ## sp9 X1             1.3218  0.3775    3.50  0.00046 ***
+    ## sp9 X2            -1.0367  0.3700   -2.80  0.00508 ** 
+    ## sp9 X3             0.7911  0.3337    2.37  0.01775 *  
+    ## sp10 (Intercept)  -0.0983  0.2091   -0.47  0.63821    
+    ## sp10 X1           -0.5550  0.3861   -1.44  0.15064    
+    ## sp10 X2           -1.2310  0.3852   -3.20  0.00139 ** 
+    ## sp10 X3           -0.5567  0.3514   -1.58  0.11315    
     ## ---
     ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 
-![](README_files/figure-gfm/unnamed-chunk-5-1.png)<!-- -->
-
-We support other distributions:
-
-- Count data with Poisson:
-
-``` r
-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"))
-```
+![](README_files/figure-gfm/unnamed-chunk-8-1.png)<!-- -->
 
-- Count data with negative Binomial (which is still experimental, if you
-  run into errors/problems, please let us know):
+We can also visualize the species-species association matrix
 
 ``` r
-model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1+X2+X3), spatial = linear(data = SP, formula = ~0+X1:X2), se = TR, family="nbinom")
+image(getCor(model))
 ```
 
-- Gaussian (normal):
+![](README_files/figure-gfm/unnamed-chunk-9-1.png)<!-- -->
 
-``` r
-model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1+X2+X3), spatial = linear(data = SP, formula = ~0+X1:X2), se = TR, family=gaussian("identity"))
-```
+### Anova / Variation partitioning
 
-### Anova
+As in other models, it can be interesting to analyze how much variation
+is explained by which parts of hte model.
 
-ANOVA can be used to partition the three components (abiotic, biotic,
-and spatial):
+![image](sjSDM/vignettes/jSDM-ANOVA.png) For the Env, Spatial,
+Covariance terms, this is implemented in
 
 ``` r
 an = anova(model)
@@ -274,19 +288,17 @@ print(an)
 
     ## Analysis of Deviance Table
     ## 
-    ## Terms added sequentially:
-    ## 
     ##           Deviance Residual deviance R2 Nagelkerke R2 McFadden
-    ## Abiotic  157.95722        1177.48500       0.79394      0.1139
-    ## Biotic   175.41278        1160.02944       0.82694      0.1265
-    ## Spatial   17.38643        1318.05579       0.15959      0.0125
-    ## Full     385.98836         949.45386       0.97893      0.2784
+    ## Abiotic 1.9347e+02        1.0452e+03    8.5554e-01      0.1413
+    ## Biotic  2.0839e+02        1.0303e+03    8.7556e-01      0.1522
+    ## Spatial 9.6113e+00        1.2291e+03    9.1639e-02      0.0070
+    ## Full    3.6694e+02        8.7175e+02    9.7451e-01      0.2680
 
 ``` r
 plot(an)
 ```
 
-![](README_files/figure-gfm/unnamed-chunk-9-1.png)<!-- -->
+![](README_files/figure-gfm/unnamed-chunk-10-1.png)<!-- -->
 
 The anova shows the relative changes in the R<sup>2</sup> of the groups
 and their intersections.
@@ -309,7 +321,7 @@ results = plotInternalStructure(an) # or plot(an, internal = TRUE)
     ##   plot.ggplot      ggplot2
     ##   print.ggplot     ggplot2
 
-![](README_files/figure-gfm/unnamed-chunk-10-1.png)<!-- -->
+![](README_files/figure-gfm/unnamed-chunk-11-1.png)<!-- -->
 
 The plot function returns the results for the internal metacommunity
 structure:
@@ -318,29 +330,17 @@ structure:
 print(results$data$Species)
 ```
 
-    ##           env         spa     codist         r2
-    ## 1  0.17677667 0.000000000 0.16810146 0.03375475
-    ## 2  0.08724636 0.026656011 0.18072040 0.02946228
-    ## 3  0.12613742 0.004529856 0.21004115 0.03407084
-    ## 4  0.16648179 0.000000000 0.15890110 0.03241345
-    ## 5  0.08585343 0.005811074 0.16168802 0.02533525
-    ## 6  0.18787936 0.012341719 0.11489709 0.03151182
-    ## 7  0.10765006 0.012898782 0.13292549 0.02534743
-    ## 8  0.12445149 0.015040332 0.06188116 0.02013730
-    ## 9  0.17762242 0.000000000 0.04357315 0.02196159
-    ## 10 0.08805858 0.017406001 0.13890574 0.02443703
-
-## Installation trouble shooting
-
-If the installation fails, check out the help of ?install_sjSDM,
-?installation_help, and vignette(“Dependencies”, package = “sjSDM”).
-
-1.  Try install_sjSDM()
-2.  New session, if no ‘PyTorch not found’ appears it should work,
-    otherwise see ?installation_help
-3.  If do not get the pkg to run, create an issue [issue
-    tracker](https://github.com/TheoreticalEcology/s-jSDM/issues) or
-    write an email to maximilian.pichler at ur.de
+    ##           env         spa     codist        r2
+    ## 1  0.19116596 0.000000000 0.15781899 0.3308529
+    ## 2  0.10087801 0.014995113 0.20387047 0.2706637
+    ## 3  0.14629446 0.016001186 0.25851625 0.3240496
+    ## 4  0.19067561 0.002588795 0.17070044 0.3312707
+    ## 5  0.09878869 0.000000000 0.19709887 0.2317703
+    ## 6  0.18667429 0.006040971 0.12480417 0.3027200
+    ## 7  0.12154125 0.007215252 0.14308883 0.2395096
+    ## 8  0.11602783 0.006021371 0.05402848 0.1866853
+    ## 9  0.14240696 0.016904452 0.04469078 0.2249502
+    ## 10 0.11681332 0.008436236 0.17090677 0.2372738
 
 ### Python Package
 
@@ -365,7 +365,7 @@ model.fit(Env, Occ, batch_size = 20, epochs = 10)
 # print(model.covariance)
 ```
 
-    ## Iter: 0/10   0%|          | [00:00, ?it/s]Iter: 0/10   0%|          | [00:00, ?it/s, loss=7.17]Iter: 1/10  10%|#         | [00:00,  3.88it/s, loss=7.17]Iter: 1/10  10%|#         | [00:00,  3.88it/s, loss=7.158]Iter: 1/10  10%|#         | [00:00,  3.88it/s, loss=7.149]Iter: 1/10  10%|#         | [00:00,  3.88it/s, loss=7.142]Iter: 1/10  10%|#         | [00:00,  3.88it/s, loss=7.124]Iter: 1/10  10%|#         | [00:00,  3.88it/s, loss=7.118]Iter: 1/10  10%|#         | [00:00,  3.88it/s, loss=7.112]Iter: 1/10  10%|#         | [00:00,  3.88it/s, loss=7.108]Iter: 1/10  10%|#         | [00:00,  3.88it/s, loss=7.118]Iter: 9/10  90%|######### | [00:00, 30.63it/s, loss=7.118]Iter: 9/10  90%|######### | [00:00, 30.63it/s, loss=7.144]Iter: 10/10 100%|##########| [00:00, 26.74it/s, loss=7.144]
+    ## Iter: 0/10   0%|          | [00:00, ?it/s]Iter: 0/10   0%|          | [00:00, ?it/s, loss=7.162]Iter: 0/10   0%|          | [00:00, ?it/s, loss=7.151]Iter: 0/10   0%|          | [00:00, ?it/s, loss=7.152]Iter: 0/10   0%|          | [00:00, ?it/s, loss=7.158]Iter: 0/10   0%|          | [00:00, ?it/s, loss=7.146]Iter: 0/10   0%|          | [00:00, ?it/s, loss=7.139]Iter: 0/10   0%|          | [00:00, ?it/s, loss=7.148]Iter: 0/10   0%|          | [00:00, ?it/s, loss=7.127]Iter: 8/10  80%|########  | [00:00, 78.42it/s, loss=7.127]Iter: 8/10  80%|########  | [00:00, 78.42it/s, loss=7.114]Iter: 8/10  80%|########  | [00:00, 78.42it/s, loss=7.138]Iter: 10/10 100%|##########| [00:00, 78.65it/s, loss=7.138]
 
 Calculate Importance:
 
@@ -377,25 +377,25 @@ covX = fa.covariance( torch.tensor(Env).t() ).data.cpu().numpy()
 fa.importance(beta=Beta, covX=covX, sigma=Sigma)
 ```
 
-    ## {'env': array([[ 1.2717709e-02,  7.8437300e-03,  6.6514793e-03, -3.3015787e-04,
-    ##          2.3806898e-04],
-    ##        [ 9.9158729e-05, -1.8891758e-06,  1.0537009e-03,  4.0511694e-04,
-    ##          1.1120385e-02],
-    ##        [ 6.1564189e-03,  5.9850062e-03,  9.2307013e-03,  5.4843356e-03,
-    ##         -3.3683516e-04],
-    ##        [ 1.3349474e-02,  4.3294221e-04,  1.8103119e-03,  1.4068705e-02,
-    ##          6.6316797e-04],
-    ##        [ 3.3953122e-05,  3.1304134e-03,  2.6658648e-03, -3.6165391e-05,
-    ##          7.3677581e-03],
-    ##        [ 2.7722977e-03,  1.9519718e-03,  4.8086399e-04,  3.0876237e-03,
-    ##          1.7828522e-04],
-    ##        [ 7.9284189e-03,  5.7881157e-04,  7.5722663e-03,  2.1802005e-06,
-    ##          4.2433664e-03],
-    ##        [-5.1329907e-06,  6.0144444e-03,  2.1059261e-05,  7.5954124e-03,
-    ##          1.5537007e-03],
-    ##        [ 7.7161851e-04,  1.7209088e-02,  4.8407568e-03,  1.8020724e-03,
-    ##          5.6920521e-04],
-    ##        [ 1.7108561e-02,  7.2742125e-04,  9.4651995e-04,  6.8342132e-03,
-    ##          1.1830850e-02]], dtype=float32), 'biotic': array([0.9728792 , 0.9873236 , 0.9734804 , 0.9696754 , 0.9868381 ,
-    ##        0.9915289 , 0.97967494, 0.9848205 , 0.9748072 , 0.9625524 ],
+    ## {'env': array([[ 7.98111968e-03,  1.13930255e-02,  1.57148903e-03,
+    ##          1.94198191e-02,  1.04639726e-03],
+    ##        [ 1.18241680e-03,  5.13114082e-03,  4.19588014e-03,
+    ##          1.49956578e-02,  1.00601939e-02],
+    ##        [ 1.30655477e-02,  3.98503980e-05,  9.63684451e-03,
+    ##          1.64120211e-04,  4.73391265e-03],
+    ##        [ 6.42327068e-04,  9.29863378e-03,  2.36115797e-04,
+    ##          3.75465327e-03,  1.57004467e-03],
+    ##        [ 3.10221454e-03,  1.19799370e-04,  2.72720726e-03,
+    ##          9.55119822e-03,  1.01295421e-02],
+    ##        [ 2.23755743e-03,  1.17474922e-03,  5.54808835e-03,
+    ##          2.66983203e-04,  1.72002669e-02],
+    ##        [ 1.64179280e-02,  1.61731802e-02,  8.98935087e-03,
+    ##          4.06762323e-04,  2.52083852e-03],
+    ##        [ 6.37685647e-03, -1.96362438e-04,  9.27473884e-03,
+    ##          8.93264223e-05,  1.40198926e-03],
+    ##        [ 1.06887252e-03,  1.00843166e-03,  1.96674443e-03,
+    ##          1.43694794e-02,  3.84828472e-03],
+    ##        [ 3.29102017e-03,  5.51448902e-03,  6.90532179e-05,
+    ##          5.58465905e-03,  4.97359270e-03]], dtype=float32), 'biotic': array([0.9585881 , 0.96443474, 0.9723598 , 0.98449826, 0.97437006,
+    ##        0.9735723 , 0.95549196, 0.98305345, 0.9777382 , 0.98056716],
     ##       dtype=float32)}
diff --git a/sjSDM/DESCRIPTION b/sjSDM/DESCRIPTION
index cee7276..cac2d0a 100644
--- a/sjSDM/DESCRIPTION
+++ b/sjSDM/DESCRIPTION
@@ -17,7 +17,7 @@ Authors@R:
              family = "Cai",
              role = "ctb",
              email = "caiwang0503@163.com"))
-Description: A scalable method to estimate joint Species Distribution Models (jSDMs) for big community datasets based on a Monte Carlo approximation of the joint likelihood.  The numerical approximation is based on 'PyTorch' and 'reticulate', and can be run on CPUs and GPUs alike. The method is described in Pichler & Hartig (2021) <doi:10.1111/2041-210X.13687>. The package contains various extensions, including support for different response families, ability to account for spatial autocorrelation, and deep neural networks instead of the linear predictor in jSDMs.
+Description: A scalable and fast method for estimating joint Species Distribution Models (jSDMs) for big community data, including eDNA data. The package estimates a full (i.e. non-latent) jSDM with different response distributions (including the traditional multivariate probit model). The package allows to perform variation partitioning (VP) / ANOVA on the fitted models to separate the contribution of environmental, spatial, and biotic associations. In addition, the total R-squared can be further partitioned per species and site to reveal the internal metacommunity structure, see Leibold et al., <doi:10.1111/oik.08618>. The internal structure can then be regressed against environmental and spatial distinctiveness, richness, and traits to analyze metacommunity assembly processes.  The package includes support for accounting for spatial autocorrelation and the option to fit responses using deep neural networks instead of a standard linear predictor. As described in Pichler & Hartig (2021) <doi:10.1111/2041-210X.13687>, scalability is achieved by using a Monte Carlo approximation of the joint likelihood implemented via 'PyTorch' and 'reticulate', which can be run on CPUs or GPUs.
 License: GPL-3
 Encoding: UTF-8
 LazyData: true
@@ -40,14 +40,17 @@ Imports:
     ggplot2,
     checkmate,
     mathjaxr,
-    ggtern
+    ggtern,
+    beeswarm,
+    qgam,
+    scales
 Suggests:
     testthat,
     knitr,
     rmarkdown,
     iml,
     fields
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.2
 URL: https://theoreticalecology.github.io/s-jSDM/
 BugReports: https://github.com/TheoreticalEcology/s-jSDM/issues
 Roxygen: list(old_usage = FALSE)
diff --git a/sjSDM/NAMESPACE b/sjSDM/NAMESPACE
index a037379..de717aa 100644
--- a/sjSDM/NAMESPACE
+++ b/sjSDM/NAMESPACE
@@ -10,6 +10,7 @@ S3method(plot,sjSDM)
 S3method(plot,sjSDM_cv)
 S3method(plot,sjSDManova)
 S3method(plot,sjSDMimportance)
+S3method(plot,sjSDMinternalStruture)
 S3method(predict,sjSDM)
 S3method(print,DNN)
 S3method(print,bioticStruct)
@@ -18,10 +19,13 @@ S3method(print,sjSDM)
 S3method(print,sjSDM_cv)
 S3method(print,sjSDManova)
 S3method(print,sjSDMimportance)
+S3method(print,sjSDMinternalStruture)
+S3method(residuals,sjSDM)
 S3method(setWeights,sjSDM)
 S3method(simulate,sjSDM)
 S3method(summary,sjSDM)
 S3method(summary,sjSDM_cv)
+S3method(summary,sjSDManova)
 S3method(update,sjSDM)
 export(AccSGD)
 export(AdaBound)
@@ -36,14 +40,13 @@ export(getCor)
 export(getCov)
 export(getSe)
 export(getWeights)
-export(importance)
 export(install_diagnostic)
 export(install_sjSDM)
+export(internalStructure)
 export(is_torch_available)
 export(linear)
 export(madgrad)
 export(plotAssemblyEffects)
-export(plotInternalStructure)
 export(setWeights)
 export(simulate_SDM)
 export(sjSDM)
diff --git a/sjSDM/R/anova.R b/sjSDM/R/anova.R
index 2539a64..7ab7624 100644
--- a/sjSDM/R/anova.R
+++ b/sjSDM/R/anova.R
@@ -1,17 +1,16 @@
-#' Anova
+#' Anova / Variation partitioning
 #' 
-#' Compute analysis of variance
+#' Compute variance explained by the three fractions env, space, associations
 #' 
 #' @param object model of object \code{\link{sjSDM}}
 #' @param samples Number of Monte Carlo samples
-#' @param fractional how to distribute the shared partitions, proportional to the fractions' R2s or equal
 #' @param ... optional arguments which are passed to the calculation of the logLikelihood
 #' 
-#' @description
+#' @details The ANOVA function removes each of the three fractions (Environment, Space, Associations) and measures the drop in variance explained, and thus the importance of the three fractions.
 #' 
-#' Calculates type II anova.
+#' Variance explained is measured by Deviance as well as the pseudo-R2 metrics of Nagelkerke and McFadden
 #' 
-#' Shared contributions (e.g. between space and environment) are also calculated (and divided proportionally) and can be optionally visualized via \code{\link{plot.sjSDManova}} with \code{add_shared=TRUE}.
+#' In downstream functions such as \code{\link{plot.sjSDManova}} or \code{\link{plot.sjSDManova}} with \code{add_shared=TRUE}.
 #' The anova can get unstable for many species and few occurrences/observations. We recommend using large numbers for 'samples'.
 #' 
 #' @return 
@@ -28,11 +27,14 @@
 #' 
 #' Implemented S3 methods are \code{\link{print.sjSDManova}} and \code{\link{plot.sjSDManova}}
 #'  
-#' @seealso \code{\link{plot.sjSDManova}}, \code{\link{print.sjSDManova}}, \code{\link{plotInternalStructure}}
+#' @seealso \code{\link{plot.sjSDManova}}, \code{\link{print.sjSDManova}},\code{\link{summary.sjSDManova}}, \code{\link{plot.sjSDMinternalStruture}}
+#' 
+#' @example /inst/examples/anova-example.R
+#' 
 #' @import stats
 #' @export
 
-anova.sjSDM = function(object, samples = 5000L, fractional = c("proportional", "equal"), ...) {
+anova.sjSDM = function(object, samples = 5000L, ...) {
   out = list()
   individual = TRUE
   samples = as.integer(samples)
@@ -69,9 +71,23 @@ anova.sjSDM = function(object, samples = 5000L, fractional = c("proportional", "
     F_AB = full_wo - A_wo-B_wo
 
     anova_rows = c("Null", "F_A", "F_B", "Full")
-    names(anova_rows) = c("Null", "Abiotic", "Biotic", "Full")
-    results = data.frame(models = c("F_A", "F_B","F_AB","Full", "Saturated", "Null"),
+    names(anova_rows) = c("Null", "Abiotic", "Assocations", "Full")
+    
+    results_discard = data.frame(models = c("F_A", "F_B","F_AB","Full", "Saturated", "Null"),
                          ll = -c(sum(null_m) + sum(F_A), sum(null_m)  + sum(F_B), sum(null_m)  + sum(F_AB), sum(full_m), sum(SAT_m), sum(null_m)))
+    F_AA = F_A + F_AB*abs(F_A)/(abs(F_A)+abs(F_B))
+    F_BB = F_B + F_AB*abs(F_B)/(abs(F_A)+abs(F_B))
+
+    results_proportional = data.frame(models = c("F_A", "F_B","F_AB","Full", "Saturated", "Null"),
+                                 ll = -c(sum(null_m) + sum(F_AA), sum(null_m)  + sum(F_BB), sum(null_m)  + sum(F_AB), sum(full_m), sum(SAT_m), sum(null_m)))
+    
+    F_AA = F_A + F_AB*0.5
+    F_BB = F_B + F_AB*0.5
+    
+    results_equal = data.frame(models = c("F_A", "F_B","F_AB","Full", "Saturated", "Null"),
+                                 ll = -c(sum(null_m) + sum(F_AA), sum(null_m)  + sum(F_BB), sum(null_m)  + sum(F_AB), sum(full_m), sum(SAT_m), sum(null_m)))
+    
+    
     results_ind = list("F_A"=-(null_m + F_A), "F_B"=-(null_m +F_B), "F_AB"=-(null_m + F_AB), "A" = -A_m, "B" = -B_m, "Full"=-full_m, "Saturated"=-SAT_m, "Null"=-null_m)
     
   } else {
@@ -108,26 +124,55 @@ anova.sjSDM = function(object, samples = 5000L, fractional = c("proportional", "
     F_AS = full_wo - B_wo - F_A - F_S
     F_AB = full_wo - S_wo - F_A - F_B
     F_BS = full_wo - A_wo - F_S - F_B
-    F_ABS = full_wo - F_BS - F_AB- F_AS- F_A- F_B - F_S        
-    results = data.frame(models = c("F_A", "F_B","F_S","F_AB","F_AS", "F_BS", "F_ABS", "Full", "Saturated", "Null"),
+    F_ABS = full_wo - F_BS - F_AB- F_AS- F_A- F_B - F_S
+    
+    ## discard
+    
+    results_discard = data.frame(models = c("F_A", "F_B","F_S","F_AB","F_AS", "F_BS", "F_ABS", "Full", "Saturated", "Null"),
                          ll = -c(sum(null_m) + sum(F_A), sum(null_m) + sum(F_B),sum(null_m) + sum(F_S), 
                                  sum(null_m) + sum(F_AB), sum(null_m) + sum(F_AS), sum(null_m) + sum(F_BS), sum(null_m) + sum(F_ABS), 
                                  sum(null_m) + sum(full_wo), sum(SAT_m), sum(null_m)))
+    ## proportional
+    F_AA = F_A + F_AB*abs(F_A)/(abs(F_A)+abs(F_B)) + F_AS*abs(F_A)/(abs(F_S)+abs(F_A))+ F_ABS*abs(F_A)/(abs(F_A)+abs(F_B)+abs(F_S))
+    F_BB = F_B + F_AB*abs(F_B)/(abs(F_A)+abs(F_B)) + F_BS*abs(F_B)/(abs(F_S)+abs(F_B))+ F_ABS*abs(F_B)/(abs(F_A)+abs(F_B)+abs(F_S))
+    F_SS = F_S + F_AS*abs(F_S)/(abs(F_S)+abs(F_A)) + F_BS*abs(F_S)/(abs(F_S)+abs(F_B))+ F_ABS*abs(F_S)/(abs(F_A)+abs(F_B)+abs(F_S))
+    results_proportional = data.frame(models = c("F_A", "F_B","F_S","F_AB","F_AS", "F_BS", "F_ABS", "Full", "Saturated", "Null"),
+                                     ll = -c(sum(null_m) + sum(F_AA, na.rm = TRUE), sum(null_m) + sum(F_BB, na.rm = TRUE),sum(null_m) + sum(F_SS, na.rm = TRUE), 
+                                           sum(null_m) + sum(F_AB), sum(null_m) + sum(F_AS), sum(null_m) + sum(F_BS), sum(null_m) + sum(F_ABS), 
+                                          sum(null_m) + sum(full_wo), sum(SAT_m), sum(null_m)))
+
+
+    
+    ## equal
+    F_AA = F_A + F_AB*0.3333333 + F_AS*0.3333333+ F_ABS*0.3333333
+    F_BB = F_B + F_AB*0.3333333 + F_BS*0.3333333+ F_ABS*0.3333333
+    F_SS = F_S + F_AB*0.3333333 + F_BS*0.3333333+ F_ABS*0.3333333
+    
+    results_equal = data.frame(models = c("F_A", "F_B","F_S","F_AB","F_AS", "F_BS", "F_ABS", "Full", "Saturated", "Null"),
+                                      ll = -c(sum(null_m) + sum(F_AA), sum(null_m) + sum(F_BB),sum(null_m) + sum(F_SS), 
+                                              sum(null_m) + sum(F_AB), sum(null_m) + sum(F_AS), sum(null_m) + sum(F_BS), sum(null_m) + sum(F_ABS), 
+                                              sum(null_m) + sum(full_wo), sum(SAT_m), sum(null_m)))
+    
     results_ind = list("F_A"=-(null_m + F_A), "F_B"=-(null_m + F_B),"F_S"=-(null_m +F_S), "F_AB"=-(null_m +F_AB),"F_AS"=-(null_m + F_AS), 
                        "F_BS"=-(null_m + F_BS), "F_ABS"=-(null_m + F_ABS), "A" = -A_m, "B" = -B_m, "S" = -S_m, AB = AB_m, AS = AS_m, BS = BS_m,
                        "Full"=-(full_m), "Saturated"= -(SAT_m), "Null"=-null_m)
     
     anova_rows = c("Null", "F_A", "F_B", "F_S", "Full")
-    names(anova_rows) = c("Null", "Abiotic", "Biotic", "Spatial", "Full")
+    names(anova_rows) = c("Null", "Abiotic", "Assocations", "Spatial", "Full")
   }
-  
-  results$`Residual deviance` = -2*(results$ll - results$ll[which(results$models == "Saturated", arr.ind = TRUE)])
-  
-  results$Deviance = results$`Residual deviance`[which(results$models == "Null", arr.ind = TRUE)] - results$`Residual deviance`
-  R21 = function(a, b) return(1-exp(2/(nrow(object$data$Y))*(-a+b)))
-  results$`R2 Nagelkerke` = R21(rep(-results$ll[which(results$models == "Null", arr.ind = TRUE)], length(results$ll)), - results$ll)
-  R22 = function(a, b) 1 - (b/a)
-  results$`R2 McFadden`= R22(rep(results$ll[which(results$models == "Null", arr.ind = TRUE)], length(results$ll)), results$ll)
+  results = 
+    lapply(list(results_discard, results_proportional, results_equal), function(res) {
+      
+      res$`Residual deviance` = -2*(res$ll - res$ll[which(res$models == "Saturated", arr.ind = TRUE)])
+      
+      res$Deviance = res$`Residual deviance`[which(res$models == "Null", arr.ind = TRUE)] - res$`Residual deviance`
+      R21 = function(a, b) return(1-exp(2/(nrow(object$data$Y))*(-a+b)))
+      res$`R2 Nagelkerke` = R21(rep(-res$ll[which(res$models == "Null", arr.ind = TRUE)], length(res$ll)), - res$ll)
+      R22 = function(a, b) 1 - (b/a)
+      res$`R2 McFadden`= R22(rep(res$ll[which(res$models == "Null", arr.ind = TRUE)], length(res$ll)), res$ll)
+      return(res)
+    
+    })
   
   # individual
   Residual_deviance_ind = lapply(results_ind, function(r) r - results_ind$Saturated)
@@ -141,23 +186,44 @@ anova.sjSDM = function(object, samples = 5000L, fractional = c("proportional", "
   R2_McFadden_ind = lapply(results_ind, function(r) R222(colSums(results_ind$Null), colSums(r)))
   R2_McFadden_sites = lapply(results_ind, function(r) R222(rowSums(results_ind$Null), rowSums(r)))
   
-  R2_McFadden_ind_shared = get_shared_anova(R2_McFadden_ind, fractional = fractional)
-  R2_McFadden_sites_shared = get_shared_anova(R2_McFadden_sites, fractional = fractional)
-  R2_Nagelkerke_ind_shared = get_shared_anova(R2_Nagelkerke_ind, fractional = fractional)
-  R2_Nagelkerke_sites_shared = get_shared_anova(R2_Nagelkerke_sites, fractional = fractional)
+  R2_McFadden_ind_shared = get_shared_anova(R2_McFadden_ind)
+  R2_McFadden_sites_shared = get_shared_anova(R2_McFadden_sites)
+  R2_Nagelkerke_ind_shared = get_shared_anova(R2_Nagelkerke_ind)
+  R2_Nagelkerke_sites_shared = get_shared_anova(R2_Nagelkerke_sites)
+  
+  #R2_McFadden_ind$Full = correct_R2(R2_McFadden_ind$Full)
+  #R2_McFadden_sites$Full = correct_R2(R2_McFadden_sites$Full)
+  
+  # precalculates reduced ANOVA tables
+  calculateResults <- function(res) {
+    rownames(res) = res$models
+    res = res[anova_rows,]
+    res$models = names(anova_rows)
+    res = res[-1,c(1, 4, 3,5,6)]
+    rownames(res) = res$models
+    res = res[,-1]
+    return(res)
+  }
+  
+  if(inherits(object, "spatial")) {
+    printFull = results[[1]][1:8,c(1, 4, 3,5,6)]
+    rownames(printFull) = printFull$models
+    printFull = printFull[,-1]
+    rownames(printFull) = c("Abiotic", "Associations","Spatial", "Shared Abiotic+Associations", "Shared Abiotic+Spatial", "Shared Spatial+Associations", "Shared Abiotic+Associations+Spatial", "Full")
+  } else {
+    printFull = results[[1]][1:4,c(1, 4, 3,5,6)]
+    rownames(printFull) = printFull$models
+    printFull = printFull[,-1]
+    rownames(printFull) = c("Abiotic", "Associations", "Shared Abiotic+Associations", "Full")
+  }
   
-  R2_McFadden_ind$Full = correct_R2(R2_McFadden_ind$Full)
-  R2_McFadden_sites$Full = correct_R2(R2_McFadden_sites$Full)
+  toPrint = list(all = printFull,
+       discard = calculateResults(results[[1]]), 
+       proportional = calculateResults(results[[2]]), 
+       equal = calculateResults(results[[3]])) 
   
-  to_print = results
-  rownames(to_print) = to_print$models
-  to_print = to_print[anova_rows,]
-  to_print$models = names(anova_rows)
-  to_print = to_print[-1,c(1, 4, 3,5,6)]
-  rownames(to_print) = to_print$models
-  to_print = to_print[,-1]
-  out$results = results
-  out$to_print = to_print
+  out$results = results[[1]] # TODO Max: check - das hier sind doch die einzigen Resultate die wir brauchen, oder? Es gibt doch eigentlich nur eine Aufteilung
+  out$to_print = toPrint
   out$N = nrow(object$data$Y)
   out$species = list(Residual_deviance = Residual_deviance_ind,
                      Deviance = Deviance_ind,
@@ -176,6 +242,7 @@ anova.sjSDM = function(object, samples = 5000L, fractional = c("proportional", "
   return(out)
 }
 
+# TODO - scale instead of 
 correct_R2 = function(R2) {
   R2 = ifelse(R2 < 0, 0, R2)
   R2 = ifelse(R2 > 1.000, 0, R2)
@@ -209,27 +276,13 @@ get_conditional_lls = function(m, null_m, ...) {
   rescaled_conditional_lls = null_m - matrix(rates, nrow = nrow(m$data$Y), ncol = ncol(m$data$Y), byrow = TRUE) * (rowSums(null_m)-joint_ll)
   
   ### Maximal/Minimal 0?
-  rescaled_conditional_lls[rescaled_conditional_lls<0] = 0
+  #rescaled_conditional_lls[rescaled_conditional_lls<0] = 0 does not work!
   
   return(rescaled_conditional_lls)
 }
 
-get_shared_anova = function(R2objt, spatial = TRUE, fractional = c("proportional", "equal")) {
-  fractional = match.arg(fractional)
+get_shared_anova = function(R2objt, spatial = TRUE) {
   if(spatial) {
-    # F_BS = R2objt$Full-R2objt$A
-    # F_AB = R2objt$Full-R2objt$S
-    # F_AS = R2objt$Full-R2objt$B
-    # F_A = R2objt$Full- R2objt$BS
-    # F_B =  R2objt$Full-R2objt$AS
-    # F_S =  R2objt$Full-R2objt$AB
-    # F_BS = R2objt$Full-R2objt$A - F_B -F_S
-    # F_AB = F_AB - F_A -F_B
-    # F_AS = F_AS - F_A -F_S
-    # F_ABS = R2objt$Full - F_BS - F_AB- F_AS- F_A- F_B - F_S
-    # A = F_A + F_AB*abs(F_A)/(abs(F_A)+abs(F_B)) + F_AS*abs(F_A)/(abs(F_S)+abs(F_A))+ F_ABS*abs(F_A)/(abs(F_A)+abs(F_B)+abs(F_S))
-    # B = F_B + F_AB*abs(F_B)/(abs(F_A)+abs(F_B)) + F_BS*abs(F_B)/(abs(F_S)+abs(F_B))+ F_ABS*abs(F_B)/(abs(F_A)+abs(F_B)+abs(F_S))
-    # S = F_S + F_AS*abs(F_S)/(abs(F_S)+abs(F_A)) + F_BS*abs(F_S)/(abs(F_S)+abs(F_B))+ F_ABS*abs(F_S)/(abs(F_A)+abs(F_B)+abs(F_S))
     F_A <- R2objt$Full - R2objt$F_BS
     F_B <- R2objt$Full - R2objt$F_AB
     F_S <- R2objt$Full - R2objt$F_AS
@@ -237,15 +290,15 @@ get_shared_anova = function(R2objt, spatial = TRUE, fractional = c("proportional
     F_AB <- R2objt$Full - R2objt$F_S - (F_A + F_B)
     F_AS <- R2objt$Full - R2objt$F_B - (F_A + F_S)
     F_ABS <- R2objt$Full - (F_A + F_B + F_S + F_BS + F_AB + F_AS)
-    if(fractional == "proportional") {
-      A = F_A + F_AB*abs(F_A)/(abs(F_A)+abs(F_B)) + F_AS*abs(F_A)/(abs(F_S)+abs(F_A))+ F_ABS*abs(F_A)/(abs(F_A)+abs(F_B)+abs(F_S))
-      B = F_B + F_AB*abs(F_B)/(abs(F_A)+abs(F_B)) + F_BS*abs(F_B)/(abs(F_S)+abs(F_B))+ F_ABS*abs(F_B)/(abs(F_A)+abs(F_B)+abs(F_S))
-      S = F_S + F_AS*abs(F_S)/(abs(F_S)+abs(F_A)) + F_BS*abs(F_S)/(abs(F_S)+abs(F_B))+ F_ABS*abs(F_S)/(abs(F_A)+abs(F_B)+abs(F_S))
-    } else {
-      A = F_A + F_AB*0.3333333 + F_AS*0.3333333+ F_ABS*0.3333333
-      B = F_B + F_AB*0.3333333 + F_BS*0.3333333+ F_ABS*0.3333333
-      S = F_S + F_AS*0.3333333 + F_BS*0.3333333+ F_ABS*0.3333333
-    }
+    A = F_A + F_AB*abs(F_A)/(abs(F_A)+abs(F_B)) + F_AS*abs(F_A)/(abs(F_S)+abs(F_A))+ F_ABS*abs(F_A)/(abs(F_A)+abs(F_B)+abs(F_S))
+    B = F_B + F_AB*abs(F_B)/(abs(F_A)+abs(F_B)) + F_BS*abs(F_B)/(abs(F_S)+abs(F_B))+ F_ABS*abs(F_B)/(abs(F_A)+abs(F_B)+abs(F_S))
+    S = F_S + F_AS*abs(F_S)/(abs(F_S)+abs(F_A)) + F_BS*abs(F_S)/(abs(F_S)+abs(F_B))+ F_ABS*abs(F_S)/(abs(F_A)+abs(F_B)+abs(F_S))
+    # R2 = correct_R2(R2objt$Full) TODO Check that this can be gone
+    proportional = list(F_A = A, F_B = B, F_S = S, R2 = R2objt$Full)
+    A = F_A + F_AB*0.3333333 + F_AS*0.3333333+ F_ABS*0.3333333
+    B = F_B + F_AB*0.3333333 + F_BS*0.3333333+ F_ABS*0.3333333
+    S = F_S + F_AB*0.3333333 + F_BS*0.3333333+ F_ABS*0.3333333
+    equal = list(F_A = A, F_B = B, F_S = S, R2 = R2objt$Full)
   } else {
     F_A = R2objt$Full-R2objt$B
     F_B =  R2objt$Full-R2objt$A
@@ -253,9 +306,13 @@ get_shared_anova = function(R2objt, spatial = TRUE, fractional = c("proportional
     A = F_A + F_AB*abs(F_A)/(abs(F_A)+abs(F_B))
     B = F_B + F_AB*abs(F_B)/(abs(F_A)+abs(F_B))
     S = 0
+    proportional = list(F_A = A, F_B = B, F_S = S, R2 = R2objt$Full)
+    A = F_A + F_AB*0.5
+    B = F_B + F_AB*0.5
+    S = 0
+    equal = list(F_A = A, F_B = B, F_S = S, R2 = R2objt$Full)
   }
-  R2 = correct_R2(R2objt$Full)
-  return(list(F_A = A, F_B = B, F_S = S, R2 = R2))
+  return(list(proportional = proportional, equal = equal))
 }
 
 get_null_ll = function(object, ...) {
@@ -293,244 +350,53 @@ get_null_ll = function(object, ...) {
 
 
 
-#' Print sjSDM anova
-#' 
-#' @param x an object of \code{\link{anova.sjSDM}}
-#' @param ... optional arguments for compatibility with the generic function, no function implemented
+#' Summary table of sjSDM anova
 #' 
-#' @return The above matrix is silently returned
+#' The function prints and returns invisible a summary table of an sjSDM ANOVA, created by \code{\link{anova.sjSDM}}
 #' 
+#' @param object an object of \code{\link{anova.sjSDM}}
+#' @param method method used to calculate the ANOVA 
+#' @param fractions how to handle the shared fractions. See details
+#' @param ... optional arguments for compatibility with the generic function, no function implemented
 #' 
-#' @export
-print.sjSDManova = function(x, ...) {
-  cat("Analysis of Deviance Table\n\n")
-  stats::printCoefmat(x$to_print)
-  return(invisible(x$to_print))
-}
-
-#' Plot internal metacommunity structure
-#' 
-#' @param object anova object from \code{\link{anova.sjSDM}}
-#' @param Rsquared which R squared should be used, McFadden or Nagelkerke (McFadden is default)
-#' @param add_shared split shared components, default is TRUE 
-#' @param env_deviance environmental deviance
-#' @param suppress_plotting should the plots be suppressed or not.
-#' 
-#' Plots and returns the internal metacommunity structure of species and sites (see Leibold et al., 2022). 
-#' Plots were heavily inspired by Leibold et al., 2022
-#' 
-#' @return 
-#' 
-#' List with the following components:
-#' 
+#' @details The function returns a ANOVA table with Deviance as well as the pseudo-R2 metrics of Nagelkerke and McFadden
 #' 
-#' \item{plots}{ggplot objects for sites and species.}
-#' \item{data}{List of data.frames with the internal metacommunity structure.}
+#' There are four options to handle shared ANOVA fractions, which is variance that can be explained, typically as a result of collinearity, by several of the fractions:
 #' 
+#' 1. "all" returns the shared fractions explicitly
+#' 2. "discard" discards the fractions, as typically in a type II Anova
+#' 3. "proportional" distributes shared fractions proportional to the unique fractions
+#' 4. "equal" distributions shared fractions equally to the unique fractions
 #' 
-#' @references 
-#' Leibold, M. A., Rudolph, F. J., Blanchet, F. G., De Meester, L., Gravel, D., Hartig, F., ... & Chase, J. M. (2022). The internal structure of metacommunities. Oikos, 2022(1).
+#' @return The matrix that is printed out is silently returned
 #' 
+#' @example /inst/examples/anova-example.R
 #' @export
-plotInternalStructure = function(object,  
-                                 Rsquared = c("McFadden", "Nagelkerke"), 
-                                 add_shared = FALSE,
-                                 env_deviance = NULL,
-                                 suppress_plotting = FALSE) {
-  Rsquared = match.arg(Rsquared)
-  out = 
-    plot.sjSDManova(x = object, 
-                    internal = TRUE, 
-                    add_shared = add_shared,
-                    type = Rsquared, 
-                    alpha = 0.15,
-                    env_deviance = env_deviance,
-                    suppress_plotting = suppress_plotting)
+summary.sjSDManova = function(object, 
+                              method = c("ANOVA"),
+                              fractions = c("all","discard", "proportional", "equal"), ...) {
+  cat("Analysis of Deviance Table\n\n")
+  method = match.arg(method)
+  fractions = match.arg(fractions)
+  
+  out = object$to_print[[fractions]]
+  stats::printCoefmat(out)
   return(invisible(out))
+  
 }
 
 
-get_eigen = function(X, double_center = TRUE, full = FALSE) {
-  D = as.matrix(dist(scale(X)))
-  if(double_center) D = D - (diag(0.0, ncol(D)) + rowMeans(D)) - t(diag(0.0, ncol(D)) + colMeans(D)) + mean(D) # Double center
-  eig = eigen(D)
-  if(!full)return(abs(eig$vectors[,which.max(abs(eig$values))]))
-  else {
-    return(list(x = eig$vectors))
-  }
-}
-
-
-#' Plot Correlations between assembly processes and predictors or traits
-#' @param object An `anova` object from the `anova.sjSDM` function.
-#' @param env Predictor variable. If `NULL`, assembly processes are plotted against environment, spatial uniqueness, and richness.
-#' @param trait Trait variable. Plotted against species R-squared for the three processes.
-#' @param Rsquared Which R-squared should be used: "McFadden" (default) or "Nagelkerke".
-#' @param cols Colors for the three assembly processes.
-#' 
-#' Correlation and plots of the three assembly processes (environment, space, and codist) against environmental and spatial uniqueness and richness. The importance of the three assembly processes is measured by the partial R-squared (shown in the internal structure plots).
-#' 
-#' Importances are available for species and sites. Custom environmental predictors or traits can be specified. Environmental predictors are plotted against site R-squared and traits are plotted against species R-squared.
+#' Print sjSDM anova object
 #' 
-#' Regression lines are estimated by 50% quantile regression models.
-#' 
-#' @return
-#' A list with the following components:
+#' This is a wrapper for \code{\link{summary.sjSDManova}}, maintained for backwards compatibility - prefer to use summary() instead
 #'
-#' \item{env}{A list of summary tables for env, space, and codist R-squared.}
-#' \item{space}{A list of summary tables for env, space, and codist R-squared.}
-#' \item{codist}{A list of summary tables for env, space, and codist R-squared.}
-#' 
-#' @references
-#' Leibold, M. A., Rudolph, F. J., Blanchet, F. G., De Meester, L., Gravel, D., Hartig, F., ... & Chase, J. M. (2022). The internal structure of metacommunities. *Oikos*, 2022(1).
+#' @param x an object of type sjSDManova created by \code{\link{anova.sjSDM}}
+#' @param ... additional arguments to \code{\link{summary.sjSDManova}}
 #' 
+#' @example /inst/examples/anova-example.R
 #' @export
-plotAssemblyEffects = function(object, env = NULL, trait = NULL, Rsquared = c("McFadden", "Nagelkerke"), cols = c("#A38310", "#B42398", "#20A382")) {
-  
-  oldpar = par(no.readonly = TRUE)
-  on.exit(par(oldpar))
-  
-  Rsquared = match.arg(Rsquared)
-  out = 
-    sjSDM:::plot.sjSDManova(x = object, 
-                    internal = TRUE, 
-                    add_shared = add_shared,
-                    type = Rsquared, 
-                    env_deviance = env_deviance,
-                    suppress_plotting = TRUE)
-  
-  
-  lwd = 2
-  
-  X = object$object$settings$env$X
-  XYcoords = object$object$settings$spatial$X
-  Y = object$object$data$Y
-  rr = out
-  out = list()
-  
-  
-  if(is.null(env) & is.null(trait)) {
-  
-    graphics::par(mfrow = c(1, 3), mar = c(4, 4, 4, 1), xaxt= "s")
-    env_eigen = get_eigen(scale(X), FALSE)
-    spatial_eigen = get_eigen(scale(XYcoords),FALSE)
-    richness = rowSums(Y)
-    
-    env_eigen_centered = scale(env_eigen, center = TRUE, scale = FALSE)
-    spatial_eigen_centered = scale(spatial_eigen, center = TRUE, scale = FALSE)
-    richness_centered = scale(richness, center = TRUE, scale = FALSE)
-    
-    
-    out$env = list()
-    out$space = list()
-    out$codist = list()
-    
-    graphics::plot(NULL, NULL, xlim = c(min(env_eigen_centered), max(env_eigen_centered)), 
-         ylim = c(0, 0.6), xlab = "Env uniqueness",main = "", ylab = "R2", las =1)
-    for(i in 1:3) {
-      graphics::points(env_eigen_centered, rr$data$Sites[,i], col = ggplot2::alpha(cols[i], 0.2), pch = 16)
-      g = qgam::qgam( Y ~ env_eigen_centered + spatial_eigen_centered + richness_centered, 
-                data = data.frame(Y = rr$data$Sites[,i], 
-                                  env_eigen_centered = env_eigen_centered, 
-                                  spatial_eigen_centered = spatial_eigen_centered, 
-                                  richness_centered = richness_centered), 
-                qu = 0.5, control = list(progress="none"))
-      out$env[[colnames(rr$data$Sites)[i]]] = g 
-      graphics::abline(a = coef(g)[c(1, 2)], col = cols[i], lwd = lwd, lty =  1*(summary(g)$p.table[2,4] > 0.05)+1 )
-      
-    }
-    graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
-    graphics::legend("topleft", legend = c("significant", "non-significant"), col = c("black", "black"),  bty = "n", lty = c(1,2))
-    
-    graphics::plot(NULL, NULL, xlim = c(min(spatial_eigen_centered), max(spatial_eigen_centered)), ylim = c(0, 0.6), xlab = "Spatial uniqueness",main = "", ylab = "R2", las =1)
-    for(i in 1:3) {
-      graphics::points(spatial_eigen_centered, rr$data$Sites[,i], col = ggplot2::alpha(cols[i], 0.2), pch = 16)
-      g = qgam::qgam( Y ~ env_eigen_centered + spatial_eigen_centered + richness_centered, 
-                data = data.frame(Y = rr$data$Sites[,i], 
-                                  env_eigen_centered = env_eigen_centered, 
-                                  spatial_eigen_centered = spatial_eigen_centered, 
-                                  richness_centered = richness_centered), 
-                qu = 0.5, control = list(progress="none"))
-      out$space[[colnames(rr$data$Sites)[i]]] = g 
-      graphics::abline(a = coef(g)[c(1, 3)], col = cols[i], lwd = lwd, lty =  1*(summary(g)$p.table[3,4] > 0.05)+1 )
-      
-    }
-    graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
-    graphics::legend("topleft", legend = c("significant", "non-significant"), col = c("black", "black"),  bty = "n", lty = c(1,2))
-    
-    
-    graphics::plot(NULL, NULL, xlim = c(min(richness_centered), max(richness_centered)), ylim = c(0, 0.6), xlab = "Richness",main = "", ylab = "R2", las =1)
-    for(i in 1:3) {
-      graphics::points(richness_centered, rr$data$Sites[,i], col = ggplot2::alpha(cols[i], 0.2), pch = 16)
-      g = qgam::qgam( Y ~ env_eigen_centered + spatial_eigen_centered + richness_centered, 
-                data = data.frame(Y = rr$data$Sites[,i], 
-                                  env_eigen_centered = env_eigen_centered, 
-                                  spatial_eigen_centered = spatial_eigen_centered, 
-                                  richness_centered = richness_centered), 
-                qu = 0.5, control = list(progress="none"))
-      out$codist[[colnames(rr$data$Sites)[i]]] = g 
-      graphics::abline(a = coef(g)[c(1, 4)], col = cols[i], lwd = lwd, lty =  1*(summary(g)$p.table[4,4] > 0.05)+1 )
-      
-    }
-    graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
-    graphics::legend("topleft", legend = c("significant", "non-significant"), col = c("black", "black"),  bty = "n", lty = c(1,2))
-  } else {
-    if(!is.null(env)) {
-      if(is.factor(env) || is.character(env)) {
-        group = env
-        df = data.frame(R2 = c(rr$data$Sites$env, rr$data$Sites$spa, rr$data$Sites$codist),
-                        part = rep(c("env", "spa", "codist"), each = nrow(rr$data$Sites)),
-                        group = rep(group, 3)
-        )
-        b = graphics::boxplot(R2~part+group, data =df,las = 2, col = alpha(cols[c(3, 1, 2)],0.7), xlab = "", main = "", notch = TRUE )
-        b = beeswarm::beeswarm(R2~part+group, data =df,las = 2, col = cols[c(3, 1, 2)], xlab = "", main = "" , add=TRUE, method = "center", spacing = 0.3)
-        graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
-      } else {
-        graphics::par(mfrow = c(1, 1), mar = c(4, 4, 4, 1), xaxt= "s")
-        pred = scale(env, center = TRUE, scale = FALSE)
-        
-        graphics::plot(NULL, NULL, xlim = c(min(pred), max(pred)), ylim = c(0, 0.6), xlab = "Predictor",main = "", ylab = "R2", las =1)
-        for(i in 1:3) {
-          graphics::points(pred, rr$data$Sites[,i], col = ggplot2::alpha(cols[i], 0.2), pch = 16)
-          g = qgam::qgam( Y ~ pred, data = data.frame(Y = rr$data$Sites[,i], pred = pred), qu = 0.5, control = list(progress="none"))
-          graphics::abline(a = coef(g)[c(1, 2)], col = cols[i], lwd = lwd, lty =  1*(summary(g)$p.table[2,4] > 0.05)+1 )
-          out$pred[[colnames(rr$data$Sites)[i]]] = g 
-        }
-        graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
-        graphics::legend("topleft", legend = c("significant", "non-significant"), col = c("black", "black"),  bty = "n", lty = c(1,2))
-      }
-    }
-    
-    if(!is.null(trait)) {
-      if(is.factor(trait) || is.character(trait)) {
-        group = trait
-        df = data.frame(R2 = c(rr$data$Species$env, rr$data$Species$spa, rr$data$Species$codist),
-                        part = rep(c("env", "spa", "codist"), each = nrow(rr$data$Species)),
-                        group = rep(group, 3)
-        )
-        b = graphics::boxplot(R2~part+group, data =df,las = 2, col = alpha(cols[c(3, 1, 2)],0.7), xlab = "", main = "", notch = TRUE )
-        b = beeswarm::beeswarm(R2~part+group, data =df,las = 2, col = cols[c(3, 1, 2)], xlab = "", main = "" , add=TRUE, method = "center", spacing = 0.3)
-        graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
-      } else {
-        graphics::par(mfrow = c(1, 1), mar = c(4, 4, 4, 1), xaxt= "s")
-        pred = scale(trait, center = TRUE, scale = FALSE)
-        
-        graphics::plot(NULL, NULL, xlim = c(min(pred), max(pred)), ylim = c(0, 0.6), xlab = "Predictor",main = "", ylab = "R2", las =1)
-        for(i in 1:3) {
-          graphics::points(pred, rr$data$Species[,i], col = ggplot2::alpha(cols[i], 0.2), pch = 16)
-          g = qgam::qgam( Y ~ pred, data = data.frame(Y = rr$data$Species[,i], pred = pred), qu = 0.5, control = list(progress="none"))
-          graphics::abline(a = coef(g)[c(1, 2)], col = cols[i], lwd = lwd, lty =  1*(summary(g)$p.table[2,4] > 0.05)+1 )
-          out$pred[[colnames(rr$data$Species)[i]]] = g 
-        }
-        graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
-        graphics::legend("topleft", legend = c("significant", "non-significant"), col = c("black", "black"),  bty = "n", lty = c(1,2))
-      }
-      
-    }
-    
-
-}
-  
+print.sjSDManova = function(x,...) {
+  out = summary(x, ...)
   return(invisible(out))
 }
 
@@ -541,26 +407,17 @@ plotAssemblyEffects = function(object, env = NULL, trait = NULL, Rsquared = c("M
 #' @param x anova object from \code{\link{anova.sjSDM}}
 #' @param y unused argument
 #' @param type deviance, Nagelkerke or McFadden R-squared
-#' @param internal logical, plot internal or total structure
-#' @param add_shared Add shared contributions when plotting the internal structure
+#' @param fractions how to handle shared fractions
 #' @param cols colors for the groups
 #' @param alpha alpha for colors
 #' @param env_deviance environmental deviance
-#' @param suppress_plotting return plots but don't plot them
 #' @param ... Additional arguments to pass to \code{plot()}
 #' 
-#' The \code{internal = TRUE} plot was heavily inspired by Leibold et al., 2022
 #' 
 #' @return 
 #' 
 #' List with the following components:
 #' 
-#' If \code{internal=TRUE}:
-#' 
-#' \item{plots}{ggplot objects for sites and species.}
-#' \item{data}{List of data.frames with the shown results.}
-#' 
-#' else:
 #' \item{VENN}{Matrix of shown results.}
 #' 
 #' @references 
@@ -570,13 +427,14 @@ plotAssemblyEffects = function(object, env = NULL, trait = NULL, Rsquared = c("M
 plot.sjSDManova = function(x, 
                            y, 
                            type = c( "McFadden", "Deviance", "Nagelkerke"), 
-                           internal = FALSE,
-                           add_shared = FALSE,
+                           fractions = c("discard", "proportional", "equal"),
                            cols = c("#7FC97F","#BEAED4","#FDC086"),
                            alpha=0.15, 
                            env_deviance = NULL,
-                           suppress_plotting = FALSE,
                            ...) {
+  
+  fractions = match.arg(fractions)
+  
   lineSeq = 0.3
   nseg = 100
   dr = 1.0
@@ -586,169 +444,49 @@ plot.sjSDManova = function(x,
   oldpar = par(no.readonly = TRUE)
   on.exit(par(oldpar))
   
-  if(!x$spatial && internal) {
-    internal=FALSE
-    warning("'internal=TRUE' currently only supported for spatial models.")
-  }
-  
-  if(!internal) {
-    
-    values = x$results
-    values$`R2 Nagelkerke` = ifelse(values$`R2 Nagelkerke`<0, 0, values$`R2 Nagelkerke`)
-    values$`R2 McFadden` = ifelse(values$`R2 McFadden`<0, 0, values$`R2 McFadden`)
-    select_rows = 
-      if(x$spatial) { 
-        sapply(c("F_A", "F_B", "F_AB","F_S", "F_AS", "F_BS", "F_ABS"), function(i) which(values$models == i, arr.ind = TRUE))
-      } else {
-        sapply(c("F_A", "F_B", "F_AB"), function(i) which(values$models == i, arr.ind = TRUE))
-      }
-    
-    values = values[select_rows,]
-    col_index = 
-      switch (type,
-        Deviance = 4,
-        Nagelkerke = 5,
-        McFadden = 6
-      )
-    
-    
-    graphics::plot(NULL, NULL, xlim = c(0,1), ylim =c(0,1),pty="s", axes = FALSE, xlab = "", ylab = "")
-    xx = 1.1*lineSeq*cos( seq(0,2*pi, length.out=nseg))
-    yy = 1.1*lineSeq*sin( seq(0,2*pi, length.out=nseg))
-    graphics::polygon(xx+lineSeq,yy+(1-lineSeq), col= addA(cols[1],alpha = alpha), border = "black", lty = 1, lwd = 1)
-    graphics::text(lineSeq-0.1, (1-lineSeq),labels = round(values[1,col_index],3))
-    graphics::text(mean(xx+lineSeq), 0.9,labels = "Environmental", pos = 3)
-    
-    graphics::polygon(xx+1-lineSeq,yy+1-lineSeq, col= addA(cols[2],alpha = alpha), border = "black", lty = 1, lwd = 1)
-    graphics::text(1-lineSeq+0.1, (1-lineSeq),labels = round(values[2,col_index],3))
-    graphics::text(1-mean(xx+lineSeq), 0.9,labels = "Associations", pos = 3)
-    graphics::text(0.5, (1-lineSeq),labels = round(values[3,col_index],3))
-    
-    if(x$spatial) {
-      graphics::polygon(xx+0.5,yy+lineSeq, col= addA(cols[3],alpha = alpha), border = "black", lty = 1, lwd = 1)
-      graphics::text(0.5, lineSeq+0.0,pos = 1,labels = round(values[4,col_index],3))
-      graphics::text(0.5, 0.1,labels = "Spatial", pos = 1)
-      graphics::text(0.3, 0.5,pos=1,labels   = round(values[5,col_index],3)) # AS
-      graphics::text(1-0.3, 0.5,pos=1,labels = round(values[6,col_index],3)) # BS
-      graphics::text(0.5, 0.5+0.05,labels    = round(values[7,col_index],3)) # ABS
-    }
-    out$VENN = values
-  } else { 
-  
-    if(type == "Deviance") {type = "R2_McFadden"
+  values = x$results
+
+  select_rows = 
+    if(x$spatial) { 
+      sapply(c("F_A", "F_B", "F_AB","F_S", "F_AS", "F_BS", "F_ABS"), function(i) which(values$models == i, arr.ind = TRUE))
     } else {
-      if(type == "McFadden") type = "R2_McFadden"
-      else type = "R2_Nagelkerke"
+      sapply(c("F_A", "F_B", "F_AB"), function(i) which(values$models == i, arr.ind = TRUE))
     }
-    internals = list()
-    
-    if(!add_shared) {
-      df = data.frame(
-          env = ifelse(x$sites[[type]]$F_A<0, 0, x$sites[[type]]$F_A),
-          spa = ifelse(x$sites[[type]]$F_S<0, 0, x$sites[[type]]$F_S),
-          codist = ifelse(x$sites[[type]]$F_B<0, 0, x$sites[[type]]$F_B),
-          r2  = ifelse(x$sites[[type]]$Full<0, 0, x$sites[[type]]$Full)#/length(x$sites[[type]]$Full)
-        )
-      internals[[1]] = df
-      names(internals)[1] = "Sites"
-      
-      df = data.frame(
-          env = ifelse(x$species[[type]]$F_A<0, 0, x$species[[type]]$F_A),
-          spa = ifelse(x$species[[type]]$F_S<0, 0, x$species[[type]]$F_S),
-          codist = ifelse(x$species[[type]]$F_B<0, 0, x$species[[type]]$F_B),
-          r2  = ifelse(x$species[[type]]$Full<0, 0, x$species[[type]]$Full)#/length(x$species[[type]]$Full)
-        )
-      
-      internals[[2]] = df
-      names(internals)[2] = "Species"
-    } else {
-      type = paste0(type, "_shared")
-      df = data.frame(
-        env = ifelse(x$sites[[type]]$F_A<0, 0, x$sites[[type]]$F_A),
-        spa = ifelse(x$sites[[type]]$F_S<0, 0, x$sites[[type]]$F_S),
-        codist = ifelse(x$sites[[type]]$F_B<0, 0, x$sites[[type]]$F_B),
-        r2  = ifelse(x$sites[[type]]$R2<0, 0, x$sites[[type]]$R2)#/length(x$sites[[type]]$R2)
-      )
-      internals[[1]] = df
-      names(internals)[1] = "Sites"
-      
-      df = data.frame(
-        env = ifelse(x$species[[type]]$F_A<0, 0, x$species[[type]]$F_A),
-        spa = ifelse(x$species[[type]]$F_S<0, 0, x$species[[type]]$F_S),
-        codist = ifelse(x$species[[type]]$F_B<0, 0, x$species[[type]]$F_B),
-        r2  = ifelse(x$species[[type]]$R2<0, 0, x$species[[type]]$R2)#/length(x$species[[type]]$R2)
-      )
-      
-      internals[[2]] = df
-      names(internals)[2] = "Species"
-    }
-    
-    
-    plots_internals = list()
-    
-    # Code taken from https://github.com/javirudolph/iStructureMetaco/blob/master/InternalStructureMetacommunities_2021_manuscript/Figures.R
-    for(i in 1:length(internals)) {
-        
-      add_grad = FALSE
-      if((i == 1) & !is.null(env_deviance)) add_grad = TRUE
-      
-      top = 7
-      if(i > 1) top = 1
-      if(is.null(env_deviance)) top = 1
-        
-        r2max = ceiling(max(internals[[i]]$r2)*1e2)/1e2
-        plt = 
-          ggtern::ggtern(internals[[i]], ggplot2::aes_string(x = "env", z = "spa", y = "codist", size = "r2"))+
-            ggtern::scale_T_continuous(limits=c(0,1),
-                                       breaks=seq(0, 1,by=0.2),
-                                       labels=seq(0,1, by= 0.2)) +
-            ggtern::scale_L_continuous(limits=c(0,1),
-                                       breaks=seq(0, 1,by=0.2),
-                                       labels=seq(0, 1,by=0.2)) +
-            ggtern::scale_R_continuous(limits=c(0,1),
-                                       breaks=seq(0, 1,by=0.2),
-                                       labels=seq(0, 1,by=0.2)) +
-            #ggplot2::scale_size_area( max_size = 3) +
-            ggplot2::labs(title = names(internals)[i],
-                          x = "E",
-                          xarrow = "Environment",
-                          y = "C",
-                          yarrow = "Species associations",
-                          z = "S", 
-                          zarrow = "Space") +
-            ggtern::theme_bw() +
-            ggtern::theme_showarrows() +
-            ggtern::theme_arrowlong() +
-            ggplot2::theme(
-              panel.grid = ggplot2::element_line(color = "darkgrey", size = 0.3),
-              plot.tag = ggplot2::element_text(size = 11),
-              plot.title = ggplot2::element_text(size = 11, hjust = 0.1 , margin = ggplot2::margin(t = 10, b = -20)),
-              tern.axis.arrow = ggplot2::element_line(size = 1),
-              tern.axis.arrow.text = ggplot2::element_text(size = 6),
-              axis.text = ggplot2::element_text(size = 4),
-              axis.title = ggplot2::element_text(size = 6),
-              legend.text = ggplot2::element_text(size = 6),
-              legend.title = ggplot2::element_text(size = 8),
-              strip.text = ggplot2::element_text(size = 8),
-              plot.margin = unit(c(top,1,1,1)*0.2, "cm"),
-            strip.background = ggplot2::element_rect(color = NA),
-          ) +
-          ggplot2::guides(size = ggplot2::guide_legend(title = expression(R^2), order = 1)) +
-          { if(!add_grad)ggplot2::geom_point(alpha = 0.7) }+
-          { if(add_grad) ggplot2::geom_point(alpha = 0.7, aes(fill=env_deviance, color = env_deviance)) }+  
-          ggplot2::scale_size_continuous(range = c(0.1,5),limits = c(0, r2max), breaks = seq(0, r2max, length.out=5)) +
-          { if(add_grad) ggplot2::scale_fill_gradient(low = "white", high = "black", guide = "none") } + 
-          { if(add_grad) ggplot2::scale_color_gradient(low = "white", high = "black", limits = c(0, max(env_deviance))) } +
-          ggplot2::theme(tern.axis.arrow.text = element_text(size = 7),legend.position = "bottom", legend.margin = margin(r = 30), legend.box="vertical") +
-          { if(!add_grad) ggplot2::guides(size = ggplot2::guide_legend(title = expression(R^2), order = 1, nrow = 1, label.position = "bottom")) } +
-          { if( add_grad) ggplot2::guides(size = ggplot2::guide_legend(title = expression(R^2), order = 1, nrow = 1, label.position = "bottom"),
-                                          color = ggplot2::guide_colorbar(title = "Environmental deviation", title.position = "top", order = 2, barheight = 0.5, barwidth = 8)) } 
-        plots_internals[[i]] = plt
-      }
-    if(!suppress_plotting) ggtern::grid.arrange(plots_internals[[1]], plots_internals[[2]], nrow=1, widths = c(5.0/10, 5/10))
-    out$plots = plots_internals
-    out$data = internals
+  
+  values = values[select_rows,]
+  col_index = 
+    switch (type,
+            Deviance = 4,
+            Nagelkerke = 5,
+            McFadden = 6
+    )
+  
+  
+  graphics::plot(NULL, NULL, xlim = c(0,1), ylim =c(0,1),pty="s", axes = FALSE, xlab = "", ylab = "")
+  xx = 1.1*lineSeq*cos( seq(0,2*pi, length.out=nseg))
+  yy = 1.1*lineSeq*sin( seq(0,2*pi, length.out=nseg))
+  graphics::polygon(xx+lineSeq,yy+(1-lineSeq), col= addA(cols[1],alpha = alpha), border = "black", lty = 1, lwd = 1)
+  graphics::text(lineSeq-0.1, (1-lineSeq),labels = round(values[1,col_index],3))
+  graphics::text(mean(xx+lineSeq), 0.9,labels = "Environmental", pos = 3)
+  
+  graphics::polygon(xx+1-lineSeq,yy+1-lineSeq, col= addA(cols[2],alpha = alpha), border = "black", lty = 1, lwd = 1)
+  graphics::text(1-lineSeq+0.1, (1-lineSeq),labels = round(values[2,col_index],3))
+  graphics::text(1-mean(xx+lineSeq), 0.9,labels = "Associations", pos = 3)
+  graphics::text(0.5, (1-lineSeq),labels = round(values[3,col_index],3))
+  
+  if(x$spatial) {
+    graphics::polygon(xx+0.5,yy+lineSeq, col= addA(cols[3],alpha = alpha), border = "black", lty = 1, lwd = 1)
+    graphics::text(0.5, lineSeq+0.0,pos = 1,labels = round(values[4,col_index],3))
+    graphics::text(0.5, 0.1,labels = "Spatial", pos = 1)
+    graphics::text(0.3, 0.5,pos=1,labels   = round(values[5,col_index],3)) # AS
+    graphics::text(1-0.3, 0.5,pos=1,labels = round(values[6,col_index],3)) # BS
+    graphics::text(0.5, 0.5+0.05,labels    = round(values[7,col_index],3)) # ABS
   }
+  out$VENN = values
   return(invisible(out))
 }
 
+
+
+
+
diff --git a/sjSDM/R/importance.R b/sjSDM/R/importance.R
index 8fa7d81..bb02fd8 100644
--- a/sjSDM/R/importance.R
+++ b/sjSDM/R/importance.R
@@ -1,4 +1,4 @@
-#' importance 
+#' Importance of environmental, spatial and association components
 #' 
 #' Computes standardized variance components with respect to abiotic, biotic, and spatial effect groups. 
 #' 
@@ -6,10 +6,11 @@
 #' @param save_memory use torch backend to calculate importance with single precision floats
 #' @param ... additional arguments
 #' 
-#' @details 
+#' @details This approach is based on Ovaskainen et al., 2017, and also used in  Leibold et al., 2021. Unlike the \code{\link{anova.sjSDM}} function in the sjSDM package, importance is not calculated by explicitly switching a particular model component of and refitting the model, but essentially by setting it ineffective. 
 #' 
-#' This variance partitioning approach is based on Ovaskainen et al., 2017. For an example how to interpret the outputs, see Leibold et al., 2021.
-#' This function will be deprecated in the future. Please use \code{plot(anova(model), internal=TRUE)} (currently only supported for spatial models).
+#' Although we have no hard reasons to discourage the use of this function, we have decided in sjSDM to measure importance maninly based on a traditional ANOVA approach. We therefore recommend users to use the \code{\link{anova.sjSDM}}.
+#' 
+#' This function is maintained hidden for comparison / benchmarking purpose, and in case there is a need to use it in the future. If you want to access it, use sjSDM:::importance. 
 #' 
 #' @return
 #' 
@@ -30,7 +31,6 @@
 #' @seealso \code{\link{print.sjSDMimportance}}, \code{\link{plot.sjSDMimportance}}
 #' @example /inst/examples/importance-example.R
 #' @author Maximilian Pichler
-#' @export
 importance = function(x, save_memory = TRUE, ...) {
   model = x
   stopifnot(
@@ -190,7 +190,6 @@ plot.sjSDMimportance= function(x, y, col.points="#24526e",cex.points=1.2, ...) {
 #' @param covSP spatial covariance matrix
 #' 
 #' @author Maximilian Pichler
-
 getImportance = function(beta, sp=NULL, association, covX, covSP=NULL) {
   nsp = ncol(beta)
   nGroups = nrow(beta)
diff --git a/sjSDM/R/installation_help.R b/sjSDM/R/installation_help.R
index 98a28f1..ef3cf22 100644
--- a/sjSDM/R/installation_help.R
+++ b/sjSDM/R/installation_help.R
@@ -1,6 +1,5 @@
 #' @title Installation help
 #' @name installation_help
-#' @docType package
 #' @description Trouble shooting guide for the installation of the sjSDM package
 #' 
 #' We provide a function \code{\link{install_sjSDM}} to install automatically 
@@ -148,4 +147,4 @@
 #' \href{https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example/}{reproducible example} 
 #' via the sjSDM \href{https://github.com/TheoreticalEcology/s-jSDM/issues/}{issue tracker} 
 #' with a copy of the \code{\link{install_diagnostic}} output as a quote. 
-NULL
\ No newline at end of file
+"_PACKAGE"
\ No newline at end of file
diff --git a/sjSDM/R/internalStructure.R b/sjSDM/R/internalStructure.R
new file mode 100644
index 0000000..fc1c329
--- /dev/null
+++ b/sjSDM/R/internalStructure.R
@@ -0,0 +1,459 @@
+
+
+#' Plot internal metacommunity structure
+#' 
+#' @param object anova object from \code{\link{anova.sjSDM}}
+#' @param Rsquared which R squared should be used, McFadden or Nagelkerke (McFadden is default)
+#' @param fractions how to handle shared fractions
+#' @param negatives how to handle negative R squareds
+#' @param plot should the plots be suppressed or not.
+#' 
+#' Plots and returns the internal metacommunity structure of species and sites (see Leibold et al., 2022). 
+#' Plots were heavily inspired by Leibold et al., 2022
+#' 
+#' @return 
+#' 
+#' List with the following components:
+#' 
+#' 
+#' \item{plots}{ggplot objects for sites and species.}
+#' \item{data}{List of data.frames with the internal metacommunity structure.}
+#' 
+#' @example /inst/examples/anova-example.R
+#' @references 
+#' Leibold, M. A., Rudolph, F. J., Blanchet, F. G., De Meester, L., Gravel, D., Hartig, F., ... & Chase, J. M. (2022). The internal structure of metacommunities. Oikos, 2022(1).
+#' 
+#' @export
+internalStructure = function(object,  
+                             Rsquared = c("McFadden", "Nagelkerke"), 
+                             fractions = c("discard", "proportional", "equal"),
+                             negatives = c("floor", "scale"), # TODO - rounding ANOVA out, here all calculations to function with option
+                             plot = FALSE) {
+  
+  fractions = match.arg(fractions)
+  Rsquared = match.arg(Rsquared)
+  negatives = match.arg(negatives)
+  
+  if(!object$spatial) stop("'internal structure' currently only supported for spatial models.")  
+  
+  if(Rsquared == "Deviance") {type = "R2_McFadden"
+  } else {
+    if(Rsquared == "McFadden") type = "R2_McFadden"
+    else type = "R2_Nagelkerke"
+  }
+  
+  internals = list()
+  
+  if(fractions == "discard") {
+    df = data.frame(
+      env = object$sites[[type]]$F_A,
+      spa = object$sites[[type]]$F_S,
+      codist = object$sites[[type]]$F_B,
+      r2  = object$sites[[type]]$Full #/length(object$sites[[type]]$Ful)
+    )
+    
+    
+    internals[[1]] = df
+    names(internals)[1] = "Sites"
+    
+    df = data.frame(
+      env = object$species[[type]]$F_A,
+      spa = object$species[[type]]$F_S,
+      codist = object$species[[type]]$F_B,
+      r2  = object$species[[type]]$Full
+    )
+    
+    internals[[2]] = df
+    names(internals)[2] = "Species"
+  } else {
+    type = paste0(type, "_shared")
+    if(fractions == "proportional") {
+      
+      df = data.frame(
+        env = object$sites[[type]]$proportional$F_A,
+        spa = object$sites[[type]]$proportional$F_S,
+        codist = object$sites[[type]]$proportional$F_B,
+        r2  = object$sites[[type]]$proportional$R2 #/length(object$sites[[type]]$R2)
+      )
+      
+      internals[[1]] = df
+      names(internals)[1] = "Sites"
+      
+      df = data.frame(
+        env = object$species[[type]]$proportional$F_A,
+        spa = object$species[[type]]$proportional$F_S,
+        codist = object$species[[type]]$proportional$F_B,
+        r2  = object$species[[type]]$proportional$R2 #/length(object$species[[type]]$R2)
+      )
+      
+      
+      internals[[2]] = df
+      names(internals)[2] = "Species"
+      
+    } else {
+      
+      df = data.frame(
+        env = object$sites[[type]]$equal$F_A ,
+        spa = object$sites[[type]]$equal$F_S ,
+        codist = object$sites[[type]]$equal$F_B ,
+        r2  = object$sites[[type]]$equal$R2 #/length(object$sites[[type]]$R2)
+      )
+      
+
+      internals[[1]] = df
+      names(internals)[1] = "Sites"
+      
+      df = data.frame(
+        env = object$species[[type]]$equal$F_A,
+        spa = object$species[[type]]$equal$F_S,
+        codist = object$species[[type]]$equal$F_B,
+        r2  = object$species[[type]]$equal$R2 #/length(object$species[[type]]$R2)
+      )
+      
+      
+      internals[[2]] = df
+      names(internals)[2] = "Species"
+    }
+      
+    }
+    
+  internals[[1]] = standardize_df(internals[[1]], standardize = negatives )
+  internals[[2]] = standardize_df(internals[[2]], standardize = negatives )
+  
+  out = list()
+  out$internals = internals
+  out$Rsquared = Rsquared
+  out$fractions = fractions
+  out$anova = object
+  class(out) = "sjSDMinternalStruture"
+  
+  if(plot == T) plot(out)
+  
+  return(out)
+}
+
+
+standardize_df = function(df, standardize) {
+  if(standardize == "floor") {
+    tmp = df[,1:3]
+    tmp[tmp<0] = 0
+    df[,1:3] = tmp
+  } else if (standardize == "abs") {
+    df[,1:3] = abs(df[,1:3])
+  } else if(standardize == "scale" ){
+    tmp = df[,1:3]
+    tmp = scales::rescale(as.matrix(tmp), to = c(0, 1))
+    df[,1:3] = tmp[,1:3]
+  }
+  return(df)
+}
+
+
+#' Print internal structure object
+#' 
+#' @param x object of class sjSDMinternalStruture
+#' @param ... no function
+#' 
+#' @export
+print.sjSDMinternalStruture <- function(x, ...){
+  return(x$internals)
+}
+
+
+#' Plot internal structure
+#' 
+#' Creates a ternary diagram of an object of class 
+#' 
+#' @param x and object of class sjSDMinternalStruture create by anova object from \code{\link{internalStructure}}
+#' @param alpha alpha of points
+#' @param env_deviance environmental deviance/gradient (points will be colored)
+#' @param ... no function
+#' 
+#' 
+#' @example /inst/examples/anova-example.R
+#' @export
+#' 
+plot.sjSDMinternalStruture <- function(x, 
+                                       alpha = 0.15,
+                                       env_deviance = NULL,
+                                       ...){
+  
+  internals = x$internals
+  
+  plots_internals = list()
+  
+  # Code taken from https://github.com/javirudolph/iStructureMetaco/blob/master/InternalStructureMetacommunities_2021_manuscript/Figures.R
+  
+  
+  if(min(internals[[1]][,1:3]) < 0 | min(internals[[2]][,1:3]) < 0) {
+    message("Negative partial R-square detected. Negative R-squareds can occur, but the ternary plots are not able to display them")
+  }
+  
+  for(i in 1:length(internals)) {
+    
+    add_grad = FALSE
+    if((i == 1) & !is.null(env_deviance)) add_grad = TRUE
+    
+    top = 7
+    if(i > 1) top = 1
+    if(is.null(env_deviance)) top = 1
+    
+    negative_r2 = FALSE
+    if(min(internals[[i]]$r2) < 0) negative_r2 = TRUE
+    
+    r2max = ceiling(max(internals[[i]]$r2)*1e2)/1e2
+    r2min = floor(min(internals[[i]]$r2)*1e2)/1e2
+    
+    color = if(!negative_r2) {NULL} else {"r2"}
+    
+    plt = 
+      ggtern::ggtern(internals[[i]], ggplot2::aes_string(x = "env", z = "spa", y = "codist", size = abs(internals[[i]]$r2), color = color) ) +
+      ggtern::scale_T_continuous(limits=c(0,1),
+                                 breaks=seq(0, 1,by=0.2),
+                                 labels=seq(0,1, by= 0.2)) +
+      ggtern::scale_L_continuous(limits=c(0,1),
+                                 breaks=seq(0, 1,by=0.2),
+                                 labels=seq(0, 1,by=0.2)) +
+      ggtern::scale_R_continuous(limits=c(0,1),
+                                 breaks=seq(0, 1,by=0.2),
+                                 labels=seq(0, 1,by=0.2)) +
+      #ggplot2::scale_size_area( max_size = 3) +
+      ggplot2::labs(title = names(internals)[i],
+                    x = "E",
+                    xarrow = "Environment",
+                    y = "C",
+                    yarrow = "Species associations",
+                    z = "S", 
+                    zarrow = "Space") +
+      ggtern::theme_bw() +
+      ggtern::theme_showarrows() +
+      ggtern::theme_arrowlong() +
+      ggplot2::theme(
+        panel.grid = ggplot2::element_line(color = "darkgrey", size = 0.3),
+        plot.tag = ggplot2::element_text(size = 11),
+        plot.title = ggplot2::element_text(size = 11, hjust = 0.1 , margin = ggplot2::margin(t = 10, b = -20)),
+        tern.axis.arrow = ggplot2::element_line(size = 1),
+        tern.axis.arrow.text = ggplot2::element_text(size = 6),
+        axis.text = ggplot2::element_text(size = 4),
+        axis.title = ggplot2::element_text(size = 6),
+        legend.text = ggplot2::element_text(size = 6),
+        legend.title = ggplot2::element_text(size = 8),
+        strip.text = ggplot2::element_text(size = 8),
+        plot.margin = unit(c(top,1,1,1)*0.2, "cm"),
+        strip.background = ggplot2::element_rect(color = NA),
+      ) +
+      { if(!negative_r2) ggplot2::guides(size = ggplot2::guide_legend(title = expression(R^2), order = 1)) } +
+      { if(!add_grad)ggplot2::geom_point(alpha = 0.7) }+
+      { if(add_grad) ggplot2::geom_point(alpha = 0.7, aes(fill=env_deviance, color = env_deviance)) }+  
+      { if(!negative_r2)  ggplot2::scale_size_continuous(range = c(0.1,5),limits = c(r2min, r2max), breaks = seq(r2min, r2max, length.out=5)) } +
+      { if(negative_r2)   ggplot2::scale_color_gradient2(low = "red", mid = "grey50", high = "black", midpoint = 0, 
+                                                  breaks = c(r2min, 0, r2max), 
+                                                  labels = c(r2min, "0", r2max), 
+                                                  limits = c(r2min, r2max),
+                                                  guide = ggplot2::guide_colorbar(title = expression(R^2))) } +
+      { if(negative_r2)  ggplot2::scale_size_continuous(range = c(0.1,5), 
+                                               breaks =  seq(0, r2max, length.out = 5), guide = "none")  } +
+      
+      { if(add_grad) ggplot2::scale_fill_gradient(low = "white", high = "black", guide = "none") } + 
+      { if(add_grad) ggplot2::scale_color_gradient(low = "white", high = "black", limits = c(0, max(env_deviance))) } +
+      ggplot2::theme(tern.axis.arrow.text = element_text(size = 7),legend.position = "bottom", legend.margin = margin(r = 30), legend.box="vertical") +
+      { if(!add_grad) { if(!negative_r2)  ggplot2::guides(size = ggplot2::guide_legend(title = expression(R^2), order = 1, nrow = 1, label.position = "bottom")) } } +
+      { if( add_grad) ggplot2::guides(size = ggplot2::guide_legend(title = expression(R^2), order = 1, nrow = 1, label.position = "bottom"),
+                                      color = ggplot2::guide_colorbar(title = "Environmental deviation", title.position = "top", order = 2, barheight = 0.5, barwidth = 8)) } 
+    plots_internals[[i]] = plt
+  }
+  
+  ggtern::grid.arrange(plots_internals[[1]], plots_internals[[2]], nrow=1, widths = c(5.0/10, 5/10))
+  
+  out = list()
+  out$plots = plots_internals
+  return(invisible(out))
+}
+
+#' Plot correlations between assembly processes and predictors or traits
+#' 
+#' @param object An \code{sjSDManova} object from the \code{\link{anova.sjSDM}} function.
+#' @param env Predictor variable. If \code{NULL}, assembly processes are plotted against environment, spatial uniqueness, and richness.
+#' @param trait Trait variable. Plotted against species R-squared for the three processes.
+#' @param cols Colors for the three assembly processes.
+#' 
+#' Correlation and plots of the three assembly processes (environment, space, and codist) against environmental and spatial uniqueness and richness. The importance of the three assembly processes is measured by the partial R-squared (shown in the internal structure plots).
+#' Importances are available for species and sites. Custom environmental predictors or traits can be specified. Environmental predictors are plotted against site R-squared and traits are plotted against species R-squared.
+#' Regression lines are estimated by 50\% quantile regression models.
+#' 
+#' @return
+#' 
+#' A list with the following components:
+#'
+#' \item{env}{A list of summary tables for env, space, and codist R-squared.}
+#' \item{space}{A list of summary tables for env, space, and codist R-squared.}
+#' \item{codist}{A list of summary tables for env, space, and codist R-squared.}
+#' 
+#' @references
+#' 
+#' Leibold, M. A., Rudolph, F. J., Blanchet, F. G., De Meester, L., Gravel, D., Hartig, F., ... & Chase, J. M. (2022). The internal structure of metacommunities. *Oikos*, 2022(1).
+#' 
+#' @example /inst/examples/anova-example.R
+#' @export
+plotAssemblyEffects = function(object, 
+                               env = NULL, 
+                               trait = NULL, 
+                               cols = c("#A38310", "#B42398", "#20A382")) {
+  
+  oldpar = par(no.readonly = TRUE)
+  on.exit(par(oldpar))
+
+  lwd = 2
+  
+  X = object$anova$object$settings$env$X
+  XYcoords = object$anova$object$settings$spatial$X
+  Y = object$anova$object$data$Y
+  rr = object
+  minR = min(rr$internals$Sites[,1:3])
+  maxR = max(rr$internals$Sites[,1:3])
+  minRS = min(rr$internals$Species[,1:3])
+  maxRS = max(rr$internals$Species[,1:3])
+  
+  out = list()
+  
+  if(is.null(env) & is.null(trait)) {
+    
+    graphics::par(mfrow = c(1, 3), mar = c(4, 4, 4, 1), xaxt= "s")
+    env_eigen = get_eigen(scale(X), FALSE)
+    spatial_eigen = get_eigen(scale(XYcoords),FALSE)
+    richness = rowSums(Y)
+    
+    env_eigen_centered = scale(env_eigen, center = TRUE, scale = FALSE)
+    spatial_eigen_centered = scale(spatial_eigen, center = TRUE, scale = FALSE)
+    richness_centered = scale(richness, center = TRUE, scale = FALSE)
+    
+    
+    out$env = list()
+    out$space = list()
+    out$codist = list()
+    
+
+    
+    graphics::plot(NULL, NULL, xlim = c(min(env_eigen_centered), max(env_eigen_centered)), 
+                   ylim = c(minR, maxR*1.5), xlab = "Env uniqueness",main = "", ylab = "R2", las =1)
+    
+    for(i in 1:3) {
+      graphics::points(env_eigen_centered, rr$internals$Sites[,i], col = ggplot2::alpha(cols[i], 0.2), pch = 16)
+      g = qgam::qgam( Y ~ env_eigen_centered + spatial_eigen_centered + richness_centered, 
+                      data = data.frame(Y = rr$internals$Sites[,i], 
+                                        env_eigen_centered = env_eigen_centered, 
+                                        spatial_eigen_centered = spatial_eigen_centered, 
+                                        richness_centered = richness_centered), 
+                      qu = 0.5, control = list(progress="none"))
+      out$env[[colnames(rr$internals$Sites)[i]]] = g 
+      graphics::abline(a = coef(g)[c(1, 2)], col = cols[i], lwd = lwd, lty =  1*(summary(g)$p.table[2,4] > 0.05)+1 )
+      
+    }
+    graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
+    graphics::legend("topleft", legend = c("significant", "non-significant"), col = c("black", "black"),  bty = "n", lty = c(1,2))
+    
+    graphics::plot(NULL, NULL, xlim = c(min(spatial_eigen_centered), max(spatial_eigen_centered)), ylim = c(minR, maxR*1.5), xlab = "Spatial uniqueness",main = "", ylab = "R2", las =1)
+    for(i in 1:3) {
+      graphics::points(spatial_eigen_centered, rr$internals$Sites[,i], col = ggplot2::alpha(cols[i], 0.2), pch = 16)
+      g = qgam::qgam( Y ~ env_eigen_centered + spatial_eigen_centered + richness_centered, 
+                      data = data.frame(Y = rr$internals$Sites[,i], 
+                                        env_eigen_centered = env_eigen_centered, 
+                                        spatial_eigen_centered = spatial_eigen_centered, 
+                                        richness_centered = richness_centered), 
+                      qu = 0.5, control = list(progress="none"))
+      out$space[[colnames(rr$internals$Sites)[i]]] = g 
+      graphics::abline(a = coef(g)[c(1, 3)], col = cols[i], lwd = lwd, lty =  1*(summary(g)$p.table[3,4] > 0.05)+1 )
+      
+    }
+    graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
+    graphics::legend("topleft", legend = c("significant", "non-significant"), col = c("black", "black"),  bty = "n", lty = c(1,2))
+    
+    
+    graphics::plot(NULL, NULL, xlim = c(min(richness_centered), max(richness_centered)), ylim = c(minR, maxR*1.5), xlab = "Richness",main = "", ylab = "R2", las =1)
+    for(i in 1:3) {
+      graphics::points(richness_centered, rr$internals$Sites[,i], col = ggplot2::alpha(cols[i], 0.2), pch = 16)
+      g = qgam::qgam( Y ~ env_eigen_centered + spatial_eigen_centered + richness_centered, 
+                      data = data.frame(Y = rr$internals$Sites[,i], 
+                                        env_eigen_centered = env_eigen_centered, 
+                                        spatial_eigen_centered = spatial_eigen_centered, 
+                                        richness_centered = richness_centered), 
+                      qu = 0.5, control = list(progress="none"))
+      out$codist[[colnames(rr$internals$Sites)[i]]] = g 
+      graphics::abline(a = coef(g)[c(1, 4)], col = cols[i], lwd = lwd, lty =  1*(summary(g)$p.table[4,4] > 0.05)+1 )
+      
+    }
+    graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
+    graphics::legend("topleft", legend = c("significant", "non-significant"), col = c("black", "black"),  bty = "n", lty = c(1,2))
+  } else {
+    if(!is.null(env)) {
+      if(is.factor(env) || is.character(env)) {
+        group = env
+        df = data.frame(R2 = c(rr$internals$Sites$env, rr$internals$Sites$spa, rr$internals$Sites$codist),
+                        part = rep(c("env", "spa", "codist"), each = nrow(rr$internals$Sites)),
+                        group = rep(group, 3)
+        )
+        b = graphics::boxplot(R2~part+group, data =df,las = 2, col = alpha(cols[c(3, 1, 2)],0.7), xlab = "", main = "", notch = TRUE )
+        b = beeswarm::beeswarm(R2~part+group, data =df,las = 2, col = cols[c(3, 1, 2)], xlab = "", main = "" , add=TRUE, method = "center", spacing = 0.3)
+        graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
+      } else {
+        graphics::par(mfrow = c(1, 1), mar = c(4, 4, 4, 1), xaxt= "s")
+        pred = scale(env, center = TRUE, scale = FALSE)
+        
+        graphics::plot(NULL, NULL, xlim = c(min(pred), max(pred)), ylim = c(minR, maxR*1.5), xlab = "Predictor",main = "", ylab = "R2", las =1)
+        for(i in 1:3) {
+          graphics::points(pred, rr$internals$Sites[,i], col = ggplot2::alpha(cols[i], 0.2), pch = 16)
+          g = qgam::qgam( Y ~ pred, data = data.frame(Y = rr$internals$Sites[,i], pred = pred), qu = 0.5, control = list(progress="none"))
+          graphics::abline(a = coef(g)[c(1, 2)], col = cols[i], lwd = lwd, lty =  1*(summary(g)$p.table[2,4] > 0.05)+1 )
+          out$pred[[colnames(rr$internals$Sites)[i]]] = g 
+        }
+        graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
+        graphics::legend("topleft", legend = c("significant", "non-significant"), col = c("black", "black"),  bty = "n", lty = c(1,2))
+      }
+    }
+    
+    if(!is.null(trait)) {
+      if(is.factor(trait) || is.character(trait)) {
+        group = trait
+        df = data.frame(R2 = c(rr$internals$Species$env, rr$internals$Species$spa, rr$internals$Species$codist),
+                        part = rep(c("env", "spa", "codist"), each = nrow(rr$internals$Species)),
+                        group = rep(group, 3)
+        )
+        b = graphics::boxplot(R2~part+group, data =df,las = 2, col = alpha(cols[c(3, 1, 2)],0.7), xlab = "", main = "", notch = TRUE )
+        b = beeswarm::beeswarm(R2~part+group, data =df,las = 2, col = cols[c(3, 1, 2)], xlab = "", main = "" , add=TRUE, method = "center", spacing = 0.3)
+        graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
+      } else {
+        graphics::par(mfrow = c(1, 1), mar = c(4, 4, 4, 1), xaxt= "s")
+        pred = scale(trait, center = TRUE, scale = FALSE)
+        
+        graphics::plot(NULL, NULL, xlim = c(min(pred), max(pred)), ylim = c(minRS, maxRS*1.5), xlab = "Predictor",main = "", ylab = "R2", las =1)
+        for(i in 1:3) {
+          graphics::points(pred, rr$internals$Species[,i], col = ggplot2::alpha(cols[i], 0.2), pch = 16)
+          g = qgam::qgam( Y ~ pred, data = data.frame(Y = rr$internals$Species[,i], pred = pred), qu = 0.5, control = list(progress="none"))
+          graphics::abline(a = coef(g)[c(1, 2)], col = cols[i], lwd = lwd, lty =  1*(summary(g)$p.table[2,4] > 0.05)+1 )
+          out$pred[[colnames(rr$internals$Species)[i]]] = g 
+        }
+        graphics::legend("topright", legend = c("env", "spa", "codist"), col = cols, pch = 15, bty = "n")
+        graphics::legend("topleft", legend = c("significant", "non-significant"), col = c("black", "black"),  bty = "n", lty = c(1,2))
+      }
+      
+    }
+    
+    
+  }
+  
+  return(invisible(out))
+}
+
+
+
+
+get_eigen = function(X, double_center = TRUE, full = FALSE) {
+  D = as.matrix(dist(scale(X)))
+  if(double_center) D = D - (diag(0.0, ncol(D)) + rowMeans(D)) - t(diag(0.0, ncol(D)) + colMeans(D)) + mean(D) # Double center
+  eig = eigen(D)
+  if(!full)return(abs(eig$vectors[,which.max(abs(eig$values))]))
+  else {
+    return(list(x = eig$vectors))
+  }
+}
+
diff --git a/sjSDM/R/plot.R b/sjSDM/R/plot.R
index 071b791..5d39360 100644
--- a/sjSDM/R/plot.R
+++ b/sjSDM/R/plot.R
@@ -43,7 +43,7 @@ plotsjSDMcoef = function(object,wrap_col=NULL,group=NULL,col=NULL,slist=NULL) {
     inherits(object$settings$env, "linear")
   )
   oldpar = par(no.readonly = TRUE)
-  on.exit(par(oldpar))
+  on.exit(suppressWarnings(par(oldpar)))
   
   if(is.null(object$se)) object=getSe(object)
   summary.se = summary(object)
diff --git a/sjSDM/R/sjSDM.R b/sjSDM/R/sjSDM.R
index 8d1e660..dccd567 100644
--- a/sjSDM/R/sjSDM.R
+++ b/sjSDM/R/sjSDM.R
@@ -69,7 +69,7 @@
 #' 
 #' \subsection{Supported distributions}{
 #' 
-#' Currently supported distributions and link functions:
+#' Currently supported distributions and link functions, which are :
 #' \itemize{
 #' \item \code{\link{binomial}}: \code{"probit"} or \code{"logit"}
 #' \item \code{\link{poisson}}: \code{"log"} 
@@ -349,6 +349,8 @@ print.sjSDM = function(x, ...) {
 #' 
 #' @return Matrix of predictions (sites by species)
 #' 
+#' @example /inst/examples/predict-example.R
+#' 
 #' @import checkmate
 #' @export
 predict.sjSDM = function(object, newdata = NULL, SP = NULL, Y = NULL, type = c("link", "raw"), dropout = FALSE,...) {
@@ -736,3 +738,21 @@ update.sjSDM = function(object, env_formula = NULL, spatial_formula = NULL, biot
   )
   return(new_model)
 }
+
+
+#' Residuals for a sjSDM model
+#'
+#' Returns residuals for a fitted sjSDM model
+#'
+#' @param object a model fitted by \code{\link{sjSDM}}
+#' @param type residual type. Currently only supports raw
+#' @param ... further arguments, not supported yet.
+#' @return residuals in the format of the provided community matrix
+#' 
+#' @export
+residuals.sjSDM <- function(object, type = "raw", ...){
+  raw = object$data$Y - predict(object, type = "raw")
+  return(raw)
+}
+
+
diff --git a/sjSDM/README.Rmd b/sjSDM/README.Rmd
deleted file mode 100644
index ca057ec..0000000
--- a/sjSDM/README.Rmd
+++ /dev/null
@@ -1,125 +0,0 @@
----
-output: github_document
----
-
-<!-- README.md is generated from README.Rmd. Please edit that file -->
-
-```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = TRUE, collapse = TRUE, fig.path = "man/figures/README-")
-```
-
-
-# sjSDM - Fast and accurate Joint Species Distribution Modeling
-
-[![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
-![R-CMD-check](https://github.com/TheoreticalEcology/s-jSDM/workflows/R-CMD-check/badge.svg?branch=master)
-[![Publication](https://img.shields.io/badge/Publication-10.1111/2041-green.svg)](https://www.doi.org/10.1111/2041-210X.13687)
-
-## Overview
-
-A scalable method to estimates joint Species Distribution Models (jSDMs) based on the multivariate probit model through Monte-Carlo approximation of the joint likelihood. The numerical approximation is based on 'PyTorch' and 'reticulate', and can be calculated on CPUs and GPUs alike. 
-
-The method is described in [Pichler & Hartig (2021)](https://www.doi.org/10.1111/2041-210X.13687) A new joint species distribution model for faster and more accurate inference of species associations from big community data.
-
-The package includes options to fit various different (j)SDM models:
-
-* jSDMs with Binomial, Poisson, negative Binomial, and Normal distributed responses
-* jSDMs based on deep neural networks
-* Spatial auto-correlation can be accounted for by spatial eigenvectors or trend surface polynomials 
-
-To get more information, install the package and run
-
-```{r, eval=FALSE}
-library(sjSDM)
-?sjSDM
-vignette("sjSDM", package="sjSDM")
-```
-
-## Installation
-
-__sjSDM__ is based on 'PyTorch', a 'python' library, and thus requires 'python' dependencies. The 'python' dependencies can be automatically installed by running:
-
-```{r,eval=FALSE}
-library(sjSDM)
-install_sjSDM()
-```
-
-If this didn't work, please check the troubleshooting guide:
-```{r,eval=FALSE}
-library(sjSDM)
-?installation_help
-```
-
-## Usage
-
-## Workflow
-
-Simulate a community and fit a sjSDM model:
-
-```{r}
-library(sjSDM)
-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)
-plot(model)
-```
-
-We support other distributions:
-
--   Count data with Poisson:
-
-    ```{r,eval=FALSE}
-    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"))
-    ```
-
--   Count data with negative Binomial (which is still experimental, if you run into errors/problems, please let us know):
-
-    ```{r,eval=FALSE}
-    model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1+X2+X3), spatial = linear(data = SP, formula = ~0+X1:X2), se = TRUE, family="nbinom")
-    ```
-
--   Gaussian (normal):
-
-    ```{r,eval=FALSE}
-    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"))
-    ```
-
-### Anova
-
-ANOVA can be used to partition the three components (abiotic, biotic, and spatial):
-
-```{r,fig.height=7, fig.width=6.3}
-an = anova(model)
-print(an)
-plot(an)
-```
-
-The anova shows the relative changes in the R^2^ of the groups and their intersections.
-
-### Internal metacommunity structure
-
-Following [Leibold et al., 2022](https://doi.org/10.1111/oik.08618) we can calculate and visualize the internal metacommunity structure (=partitioning of the three components for species and sites). The internal structure is already calculated by the ANOVA and we can visualize it with the plot method:
-
-```{r,fig.height=7, fig.width=8, warning=FALSE}
-results = plotInternalStructure(an) # or plot(an, internal = TRUE)
-```
-
-The plot function returns the results for the internal metacommunity structure:
-
-```{r}
-print(results$data$Species)
-```
-
-## Deep neural networks
-
-Change linear part of model to a deep neural network:
-
-```{r}
-DNN <- sjSDM(Y = Occ, env = DNN(data = Env, formula = ~.), spatial = linear(data = SP, formula = ~0+X1:X2), se = TRUE, family=binomial("probit"), sampling = 100L)
-summary(DNN)
-```
diff --git a/sjSDM/README.md b/sjSDM/README.md
index 3734538..484e56e 100644
--- a/sjSDM/README.md
+++ b/sjSDM/README.md
@@ -1,355 +1,6 @@
-
-<!-- README.md is generated from README.Rmd. Please edit that file -->
-
 # sjSDM - Fast and accurate Joint Species Distribution Modeling
 
-[![License: GPL
-v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
-![R-CMD-check](https://github.com/TheoreticalEcology/s-jSDM/workflows/R-CMD-check/badge.svg?branch=master)
-[![Publication](https://img.shields.io/badge/Publication-10.1111/2041-green.svg)](https://www.doi.org/10.1111/2041-210X.13687)
-
-## Overview
-
-A scalable method to estimates joint Species Distribution Models (jSDMs)
-based on the multivariate probit model through Monte-Carlo approximation
-of the joint likelihood. The numerical approximation is based on
-‘PyTorch’ and ‘reticulate’, and can be calculated on CPUs and GPUs
-alike.
-
-The method is described in [Pichler & Hartig
-(2021)](https://www.doi.org/10.1111/2041-210X.13687)
-A new joint species distribution model for faster and more accurate
-inference of species associations from big community data.
-
-The package includes options to fit various different (j)SDM models:
-
-- jSDMs with Binomial, Poisson, negative Binomial, and Normal
-  distributed responses
-- jSDMs based on deep neural networks
-- Spatial auto-correlation can be accounted for by spatial eigenvectors
-  or trend surface polynomials
-
-To get more information, install the package and run
-
-``` r
-library(sjSDM)
-?sjSDM
-vignette("sjSDM", package="sjSDM")
-```
-
-## Installation
-
-**sjSDM** is based on ‘PyTorch’, a ‘python’ library, and thus requires
-‘python’ dependencies. The ‘python’ dependencies can be automatically
-installed by running:
-
-``` r
-library(sjSDM)
-install_sjSDM()
-```
-
-If this didn’t work, please check the troubleshooting guide:
-
-``` r
-library(sjSDM)
-?installation_help
-```
-
-## Usage
-
-## Workflow
-
-Simulate a community and fit a sjSDM model:
-
-``` r
-library(sjSDM)
-## ── Attaching sjSDM ──────────────────────────────────────────────────── 1.0.4 ──
-## ✔ 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)
-## Family:  binomial 
-## 
-## LogLik:  -508.7274 
-## Regularization loss:  0 
-## 
-## Species-species correlation matrix: 
-## 
-##  sp1  1.0000                                 
-##  sp2 -0.3650  1.0000                             
-##  sp3 -0.1910 -0.4300  1.0000                         
-##  sp4 -0.1800 -0.3670  0.8280  1.0000                     
-##  sp5  0.6860 -0.3800 -0.1050 -0.0890  1.0000                 
-##  sp6 -0.2860  0.4900  0.1730  0.1960 -0.0970  1.0000             
-##  sp7  0.5640 -0.1080  0.1310  0.1630  0.5580  0.2810  1.0000         
-##  sp8  0.2920  0.1920 -0.5070 -0.5220  0.2050 -0.0380  0.1150  1.0000     
-##  sp9 -0.0670 -0.0540  0.0480  0.0560 -0.3910 -0.3570 -0.2350 -0.1250  1.0000 
-##  sp10     0.2110  0.4840 -0.7090 -0.6490  0.2550  0.1420  0.1330  0.4510 -0.2690  1.0000
-## 
-## 
-## 
-## Spatial: 
-##            sp1       sp2      sp3       sp4      sp5     sp6      sp7      sp8
-## X1:X2 1.807111 -3.891649 3.616006 0.2416204 2.371867 1.23862 3.087861 1.890661
-##            sp9     sp10
-## X1:X2 1.237618 1.302866
-## 
-## 
-## 
-##                  Estimate  Std.Err Z value Pr(>|z|)    
-## sp1 (Intercept)  -0.05573  0.27360   -0.20  0.83858    
-## sp1 X1            1.34006  0.56886    2.36  0.01849 *  
-## sp1 X2           -2.42857  0.51421   -4.72  2.3e-06 ***
-## sp1 X3           -0.27071  0.43979   -0.62  0.53819    
-## sp2 (Intercept)   0.00232  0.27897    0.01  0.99336    
-## sp2 X1            1.38987  0.59367    2.34  0.01922 *  
-## sp2 X2            0.36444  0.50596    0.72  0.47135    
-## sp2 X3            0.68465  0.43890    1.56  0.11878    
-## sp3 (Intercept)  -0.55909  0.27776   -2.01  0.04413 *  
-## sp3 X1            1.49674  0.51045    2.93  0.00337 ** 
-## sp3 X2           -0.47064  0.49361   -0.95  0.34036    
-## sp3 X3           -1.12252  0.49705   -2.26  0.02392 *  
-## sp4 (Intercept)  -0.09188  0.23848   -0.39  0.70003    
-## sp4 X1           -1.55543  0.47517   -3.27  0.00106 ** 
-## sp4 X2           -1.91234  0.45110   -4.24  2.2e-05 ***
-## sp4 X3           -0.37758  0.40194   -0.94  0.34752    
-## sp5 (Intercept)  -0.22248  0.26192   -0.85  0.39564    
-## sp5 X1            0.76319  0.51134    1.49  0.13556    
-## sp5 X2            0.56888  0.50411    1.13  0.25911    
-## sp5 X3           -0.73548  0.45295   -1.62  0.10443    
-## sp6 (Intercept)   0.30172  0.26113    1.16  0.24792    
-## sp6 X1            2.66774  0.53908    4.95  7.5e-07 ***
-## sp6 X2           -1.09261  0.49401   -2.21  0.02699 *  
-## sp6 X3            0.19688  0.42018    0.47  0.63938    
-## sp7 (Intercept)  -0.02088  0.23800   -0.09  0.93009    
-## sp7 X1           -0.30568  0.47481   -0.64  0.51970    
-## sp7 X2            0.32681  0.43209    0.76  0.44945    
-## sp7 X3           -1.51884  0.40833   -3.72  0.00020 ***
-## sp8 (Intercept)   0.16312  0.16370    1.00  0.31905    
-## sp8 X1            0.34875  0.31925    1.09  0.27466    
-## sp8 X2            0.31206  0.30730    1.02  0.30988    
-## sp8 X3           -1.20898  0.28881   -4.19  2.8e-05 ***
-## sp9 (Intercept)   0.02942  0.19061    0.15  0.87733    
-## sp9 X1            1.39647  0.37376    3.74  0.00019 ***
-## sp9 X2           -1.09693  0.36395   -3.01  0.00258 ** 
-## sp9 X3            0.76818  0.31463    2.44  0.01463 *  
-## sp10 (Intercept) -0.10078  0.20713   -0.49  0.62656    
-## sp10 X1          -0.51072  0.37896   -1.35  0.17776    
-## sp10 X2          -1.23867  0.37494   -3.30  0.00095 ***
-## sp10 X3          -0.55472  0.35597   -1.56  0.11916    
-## ---
-## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-plot(model)
-## Family:  binomial 
-## 
-## LogLik:  -508.7274 
-## Regularization loss:  0 
-## 
-## Species-species correlation matrix: 
-## 
-##  sp1  1.0000                                 
-##  sp2 -0.3650  1.0000                             
-##  sp3 -0.1910 -0.4300  1.0000                         
-##  sp4 -0.1800 -0.3670  0.8280  1.0000                     
-##  sp5  0.6860 -0.3800 -0.1050 -0.0890  1.0000                 
-##  sp6 -0.2860  0.4900  0.1730  0.1960 -0.0970  1.0000             
-##  sp7  0.5640 -0.1080  0.1310  0.1630  0.5580  0.2810  1.0000         
-##  sp8  0.2920  0.1920 -0.5070 -0.5220  0.2050 -0.0380  0.1150  1.0000     
-##  sp9 -0.0670 -0.0540  0.0480  0.0560 -0.3910 -0.3570 -0.2350 -0.1250  1.0000 
-##  sp10     0.2110  0.4840 -0.7090 -0.6490  0.2550  0.1420  0.1330  0.4510 -0.2690  1.0000
-## 
-## 
-## 
-## Spatial: 
-##            sp1       sp2      sp3       sp4      sp5     sp6      sp7      sp8
-## X1:X2 1.807111 -3.891649 3.616006 0.2416204 2.371867 1.23862 3.087861 1.890661
-##            sp9     sp10
-## X1:X2 1.237618 1.302866
-## 
-## 
-## 
-##                  Estimate  Std.Err Z value Pr(>|z|)    
-## sp1 (Intercept)  -0.05573  0.27360   -0.20  0.83858    
-## sp1 X1            1.34006  0.56886    2.36  0.01849 *  
-## sp1 X2           -2.42857  0.51421   -4.72  2.3e-06 ***
-## sp1 X3           -0.27071  0.43979   -0.62  0.53819    
-## sp2 (Intercept)   0.00232  0.27897    0.01  0.99336    
-## sp2 X1            1.38987  0.59367    2.34  0.01922 *  
-## sp2 X2            0.36444  0.50596    0.72  0.47135    
-## sp2 X3            0.68465  0.43890    1.56  0.11878    
-## sp3 (Intercept)  -0.55909  0.27776   -2.01  0.04413 *  
-## sp3 X1            1.49674  0.51045    2.93  0.00337 ** 
-## sp3 X2           -0.47064  0.49361   -0.95  0.34036    
-## sp3 X3           -1.12252  0.49705   -2.26  0.02392 *  
-## sp4 (Intercept)  -0.09188  0.23848   -0.39  0.70003    
-## sp4 X1           -1.55543  0.47517   -3.27  0.00106 ** 
-## sp4 X2           -1.91234  0.45110   -4.24  2.2e-05 ***
-## sp4 X3           -0.37758  0.40194   -0.94  0.34752    
-## sp5 (Intercept)  -0.22248  0.26192   -0.85  0.39564    
-## sp5 X1            0.76319  0.51134    1.49  0.13556    
-## sp5 X2            0.56888  0.50411    1.13  0.25911    
-## sp5 X3           -0.73548  0.45295   -1.62  0.10443    
-## sp6 (Intercept)   0.30172  0.26113    1.16  0.24792    
-## sp6 X1            2.66774  0.53908    4.95  7.5e-07 ***
-## sp6 X2           -1.09261  0.49401   -2.21  0.02699 *  
-## sp6 X3            0.19688  0.42018    0.47  0.63938    
-## sp7 (Intercept)  -0.02088  0.23800   -0.09  0.93009    
-## sp7 X1           -0.30568  0.47481   -0.64  0.51970    
-## sp7 X2            0.32681  0.43209    0.76  0.44945    
-## sp7 X3           -1.51884  0.40833   -3.72  0.00020 ***
-## sp8 (Intercept)   0.16312  0.16370    1.00  0.31905    
-## sp8 X1            0.34875  0.31925    1.09  0.27466    
-## sp8 X2            0.31206  0.30730    1.02  0.30988    
-## sp8 X3           -1.20898  0.28881   -4.19  2.8e-05 ***
-## sp9 (Intercept)   0.02942  0.19061    0.15  0.87733    
-## sp9 X1            1.39647  0.37376    3.74  0.00019 ***
-## sp9 X2           -1.09693  0.36395   -3.01  0.00258 ** 
-## sp9 X3            0.76818  0.31463    2.44  0.01463 *  
-## sp10 (Intercept) -0.10078  0.20713   -0.49  0.62656    
-## sp10 X1          -0.51072  0.37896   -1.35  0.17776    
-## sp10 X2          -1.23867  0.37494   -3.30  0.00095 ***
-## sp10 X3          -0.55472  0.35597   -1.56  0.11916    
-## ---
-## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-```
-
-![](man/figures/README-unnamed-chunk-4-1.png)<!-- -->
-
-We support other distributions:
-
-- Count data with Poisson:
-
-  ``` r
-  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"))
-  ```
-
-- Count data with negative Binomial (which is still experimental, if you
-  run into errors/problems, please let us know):
-
-  ``` r
-  model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1+X2+X3), spatial = linear(data = SP, formula = ~0+X1:X2), se = TRUE, family="nbinom")
-  ```
-
-- Gaussian (normal):
-
-  ``` r
-  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"))
-  ```
-
-### Anova
-
-ANOVA can be used to partition the three components (abiotic, biotic,
-and spatial):
-
-``` r
-an = anova(model)
-print(an)
-## Analysis of Deviance Table
-## 
-## Terms added sequentially:
-## 
-##           Deviance Residual deviance R2 Nagelkerke R2 McFadden
-## Abiotic  157.92268        1177.37216       0.79387      0.1139
-## Biotic   175.40897        1159.88588       0.82694      0.1265
-## Spatial   16.28744        1319.00741       0.15030      0.0117
-## Full     387.13884         948.15601       0.97917      0.2793
-plot(an)
-```
-
-![](man/figures/README-unnamed-chunk-8-1.png)<!-- -->
-
-The anova shows the relative changes in the R<sup>2</sup> of the groups
-and their intersections.
-
-### Internal metacommunity structure
-
-Following [Leibold et al., 2022](https://doi.org/10.1111/oik.08618) we
-can calculate and visualize the internal metacommunity structure
-(=partitioning of the three components for species and sites). The
-internal structure is already calculated by the ANOVA and we can
-visualize it with the plot method:
-
-``` r
-results = plotInternalStructure(an) # or plot(an, internal = TRUE)
-## Registered S3 methods overwritten by 'ggtern':
-##   method           from   
-##   grid.draw.ggplot ggplot2
-##   plot.ggplot      ggplot2
-##   print.ggplot     ggplot2
-```
-
-![](man/figures/README-unnamed-chunk-9-1.png)<!-- -->
-
-The plot function returns the results for the internal metacommunity
-structure:
-
-``` r
-print(results$data$Species)
-##           env         spa     codist         r2
-## 1  0.18030741 0.000000000 0.16283743 0.03411591
-## 2  0.08715478 0.019089232 0.18394478 0.02901888
-## 3  0.11717770 0.013074584 0.20191543 0.03321677
-## 4  0.16619488 0.003832403 0.16366855 0.03336958
-## 5  0.08432392 0.000000000 0.16981670 0.02522153
-## 6  0.18586884 0.000000000 0.11975589 0.03054916
-## 7  0.11465238 0.016740799 0.13184580 0.02632390
-## 8  0.13783609 0.008329915 0.05558837 0.02017544
-## 9  0.16799688 0.014057668 0.04716048 0.02292150
-## 10 0.09935810 0.007978024 0.13615356 0.02434897
-```
+sjSDM - an R package to estimates joint Species Distribution Models (jSDMs) in a scalable way. For more info on the package, go to https://theoreticalecology.github.io/s-jSDM/
 
-## Deep neural networks
 
-Change linear part of model to a deep neural network:
 
-``` r
-DNN <- sjSDM(Y = Occ, env = DNN(data = Env, formula = ~.), spatial = linear(data = SP, formula = ~0+X1:X2), se = TRUE, family=binomial("probit"), sampling = 100L)
-summary(DNN)
-## Family:  binomial 
-## 
-## LogLik:  -466.828 
-## Regularization loss:  0 
-## 
-## Species-species correlation matrix: 
-## 
-##  sp1  1.0000                                 
-##  sp2 -0.4510  1.0000                             
-##  sp3 -0.1620 -0.3370  1.0000                         
-##  sp4 -0.0730 -0.4010  0.8670  1.0000                     
-##  sp5  0.6450 -0.3200 -0.1800 -0.1180  1.0000                 
-##  sp6 -0.3620  0.4120  0.3070  0.1980 -0.0680  1.0000             
-##  sp7  0.5680 -0.0980  0.1680  0.2080  0.5430  0.2750  1.0000         
-##  sp8  0.2310  0.1720 -0.5000 -0.5500  0.1380 -0.0170  0.0480  1.0000     
-##  sp9 -0.0190  0.0530  0.0130  0.0650 -0.4540 -0.4090 -0.2020 -0.0820  1.0000 
-##  sp10     0.1290  0.4600 -0.7080 -0.7200  0.3200  0.0790  0.1160  0.4290 -0.2620  1.0000
-## 
-## 
-## 
-## Spatial: 
-##            sp1       sp2      sp3       sp4      sp5      sp6      sp7      sp8
-## X1:X2 1.763654 -3.796902 3.750247 0.6460114 3.022765 1.352927 3.291612 2.691203
-##            sp9     sp10
-## X1:X2 1.119004 1.262632
-## 
-## 
-## 
-## Env architecture:
-## ===================================
-## Layer_1:  (4, 10)
-## Layer_2:  SELU
-## Layer_3:  (10, 10)
-## Layer_4:  SELU
-## Layer_5:  (10, 10)
-## Layer_6:  SELU
-## Layer_7:  (10, 10)
-## ===================================
-## Weights :     340
-```
diff --git a/sjSDM/inst/examples/Anova-example.R b/sjSDM/inst/examples/Anova-example.R
new file mode 100644
index 0000000..6c50488
--- /dev/null
+++ b/sjSDM/inst/examples/Anova-example.R
@@ -0,0 +1,41 @@
+library(sjSDM)
+# simulate community:
+community = simulate_SDM(env = 3L, species = 10L, sites = 100L)
+
+Occ <- community$response
+Env <- community$env_weights
+SP <- data.frame(matrix(rnorm(200, 0, 0.3), 100, 2)) # spatial coordinates
+
+
+# fit model:
+model <- sjSDM(Y = Occ, 
+               env = linear(data = Env, formula = ~X1+X2+X3), 
+               spatial = linear(data = SP, formula = ~0+X1*X2), 
+               family=binomial("probit"),
+               iter = 20) # increase iter for real analysis
+
+# Calculate ANOVA for env, space, associations, for details see ?anova.sjSDM
+an = anova(model, samples = 10) # increase iter for real analysis
+
+# Show anova fractions
+plot(an)
+
+# ANOVA tables with different way to handle fractions
+summary(an)
+summary(an, fractions = "discard")
+summary(an, fractions = "proportional")
+summary(an, fractions = "equal")
+
+# Internal structure
+int = internalStructure(an, fractions = "equal")
+
+print(int)
+
+plot(int)
+
+plotAssemblyEffects(int)
+plotAssemblyEffects(int, env = runif(100))
+plotAssemblyEffects(int, env = as.factor(c(rep(1, 50), rep(2, 50))))
+plotAssemblyEffects(int, trait = runif(10))
+plotAssemblyEffects(int, trait = as.factor(c(rep(1, 5), rep(2, 5))))
+
diff --git a/sjSDM/inst/examples/plot.sjSDM-example.R b/sjSDM/inst/examples/plot.sjSDM-example.R
index 236088e..d1155cf 100644
--- a/sjSDM/inst/examples/plot.sjSDM-example.R
+++ b/sjSDM/inst/examples/plot.sjSDM-example.R
@@ -1,15 +1,16 @@
-\dontrun{
 library(sjSDM)
 # simulate community:
 com = simulate_SDM(env = 6L, species = 7L, sites = 100L)
 
 # fit model:
-model = sjSDM(Y = com$response,env = com$env_weights, iter = 2L,se = TRUE) 
+model = sjSDM(Y = com$response,env = com$env_weights, iter = 2L, se = TRUE) 
 
-#create a group dataframe for plot
+# normal plot
+plot(model)
+
+# colored by groups
 species=c("sp1","sp2","sp3","sp4","sp5","sp6","sp7")
 group=c("mammal","bird","fish","fish","mammal","amphibian","amphibian")
 group = data.frame(species=species,group=group)
 
 plot(model,group=group)
-}
\ No newline at end of file
diff --git a/sjSDM/inst/examples/predict-example.R b/sjSDM/inst/examples/predict-example.R
new file mode 100644
index 0000000..c5051d3
--- /dev/null
+++ b/sjSDM/inst/examples/predict-example.R
@@ -0,0 +1,22 @@
+\dontrun{
+
+## Conditional predictions based on focal species
+com = simulate_SDM(sites = 200L)
+## first 100 observations are the training data
+model = sjSDM(com$response[1:100, ], com$env_weights[1:100,])
+## Assume that for the other 100 observations, only the first species is missing 
+## and we want to use the other 4 species to improve the predictions:
+Y_focal = com$response[101:200, ]
+Y_focal[,1] = NA # set to NA because occurrences are unknown
+
+pred_conditional = predict(model, newdata = com$env_weights[101:200,], Y = Y_focal)
+pred_unconditional = predict(model, newdata = com$env_weights[101:200,])[,1]
+
+## Compare performance:
+Metrics::auc(com$response[101:200, 1], pred_conditional)
+Metrics::auc(com$response[101:200, 1], pred_unconditional)
+
+## Conditional predictions are better, however, it only works if occurrences of
+## other species for new sites are known!
+
+}
\ No newline at end of file
diff --git a/sjSDM/inst/examples/sjSDM-example.R b/sjSDM/inst/examples/sjSDM-example.R
index 780114a..23810cd 100644
--- a/sjSDM/inst/examples/sjSDM-example.R
+++ b/sjSDM/inst/examples/sjSDM-example.R
@@ -8,6 +8,10 @@ com = simulate_SDM(env = 3L, species = 7L, sites = 100L)
 model = sjSDM(Y = com$response,env = com$env_weights, iter = 50L) 
 # increase iter for your own data 
 
+# Default distribution is binomial("probit"). Alternatively, you can use
+# binomial(logit), poisson("log"), "nbinom" (with log, still somewhat 
+# experimental) and gaussian("identity")
+
 coef(model)
 summary(model)
 getCov(model)
@@ -69,7 +73,16 @@ predict(model, newdata = com$env_weights, SP = SPV)
 
 ## Visualize internal meta-community structure
 an = anova(model)
-plot(an, internal=TRUE)
+
+internal = internalStructure(an)
+plot(internal)
+
+## Visualize community assemlby effects 
+
+plotAssemblyEffects(internal)
+
+### see ?anova.sjSDM for mroe details
+
 
 ## non-linear(deep neural network) model
 model = sjSDM(Y = com$response, env = linear(com$env_weights), 
@@ -108,30 +121,10 @@ result = anova(model)
 print(result)
 plot(result)
 
-## visualize meta-community structure
-plot(result, internal=TRUE)
-
-
-
-# Conditional predictions
-## Conditional predictions based on focal species
-com = simulate_SDM(sites = 200L)
-## first 100 observations are the training data
-model = sjSDM(com$response[1:100, ], com$env_weights[1:100,])
-## Assume that for the other 100 observations, only the first species is missing 
-## and we want to use the other 4 species to improve the predictions:
-Y_focal = com$response[101:200, ]
-Y_focal[,1] = NA # set to NA because occurrences are unknown
-
-pred_conditional = predict(model, newdata = com$env_weights[101:200,], Y = Y_focal)
-pred_unconditional = predict(model, newdata = com$env_weights[101:200,])[,1]
-
-## Compare performance:
-Metrics::auc(com$response[101:200, 1], pred_conditional)
-Metrics::auc(com$response[101:200, 1], pred_unconditional)
+## visualize internal meta-community structure
+internal = internalStructure(an)
+plot(internal)
 
-## Conditional predictions are better, however, it only works if occurrences of
-## other species for new sites are known!
 
 
 # Deep neural networks
diff --git a/sjSDM/man/DNN.Rd b/sjSDM/man/DNN.Rd
index bf2f2de..a17b098 100644
--- a/sjSDM/man/DNN.Rd
+++ b/sjSDM/man/DNN.Rd
@@ -60,6 +60,10 @@ com = simulate_SDM(env = 3L, species = 7L, sites = 100L)
 model = sjSDM(Y = com$response,env = com$env_weights, iter = 50L) 
 # increase iter for your own data 
 
+# Default distribution is binomial("probit"). Alternatively, you can use
+# binomial(logit), poisson("log"), "nbinom" (with log, still somewhat 
+# experimental) and gaussian("identity")
+
 coef(model)
 summary(model)
 getCov(model)
@@ -121,7 +125,16 @@ predict(model, newdata = com$env_weights, SP = SPV)
 
 ## Visualize internal meta-community structure
 an = anova(model)
-plot(an, internal=TRUE)
+
+internal = internalStructure(an)
+plot(internal)
+
+## Visualize community assemlby effects 
+
+plotAssemblyEffects(internal)
+
+### see ?anova.sjSDM for mroe details
+
 
 ## non-linear(deep neural network) model
 model = sjSDM(Y = com$response, env = linear(com$env_weights), 
@@ -160,30 +173,10 @@ result = anova(model)
 print(result)
 plot(result)
 
-## visualize meta-community structure
-plot(result, internal=TRUE)
-
-
-
-# Conditional predictions
-## Conditional predictions based on focal species
-com = simulate_SDM(sites = 200L)
-## first 100 observations are the training data
-model = sjSDM(com$response[1:100, ], com$env_weights[1:100,])
-## Assume that for the other 100 observations, only the first species is missing 
-## and we want to use the other 4 species to improve the predictions:
-Y_focal = com$response[101:200, ]
-Y_focal[,1] = NA # set to NA because occurrences are unknown
-
-pred_conditional = predict(model, newdata = com$env_weights[101:200,], Y = Y_focal)
-pred_unconditional = predict(model, newdata = com$env_weights[101:200,])[,1]
-
-## Compare performance:
-Metrics::auc(com$response[101:200, 1], pred_conditional)
-Metrics::auc(com$response[101:200, 1], pred_unconditional)
+## visualize internal meta-community structure
+internal = internalStructure(an)
+plot(internal)
 
-## Conditional predictions are better, however, it only works if occurrences of
-## other species for new sites are known!
 
 
 # Deep neural networks
diff --git a/sjSDM/man/anova.sjSDM.Rd b/sjSDM/man/anova.sjSDM.Rd
index 09c1011..a4e9bce 100644
--- a/sjSDM/man/anova.sjSDM.Rd
+++ b/sjSDM/man/anova.sjSDM.Rd
@@ -2,17 +2,15 @@
 % Please edit documentation in R/anova.R
 \name{anova.sjSDM}
 \alias{anova.sjSDM}
-\title{Anova}
+\title{Anova / Variation partitioning}
 \usage{
-\method{anova}{sjSDM}(object, samples = 5000L, fractional = c("proportional", "equal"), ...)
+\method{anova}{sjSDM}(object, samples = 5000L, ...)
 }
 \arguments{
 \item{object}{model of object \code{\link{sjSDM}}}
 
 \item{samples}{Number of Monte Carlo samples}
 
-\item{fractional}{how to distribute the shared partitions, proportional to the fractions' R2s or equal}
-
 \item{...}{optional arguments which are passed to the calculation of the logLikelihood}
 }
 \value{
@@ -30,14 +28,59 @@ An S3 class of type 'sjSDManova' including the following components:
 Implemented S3 methods are \code{\link{print.sjSDManova}} and \code{\link{plot.sjSDManova}}
 }
 \description{
-Calculates type II anova.
+Compute variance explained by the three fractions env, space, associations
+}
+\details{
+The ANOVA function removes each of the three fractions (Environment, Space, Associations) and measures the drop in variance explained, and thus the importance of the three fractions.
+
+Variance explained is measured by Deviance as well as the pseudo-R2 metrics of Nagelkerke and McFadden
 
-Shared contributions (e.g. between space and environment) are also calculated (and divided proportionally) and can be optionally visualized via \code{\link{plot.sjSDManova}} with \code{add_shared=TRUE}.
+In downstream functions such as \code{\link{plot.sjSDManova}} or \code{\link{plot.sjSDManova}} with \code{add_shared=TRUE}.
 The anova can get unstable for many species and few occurrences/observations. We recommend using large numbers for 'samples'.
 }
-\details{
-Compute analysis of variance
+\examples{
+library(sjSDM)
+# simulate community:
+community = simulate_SDM(env = 3L, species = 10L, sites = 100L)
+
+Occ <- community$response
+Env <- community$env_weights
+SP <- data.frame(matrix(rnorm(200, 0, 0.3), 100, 2)) # spatial coordinates
+
+
+# fit model:
+model <- sjSDM(Y = Occ, 
+               env = linear(data = Env, formula = ~X1+X2+X3), 
+               spatial = linear(data = SP, formula = ~0+X1*X2), 
+               family=binomial("probit"),
+               iter = 20) # increase iter for real analysis
+
+# Calculate ANOVA for env, space, associations, for details see ?anova.sjSDM
+an = anova(model, samples = 10) # increase iter for real analysis
+
+# Show anova fractions
+plot(an)
+
+# ANOVA tables with different way to handle fractions
+summary(an)
+summary(an, fractions = "discard")
+summary(an, fractions = "proportional")
+summary(an, fractions = "equal")
+
+# Internal structure
+int = internalStructure(an, fractions = "equal")
+
+print(int)
+
+plot(int)
+
+plotAssemblyEffects(int)
+plotAssemblyEffects(int, env = runif(100))
+plotAssemblyEffects(int, env = as.factor(c(rep(1, 50), rep(2, 50))))
+plotAssemblyEffects(int, trait = runif(10))
+plotAssemblyEffects(int, trait = as.factor(c(rep(1, 5), rep(2, 5))))
+
 }
 \seealso{
-\code{\link{plot.sjSDManova}}, \code{\link{print.sjSDManova}}, \code{\link{plotInternalStructure}}
+\code{\link{plot.sjSDManova}}, \code{\link{print.sjSDManova}},\code{\link{summary.sjSDManova}}, \code{\link{plot.sjSDMinternalStruture}}
 }
diff --git a/sjSDM/man/bioticStruct.Rd b/sjSDM/man/bioticStruct.Rd
index 63a6c27..7f3c292 100644
--- a/sjSDM/man/bioticStruct.Rd
+++ b/sjSDM/man/bioticStruct.Rd
@@ -55,6 +55,10 @@ com = simulate_SDM(env = 3L, species = 7L, sites = 100L)
 model = sjSDM(Y = com$response,env = com$env_weights, iter = 50L) 
 # increase iter for your own data 
 
+# Default distribution is binomial("probit"). Alternatively, you can use
+# binomial(logit), poisson("log"), "nbinom" (with log, still somewhat 
+# experimental) and gaussian("identity")
+
 coef(model)
 summary(model)
 getCov(model)
@@ -116,7 +120,16 @@ predict(model, newdata = com$env_weights, SP = SPV)
 
 ## Visualize internal meta-community structure
 an = anova(model)
-plot(an, internal=TRUE)
+
+internal = internalStructure(an)
+plot(internal)
+
+## Visualize community assemlby effects 
+
+plotAssemblyEffects(internal)
+
+### see ?anova.sjSDM for mroe details
+
 
 ## non-linear(deep neural network) model
 model = sjSDM(Y = com$response, env = linear(com$env_weights), 
@@ -155,30 +168,10 @@ result = anova(model)
 print(result)
 plot(result)
 
-## visualize meta-community structure
-plot(result, internal=TRUE)
-
-
-
-# Conditional predictions
-## Conditional predictions based on focal species
-com = simulate_SDM(sites = 200L)
-## first 100 observations are the training data
-model = sjSDM(com$response[1:100, ], com$env_weights[1:100,])
-## Assume that for the other 100 observations, only the first species is missing 
-## and we want to use the other 4 species to improve the predictions:
-Y_focal = com$response[101:200, ]
-Y_focal[,1] = NA # set to NA because occurrences are unknown
-
-pred_conditional = predict(model, newdata = com$env_weights[101:200,], Y = Y_focal)
-pred_unconditional = predict(model, newdata = com$env_weights[101:200,])[,1]
-
-## Compare performance:
-Metrics::auc(com$response[101:200, 1], pred_conditional)
-Metrics::auc(com$response[101:200, 1], pred_unconditional)
+## visualize internal meta-community structure
+internal = internalStructure(an)
+plot(internal)
 
-## Conditional predictions are better, however, it only works if occurrences of
-## other species for new sites are known!
 
 
 # Deep neural networks
diff --git a/sjSDM/man/importance.Rd b/sjSDM/man/importance.Rd
index 9a6edfe..c3e548e 100644
--- a/sjSDM/man/importance.Rd
+++ b/sjSDM/man/importance.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/importance.R
 \name{importance}
 \alias{importance}
-\title{importance}
+\title{Importance of environmental, spatial and association components}
 \usage{
 importance(x, save_memory = TRUE, ...)
 }
@@ -26,8 +26,11 @@ Implemented S3 methods include \code{\link{print.sjSDMimportance}} and \code{\li
 Computes standardized variance components with respect to abiotic, biotic, and spatial effect groups.
 }
 \details{
-This variance partitioning approach is based on Ovaskainen et al., 2017. For an example how to interpret the outputs, see Leibold et al., 2021.
-This function will be deprecated in the future. Please use \code{plot(anova(model), internal=TRUE)} (currently only supported for spatial models).
+This approach is based on Ovaskainen et al., 2017, and also used in  Leibold et al., 2021. Unlike the \code{\link{anova.sjSDM}} function in the sjSDM package, importance is not calculated by explicitly switching a particular model component of and refitting the model, but essentially by setting it ineffective. 
+
+Although we have no hard reasons to discourage the use of this function, we have decided in sjSDM to measure importance maninly based on a traditional ANOVA approach. We therefore recommend users to use the \code{\link{anova.sjSDM}}.
+
+This function is maintained hidden for comparison / benchmarking purpose, and in case there is a need to use it in the future. If you want to access it, use sjSDM:::importance.
 }
 \examples{
 \dontrun{
diff --git a/sjSDM/man/installation_help.Rd b/sjSDM/man/installation_help.Rd
index 7f14e0b..cc37c16 100644
--- a/sjSDM/man/installation_help.Rd
+++ b/sjSDM/man/installation_help.Rd
@@ -2,6 +2,7 @@
 % Please edit documentation in R/installation_help.R
 \docType{package}
 \name{installation_help}
+\alias{sjSDM-package}
 \alias{installation_help}
 \title{Installation help}
 \description{
@@ -163,3 +164,25 @@ via the sjSDM \href{https://github.com/TheoreticalEcology/s-jSDM/issues/}{issue
 with a copy of the \code{\link{install_diagnostic}} output as a quote.
 }
 
+\seealso{
+Useful links:
+\itemize{
+  \item \url{https://theoreticalecology.github.io/s-jSDM/}
+  \item Report bugs at \url{https://github.com/TheoreticalEcology/s-jSDM/issues}
+}
+
+}
+\author{
+\strong{Maintainer}: Maximilian Pichler \email{maximilian.pichler@biologie.uni-regensburg.de} (\href{https://orcid.org/0000-0003-2252-8327}{ORCID})
+
+Authors:
+\itemize{
+  \item Florian Hartig \email{florian.hartig@biologie.uni-regensburg.de} (\href{https://orcid.org/0000-0002-6255-9059}{ORCID})
+}
+
+Other contributors:
+\itemize{
+  \item Wang Cai \email{caiwang0503@163.com} [contributor]
+}
+
+}
diff --git a/sjSDM/man/internalStructure.Rd b/sjSDM/man/internalStructure.Rd
new file mode 100644
index 0000000..460b457
--- /dev/null
+++ b/sjSDM/man/internalStructure.Rd
@@ -0,0 +1,84 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/internalStructure.R
+\name{internalStructure}
+\alias{internalStructure}
+\title{Plot internal metacommunity structure}
+\usage{
+internalStructure(
+  object,
+  Rsquared = c("McFadden", "Nagelkerke"),
+  fractions = c("discard", "proportional", "equal"),
+  negatives = c("floor", "scale"),
+  plot = FALSE
+)
+}
+\arguments{
+\item{object}{anova object from \code{\link{anova.sjSDM}}}
+
+\item{Rsquared}{which R squared should be used, McFadden or Nagelkerke (McFadden is default)}
+
+\item{fractions}{how to handle shared fractions}
+
+\item{negatives}{how to handle negative R squareds}
+
+\item{plot}{should the plots be suppressed or not.
+
+Plots and returns the internal metacommunity structure of species and sites (see Leibold et al., 2022). 
+Plots were heavily inspired by Leibold et al., 2022}
+}
+\value{
+List with the following components:
+
+
+\item{plots}{ggplot objects for sites and species.}
+\item{data}{List of data.frames with the internal metacommunity structure.}
+}
+\description{
+Plot internal metacommunity structure
+}
+\examples{
+library(sjSDM)
+# simulate community:
+community = simulate_SDM(env = 3L, species = 10L, sites = 100L)
+
+Occ <- community$response
+Env <- community$env_weights
+SP <- data.frame(matrix(rnorm(200, 0, 0.3), 100, 2)) # spatial coordinates
+
+
+# fit model:
+model <- sjSDM(Y = Occ, 
+               env = linear(data = Env, formula = ~X1+X2+X3), 
+               spatial = linear(data = SP, formula = ~0+X1*X2), 
+               family=binomial("probit"),
+               iter = 20) # increase iter for real analysis
+
+# Calculate ANOVA for env, space, associations, for details see ?anova.sjSDM
+an = anova(model, samples = 10) # increase iter for real analysis
+
+# Show anova fractions
+plot(an)
+
+# ANOVA tables with different way to handle fractions
+summary(an)
+summary(an, fractions = "discard")
+summary(an, fractions = "proportional")
+summary(an, fractions = "equal")
+
+# Internal structure
+int = internalStructure(an, fractions = "equal")
+
+print(int)
+
+plot(int)
+
+plotAssemblyEffects(int)
+plotAssemblyEffects(int, env = runif(100))
+plotAssemblyEffects(int, env = as.factor(c(rep(1, 50), rep(2, 50))))
+plotAssemblyEffects(int, trait = runif(10))
+plotAssemblyEffects(int, trait = as.factor(c(rep(1, 5), rep(2, 5))))
+
+}
+\references{
+Leibold, M. A., Rudolph, F. J., Blanchet, F. G., De Meester, L., Gravel, D., Hartig, F., ... & Chase, J. M. (2022). The internal structure of metacommunities. Oikos, 2022(1).
+}
diff --git a/sjSDM/man/linear.Rd b/sjSDM/man/linear.Rd
index eae1a6b..c2802d5 100644
--- a/sjSDM/man/linear.Rd
+++ b/sjSDM/man/linear.Rd
@@ -40,6 +40,10 @@ com = simulate_SDM(env = 3L, species = 7L, sites = 100L)
 model = sjSDM(Y = com$response,env = com$env_weights, iter = 50L) 
 # increase iter for your own data 
 
+# Default distribution is binomial("probit"). Alternatively, you can use
+# binomial(logit), poisson("log"), "nbinom" (with log, still somewhat 
+# experimental) and gaussian("identity")
+
 coef(model)
 summary(model)
 getCov(model)
@@ -101,7 +105,16 @@ predict(model, newdata = com$env_weights, SP = SPV)
 
 ## Visualize internal meta-community structure
 an = anova(model)
-plot(an, internal=TRUE)
+
+internal = internalStructure(an)
+plot(internal)
+
+## Visualize community assemlby effects 
+
+plotAssemblyEffects(internal)
+
+### see ?anova.sjSDM for mroe details
+
 
 ## non-linear(deep neural network) model
 model = sjSDM(Y = com$response, env = linear(com$env_weights), 
@@ -140,30 +153,10 @@ result = anova(model)
 print(result)
 plot(result)
 
-## visualize meta-community structure
-plot(result, internal=TRUE)
-
-
-
-# Conditional predictions
-## Conditional predictions based on focal species
-com = simulate_SDM(sites = 200L)
-## first 100 observations are the training data
-model = sjSDM(com$response[1:100, ], com$env_weights[1:100,])
-## Assume that for the other 100 observations, only the first species is missing 
-## and we want to use the other 4 species to improve the predictions:
-Y_focal = com$response[101:200, ]
-Y_focal[,1] = NA # set to NA because occurrences are unknown
-
-pred_conditional = predict(model, newdata = com$env_weights[101:200,], Y = Y_focal)
-pred_unconditional = predict(model, newdata = com$env_weights[101:200,])[,1]
-
-## Compare performance:
-Metrics::auc(com$response[101:200, 1], pred_conditional)
-Metrics::auc(com$response[101:200, 1], pred_unconditional)
+## visualize internal meta-community structure
+internal = internalStructure(an)
+plot(internal)
 
-## Conditional predictions are better, however, it only works if occurrences of
-## other species for new sites are known!
 
 
 # Deep neural networks
diff --git a/sjSDM/man/plot.sjSDM.DNN.Rd b/sjSDM/man/plot.sjSDM.DNN.Rd
index e4400a6..652e8ba 100644
--- a/sjSDM/man/plot.sjSDM.DNN.Rd
+++ b/sjSDM/man/plot.sjSDM.DNN.Rd
@@ -18,19 +18,20 @@ No return value, called for side effects.
 Plot training loss history
 }
 \examples{
-\dontrun{
 library(sjSDM)
 # simulate community:
 com = simulate_SDM(env = 6L, species = 7L, sites = 100L)
 
 # fit model:
-model = sjSDM(Y = com$response,env = com$env_weights, iter = 2L,se = TRUE) 
+model = sjSDM(Y = com$response,env = com$env_weights, iter = 2L, se = TRUE) 
 
-#create a group dataframe for plot
+# normal plot
+plot(model)
+
+# colored by groups
 species=c("sp1","sp2","sp3","sp4","sp5","sp6","sp7")
 group=c("mammal","bird","fish","fish","mammal","amphibian","amphibian")
 group = data.frame(species=species,group=group)
 
 plot(model,group=group)
 }
-}
diff --git a/sjSDM/man/plot.sjSDM.Rd b/sjSDM/man/plot.sjSDM.Rd
index 192f972..1112785 100644
--- a/sjSDM/man/plot.sjSDM.Rd
+++ b/sjSDM/man/plot.sjSDM.Rd
@@ -19,22 +19,23 @@ Plotting coefficients returned by sjSDM model.
 This function only for model fitted by linear, fitted by DNN is not yet supported.
 }
 \examples{
-\dontrun{
 library(sjSDM)
 # simulate community:
 com = simulate_SDM(env = 6L, species = 7L, sites = 100L)
 
 # fit model:
-model = sjSDM(Y = com$response,env = com$env_weights, iter = 2L,se = TRUE) 
+model = sjSDM(Y = com$response,env = com$env_weights, iter = 2L, se = TRUE) 
 
-#create a group dataframe for plot
+# normal plot
+plot(model)
+
+# colored by groups
 species=c("sp1","sp2","sp3","sp4","sp5","sp6","sp7")
 group=c("mammal","bird","fish","fish","mammal","amphibian","amphibian")
 group = data.frame(species=species,group=group)
 
 plot(model,group=group)
 }
-}
 \seealso{
 \code{\link{plotsjSDMcoef}}
 }
diff --git a/sjSDM/man/plot.sjSDManova.Rd b/sjSDM/man/plot.sjSDManova.Rd
index 5fb7906..14ca4ed 100644
--- a/sjSDM/man/plot.sjSDManova.Rd
+++ b/sjSDM/man/plot.sjSDManova.Rd
@@ -8,12 +8,10 @@
   x,
   y,
   type = c("McFadden", "Deviance", "Nagelkerke"),
-  internal = FALSE,
-  add_shared = FALSE,
+  fractions = c("discard", "proportional", "equal"),
   cols = c("#7FC97F", "#BEAED4", "#FDC086"),
   alpha = 0.15,
   env_deviance = NULL,
-  suppress_plotting = FALSE,
   ...
 )
 }
@@ -24,9 +22,7 @@
 
 \item{type}{deviance, Nagelkerke or McFadden R-squared}
 
-\item{internal}{logical, plot internal or total structure}
-
-\item{add_shared}{Add shared contributions when plotting the internal structure}
+\item{fractions}{how to handle shared fractions}
 
 \item{cols}{colors for the groups}
 
@@ -34,21 +30,11 @@
 
 \item{env_deviance}{environmental deviance}
 
-\item{suppress_plotting}{return plots but don't plot them}
-
-\item{...}{Additional arguments to pass to \code{plot()}
-
-The \code{internal = TRUE} plot was heavily inspired by Leibold et al., 2022}
+\item{...}{Additional arguments to pass to \code{plot()}}
 }
 \value{
 List with the following components:
 
-If \code{internal=TRUE}:
-
-\item{plots}{ggplot objects for sites and species.}
-\item{data}{List of data.frames with the shown results.}
-
-else:
 \item{VENN}{Matrix of shown results.}
 }
 \description{
diff --git a/sjSDM/man/plot.sjSDMinternalStruture.Rd b/sjSDM/man/plot.sjSDMinternalStruture.Rd
new file mode 100644
index 0000000..33a7da3
--- /dev/null
+++ b/sjSDM/man/plot.sjSDMinternalStruture.Rd
@@ -0,0 +1,63 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/internalStructure.R
+\name{plot.sjSDMinternalStruture}
+\alias{plot.sjSDMinternalStruture}
+\title{Plot internal structure}
+\usage{
+\method{plot}{sjSDMinternalStruture}(x, alpha = 0.15, env_deviance = NULL, ...)
+}
+\arguments{
+\item{x}{and object of class sjSDMinternalStruture create by anova object from \code{\link{internalStructure}}}
+
+\item{alpha}{alpha of points}
+
+\item{env_deviance}{environmental deviance/gradient (points will be colored)}
+
+\item{...}{no function}
+}
+\description{
+Creates a ternary diagram of an object of class
+}
+\examples{
+library(sjSDM)
+# simulate community:
+community = simulate_SDM(env = 3L, species = 10L, sites = 100L)
+
+Occ <- community$response
+Env <- community$env_weights
+SP <- data.frame(matrix(rnorm(200, 0, 0.3), 100, 2)) # spatial coordinates
+
+
+# fit model:
+model <- sjSDM(Y = Occ, 
+               env = linear(data = Env, formula = ~X1+X2+X3), 
+               spatial = linear(data = SP, formula = ~0+X1*X2), 
+               family=binomial("probit"),
+               iter = 20) # increase iter for real analysis
+
+# Calculate ANOVA for env, space, associations, for details see ?anova.sjSDM
+an = anova(model, samples = 10) # increase iter for real analysis
+
+# Show anova fractions
+plot(an)
+
+# ANOVA tables with different way to handle fractions
+summary(an)
+summary(an, fractions = "discard")
+summary(an, fractions = "proportional")
+summary(an, fractions = "equal")
+
+# Internal structure
+int = internalStructure(an, fractions = "equal")
+
+print(int)
+
+plot(int)
+
+plotAssemblyEffects(int)
+plotAssemblyEffects(int, env = runif(100))
+plotAssemblyEffects(int, env = as.factor(c(rep(1, 50), rep(2, 50))))
+plotAssemblyEffects(int, trait = runif(10))
+plotAssemblyEffects(int, trait = as.factor(c(rep(1, 5), rep(2, 5))))
+
+}
diff --git a/sjSDM/man/plotAssemblyEffects.Rd b/sjSDM/man/plotAssemblyEffects.Rd
index 9edb7d8..b1fa0cd 100644
--- a/sjSDM/man/plotAssemblyEffects.Rd
+++ b/sjSDM/man/plotAssemblyEffects.Rd
@@ -1,33 +1,28 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/anova.R
+% Please edit documentation in R/internalStructure.R
 \name{plotAssemblyEffects}
 \alias{plotAssemblyEffects}
-\title{Plot Correlations between assembly processes and predictors or traits}
+\title{Plot correlations between assembly processes and predictors or traits}
 \usage{
 plotAssemblyEffects(
   object,
   env = NULL,
   trait = NULL,
-  Rsquared = c("McFadden", "Nagelkerke"),
   cols = c("#A38310", "#B42398", "#20A382")
 )
 }
 \arguments{
-\item{object}{An `anova` object from the `anova.sjSDM` function.}
+\item{object}{An \code{sjSDManova} object from the \code{\link{anova.sjSDM}} function.}
 
-\item{env}{Predictor variable. If `NULL`, assembly processes are plotted against environment, spatial uniqueness, and richness.}
+\item{env}{Predictor variable. If \code{NULL}, assembly processes are plotted against environment, spatial uniqueness, and richness.}
 
 \item{trait}{Trait variable. Plotted against species R-squared for the three processes.}
 
-\item{Rsquared}{Which R-squared should be used: "McFadden" (default) or "Nagelkerke".}
-
 \item{cols}{Colors for the three assembly processes.
 
 Correlation and plots of the three assembly processes (environment, space, and codist) against environmental and spatial uniqueness and richness. The importance of the three assembly processes is measured by the partial R-squared (shown in the internal structure plots).
-
 Importances are available for species and sites. Custom environmental predictors or traits can be specified. Environmental predictors are plotted against site R-squared and traits are plotted against species R-squared.
-
-Regression lines are estimated by 50% quantile regression models.}
+Regression lines are estimated by 50\% quantile regression models.}
 }
 \value{
 A list with the following components:
@@ -37,7 +32,50 @@ A list with the following components:
 \item{codist}{A list of summary tables for env, space, and codist R-squared.}
 }
 \description{
-Plot Correlations between assembly processes and predictors or traits
+Plot correlations between assembly processes and predictors or traits
+}
+\examples{
+library(sjSDM)
+# simulate community:
+community = simulate_SDM(env = 3L, species = 10L, sites = 100L)
+
+Occ <- community$response
+Env <- community$env_weights
+SP <- data.frame(matrix(rnorm(200, 0, 0.3), 100, 2)) # spatial coordinates
+
+
+# fit model:
+model <- sjSDM(Y = Occ, 
+               env = linear(data = Env, formula = ~X1+X2+X3), 
+               spatial = linear(data = SP, formula = ~0+X1*X2), 
+               family=binomial("probit"),
+               iter = 20) # increase iter for real analysis
+
+# Calculate ANOVA for env, space, associations, for details see ?anova.sjSDM
+an = anova(model, samples = 10) # increase iter for real analysis
+
+# Show anova fractions
+plot(an)
+
+# ANOVA tables with different way to handle fractions
+summary(an)
+summary(an, fractions = "discard")
+summary(an, fractions = "proportional")
+summary(an, fractions = "equal")
+
+# Internal structure
+int = internalStructure(an, fractions = "equal")
+
+print(int)
+
+plot(int)
+
+plotAssemblyEffects(int)
+plotAssemblyEffects(int, env = runif(100))
+plotAssemblyEffects(int, env = as.factor(c(rep(1, 50), rep(2, 50))))
+plotAssemblyEffects(int, trait = runif(10))
+plotAssemblyEffects(int, trait = as.factor(c(rep(1, 5), rep(2, 5))))
+
 }
 \references{
 Leibold, M. A., Rudolph, F. J., Blanchet, F. G., De Meester, L., Gravel, D., Hartig, F., ... & Chase, J. M. (2022). The internal structure of metacommunities. *Oikos*, 2022(1).
diff --git a/sjSDM/man/plotInternalStructure.Rd b/sjSDM/man/plotInternalStructure.Rd
deleted file mode 100644
index 643b10d..0000000
--- a/sjSDM/man/plotInternalStructure.Rd
+++ /dev/null
@@ -1,41 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/anova.R
-\name{plotInternalStructure}
-\alias{plotInternalStructure}
-\title{Plot internal metacommunity structure}
-\usage{
-plotInternalStructure(
-  object,
-  Rsquared = c("McFadden", "Nagelkerke"),
-  add_shared = FALSE,
-  env_deviance = NULL,
-  suppress_plotting = FALSE
-)
-}
-\arguments{
-\item{object}{anova object from \code{\link{anova.sjSDM}}}
-
-\item{Rsquared}{which R squared should be used, McFadden or Nagelkerke (McFadden is default)}
-
-\item{add_shared}{split shared components, default is TRUE}
-
-\item{env_deviance}{environmental deviance}
-
-\item{suppress_plotting}{should the plots be suppressed or not.
-
-Plots and returns the internal metacommunity structure of species and sites (see Leibold et al., 2022). 
-Plots were heavily inspired by Leibold et al., 2022}
-}
-\value{
-List with the following components:
-
-
-\item{plots}{ggplot objects for sites and species.}
-\item{data}{List of data.frames with the internal metacommunity structure.}
-}
-\description{
-Plot internal metacommunity structure
-}
-\references{
-Leibold, M. A., Rudolph, F. J., Blanchet, F. G., De Meester, L., Gravel, D., Hartig, F., ... & Chase, J. M. (2022). The internal structure of metacommunities. Oikos, 2022(1).
-}
diff --git a/sjSDM/man/plotsjSDMcoef.Rd b/sjSDM/man/plotsjSDMcoef.Rd
index 3b758f0..ee20969 100644
--- a/sjSDM/man/plotsjSDMcoef.Rd
+++ b/sjSDM/man/plotsjSDMcoef.Rd
@@ -25,22 +25,23 @@ Plotting coefficients returned by sjSDM model.
 This function only for model fitted by linear, fitted by DNN is not yet supported.
 }
 \examples{
-\dontrun{
 library(sjSDM)
 # simulate community:
 com = simulate_SDM(env = 6L, species = 7L, sites = 100L)
 
 # fit model:
-model = sjSDM(Y = com$response,env = com$env_weights, iter = 2L,se = TRUE) 
+model = sjSDM(Y = com$response,env = com$env_weights, iter = 2L, se = TRUE) 
 
-#create a group dataframe for plot
+# normal plot
+plot(model)
+
+# colored by groups
 species=c("sp1","sp2","sp3","sp4","sp5","sp6","sp7")
 group=c("mammal","bird","fish","fish","mammal","amphibian","amphibian")
 group = data.frame(species=species,group=group)
 
 plot(model,group=group)
 }
-}
 \author{
 CAI Wang
 }
diff --git a/sjSDM/man/predict.sjSDM.Rd b/sjSDM/man/predict.sjSDM.Rd
index f73004a..a1f4b94 100644
--- a/sjSDM/man/predict.sjSDM.Rd
+++ b/sjSDM/man/predict.sjSDM.Rd
@@ -35,3 +35,27 @@ Matrix of predictions (sites by species)
 \description{
 Predict from a fitted sjSDM model
 }
+\examples{
+\dontrun{
+
+## Conditional predictions based on focal species
+com = simulate_SDM(sites = 200L)
+## first 100 observations are the training data
+model = sjSDM(com$response[1:100, ], com$env_weights[1:100,])
+## Assume that for the other 100 observations, only the first species is missing 
+## and we want to use the other 4 species to improve the predictions:
+Y_focal = com$response[101:200, ]
+Y_focal[,1] = NA # set to NA because occurrences are unknown
+
+pred_conditional = predict(model, newdata = com$env_weights[101:200,], Y = Y_focal)
+pred_unconditional = predict(model, newdata = com$env_weights[101:200,])[,1]
+
+## Compare performance:
+Metrics::auc(com$response[101:200, 1], pred_conditional)
+Metrics::auc(com$response[101:200, 1], pred_unconditional)
+
+## Conditional predictions are better, however, it only works if occurrences of
+## other species for new sites are known!
+
+}
+}
diff --git a/sjSDM/man/print.sjSDManova.Rd b/sjSDM/man/print.sjSDManova.Rd
index fa6bb25..5fda8e1 100644
--- a/sjSDM/man/print.sjSDManova.Rd
+++ b/sjSDM/man/print.sjSDManova.Rd
@@ -2,18 +2,58 @@
 % Please edit documentation in R/anova.R
 \name{print.sjSDManova}
 \alias{print.sjSDManova}
-\title{Print sjSDM anova}
+\title{Print sjSDM anova object}
 \usage{
 \method{print}{sjSDManova}(x, ...)
 }
 \arguments{
-\item{x}{an object of \code{\link{anova.sjSDM}}}
+\item{x}{an object of type sjSDManova created by \code{\link{anova.sjSDM}}}
 
-\item{...}{optional arguments for compatibility with the generic function, no function implemented}
-}
-\value{
-The above matrix is silently returned
+\item{...}{additional arguments to \code{\link{summary.sjSDManova}}}
 }
 \description{
-Print sjSDM anova
+This is a wrapper for \code{\link{summary.sjSDManova}}, maintained for backwards compatibility - prefer to use summary() instead
+}
+\examples{
+library(sjSDM)
+# simulate community:
+community = simulate_SDM(env = 3L, species = 10L, sites = 100L)
+
+Occ <- community$response
+Env <- community$env_weights
+SP <- data.frame(matrix(rnorm(200, 0, 0.3), 100, 2)) # spatial coordinates
+
+
+# fit model:
+model <- sjSDM(Y = Occ, 
+               env = linear(data = Env, formula = ~X1+X2+X3), 
+               spatial = linear(data = SP, formula = ~0+X1*X2), 
+               family=binomial("probit"),
+               iter = 20) # increase iter for real analysis
+
+# Calculate ANOVA for env, space, associations, for details see ?anova.sjSDM
+an = anova(model, samples = 10) # increase iter for real analysis
+
+# Show anova fractions
+plot(an)
+
+# ANOVA tables with different way to handle fractions
+summary(an)
+summary(an, fractions = "discard")
+summary(an, fractions = "proportional")
+summary(an, fractions = "equal")
+
+# Internal structure
+int = internalStructure(an, fractions = "equal")
+
+print(int)
+
+plot(int)
+
+plotAssemblyEffects(int)
+plotAssemblyEffects(int, env = runif(100))
+plotAssemblyEffects(int, env = as.factor(c(rep(1, 50), rep(2, 50))))
+plotAssemblyEffects(int, trait = runif(10))
+plotAssemblyEffects(int, trait = as.factor(c(rep(1, 5), rep(2, 5))))
+
 }
diff --git a/sjSDM/man/print.sjSDMinternalStruture.Rd b/sjSDM/man/print.sjSDMinternalStruture.Rd
new file mode 100644
index 0000000..280149c
--- /dev/null
+++ b/sjSDM/man/print.sjSDMinternalStruture.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/internalStructure.R
+\name{print.sjSDMinternalStruture}
+\alias{print.sjSDMinternalStruture}
+\title{Print internal structure object}
+\usage{
+\method{print}{sjSDMinternalStruture}(x, ...)
+}
+\arguments{
+\item{x}{object of class sjSDMinternalStruture}
+
+\item{...}{no function}
+}
+\description{
+Print internal structure object
+}
diff --git a/sjSDM/man/residuals.sjSDM.Rd b/sjSDM/man/residuals.sjSDM.Rd
new file mode 100644
index 0000000..29c9401
--- /dev/null
+++ b/sjSDM/man/residuals.sjSDM.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sjSDM.R
+\name{residuals.sjSDM}
+\alias{residuals.sjSDM}
+\title{Residuals for a sjSDM model}
+\usage{
+\method{residuals}{sjSDM}(object, type = "raw", ...)
+}
+\arguments{
+\item{object}{a model fitted by \code{\link{sjSDM}}}
+
+\item{type}{residual type. Currently only supports raw}
+
+\item{...}{further arguments, not supported yet.}
+}
+\value{
+residuals in the format of the provided community matrix
+}
+\description{
+Returns residuals for a fitted sjSDM model
+}
diff --git a/sjSDM/man/sjSDM.Rd b/sjSDM/man/sjSDM.Rd
index aea6d03..fe2d749 100644
--- a/sjSDM/man/sjSDM.Rd
+++ b/sjSDM/man/sjSDM.Rd
@@ -134,7 +134,7 @@ See Pichler and Hartig, 2020 for benchmark results.
 
 \subsection{Supported distributions}{
 
-Currently supported distributions and link functions:
+Currently supported distributions and link functions, which are :
 \itemize{
 \item \code{\link{binomial}}: \code{"probit"} or \code{"logit"}
 \item \code{\link{poisson}}: \code{"log"} 
@@ -187,6 +187,10 @@ com = simulate_SDM(env = 3L, species = 7L, sites = 100L)
 model = sjSDM(Y = com$response,env = com$env_weights, iter = 50L) 
 # increase iter for your own data 
 
+# Default distribution is binomial("probit"). Alternatively, you can use
+# binomial(logit), poisson("log"), "nbinom" (with log, still somewhat 
+# experimental) and gaussian("identity")
+
 coef(model)
 summary(model)
 getCov(model)
@@ -248,7 +252,16 @@ predict(model, newdata = com$env_weights, SP = SPV)
 
 ## Visualize internal meta-community structure
 an = anova(model)
-plot(an, internal=TRUE)
+
+internal = internalStructure(an)
+plot(internal)
+
+## Visualize community assemlby effects 
+
+plotAssemblyEffects(internal)
+
+### see ?anova.sjSDM for mroe details
+
 
 ## non-linear(deep neural network) model
 model = sjSDM(Y = com$response, env = linear(com$env_weights), 
@@ -287,30 +300,10 @@ result = anova(model)
 print(result)
 plot(result)
 
-## visualize meta-community structure
-plot(result, internal=TRUE)
-
-
-
-# Conditional predictions
-## Conditional predictions based on focal species
-com = simulate_SDM(sites = 200L)
-## first 100 observations are the training data
-model = sjSDM(com$response[1:100, ], com$env_weights[1:100,])
-## Assume that for the other 100 observations, only the first species is missing 
-## and we want to use the other 4 species to improve the predictions:
-Y_focal = com$response[101:200, ]
-Y_focal[,1] = NA # set to NA because occurrences are unknown
-
-pred_conditional = predict(model, newdata = com$env_weights[101:200,], Y = Y_focal)
-pred_unconditional = predict(model, newdata = com$env_weights[101:200,])[,1]
-
-## Compare performance:
-Metrics::auc(com$response[101:200, 1], pred_conditional)
-Metrics::auc(com$response[101:200, 1], pred_unconditional)
+## visualize internal meta-community structure
+internal = internalStructure(an)
+plot(internal)
 
-## Conditional predictions are better, however, it only works if occurrences of
-## other species for new sites are known!
 
 
 # Deep neural networks
diff --git a/sjSDM/man/summary.sjSDManova.Rd b/sjSDM/man/summary.sjSDManova.Rd
new file mode 100644
index 0000000..dda2fb4
--- /dev/null
+++ b/sjSDM/man/summary.sjSDManova.Rd
@@ -0,0 +1,81 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/anova.R
+\name{summary.sjSDManova}
+\alias{summary.sjSDManova}
+\title{Summary table of sjSDM anova}
+\usage{
+\method{summary}{sjSDManova}(
+  object,
+  method = c("ANOVA"),
+  fractions = c("all", "discard", "proportional", "equal"),
+  ...
+)
+}
+\arguments{
+\item{object}{an object of \code{\link{anova.sjSDM}}}
+
+\item{method}{method used to calculate the ANOVA}
+
+\item{fractions}{how to handle the shared fractions. See details}
+
+\item{...}{optional arguments for compatibility with the generic function, no function implemented}
+}
+\value{
+The matrix that is printed out is silently returned
+}
+\description{
+The function prints and returns invisible a summary table of an sjSDM ANOVA, created by \code{\link{anova.sjSDM}}
+}
+\details{
+The function returns a ANOVA table with Deviance as well as the pseudo-R2 metrics of Nagelkerke and McFadden
+
+There are four options to handle shared ANOVA fractions, which is variance that can be explained, typically as a result of collinearity, by several of the fractions:
+
+1. "all" returns the shared fractions explicitly
+2. "discard" discards the fractions, as typically in a type II Anova
+3. "proportional" distributes shared fractions proportional to the unique fractions
+4. "equal" distributions shared fractions equally to the unique fractions
+}
+\examples{
+library(sjSDM)
+# simulate community:
+community = simulate_SDM(env = 3L, species = 10L, sites = 100L)
+
+Occ <- community$response
+Env <- community$env_weights
+SP <- data.frame(matrix(rnorm(200, 0, 0.3), 100, 2)) # spatial coordinates
+
+
+# fit model:
+model <- sjSDM(Y = Occ, 
+               env = linear(data = Env, formula = ~X1+X2+X3), 
+               spatial = linear(data = SP, formula = ~0+X1*X2), 
+               family=binomial("probit"),
+               iter = 20) # increase iter for real analysis
+
+# Calculate ANOVA for env, space, associations, for details see ?anova.sjSDM
+an = anova(model, samples = 10) # increase iter for real analysis
+
+# Show anova fractions
+plot(an)
+
+# ANOVA tables with different way to handle fractions
+summary(an)
+summary(an, fractions = "discard")
+summary(an, fractions = "proportional")
+summary(an, fractions = "equal")
+
+# Internal structure
+int = internalStructure(an, fractions = "equal")
+
+print(int)
+
+plot(int)
+
+plotAssemblyEffects(int)
+plotAssemblyEffects(int, env = runif(100))
+plotAssemblyEffects(int, env = as.factor(c(rep(1, 50), rep(2, 50))))
+plotAssemblyEffects(int, trait = runif(10))
+plotAssemblyEffects(int, trait = as.factor(c(rep(1, 5), rep(2, 5))))
+
+}
diff --git a/sjSDM/vignettes/Dependencies.Rmd b/sjSDM/vignettes/Dependencies.Rmd
deleted file mode 100644
index ae19a53..0000000
--- a/sjSDM/vignettes/Dependencies.Rmd
+++ /dev/null
@@ -1,168 +0,0 @@
----
-title: "sjSDM: Help with the installation of dependencies"
-author: "Maximilian Pichler & Florian Hartig, Theoretical Ecology, University of Regensburg"
-date: "`r Sys.Date()`"
-output: 
-  rmarkdown::html_vignette:
-    toc: true
-    toc_depth: 3
-vignette: >
-  %\VignetteIndexEntry{sjSDM: Help with the installation of dependencies}
-  \usepackage[utf8]{inputenc}
-  %\VignetteEngine{knitr::rmarkdown}
-abstract: "This document explains how to install dependencies for the sjSDM package. \n \n \n"
-editor_options: 
-  
-  chunk_output_type: console
----
-
-```{r, echo = F, message = F, purl=FALSE}
-set.seed(123)
-```
-
-
-```{r global_options, include=FALSE, purl=FALSE}
-knitr::opts_chunk$set(fig.width=7, fig.height=4.5, fig.align='center', warning=FALSE, message=FALSE, cache = FALSE, eval = FALSE)
-```
-
-The `r knitr::inline_expr('sjSDM::install_sjSDM()')` function can install automatically all necessary 'python' dependencies but it can fail sometimes because of individual system settings or if other 'python'/'conda' installations get into the way. 
-
-## PyTorch Installation - Before you start:
-
-A few notes before you start with the installation (skip this point if you 
-do not know conda):
-
-* existing 'conda' installations: make sure you have the latest 'conda' or 'miniconda' version and remove unnecessary 'conda' installations.
-* existing 'conda' or 'virtualenv' environments (skip this point if you do not know 'conda'): we currently enforce the usage of a specific environment called 'r-sjsdm',so if you want use a custom environment it should be named 'r-sjsdm'
- 
- 
-## Windows - automatic installation:
- 
-Sometimes the automatic 'miniconda' installation (via 
-`r knitr::inline_expr('sjSDM::install_sjSDM()')` ).doesn't work because of white spaces in the user's name. But you can easily download and install 'conda' on your own:
-
-Download and install the latest ['conda' version](https://www.anaconda.com/download/)
-
-Afterwards run:
-```{r,eval=FALSE, purl=FALSE}
-install_sjSDM(version = c("gpu")) # or "cpu" if you do not have a proper gpu device
-```
-
- 
-Reload the package and run the example, if this doesn't work:
-
-* Restart RStudio
-* Install manually pytorch, see the following section
- 
-  
-## Windows - manual installation:
- 
-Download and install the latest ['conda' version](https://www.anaconda.com/download/)
-
-Open the command window (cmd.exe - hit windows key + r and write cmd)
-
-Run in cmd.exe:
-
-```{bash, eval=FALSE, purl=FALSE}
-conda create --name r-sjsdm python=3.7
-conda activate r-sjsdm
-conda install pytorch torchvision cpuonly -c pytorch # cpu
-conda install pytorch torchvision cudatoolkit=11.3 -c pytorch #gpu
-python -m pip install pyro-ppl torch_optimizer madgrad
-```
- 
- 
-Restart R, try to run the example, and if this doesn't work:
-
-* Restart RStudio
-* See the 'Help and bugs' section
- 
-    
-## Linux - automatic installation:
- 
-Run in R: 
-```{r,eval=FALSE, purl=FALSE}
-install_sjSDM(version = c("gpu")) # or "cpu" if you do not have a proper gpu device
-```
-
-Restart R try to run the example, if this doesn't work:
-
-* Restart RStudio
-* Install manually 'PyTorch', see the following section
-
-  
-  
-## Linux - manual installation:
- 
-We strongly advise to use a 'conda' environment but a virtual environment should also work. The only requirement is that it is named 'r-sjsdm'
-
-Download and install the latest ['conda' version](https://www.anaconda.com/download/)
-
-Open your terminal and run:
-```{bash, eval=FALSE, purl=FALSE}
-conda create --name r-sjsdm python=3.7
-conda activate r-sjsdm
-conda install pytorch torchvision cpuonly -c pytorch # cpu
-conda install pytorch torchvision cudatoolkit=11.3 -c pytorch #gpu
-python -m pip install pyro-ppl torch_optimizer madgrad
-```
- 
- 
-Restart R try to run the example, if this doesn't work:
-
-* Restart RStudio
-* See the 'Help and bugs' section
- 
- 
-## MacOS - automatic installation:
- 
-Run in R: 
-```{r,eval=FALSE, purl=FALSE}
-install_sjSDM()
-```
-
-Restart R try to run the example, if this doesn't work:
-
-* Restart RStudio
-* Install manually 'PyTorch', see the following section
-
- 
- 
-## MacOS - manual installation:
- 
-We strongly advise to use a 'conda' environment but a virtual environment should also work. The only requirement is that it is named 'r-sjsdm'
-
-Download and install the latest conda [conda version](https://www.anaconda.com/download/)
-
-Open your terminal and run:
-```{bash, eval=FALSE, purl=FALSE}
-conda create --name r-sjsdm python=3.7
-conda activate r-sjsdm
-python -m pip install torch torchvision torchaudio 
-python -m pip install pyro-ppl torch_optimizer madgrad
-```
- 
- 
-Restart R try to run the example, if this doesn't work:
-
-* Restart RStudio
-* See the 'Help and bugs' section
-
-## Help and bugs
-
-* ask for help via the [issue tracker](https://github.com/TheoreticalEcology/s-jSDM/issues). Please post the output of  `r knitr::inline_expr('sjSDM::install_diagnostic()')`as a quote. 
-* report bugs with a [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-emple)
- via the [issue tracker](https://github.com/TheoreticalEcology/s-jSDM/issues)
-
-
-
-## Example
-```{r,eval=FALSE}
-library(sjSDM)
-community <- simulate_SDM(sites = 100, species = 10, env = 5)
-Env <- community$env_weights
-Occ <- community$response
-
-model <- sjSDM(Y = Occ, env = linear(Env, ~0+X1*X2 + X3 + X4))
-summary(model)
-```
\ No newline at end of file
diff --git a/sjSDM/vignettes/jSDM-ANOVA.png b/sjSDM/vignettes/jSDM-ANOVA.png
new file mode 100644
index 0000000..b5764a5
Binary files /dev/null and b/sjSDM/vignettes/jSDM-ANOVA.png differ
diff --git a/sjSDM/vignettes/jSDM-structure.png b/sjSDM/vignettes/jSDM-structure.png
new file mode 100644
index 0000000..4611d97
Binary files /dev/null and b/sjSDM/vignettes/jSDM-structure.png differ
diff --git a/sjSDM/vignettes/sjSDM_Introduction.Rmd b/sjSDM/vignettes/sjSDM_Introduction.Rmd
index da17bf9..6ed7125 100644
--- a/sjSDM/vignettes/sjSDM_Introduction.Rmd
+++ b/sjSDM/vignettes/sjSDM_Introduction.Rmd
@@ -3,14 +3,13 @@ title: "sjSDM: Getting started with sjSDM - a scalable joint Species Distributio
 author: "Maximilian Pichler & Florian Hartig, Theoretical Ecology, University of Regensburg"
 date: "`r Sys.Date()`"
 output: 
-
   rmarkdown::html_vignette:
     toc: true
 vignette: >
-  %\VignetteIndexEntry{sjSDM: Getting started with sjSDM - a scalable joint Species Distribution Model}
   %\VignetteEncoding{UTF-8}
+  %\VignetteIndexEntry{sjSDM: Getting started with the sjSDM package}
   %\VignetteEngine{knitr::rmarkdown}
-abstract: "sjSDM is a scalable joint species distribution model based on Monte-Carlo intergration of the joint likelihood. \n \n \n"
+abstract: "A scalable and fast method for estimating joint Species Distribution Models (jSDMs) for big community data, including eDNA data. The package estimates a full (i.e. non-latent) jSDM with different response functions (including the traditional multivariate probit model). As described in Pichler & Hartig (2021) <doi:10.1111/2041-210X.13687>, scalability is achieved by using a Monte Carlo approximation of the joint likelihood implemented via 'PyTorch' and 'reticulate', which can be run on CPUs or GPUs. The package includes support for different response families, the ability to account for spatial autocorrelation, and the option to fit responses via deep neural networks instead of a standard linear predictor. In addition, you can perform variation partitioning (VP) / ANOVA on the fitted models for environmental, spatial, and co-association model components, and further subdivide these three components per species and sites (internal metacommunity structure, see Leibold et al., <doi:10.1111/oik.08618>). \n \n \n"
 editor_options: 
   chunk_output_type: console
 ---
@@ -33,291 +32,200 @@ knitr::opts_chunk$set(
 
 # Getting started
 
-Load the package via
+## About sjSDM
 
-```{r,eval=FALSE}
-library(sjSDM)
-```
+The sjSDM package is an R package for estimating so-called joint species distribution models. A jSDM is a GLMM that models a multivariate (i.e. a many-species) response to the environment, space and a covariance term that models conditional (on the other terms) correlations between the outputs (i.e. species). 
 
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
-── Attaching sjSDM ──────────────────────────────────────── 1.0.4 ──
-✔ torch <environment> 
-✔ torch_optimizer  
-✔ pyro  
-✔ madgrad"
-)
-```
+![image](jSDM-structure.png)
 
-## Installing dependencies
+A big challenge in jSDM implementation is computational speed. The goal of the sjSDM (which stands for "scalable joint species distribution models") is to make jSDM computations fast and scalable. Unlike many other packages, which use a latent-variable approximation to make estimating jSDMs faster, sjSDM fits a full covariance matrix in the likelihood, which is, however, numerically approximated via simulations. 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://www.doi.org/10.1111/2041-210X.13687. 
 
-sjSDM depends on Anaconda which will be installed automatically if it is not available.
+The core code of sjSDM is implemented in Python / PyTorch, which is then wrapped into an R package. In principle, you can also use it stand-alone under Python (see instructions below). Note: for both the R and the python package, python \>= 3.7 and pytorch must be installed (more details below). However, for most users, it will be more convenient to use sjSDM via the sjSDM R package, which also provides a large number of downstream functionalities. 
 
-To install the CPU version, run:
+To get citation info for sjSDM when you use it for your reseach, type  
 
-```{r, eval = F, purl=FALSE}
-install_sjSDM()
+```{r,eval=FALSE}
+citation("sjSDM")
 ```
 
-For the GPU version (only NVIDIA GPUs are supported), run:
+## Installing the R package
 
-```{r, eval = F, purl=FALSE}
-install_sjSDM(method = "gpu")
+sjSDM is distributed via [CRAN](https://cran.rstudio.com/web/packages/sjSDM/index.html). For most users, it will be best to install the package from CRAN
+
+```{r,eval=FALSE}
+install.packages("sjSDM")
 ```
 
-More details on and trouble-shooting for installing the necessary dependencies is explained in the vignette dependencies, call
+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](https://github.com/TheoreticalEcology/s-jSDM/blob/master/sjSDM/vignettes/Dependencies.Rmd). Very briefly, the dependencies can be automatically installed from within R:
 
-```{r,eval = F, purl=FALSE}
-vignette("Dependencies", package = "sjSDM")
+```{r,eval=FALSE}
+sjSDM::install_sjSDM(version = "gpu") # or
+sjSDM::install_sjSDM(version = "cpu")
 ```
 
-## Citing sjSDM
-
-To cite sjSDM in a publication, use
+For advanced users: if you want to install the current (development) version from this repository, run
 
 ```{r,eval=FALSE}
-citation("sjSDM")
+devtools::install_github("https://github.com/TheoreticalEcology/s-jSDM", subdir = "sjSDM", ref = "master")
 ```
 
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
+dependencies should be installed as above. If the installation fails, check out the help of ?install_sjSDM, ?installation_help, and vignette("Dependencies", package = "sjSDM").
 
-To cite sjSDM in publications use:
+1.  Try install_sjSDM()
+2.  New session, if no 'PyTorch not found' appears it should work, otherwise see ?installation_help
+3.  If do not get the pkg to run, create an issue [issue tracker](https://github.com/TheoreticalEcology/s-jSDM/issues) or write an email to maximilian.pichler at ur.de
 
-  Pichler, M. and Hartig, F. (2021), A new joint species distribution model for faster and more accurate inference of species associations
-  from big community data. Methods in Ecology and Evolution. Accepted Author Manuscript. https://doi.org/10.1111/2041-210X.13687
-
-A BibTeX entry for LaTeX users is
+# Working with sjSDM
 
-  @Article{,
-    title = {A new joint species distribution model for faster and more accurate inference of species associations from big community data},
-    author = {Maximilian Pichler and Florian Hartig},
-    journal = {Methods in Ecology and Evolution},
-    year = {2021},
-    doi = {10.1111/2041-210X.13687},
-  }
+sjSDM has a function to create test data: We create a dataset with 3 environmental predictors, 100 species and 100 sites.
 
-"
-)
+```{r}
+library(sjSDM)
+set.seed(42)
+community <- simulate_SDM(sites = 100, species = 10, env = 3, se = TRUE)
 ```
 
-# Working with sjSDM
-
-sjSDM has a function to create test data: We create a dataset with 3 environmental predictors, 5 species and 100 sites.
+For convenience, we extract the community matrix, the environmental matrix and the spatial position
 
-```{r,eval=FALSE}
-com = simulate_SDM(env = 3L, species = 5L, sites = 100L)
+```{r}
+Occ <- community$response
+Env <- community$env_weights
+SP <- matrix(rnorm(200, 0, 0.3), 100, 2) # spatial coordinates
 ```
 
 ## Fit model
 
-The model is then fit with the function sjSDM(). You have to provide predictors (can be also a data.frame) and response as matrices.
+The model is then fit with the function sjSDM(). You have to provide predictors (can be also a data.frame) and response as matrices. If you want/need p-values, you have to set se = TRUE.
 
-```{r,eval=FALSE}
-model = sjSDM(Y = com$response, env = com$env_weights, iter = 100L, se=TRUE)
+```{r}
+model <- sjSDM(Y = Occ, 
+               env = linear(data = Env, formula = ~X1+X2+X3), 
+               spatial = linear(data = SP, formula = ~0+X1*X2), 
+               family=binomial("probit"),
+               se = TRUE)
+summary(model)
 ```
 
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
-Iter: 100/100 100%|██████████| [00:01, 55.40it/s, loss=2.845]
+As you see, for environmental and spatial predictors, you can provide the dependency you want via the standard R formular syntax. E.g. for an interaction between X1 and X2, write 
 
-Calculating standard errors...
-Species: 5/5 100%|██████████| [00:00,  7.02it/s]
-"
-)
+```{r, eval = F}
+env = linear(data = Env, formula = ~ X1*X2,)
 ```
 
-## Interpreting model output
+or for a quadratic effect of X1, but without an intercept, write
 
-Things you can do with the model output
+```{r, eval = F}
+env = linear(data = Env, formula = ~0+ I(X1^2)+X1)
+```
 
-```{r,eval=FALSE}
-coef(model)
+The linear means that you essentially fit a linear regression structure for the respective component. Alternatively, you could also fit a neural network (see later). 
+
+## Interpreting the estimated model coefficients
+
+To get an overall summary of the model, type 
+
+```{r}
 summary(model)
-Rsquared(model)
 ```
 
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
-[[1]]
-           [,1]       [,2]         [,3]       [,4]
-[1,]  0.2123860  0.1510233 -0.921597779 -0.8924292
-[2,]  0.1897056  0.9265483 -0.004749471  1.1936187
-[3,]  0.1862375 -0.4717873  1.039884210  1.9334294
-[4,] -0.2361571  1.3922049  0.119922511 -0.7957588
-[5,]  0.1032345  0.4762956  1.235011458 -0.7303844
-"
-)
-```
+This output will show you the fitted environmental, spatial and covariance parameters (provided that those components were specified). Implemented S3 functions for a model object include
 
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
-Family:  binomial 
+```{r, echo=F}
+coef(model)
+residuals(model)
+Rsquared(model)
+```
 
-LogLik:  -283.4419 
-Regularization loss:  0 
+### Environmental effects
 
-Species-species correlation matrix: 
+The environmental effects are displayed in the summary() table. To get a visual plot, you can use 
 
-	sp1	 1.0000				
-	sp2	 0.2110	 1.0000			
-	sp3	 0.1020	 0.6660	 1.0000		
-	sp4	-0.2530	 0.2630	 0.4920	 1.0000	
-	sp5	 0.0130	 0.3430	 0.2780	-0.1490	 1.0000
-
-
-
-                Estimate  Std.Err Z value Pr(>|z|)    
-sp1 (Intercept)  0.21239  0.19224    1.10  0.26924    
-sp1 X1           0.15102  0.35645    0.42  0.67180    
-sp1 X2          -0.92160  0.33932   -2.72  0.00661 ** 
-sp1 X3          -0.89243  0.32367   -2.76  0.00583 ** 
-sp2 (Intercept)  0.18971  0.19152    0.99  0.32191    
-sp2 X1           0.92655  0.37360    2.48  0.01314 *  
-sp2 X2          -0.00475  0.34082   -0.01  0.98888    
-sp2 X3           1.19362  0.32748    3.64  0.00027 ***
-sp3 (Intercept)  0.18624  0.25774    0.72  0.46995    
-sp3 X1          -0.47179  0.49147   -0.96  0.33707    
-sp3 X2           1.03988  0.46161    2.25  0.02428 *  
-sp3 X3           1.93343  0.43018    4.49    7e-06 ***
-sp4 (Intercept) -0.23616  0.24017   -0.98  0.32547    
-sp4 X1           1.39220  0.45671    3.05  0.00230 ** 
-sp4 X2           0.11992  0.43325    0.28  0.78193    
-sp4 X3          -0.79576  0.39868   -2.00  0.04594 *  
-sp5 (Intercept)  0.10323  0.19424    0.53  0.59510    
-sp5 X1           0.47630  0.35038    1.36  0.17403    
-sp5 X2           1.23501  0.36023    3.43  0.00061 ***
-sp5 X3          -0.73038  0.33345   -2.19  0.02850 *  
----
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-"
-)
+```{r}
+plot(model)
 ```
 
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
-[1] 0.1833771
-"
-)
+If you have species groups in your data, a bit nicer plot is produced by
+
+```{r}
+species=c("sp1","sp2","sp3","sp4","sp5","sp6","sp7","sp8","sp9","sp10")
+group=c("mammal","bird","fish","fish","mammal","amphibian","amphibian", "fish","mammal","amphibian")
+group = data.frame(species=species,group=group)
+plot(model,group=group)
 ```
 
-### Get species-species association matrix
 
-```{r,eval=FALSE}
+### Spatial effects 
+
+The spatial effects are displayed in the summary() table. Currently, there are not additional options implemented to visualize the spatial effects
+
+### Species-species covariance (associations)
+
+The species-species associations are displayed in the summary() table. You can extract them as follows:
+
+```{r}
 association = getCor(model)
-print(association)
+sp = x = 1:length(species)
+image(x = sp, y = sp, z= association, )
 ```
 
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
-            [,1]      [,2]      [,3]       [,4]        [,5]
-[1,]  1.00000000 0.2107329 0.1017986 -0.2527094  0.01321193
-[2,]  0.21073294 1.0000000 0.6657167  0.2625644  0.34263917
-[3,]  0.10179862 0.6657167 1.0000000  0.4915384  0.27815582
-[4,] -0.25270941 0.2625644 0.4915384  1.0000000 -0.14870448
-[5,]  0.01321191 0.3426392 0.2781559 -0.1487045  1.00000000
-"
-)
-```
+## ANOVA and variation partitioning
 
-## Adding quadratic predictors and interactions
+### Global ANOVA
 
-sjSDM supports the R formula syntax.
+![image](jSDM-ANOVA.png)
 
-E.g. interaction with intercept (if you want/need p-values, you have to set se = TRUE)
+ANOVA (Type II) will calculate the amount of variation explained by the three components (environment, associations, and space):
 
-```{r,eval=FALSE}
-model = sjSDM(Y = com$response, 
-              env = linear(data = com$env_weights, formula = ~ X1*X2,), 
-              iter = 50L, se = TRUE)
-summary(model)
+```{r}
+an = anova(model)
 ```
 
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
-Family:  binomial 
+Amount of variance explained is measured by 
 
-LogLik:  -302.3292 
-Regularization loss:  0 
 
-Species-species correlation matrix: 
 
-	sp1	 1.0000				
-	sp2	 0.0290	 1.0000			
-	sp3	-0.0630	 0.7200	 1.0000		
-	sp4	-0.1850	 0.1670	 0.3530	 1.0000	
-	sp5	 0.0730	 0.2410	 0.1370	-0.1400	 1.0000
+```{r}
+plot(an)
+```
 
 
+print(an)
+
 
-                Estimate Std.Err Z value Pr(>|z|)    
-sp1 (Intercept)   0.1986  0.1555    1.28  0.20144    
-sp1 X1           -0.0040  0.2802   -0.01  0.98861    
-sp1 X2           -0.8405  0.2730   -3.08  0.00208 ** 
-sp1 X1:X2         0.1894  0.4785    0.40  0.69227    
-sp2 (Intercept)   0.2243  0.2261    0.99  0.32113    
-sp2 X1            1.0708  0.4270    2.51  0.01215 *  
-sp2 X2            0.1217  0.3986    0.31  0.76008    
-sp2 X1:X2         0.8206  0.7316    1.12  0.26200    
-sp3 (Intercept)   0.1152  0.2414    0.48  0.63333    
-sp3 X1           -0.2812  0.4482   -0.63  0.53038    
-sp3 X2            1.1449  0.4289    2.67  0.00761 ** 
-sp3 X1:X2        -0.1565  0.7788   -0.20  0.84073    
-sp4 (Intercept)  -0.1829  0.1686   -1.09  0.27783    
-sp4 X1            0.8784  0.3126    2.81  0.00495 ** 
-sp4 X2            0.0364  0.2886    0.13  0.89955    
-sp4 X1:X2        -0.5055  0.5256   -0.96  0.33618    
-sp5 (Intercept)   0.0433  0.1542    0.28  0.77874    
-sp5 X1            0.3540  0.2772    1.28  0.20148    
-sp5 X2            0.9742  0.2885    3.38  0.00074 ***
-sp5 X1:X2        -1.0459  0.5089   -2.06  0.03987 *  
----
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-"
-)
-```
 
-E.g. quadratic effect without intercept:
+### Internal metacommunity structure
+
+We can calculate and visualize the internal metacommunity structure for species and sites (only for models with space):
 
 ```{r,eval=FALSE}
-model = update(model, env_formula = ~0+ I(X1^2)+X1)
-summary(model)
+results = plotInternalStructure(an)
+print(results$data$Species)
 ```
 
+![](internal.png)
+
 ```{r,  eval = TRUE,echo = FALSE, results = TRUE}
 cat("
-Family:  binomial 
-
-LogLik:  -318.0733 
-Regularization loss:  0 
-
-Species-species correlation matrix: 
-
-	sp1	 1.0000				
-	sp2	-0.0230	 1.0000			
-	sp3	-0.1680	 0.7150	 1.0000		
-	sp4	-0.1660	 0.1260	 0.3010	 1.0000	
-	sp5	-0.1470	 0.2070	 0.2380	-0.0670	 1.0000
-
-
-
-            Estimate Std.Err Z value Pr(>|z|)   
-sp1 I(X1^2)   0.4843  0.3651    1.33   0.1846   
-sp1 X1        0.1387  0.2860    0.49   0.6276   
-sp2 I(X1^2)   0.0865  0.5380    0.16   0.8723   
-sp2 X1        1.0220  0.4157    2.46   0.0139 * 
-sp3 I(X1^2)  -0.3807  0.5516   -0.69   0.4901   
-sp3 X1       -0.4502  0.4267   -1.05   0.2915   
-sp4 I(X1^2)  -0.2723  0.4015   -0.68   0.4976   
-sp4 X1        0.8490  0.3105    2.73   0.0063 **
-sp5 I(X1^2)   0.6131  0.3653    1.68   0.0933 . 
-sp5 X1        0.1927  0.2886    0.67   0.5044   
----
-Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
+         env        spa     codist         r2
+1 0.20700207 0.06585441 0.06470110 0.06751152
+2 0.08447389 0.04961366 0.04666684 0.03615088
+3 0.12547761 0.05222821 0.00932909 0.03740698
+4 0.25344941 0.22698509 0.00000000 0.09249987
+5 0.14541311 0.05033103 0.00000000 0.03742124
 "
 )
 ```
 
+
+
+## Predictions 
+
+
+
+
+
+# Advanced topics
+
 ## Fitting other distributions (e.g. species frequencies)
 
 sjSDM supports other responses than presence-absence data: Simulate non-presence-absence data:
@@ -329,7 +237,7 @@ X = com$env_weights
 Y = com$response
 ```
 
-## Poisson
+### Poisson
 
 ```{r,eval=FALSE}
 model = sjSDM(Y, env = linear(X, ~.), se = TRUE, 
@@ -381,7 +289,7 @@ Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’
 )
 ```
 
-## Negative Binomial
+### Negative Binomial
 
 ```{r,eval=FALSE}
 model = sjSDM(Y, env = linear(X, ~.), se = TRUE, iter = 50L, family = "nbinom")
@@ -434,7 +342,7 @@ Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’
 )
 ```
 
-## Normal (gaussian)
+### Normal (gaussian)
 
 ```{r,eval=FALSE}
 model = sjSDM(log(Y+0.01), env = linear(X, ~.), se = TRUE, 
@@ -486,9 +394,7 @@ Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’
 )
 ```
 
-# Handling spatial autocorrelation
-
-## Simulating spatial data
+## Modifying the spatial component 
 
 jSDMs account for correlation between species within communities (sites), in real datasets, however, communities (sites) are often also correlated (== spatial autocorrelation). Usually conditional autoregressive (CAR) models are used to account for the spatial autocorrelation in the residuals, which we, however, do not support yet. A similar approach is to condition the model on space, which we can do by using space as predictors.
 
@@ -519,8 +425,6 @@ Ysp = Y + t(spatialResiduals)
 Y = ifelse(Ysp < 0, 0, 1) # multivariate probit model
 ```
 
-## Fitting spatial jSDM
-
 There are three options to condition our model on space:
 
 ### Using Moran's eigenvector map predictors
@@ -702,88 +606,10 @@ X3          -0.1201278 -0.2378207  0.6429102 -0.7199039  0.08894651
 )
 ```
 
-# VENN diagram and community analyses
 
-## Anova
+## Adjusting regularization parameters
 
-ANOVA (Type II) will separate the three components (environment, associations, and space):
-
-```{r,eval=FALSE}
-model = sjSDM(Y, 
-              env = linear(X, ~.), 
-              spatial = linear(XYcoords, ~0+XX+YY+XX:YY+I(XX^2)+I(YY^2)), 
-              iter = 100L)
-
-an = anova(model)
-print(an)
-plot(an)
-```
-
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
-Iter: 100/100 100%|██████████| [00:01, 64.31it/s, loss=2.819]
-Iter: 100/100 100%|██████████| [00:01, 58.37it/s, loss=3.373]
-Iter: 100/100 100%|██████████| [00:01, 64.17it/s, loss=3.099]
-Iter: 100/100 100%|██████████| [00:01, 56.71it/s, loss=2.731]
-Iter: 100/100 100%|██████████| [00:01, 58.89it/s, loss=2.646]
-Iter: 100/100 100%|██████████| [00:01, 51.46it/s, loss=2.984]
-Iter: 100/100 100%|██████████| [00:02, 48.46it/s, loss=0.121]
-"
-)
-```
-
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
-Analysis of Deviance Table
-
-         Deviance Residual deviance R2 Nagelkerke R2 McFadden
-Abiotic 135.46469         534.13894       0.74196      0.1954
-Biotic   24.48839         645.11523       0.21720      0.0353
-Spatial  91.44113         578.16250       0.59925      0.1319
-Full    187.83630         481.76733       0.84716      0.2710
-"
-)
-```
-
-![](VENN.png)
-
-## Internal metacommunity structure
-
-We can calculate and visualize the internal metacommunity structure for species and sites (only for models with space):
-
-```{r,eval=FALSE}
-results = plotInternalStructure(an)
-print(results$data$Species)
-```
-
-![](internal.png)
-
-```{r,  eval = TRUE,echo = FALSE, results = TRUE}
-cat("
-         env        spa     codist         r2
-1 0.20700207 0.06585441 0.06470110 0.06751152
-2 0.08447389 0.04961366 0.04666684 0.03615088
-3 0.12547761 0.05222821 0.00932909 0.03740698
-4 0.25344941 0.22698509 0.00000000 0.09249987
-5 0.14541311 0.05033103 0.00000000 0.03742124
-"
-)
-```
-
-## Importance
-
-We can also partition the effects of the different components. This function will be deprecated in future. Please use `plot(anova(model), internal = TRUE)`
-
-```{r,eval=FALSE}
-imp = importance(model)
-plot(imp)
-```
-
-![](importance.png)
-
-# Predictive modeling 
-
-## Regularization on abiotic coefficients
+### Regularization on abiotic coefficients
 
 sjSDM supports l1 (lasso) and l2 (ridge) regularization: \* alpha is the weighting between lasso and ridge \* alpha = 0.0 corresponds to pure lasso \* alpha = 1.0 corresponds to pure ridge
 
@@ -822,7 +648,7 @@ I(X1^2) 0.004259373 -0.01279781 -0.00810579 0.001065778 -0.007236504
 ```
 
 
-## Regularization on species-species associations
+### Regularization on species-species associations
 
 We can do the same for the species associations:
 
@@ -862,7 +688,7 @@ I(X1^2) -0.002930772 -0.01018739 0.0004073409 -0.01492386 -0.0006310232
 ```
 
 
-## Regularization on the spatial model:
+### Regularization on the spatial model:
 
 ```{r,eval=FALSE}
 
@@ -975,3 +801,136 @@ Plot the training history:
 ```{r,eval=FALSE}
 plot(model)
 ```
+
+# Installation problems
+
+The `r knitr::inline_expr('sjSDM::install_sjSDM()')` function can install automatically all necessary 'python' dependencies but it can fail sometimes because of individual system settings or if other 'python'/'conda' installations get into the way. 
+
+## PyTorch Installation - Before you start:
+
+A few notes before you start with the installation (skip this point if you 
+do not know conda):
+
+* existing 'conda' installations: make sure you have the latest 'conda' or 'miniconda' version and remove unnecessary 'conda' installations.
+* existing 'conda' or 'virtualenv' environments (skip this point if you do not know 'conda'): we currently enforce the usage of a specific environment called 'r-sjsdm',so if you want use a custom environment it should be named 'r-sjsdm'
+ 
+ 
+## Windows - automatic installation:
+ 
+Sometimes the automatic 'miniconda' installation (via 
+`r knitr::inline_expr('sjSDM::install_sjSDM()')` ).doesn't work because of white spaces in the user's name. But you can easily download and install 'conda' on your own:
+
+Download and install the latest ['conda' version](https://www.anaconda.com/download/)
+
+Afterwards run:
+```{r,eval=FALSE, purl=FALSE}
+install_sjSDM(version = c("gpu")) # or "cpu" if you do not have a proper gpu device
+```
+
+ 
+Reload the package and run the example, if this doesn't work:
+
+* Restart RStudio
+* Install manually pytorch, see the following section
+ 
+  
+## Windows - manual installation:
+ 
+Download and install the latest ['conda' version](https://www.anaconda.com/download/)
+
+Open the command window (cmd.exe - hit windows key + r and write cmd)
+
+Run in cmd.exe:
+
+```{bash, eval=FALSE, purl=FALSE}
+conda create --name r-sjsdm python=3.7
+conda activate r-sjsdm
+conda install pytorch torchvision cpuonly -c pytorch # cpu
+conda install pytorch torchvision cudatoolkit=11.3 -c pytorch #gpu
+python -m pip install pyro-ppl torch_optimizer madgrad
+```
+ 
+ 
+Restart R, try to run the example, and if this doesn't work:
+
+* Restart RStudio
+* See the 'Help and bugs' section
+ 
+    
+## Linux - automatic installation:
+ 
+Run in R: 
+```{r,eval=FALSE, purl=FALSE}
+install_sjSDM(version = c("gpu")) # or "cpu" if you do not have a proper gpu device
+```
+
+Restart R try to run the example, if this doesn't work:
+
+* Restart RStudio
+* Install manually 'PyTorch', see the following section
+
+  
+  
+## Linux - manual installation:
+ 
+We strongly advise to use a 'conda' environment but a virtual environment should also work. The only requirement is that it is named 'r-sjsdm'
+
+Download and install the latest ['conda' version](https://www.anaconda.com/download/)
+
+Open your terminal and run:
+```{bash, eval=FALSE, purl=FALSE}
+conda create --name r-sjsdm python=3.7
+conda activate r-sjsdm
+conda install pytorch torchvision cpuonly -c pytorch # cpu
+conda install pytorch torchvision cudatoolkit=11.3 -c pytorch #gpu
+python -m pip install pyro-ppl torch_optimizer madgrad
+```
+ 
+ 
+Restart R try to run the example, if this doesn't work:
+
+* Restart RStudio
+* See the 'Help and bugs' section
+ 
+ 
+## MacOS - automatic installation:
+ 
+Run in R: 
+```{r,eval=FALSE, purl=FALSE}
+install_sjSDM()
+```
+
+Restart R try to run the example, if this doesn't work:
+
+* Restart RStudio
+* Install manually 'PyTorch', see the following section
+
+ 
+ 
+## MacOS - manual installation:
+ 
+We strongly advise to use a 'conda' environment but a virtual environment should also work. The only requirement is that it is named 'r-sjsdm'
+
+Download and install the latest conda [conda version](https://www.anaconda.com/download/)
+
+Open your terminal and run:
+```{bash, eval=FALSE, purl=FALSE}
+conda create --name r-sjsdm python=3.7
+conda activate r-sjsdm
+python -m pip install torch torchvision torchaudio 
+python -m pip install pyro-ppl torch_optimizer madgrad
+```
+ 
+ 
+Restart R try to run the example, if this doesn't work:
+
+* Restart RStudio
+* See the 'Help and bugs' section
+
+## Help and bugs
+
+* ask for help via the [issue tracker](https://github.com/TheoreticalEcology/s-jSDM/issues). Please post the output of  `r knitr::inline_expr('sjSDM::install_diagnostic()')`as a quote. 
+* report bugs with a [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-emple)
+ via the [issue tracker](https://github.com/TheoreticalEcology/s-jSDM/issues)
+
+
diff --git a/sjSDM/vignettes/sjSDM_MSDM.Rmd b/sjSDM/vignettes/sjSDM_MSDM.Rmd
index 6f65b04..819e7e7 100644
--- a/sjSDM/vignettes/sjSDM_MSDM.Rmd
+++ b/sjSDM/vignettes/sjSDM_MSDM.Rmd
@@ -6,8 +6,8 @@ output:
   rmarkdown::html_vignette:
     toc: true
 vignette: >
-  %\VignetteIndexEntry{sjSDM: Multispecies distribution modeling, Deep Neural Networks, and explainable AI}
   %\VignetteEncoding{UTF-8}
+  %\VignetteIndexEntry{sjSDM: Multispecies distribution modeling, Deep Neural Networks, and explainable AI}
   %\VignetteEngine{knitr::rmarkdown}
 abstract: "sjSDM can be used to estimate multispecies distribution models based on DNN and we can use explainable AI (xAI) to study these black-box models  \n \n \n"
 editor_options: 


Test patch:
diff --git a/sjSDM/tests/testthat/test-1-anova.R b/sjSDM/tests/testthat/test-1-anova.R
index 1467e43..ae323c5 100644
--- a/sjSDM/tests/testthat/test-1-anova.R
+++ b/sjSDM/tests/testthat/test-1-anova.R
@@ -146,9 +146,9 @@ test_model = function(occ = NULL, env, spatial=NULL, biotic = bioticStruct(),
     gl2 = stats::glm(Y~1, family = binomial("logit"))
     R0 = 1 - stats::logLik(gl1)/stats::logLik(gl2)
     m = sjSDM(matrix(Y, ncol = 1L), data.frame(X = X), biotic = bioticStruct(diag = TRUE),sampling = 500L, family = binomial("logit"))
-    R1 = mean(replicate(500, {Rsquared(m)}))
+    R1 = mean(replicate(10, {Rsquared(m)}))
     
-    testthat::expect_true(abs(R0 - R1) < 0.03)
+    testthat::expect_true(abs(R0 - R1) < 0.1)
     
     ## Poisson ##
     X = runif(500)
@@ -159,8 +159,8 @@ test_model = function(occ = NULL, env, spatial=NULL, biotic = bioticStruct(),
     m = sjSDM(matrix(Y, ncol = 1L), linear(data.frame(X = X), ~1+X), 
               biotic = bioticStruct(diag = TRUE),
               learning_rate = 0.1,sampling = 10L, family = poisson(), control = sjSDMControl(RMSprop(weight_decay = 0.0)))
-    R1 = mean(replicate(500, {Rsquared(m)}))
-    testthat::expect_true(abs(R0 - R1) < 0.08)
+    R1 = mean(replicate(10, {Rsquared(m)}))
+    testthat::expect_true(abs(R0 - R1) < 0.1)
     
   })