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

Issue:
error in function sjSDM.tune
Hi,

Thanks for doing some improvements and updates! I just tried this new function but it seems something is not working properly. Could you please check this?

**This is what I try to run:**
tune2 = sjSDM_cv(env = scale(as.matrix(Env)), Y = as.matrix(Occ), 
                learning_rate = 0.01, iter = 100L, CV = 5,
                tune_steps = 40,
                lambda_cov = seq(0, 0.1, 0.01),
                lambda_coef = seq(0, 0.1, 0.01),
                alpha_cov = seq(0, 1, 0.05),
                alpha_coef = seq(0, 1, 0.05),
                alpha_spatial = seq(0, 1, 0.05),
                sampling = 100L,
                biotic =  bioticStruct(df=dim(Occ)[2]),
                spatial = linear(scale(SPV), ~0+.),
                lambda_spatial = 2^seq(-10, -0.5, length.out = 20),
                step_size = 5L,
                blocks = 3L,
                family=binomial("probit"),
                n_cores = 8)
#Done, tune2 object was plotted.

model2=sjSDM.tune(tune2)

**Message:**
> model2=sjSDM.tune(tune2)
Error in (function (Y = NULL, env = NULL, biotic = bioticStruct(), spatial = NULL,  : 
  unused arguments (family.family = "binomial", family.link = "probit", family.linkfun = function (mu) 
qnorm(mu), family.linkinv = function (eta) 
{
    thresh <- -qnorm(.Machine$double.eps)
    eta <- pmin(pmax(eta, -thresh), thresh)
    pnorm(eta)
}, family.variance = function (mu) 
mu * (1 - mu), family.dev.resids = function (y, mu, wt) 
.Call(C_binomial_dev_resids, y, mu, wt), family.aic = function (y, n, mu, wt, dev) 
{
    m <- if (any(n > 1)) n else wt
    -2 * sum(ifelse(m > 0, (wt/m), 0) * dbinom(round(m * y), round(m), mu, log = TRUE))
}, family.mu.eta = function (eta) 
pmax(dnorm(eta), .Machine$double.eps), family.initialize = {
    if (NCOL(y) == 1) {
        if (is.factor(y)) y <- y != levels(y)[1]
        n <- rep.int(1, nobs)
        y[weights == 0] <- 0
        if (any(y < 0 | y > 1)) stop("y values must be 0 <= y <= 1")
        mustart <- (weights * y + 0.5)/(weights + 1)
        m <- weights * y
        if ("binomial" == "binomial" && any(abs(m - round(m)) >


Many thanks!
Best,

Máté
Rsquared()
Hi,

I try to use the Rsquared function but it somehow fails. How could I make it work?

> model2$spatial
~0 + .
<environment: 0x7fbd9b72ecf0>
> R2 = Rsquared(model2)
Error in eval(mf[4]$spatial_formula, envir = parent.env(environment())) : 
  object 'spatial_formula' not found


Many thanks for your help in advance!
/Máté


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/automatic-installation.yaml b/.github/workflows/automatic-installation.yaml
new file mode 100644
index 0000000..1072b43
--- /dev/null
+++ b/.github/workflows/automatic-installation.yaml
@@ -0,0 +1,109 @@
+# For help debugging build failures open an issue on the RStudio community with the 'github-actions' tag.
+# https://community.rstudio.com/new-topic?category=Package%20development&tags=github-actions
+on:
+  schedule:
+  - cron: '0 0 * * *'
+  push:
+    branches:
+      - main
+      - master
+  pull_request:
+    branches:
+      - main
+      - master
+
+name: Automatic-installation
+
+jobs:
+  Automatic-installation:
+    runs-on: ${{ matrix.config.os }}
+
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: windows-latest, r: 'release'}
+          - {os: ubuntu-20.04, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+
+    env:
+      R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
+      RSPM: ${{ matrix.config.rspm }}
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      KMP_DUPLICATE_LIB_OK: TRUE
+    steps:
+      - uses: actions/checkout@v2
+
+      - uses: r-lib/actions/setup-r@v1
+        with:
+          r-version: ${{ matrix.config.r }}
+
+      - uses: r-lib/actions/setup-pandoc@v1
+      
+      - uses: r-lib/actions/setup-tinytex@master
+      
+      - name: Install makeindex
+        run: tlmgr install makeindex
+
+      - name: Query dependencies
+        run: |
+          setwd("sjSDM")
+          install.packages("pak", repos = "https://r-lib.github.io/p/pak/dev/")
+          saveRDS(pak::pkg_deps("local::.", dependencies = TRUE), "../.github/r-depends.rds")
+        shell: Rscript {0}
+
+      - name: Cache R packages
+        if: runner.os != 'Windows'
+        uses: actions/cache@v2
+        with:
+          path: ${{ env.R_LIBS_USER }}
+          key: ${{ matrix.config.os }}-${{ steps.install-r.outputs.installed-r-version }}-1-${{ hashFiles('.github/r-depends.rds') }}
+          restore-keys: ${{ matrix.config.os }}-${{ steps.install-r.outputs.installed-r-version }}-1-
+
+      - name: Install system dependencies
+        if: runner.os == 'Linux'
+        run: |
+          setwd("sjSDM")
+          pak::local_system_requirements(execute = TRUE)
+          pak::pkg_system_requirements("rcmdcheck", execute = TRUE)
+        shell: Rscript {0}
+        
+      - name: Install dependencies
+        run: |
+          setwd("sjSDM")
+          pak::local_install_dev_deps(upgrade = TRUE)
+          pak::pkg_install("devtools")
+        shell: Rscript {0}
+        
+      - name: Install Pkg
+        run: |
+          cd sjSDM
+          Rscript -e "pak::local_install()"
+          Rscript -e "pak::pkg_install('rstudio/reticulate')"
+          Rscript -e "sjSDM::install_sjSDM(version='cpu')"
+          
+      - name: Run-sjSDM
+        env:
+          _R_CHECK_CRAN_INCOMING_REMOTE_: false
+          CI: true
+        run: |
+          setwd("sjSDM")
+          options(crayon.enabled = TRUE)
+          library(sjSDM)
+          com = simulate_SDM(env = 3L, species = 7L, sites = 10L)
+          model = sjSDM(Y = com$response,env = com$env_weights, iter = 2L)
+        shell: Rscript {0}
+        
+      - name: testthat output
+        if: always()
+        run: find check -name 'testthat.Rout*' -exec cat '{}' \; || true
+        shell: bash
+        
+        
+      - name: Upload check results
+        if: failure()
+        uses: actions/upload-artifact@main
+        with:
+          name: ${{ matrix.config.os }}-r${{ matrix.config.r }}-results
+          path: check
diff --git a/.gitignore b/.gitignore
index 2ae58a3..28bf461 100644
--- a/.gitignore
+++ b/.gitignore
@@ -24,3 +24,4 @@ sjSDM/vignettes/sjSDM.R
 sjSDM.Rproj
 .DS_Store
 .DS_Store
+.DS_Store
diff --git a/Code/Theory/residulas.R b/Code/Theory/residulas.R
new file mode 100644
index 0000000..9015171
--- /dev/null
+++ b/Code/Theory/residulas.R
@@ -0,0 +1,50 @@
+library(sjSDM)
+library(DHARMa)
+library(mltools)
+library(data.table)
+
+
+## simulate community:
+com = simulate_SDM(env = 2L, species = 7L, sites = 100L)
+
+## fit model:
+model = sjSDM(Y = com$response,
+              env = linear(com$env_weights, ~.), 
+              iter = 50L,
+              family = binomial("probit"))
+# increase iter for your own data 
+sims = simulate(model, nsim = 200L)
+dim(sims)
+
+sp = 1
+dim(sims)
+sim = 
+  createDHARMa(t(array(sims, c(200, 100*7))), 
+              c(com$response), 
+              c(predict(model)),
+              integerResponse = TRUE, method = 'traditional')
+plot(sim)
+
+
+
+vs = 
+sapply(1:100, function(i) {
+  empirical_cdf(as.data.table(sims[,i,]) + matrix(runif(200*7, -0.5, 0.5), 200, 7),
+              ubounds = as.data.table(com$response[i,]+runif(7,-0.5,0.5)))$CDF
+})
+
+
+sim$scaledResiduals = c(vs)
+sim$observedResponse = c(com$response)
+sim$fittedPredictedResponse = c(predict(model))
+sim$nObs = 700
+plot(sim)
+
+
+# for (i in 1:n) {
+#   minSim <- mean(simulations[i, ] < observed[i])
+#   maxSim <- mean(simulations[i, ] <= observed[i])
+#   if (minSim == maxSim) 
+#     scaledResiduals[i] = minSim
+#   else scaledResiduals[i] = runif(1, minSim, maxSim)
+# }
\ No newline at end of file
diff --git a/Code/sjSDM with random effects/EM sjSDM.ipynb b/Code/sjSDM with random effects/EM sjSDM.ipynb
index 6376cbd..af781d4 100644
--- a/Code/sjSDM with random effects/EM sjSDM.ipynb
+++ b/Code/sjSDM with random effects/EM sjSDM.ipynb
@@ -1,14 +1,28 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "1a375515",
+   "metadata": {},
+   "source": [
+    "## Playground for random effect integration in sjSDM\n",
+    "Dependencies + sjSDM simulation including:\n",
+    "\n",
+    "* random intercepts per sites\n",
+    "* random intercepts per sites and species\n",
+    "* random intercepts following CAR model"
+   ]
+  },
   {
    "cell_type": "code",
-   "execution_count": 697,
+   "execution_count": 1,
    "id": "8f34bc58-a745-4b0a-b6ac-ff46115fe65e",
    "metadata": {},
    "outputs": [],
    "source": [
     "import torch\n",
     "import numpy as np\n",
+    "import numpy as np, numpy.linalg\n",
     "import matplotlib.pyplot as plt\n",
     "from dataclasses import dataclass\n",
     "import statsmodels.api as sm\n",
@@ -17,6 +31,7 @@
     "from scipy.optimize import minimize\n",
     "from dataclasses import dataclass\n",
     "import numpy as np\n",
+    "\n",
     "@dataclass\n",
     "class Simulation:\n",
     "    X: np.ndarray\n",
@@ -37,8 +52,6 @@
     "    correlation[covariance == 0] = 0\n",
     "    return correlation\n",
     "\n",
-    "import numpy as np,numpy.linalg\n",
-    "\n",
     "def _getAplus(A):\n",
     "    eigval, eigvec = np.linalg.eig(A)\n",
     "    Q = np.matrix(eigvec)\n",
@@ -68,6 +81,15 @@
     "        Yk = _getPu(Xk, W=W)\n",
     "    return Yk\n",
     "\n",
+    "def simulate_CAR(sides = 10, l = 0.7):\n",
+    "    row_coords = np.tile(np.arange(0, sides), sides).reshape([-1, 1])\n",
+    "    col_coords = np.tile(np.arange(0, sides), sides).reshape([-1, 1])\n",
+    "    d = np.concatenate([row_coords, col_coords], axis= 1)\n",
+    "    D = np.sum((d[:, np.newaxis, :] - d[np.newaxis, :, :]) ** 2, axis = -1)\n",
+    "    Sigma = nearPD(np.exp(-l*D))+np.eye(sides**2)*0.01\n",
+    "    re= np.random.multivariate_normal(np.zeros([sides**2]), cov=Sigma,size=[1])\n",
+    "    return D, Sigma, re\n",
+    "\n",
     "def simulate(N=100, E=2, SP=3, G = 20, scale_re = 0.5, CAR=False):\n",
     "    if CAR is not True:\n",
     "        re = np.random.normal(0, scale_re, [G,1])\n",
@@ -87,35 +109,20 @@
     "    YY = np.concatenate([np.random.multivariate_normal(Y[i,:], cov = SS, size=[1,]) for i in range(N)], 0)\n",
     "    YY = YY > 0\n",
     "    YY = YY.astype(np.float64)\n",
-    "    return Simulation(X, YY, SS, W, scale_re, g, G, re, D)\n",
-    "\n",
-    "def simulate_CAR(sides = 10, l = 0.7):\n",
-    "    row_coords = np.tile(np.arange(0, sides), sides).reshape([-1, 1])\n",
-    "    col_coords = np.tile(np.arange(0, sides), sides).reshape([-1, 1])\n",
-    "    d = np.concatenate([row_coords, col_coords], axis= 1)\n",
-    "    D = np.sum((d[:, np.newaxis, :] - d[np.newaxis, :, :]) ** 2, axis = -1)\n",
-    "    Sigma = nearPD(np.exp(-l*D))+np.eye(sides**2)*0.01\n",
-    "    re= np.random.multivariate_normal(np.zeros([sides**2]), cov=Sigma,size=[1])\n",
-    "    return D, Sigma, re"
+    "    return Simulation(X, YY, SS, W, scale_re, g, G, re, D)"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": 375,
-   "id": "2bb720b4-535b-428a-a406-1079548d61a8",
+   "cell_type": "markdown",
+   "id": "dabd2544",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "data = simulate(SP=4,N=1000,scale_re=0.01,G = 10)\n",
-    "E = data.X.shape[1]\n",
-    "N = data.X.shape[0]\n",
-    "SP = data.Y.shape[1]\n",
-    "\n"
+    "`fit_model` estimates a MVP jSDM based on the sjSDM approach and returns the standard deviation of the random effects and the covariance accuracy"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 716,
+   "execution_count": 12,
    "id": "b7bbb930-6320-4ab0-9f33-61f6bec25d2e",
    "metadata": {},
    "outputs": [],
@@ -146,18 +153,17 @@
     "    high = np.sqrt(6.0/(r_dim+df)) # type: ignore      \n",
     "    sigma = torch.tensor(np.random.uniform(low, high, [r_dim, df]), requires_grad = True, dtype=torch.float32, device=dev) # type: ignore\n",
     "\n",
-    "    #@torch.jit.script\n",
+    "    @torch.jit.script\n",
     "    def likelihood(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):\n",
     "        noise = torch.randn(size = [sampling, batch_size, df], device=torch.device(device), dtype=dtype)\n",
     "        E = torch.sigmoid(   torch.einsum(\"ijk, lk -> ijl\", [noise, sigma]).add(mu).mul(alpha)   ).mul(0.999999).add(0.0000005)\n",
     "        logprob = E.log().mul(Ys).add((1.0 - E).log().mul(1.0 - Ys)).neg().sum(dim = 2).neg()\n",
     "        maxlogprob = logprob.max(dim = 0).values\n",
     "        Eprob = logprob.sub(maxlogprob).exp().mean(dim = 0)\n",
-    "        loss = Eprob.log().neg().sub(maxlogprob)\n",
-    "        return loss\n",
+    "        return Eprob.log().neg().sub(maxlogprob)\n",
     "\n",
-    "    res = torch.tensor(np.random.normal(0, 0.01, size=(G, 1)), dtype=torch.float32, requires_grad=True, device=dev)\n",
     "    scale_log = torch.tensor(np.random.normal(0.0,0.001, [1]), dtype=torch.float32, requires_grad=True, device=dev)\n",
+    "    res = torch.tensor(np.random.normal(0.0,0.001, [G, 1]), dtype=torch.float32, requires_grad=True, device=dev)\n",
     "    zero_intercept = torch.zeros([1], dtype=torch.float32, device=dev)\n",
     "    zero_CAR = torch.zeros([G], dtype=torch.float32, device=dev)\n",
     "\n",
@@ -176,10 +182,9 @@
     "    \n",
     "    def re_loss():\n",
     "        if CAR is not True:\n",
-    "            loss = -torch.distributions.Normal(zero_intercept, soft(scale_log)).log_prob(res).sum()\n",
+    "            return -torch.distributions.Normal(zero_intercept, soft(scale_log)).log_prob(res).sum()\n",
     "        else:\n",
-    "            loss = -torch.distributions.MultivariateNormal(zero_CAR, (-soft(scale_log)*D).exp()+const_cov).log_prob(res.reshape([1,-1])).sum()\n",
-    "        return loss\n",
+    "            return -torch.distributions.MultivariateNormal(zero_CAR, (-soft(scale_log)*D).exp()+const_cov).log_prob(res.reshape([1,-1])).sum()\n",
     "\n",
     "    \n",
     "    def ll(res, W,sigma, XT, YT, indices_T):\n",
@@ -194,7 +199,7 @@
     "    \n",
     "    def torch_optim(epoch):\n",
     "        if epoch % inner_epochs == 0:\n",
-    "            for i in range(20):\n",
+    "            for _ in range(20):\n",
     "                opt2.zero_grad()\n",
     "                loss = ll(res, W,sigma, XT, YT, indices_T)\n",
     "                loss.backward()\n",
@@ -227,7 +232,6 @@
     "    \n",
     "    if det is not True:\n",
     "        const_val = torch.tensor(0.0, dtype=torch.float32, device=dev)\n",
-    "\n",
     "    for epoch in range(outer_epochs):\n",
     "        if optim == \"torch\": \n",
     "            torch_optim(epoch)\n",
@@ -250,81 +254,379 @@
     "    del XT\n",
     "    del YT\n",
     "    del sigma\n",
-    "    #torch.cuda.empty_cache()\n",
     "    return avg_scale.tolist(), avg_acc"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 723,
-   "id": "4ea51833-7de6-46a3-a2b2-caa8b64d247c",
+   "execution_count": 11,
+   "id": "79420b4e",
    "metadata": {},
    "outputs": [
     {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/var/folders/5c/9lp9691n1g15x6vj_bm1fjbm0000gn/T/ipykernel_72226/3649358109.py:89: ComplexWarning: Casting complex values to real discards the imaginary part\n",
-      "  re= np.random.multivariate_normal(np.zeros([sides**2]), cov=Sigma,size=[1])\n"
-     ]
+     "data": {
+      "text/plain": [
+       "array([[-0.40456842],\n",
+       "       [-0.47768778],\n",
+       "       [-0.57108772],\n",
+       "       [-0.21552874],\n",
+       "       [ 0.71881768],\n",
+       "       [ 0.31674885],\n",
+       "       [ 1.42935109],\n",
+       "       [-0.69193213],\n",
+       "       [ 0.18393983],\n",
+       "       [ 0.34073885],\n",
+       "       [ 0.20831082],\n",
+       "       [ 0.05193391],\n",
+       "       [-0.06783768],\n",
+       "       [-0.24047788],\n",
+       "       [ 0.3399139 ],\n",
+       "       [-0.52989648],\n",
+       "       [-0.02116909],\n",
+       "       [ 0.23199976],\n",
+       "       [ 0.02096044],\n",
+       "       [ 0.3455103 ],\n",
+       "       [ 0.06003734],\n",
+       "       [-0.93551233],\n",
+       "       [-0.39061279],\n",
+       "       [-0.34171041],\n",
+       "       [ 0.65986416],\n",
+       "       [ 0.7608267 ],\n",
+       "       [-0.36777635],\n",
+       "       [-1.2574748 ],\n",
+       "       [-0.25170756],\n",
+       "       [-0.3295159 ],\n",
+       "       [-0.31556222],\n",
+       "       [-0.77594736],\n",
+       "       [-0.14177741],\n",
+       "       [-0.52532314],\n",
+       "       [-0.09569401],\n",
+       "       [ 1.15591509],\n",
+       "       [-0.08016586],\n",
+       "       [ 0.82441512],\n",
+       "       [ 0.17998155],\n",
+       "       [ 1.08525601],\n",
+       "       [ 0.14132331],\n",
+       "       [ 0.65178794],\n",
+       "       [-0.41840465],\n",
+       "       [-0.45944114],\n",
+       "       [ 0.04853268],\n",
+       "       [ 0.90387567],\n",
+       "       [-0.30684973],\n",
+       "       [-0.0371868 ],\n",
+       "       [ 0.03796143],\n",
+       "       [-0.40112712],\n",
+       "       [ 0.47779864],\n",
+       "       [ 0.44921315],\n",
+       "       [ 0.40080424],\n",
+       "       [-0.07006795],\n",
+       "       [ 0.5423653 ],\n",
+       "       [ 0.36965576],\n",
+       "       [-0.14301217],\n",
+       "       [-0.22619527],\n",
+       "       [ 0.06485423],\n",
+       "       [ 0.29085903],\n",
+       "       [ 0.28651125],\n",
+       "       [ 0.24376217],\n",
+       "       [ 0.1374285 ],\n",
+       "       [ 0.71173097],\n",
+       "       [-0.27467802],\n",
+       "       [ 0.35474841],\n",
+       "       [ 0.19161101],\n",
+       "       [ 1.05510804],\n",
+       "       [-0.36710234],\n",
+       "       [ 0.32709943],\n",
+       "       [ 0.06457097],\n",
+       "       [-0.02826677],\n",
+       "       [-0.18084159],\n",
+       "       [-0.21584778],\n",
+       "       [-0.12073306],\n",
+       "       [ 0.71320465],\n",
+       "       [ 0.32572252],\n",
+       "       [ 0.758707  ],\n",
+       "       [-0.14699749],\n",
+       "       [ 0.54849441],\n",
+       "       [ 0.6140983 ],\n",
+       "       [ 0.26618183],\n",
+       "       [-0.58234565],\n",
+       "       [-0.60209572],\n",
+       "       [-0.34503113],\n",
+       "       [ 0.37050768],\n",
+       "       [-0.21642521],\n",
+       "       [-0.55042142],\n",
+       "       [-0.00415465],\n",
+       "       [ 0.61650127],\n",
+       "       [ 1.0966027 ],\n",
+       "       [-0.20510165],\n",
+       "       [ 0.20709255],\n",
+       "       [ 0.02846783],\n",
+       "       [ 0.10907124],\n",
+       "       [ 0.21588155],\n",
+       "       [-0.1023169 ],\n",
+       "       [ 0.78601463],\n",
+       "       [-0.02174354],\n",
+       "       [ 0.14260426]])"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
     }
    ],
    "source": [
-    "data = simulate(SP=2,N=400,scale_re=0.3,G = 20, CAR=True)"
+    "data = simulate(SP=2,N=5000,scale_re=0.5,G = 100, CAR=False)\n",
+    "#fit_model(data, CAR=True)\n",
+    "data.re"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 725,
-   "id": "35ec11ba-4e1f-4d31-8910-fbd995505155",
-   "metadata": {
-    "tags": []
-   },
+   "execution_count": 12,
+   "id": "088c1413",
+   "metadata": {},
    "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[55.80103302001953, 21.666128158569336]\n",
+      "[59.56851577758789, 26.515663146972656]\n",
+      "[53.23011016845703, 21.068098068237305]\n",
+      "[48.70780944824219, 18.318222045898438]\n",
+      "[49.03778076171875, 20.532358169555664]\n",
+      "[47.15538787841797, 19.405576705932617]\n",
+      "[47.09665298461914, 19.003772735595703]\n",
+      "[48.80914306640625, 20.8721923828125]\n",
+      "[51.730995178222656, 24.477840423583984]\n",
+      "[38.06648635864258, 12.463699340820312]\n"
+     ]
+    },
     {
      "data": {
       "text/plain": [
-       "(0.6581300497055054, 0.6666666666666666)"
+       "[0.20226380921304, 0.458190164843836, 0.4199224385975236, 1.0]"
       ]
      },
-     "execution_count": 725,
+     "execution_count": 12,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "fit_model(data, 550, 20, optim=\"scipy\", likelihood_type=\"mvp\", CAR=True)"
+    "fit_model(data, CAR=False, det=True)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 580,
-   "id": "bd167e91-bea3-4ba9-bf79-8ba3004442e5",
+   "execution_count": 75,
+   "id": "e7f56df6",
    "metadata": {},
    "outputs": [],
    "source": [
-    "results_y = [fit_model(data, 200, 12, optim=\"scipy\", likelihood_type=\"mvp\", det=True)  for _ in range(15)]\n",
-    "#results_n = [fit_model(data, 60, 5, optim=\"scipy\", likelihood_type=\"mvp\", det=False) for _ in range(60)]"
+    "def fit_model(data: Simulation, det=True, CAR=False, device = \"cpu:0\", batch_size=50, ll=\"MVP\", intercept=\"sites\", df=None) -> list:\n",
+    "     N, E = data.X.shape\n",
+    "     SP = data.Y.shape[1]\n",
+    "     X, Y, G, indices, D = data.X, data.Y, data.G, data.g, data.D\n",
+    "     dev = torch.device(device)\n",
+    "     r_dim = Y.shape[1]\n",
+    "     if df is None:\n",
+    "          df = int(np.rint(Y.shape[1]/2))\n",
+    "     low = -np.sqrt(6.0/(r_dim+df)) \n",
+    "     high = np.sqrt(6.0/(r_dim+df))     \n",
+    "     XT = torch.tensor(X, dtype=torch.float32, device=torch.device(\"cpu:0\"))\n",
+    "     YT = torch.tensor(Y, dtype=torch.float32, device=torch.device(\"cpu:0\"))\n",
+    "     indices_T = torch.tensor(indices, dtype=torch.long, device=torch.device(\"cpu:0\"))\n",
+    "     init_scale = 10.0\n",
+    "     if CAR is True:\n",
+    "          G = G**2\n",
+    "          D = torch.tensor(D, dtype=torch.float32, device=dev)\n",
+    "          init_scale = 0.0\n",
+    "     # Variables\n",
+    "     W = torch.tensor(np.random.normal(0.0,0.001, [XT.shape[1], YT.shape[1]]), dtype=torch.float32, device=dev, requires_grad=True)\n",
+    "     \n",
+    "     if intercept is not \"species\":\n",
+    "          scale_log = torch.tensor(1.0, dtype=torch.float32, requires_grad=True, device=dev)\n",
+    "          res = torch.tensor(np.random.normal(0.0,0.001, [G, 1]), dtype=torch.float32, requires_grad=True, device=dev)\n",
+    "     if intercept is \"species\":\n",
+    "          res = torch.tensor(np.random.normal(0.0,0.001, [1, SP]), dtype=torch.float32, requires_grad=True, device=dev)\n",
+    "          scale_log = torch.tensor(np.random.normal(0.0,0.001, [SP]), dtype=torch.float32, requires_grad=True, device=dev)\n",
+    "          init_scale = torch.ones_like(scale_log)\n",
+    "     sigma = torch.tensor(np.random.uniform(low, high, [r_dim, df]), requires_grad = True, dtype=torch.float32, device=dev) # type: ignore\n",
+    "\n",
+    "     soft = lambda t: torch.nn.functional.softplus(t)+0.0001\n",
+    "     zero_intercept = torch.zeros([1], dtype=torch.float32, device=dev)\n",
+    "     zero_CAR = torch.zeros([G], dtype=torch.float32, device=dev)\n",
+    "     loss2 = torch.zeros([1], dtype=torch.float32, device=dev)\n",
+    "     adapt = torch.tensor(np.rint(XT.shape[0]/batch_size).tolist(), dtype=torch.float32, device=dev)\n",
+    "\n",
+    "\n",
+    "     @torch.jit.script\n",
+    "     def MVP(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):\n",
+    "         noise = torch.randn(size = [sampling, batch_size, df], device=torch.device(device), dtype=dtype)\n",
+    "         E = torch.sigmoid(   torch.einsum(\"ijk, lk -> ijl\", [noise, sigma]).add(mu).mul(alpha)   ).mul(0.999999).add(0.0000005)\n",
+    "         logprob = E.log().mul(Ys).add((1.0 - E).log().mul(1.0 - Ys)).neg().sum(dim = 2).neg()\n",
+    "         maxlogprob = logprob.max(dim = 0).values\n",
+    "         Eprob = logprob.sub(maxlogprob).exp().mean(dim = 0)\n",
+    "         return Eprob.log().neg().sub(maxlogprob)\n",
+    "    \n",
+    "     def Binomial(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):\n",
+    "          return - torch.distributions.Binomial(1, mu.sigmoid()).log_prob(Ys)\n",
+    "     \n",
+    "     if ll is \"MVP\":\n",
+    "          likelihood = MVP\n",
+    "     else:\n",
+    "          likelihood = Binomial\n",
+    "\n",
+    "     def ll(res, W, sigma, XT, YT, indices_T, scale_log):\n",
+    "          if intercept is not \"species\":\n",
+    "               pred = XT@W+res[indices_T,:]\n",
+    "          else:\n",
+    "               pred = XT@W+res\n",
+    "          loss = likelihood(pred, YT, sigma, XT.shape[0], 100, df, 1.7012, device, torch.float32).sum()/XT.shape[0]\n",
+    "          if intercept is \"sites\":\n",
+    "               loss += -torch.distributions.Normal(zero_intercept, (scale_log.exp())).log_prob(res[indices_T.unique()]).sum()/adapt/XT.shape[0]\n",
+    "          if intercept is \"species\":\n",
+    "               loss += -torch.distributions.LowRankMultivariateNormal(torch.zeros([YT.shape[1]], dtype=torch.float32, device=dev), sigma, scale_log.exp() ).log_prob(res).sum()/adapt/XT.shape[0]\n",
+    "          if CAR is True:\n",
+    "               ind2 = indices_T.unique()\n",
+    "               D_tmp = D.index_select(0, ind2).index_select(1, ind2)\n",
+    "               const_val = torch.eye(ind2.shape[0], device=dev, dtype=torch.float32)*0.001\n",
+    "               loss += -torch.distributions.MultivariateNormal(zero_CAR[ind2], (-(scale_log.exp())*D_tmp).exp()+const_val).log_prob(res[indices_T.unique()].reshape([1,-1])).sum()/adapt/XT.shape[0]\n",
+    "          return loss\n",
+    "\n",
+    "     \n",
+    "     optimizer = torch.optim.Adamax([W, scale_log, sigma], lr = 0.03)\n",
+    "     optimizer_re = torch.optim.Adamax([res], lr = 0.03)\n",
+    "\n",
+    "     dataset = torch.utils.data.TensorDataset(XT, YT, indices_T)\n",
+    "     dataLoader = torch.utils.data.DataLoader(dataset, batch_size=batch_size, shuffle=True)\n",
+    "     for _ in range(50):\n",
+    "          for x, y, inds in dataLoader:\n",
+    "               optimizer_re.zero_grad()\n",
+    "               loss = ll(res, W.detach(),sigma.detach(), x.to(dev), y.to(dev), inds.to(dev),torch.tensor(init_scale, dtype=torch.float32, device=dev) )\n",
+    "               loss.backward(  )\n",
+    "               optimizer_re.step()\n",
+    "     optimizer_re.zero_grad()\n",
+    "     print(torch.std(res))\n",
+    "\n",
+    "     for i in range(100):\n",
+    "          if i > 0:\n",
+    "               #print(torch.std(res))\n",
+    "               for x, y, inds in dataLoader:\n",
+    "                    optimizer_re.zero_grad()\n",
+    "                    loss = ll(res, W.detach(),sigma, x.to(dev), y.to(dev), inds.to(dev), scale_log.detach())#/x.shape[0]\n",
+    "                    loss.backward()\n",
+    "                    optimizer_re.step()\n",
+    "          optimizer_re.zero_grad()\n",
+    "\n",
+    "          for x, y, inds in dataLoader:\n",
+    "               optimizer.zero_grad()\n",
+    "               loss = ll(res, W, sigma, x.to(dev), y.to(dev), inds.to(dev), scale_log)#/x.shape[0]\n",
+    "               if det is True:\n",
+    "                    loss.backward(  retain_graph=True )\n",
+    "                    gg = torch.autograd.grad(loss, res, retain_graph=True, create_graph=True)[0]\n",
+    "                    gg=gg[gg.nonzero(as_tuple=True)].reshape([-1,1])\n",
+    "                    logDA = torch.reciprocal(gg**2*gg.shape[0]).sqrt().reshape([-1]).diag().inverse().logdet()\n",
+    "                    #hess = torch.autograd.functional.hessian(lambda res: ll(res, W, sigma, x.to(dev), y.to(dev), inds.to(dev),scale_log), res, create_graph=True).squeeze()\n",
+    "                    #ind2 = inds.to(dev).unique()\n",
+    "                    #D_tmp = hess.index_select(0, ind2).index_select(1, ind2)\n",
+    "                    #const_val = torch.eye(ind2.shape[0], device=dev, dtype=torch.float32)*0.01\n",
+    "                    #logDA=(D_tmp+const_val).inverse().logdet()\n",
+    "                    loss2 = ((gg.shape[0]*0.5*torch.log((2*torch.tensor(3.14, dtype=torch.float32, device=dev))) - 0.5*logDA))/adapt/x.shape[0]\n",
+    "                    #loss2 = gg.shape[0]*0.5*torch.log((2*torch.tensor(3.14, dtype=torch.float32))) - 0.5*logDA\n",
+    "                    loss+=loss2\n",
+    "               loss = loss\n",
+    "               loss.backward()\n",
+    "               optimizer.step()\n",
+    "\n",
+    "          optimizer.zero_grad()\n",
+    "\n",
+    "          #if i % 10 == 0:\n",
+    "          #     print([loss.item(), loss2.item()])\n",
+    "               \n",
+    "     return [(scale_log.exp().mean()).cpu().data.numpy().tolist(), \n",
+    "             np.mean((np.sign(correlation_from_covariance((sigma @ sigma.t()).data.cpu().numpy())) == np.sign(data.Sigma) )[np.triu_indices(SP)]), \n",
+    "             correlation_from_covariance((sigma @ sigma.t() + scale_log.exp()).data.cpu().numpy()) ]"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 581,
-   "id": "19adb7ac-a9d5-4f97-8621-e47ee8834fc5",
+   "execution_count": 69,
+   "id": "83e59b4d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data = simulate(SP = 100, N = 1000  )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 79,
+   "id": "14e7144d",
    "metadata": {},
    "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/maximilianpichler/Library/r-miniconda/envs/r-sjsdm/lib/python3.7/site-packages/ipykernel_launcher.py:81: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n"
+     ]
+    },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[4.13961991 0.51632461]\n",
-      "[0.3287191 0.0027272]\n"
+      "tensor(0.0823, grad_fn=<StdBackward0>)\n"
      ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "[8.509558028890751e-06,\n",
+       " 0.5134653465346535,\n",
+       " array([[1.        , 0.27305794, 0.3846035 , ..., 0.5822861 , 0.5183367 ,\n",
+       "         0.31661004],\n",
+       "        [0.2729554 , 1.        , 0.69973063, ..., 0.5865577 , 0.51317865,\n",
+       "         0.41813707],\n",
+       "        [0.38451353, 0.6997175 , 0.9999999 , ..., 0.7276242 , 0.624743  ,\n",
+       "         0.5473529 ],\n",
+       "        ...,\n",
+       "        [0.58230746, 0.58664626, 0.7276993 , ..., 1.0000001 , 0.58745205,\n",
+       "         0.64994854],\n",
+       "        [0.5183621 , 0.5132928 , 0.62483996, ..., 0.58745027, 1.        ,\n",
+       "         0.6510153 ],\n",
+       "        [0.316561  , 0.41817302, 0.5473907 , ..., 0.64989865, 0.65095186,\n",
+       "         1.0000001 ]], dtype=float32)]"
+      ]
+     },
+     "execution_count": 79,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fit_model(data, ll=\"binomial\", intercept=\"species\", df = 40, batch_size = 200, det = False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "id": "ee74000a",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "tensor([[ 0.1462, -0.2108,  1.5039, -0.3509,  0.9086]])"
+      ]
+     },
+     "execution_count": 39,
+     "metadata": {},
+     "output_type": "execute_result"
     }
    ],
    "source": [
-    "print( np.mean(np.asarray( results_y ), axis = 0) )\n",
-    "print( np.std(np.asarray( results_y ), axis = 0) )\n"
+    "torch.distributions.LowRankMultivariateNormal(torch.zeros(5), torch.ones([5,2]), torch.ones([5])).sample([1])"
    ]
   }
  ],
@@ -344,7 +646,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.12"
+   "version": "3.7.1"
   }
  },
  "nbformat": 4,
diff --git a/Code/sjSDM with random effects/LA.ipynb b/Code/sjSDM with random effects/LA.ipynb
new file mode 100644
index 0000000..f186e84
--- /dev/null
+++ b/Code/sjSDM with random effects/LA.ipynb
@@ -0,0 +1,378 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import torch\n",
+    "import numpy as np\n",
+    "import numpy as np, numpy.linalg\n",
+    "import matplotlib.pyplot as plt\n",
+    "from dataclasses import dataclass\n",
+    "import statsmodels.api as sm\n",
+    "import statsmodels.formula.api as smf\n",
+    "import pandas as pd\n",
+    "from scipy.optimize import minimize\n",
+    "from dataclasses import dataclass\n",
+    "import numpy as np\n",
+    "\n",
+    "@dataclass\n",
+    "class Simulation:\n",
+    "    X: np.ndarray\n",
+    "    Y: np.ndarray\n",
+    "    W: np.ndarray\n",
+    "    g: np.ndarray\n",
+    "    G: np.ndarray\n",
+    "    re: np.ndarray\n",
+    "    sd_re: float\n",
+    "    \n",
+    "def simulate(N=100,E=2,G=10, sd_re=1.0):\n",
+    "    X = np.random.uniform(-1, 1, [N, E])\n",
+    "    W = np.random.normal(0, 1.0, [E, 1])\n",
+    "    Y = np.matmul(X, W) + np.random.normal(0.0, 0.8, [N,1])\n",
+    "    re = np.random.normal(0, sd_re, [G,1])\n",
+    "    g = np.repeat(np.arange(0,G), np.round(N/G))\n",
+    "    Y = Y+re[g,:]\n",
+    "    return Simulation(X, Y,W, g, G, re, sd_re)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 429,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#@ray.remote(num_gpus=1,num_cpus=2)\n",
+    "def fit_model(data: Simulation, det=True, CAR=False, device = \"cuda:0\", batch_size = 100) -> list:\n",
+    "     N, E = data.X.shape\n",
+    "     SP = data.Y.shape[1]\n",
+    "     X, Y, G, indices = data.X, data.Y, data.G, data.g\n",
+    "     dev = torch.device(device)\n",
+    "    \n",
+    "     XT = torch.tensor(X, dtype=torch.float32, device=torch.device(\"cpu:0\"))\n",
+    "     YT = torch.tensor(Y, dtype=torch.float32, device=torch.device(\"cpu:0\"))\n",
+    "     indices_T = torch.tensor(indices, dtype=torch.long, device=torch.device(\"cpu:0\"))\n",
+    "\n",
+    "     # Variables\n",
+    "     W = torch.tensor(np.random.normal(0.0,0.001, [XT.shape[1], YT.shape[1]]), dtype=torch.float32, device=dev, requires_grad=True)\n",
+    "     scale_log = torch.tensor(1.0, dtype=torch.float32, requires_grad=True, device=dev)\n",
+    "     scale_log_normal = torch.tensor(0.0, dtype=torch.float32, requires_grad=True, device=dev)\n",
+    "     res = torch.tensor(np.random.normal(0.0,0.001, [G, 1]), dtype=torch.float32, requires_grad=True, device=dev)\n",
+    "     soft = lambda t: torch.nn.functional.softplus(t)+0.0001\n",
+    "     zero_intercept = torch.zeros([1], dtype=torch.float32, device=dev)\n",
+    "     loss2 = torch.zeros([1], dtype=torch.float32, device=dev)\n",
+    "\n",
+    "     mse = torch.nn.MSELoss(reduction=\"mean\")\n",
+    "     \n",
+    "     adapt = torch.tensor(np.rint(XT.shape[0]/batch_size).tolist(), dtype=torch.float32, device=dev)\n",
+    "     \n",
+    "     def ll(res, W, XT, YT, indices_T, scale_log, scale_log_normal):\n",
+    "        pred = XT@W+res[indices_T,:]\n",
+    "        loss = -torch.distributions.Normal(pred, torch.exp(scale_log_normal)).log_prob(YT).sum()/XT.shape[0]\n",
+    "        #loss = mse(YT, pred).mean()\n",
+    "        #sigma_noise = torch.exp(scale_log_normal)\n",
+    "        #sigma_res = 1./torch.square(torch.exp(scale_log))\n",
+    "        # crit_factor = 1.0 / (2 * sigma_noise.square())\n",
+    "        loss += -torch.distributions.Normal(zero_intercept, torch.exp(scale_log)).log_prob(res[indices_T.unique()]).sum()/adapt/XT.shape[0]#/XT.shape[0]\n",
+    "        #loss += ((res[indices_T.unique()].pow(2.0))*(sigma_res)*0.5).sum()/res[indices_T.unique()].shape[0]#/crit_factor\n",
+    "        return loss\n",
+    "\n",
+    "\n",
+    "     optimizer = torch.optim.Adamax([W, scale_log,scale_log_normal], lr = 0.1)\n",
+    "     optimizer_re = torch.optim.Adamax([res], lr = 0.1)\n",
+    "\n",
+    "     dataset = torch.utils.data.TensorDataset(XT, YT, indices_T)\n",
+    "     dataLoader = torch.utils.data.DataLoader(dataset, batch_size=batch_size, shuffle=True)\n",
+    "     \n",
+    "     for _ in range(20):\n",
+    "          for x, y, inds in dataLoader:\n",
+    "               optimizer_re.zero_grad()\n",
+    "               loss = ll(res, W.detach(), x.to(dev), y.to(dev), inds.to(dev),torch.tensor(100.0, dtype=torch.float32), scale_log_normal.detach() )\n",
+    "               loss.backward(  )\n",
+    "               optimizer_re.step()\n",
+    "     optimizer_re.zero_grad()\n",
+    "     print(torch.std(res))\n",
+    "\n",
+    "     for i in range(40):\n",
+    "          if i > 1:\n",
+    "               #print(torch.std(res))\n",
+    "               for x, y, inds in dataLoader:\n",
+    "                    optimizer_re.zero_grad()\n",
+    "                    loss = ll(res, W, x.to(dev), y.to(dev), inds.to(dev), scale_log.detach(), scale_log_normal.detach())#/x.shape[0]\n",
+    "                    loss.backward()\n",
+    "                    optimizer_re.step()\n",
+    "               optimizer_re.zero_grad()\n",
+    "          \n",
+    "          if i > 0:\n",
+    "               for _ in range(3):\n",
+    "                    for x, y, inds in dataLoader:\n",
+    "                         optimizer.zero_grad()\n",
+    "                         loss = ll(res, W, x.to(dev), y.to(dev), inds.to(dev),scale_log, scale_log_normal)#/x.shape[0]\n",
+    "                         if det is True:\n",
+    "                              #loss.backward(  retain_graph=True )\n",
+    "                              #gg = torch.autograd.grad(loss, res, retain_graph=True, create_graph=True)[0]\n",
+    "                              #gg=gg[gg.nonzero(as_tuple=True)].reshape([-1,1])\n",
+    "                              #logDA = torch.reciprocal(gg**2).sqrt().reshape([-1]).diag().inverse().logdet()\n",
+    "                              #loss2 = (-0.5*logDA)#/x.shape[0]\n",
+    "                              hess = torch.autograd.functional.hessian(lambda res: ll(res, W, x.to(dev), y.to(dev), inds.to(dev),scale_log, scale_log_normal), res, create_graph=True).squeeze()\n",
+    "                              ind2 = inds.unique()\n",
+    "                              D_tmp = hess.index_select(0, ind2).index_select(1, ind2)\n",
+    "                              const_val = torch.eye(ind2.shape[0], device=dev, dtype=torch.float32)*0.01\n",
+    "                              logDA=(D_tmp+const_val).inverse().logdet()\n",
+    "                              loss2 = ((ind2.shape[0]*0.5*torch.log((2*torch.tensor(3.14, dtype=torch.float32))) - 0.5*logDA))/adapt/x.shape[0]\n",
+    "                              loss+=loss2\n",
+    "                         loss = loss\n",
+    "                         loss.backward()\n",
+    "                         optimizer.step()\n",
+    "          \n",
+    "          optimizer.zero_grad()\n",
+    "          \n",
+    "\n",
+    "\n",
+    "          if i % 20 == 0:\n",
+    "               print([loss.item(), loss2])\n",
+    "               \n",
+    "     return [(torch.exp(scale_log)).cpu().data.numpy().tolist(), \n",
+    "             (torch.exp(scale_log_normal)).cpu().data.numpy().tolist(), \n",
+    "             W.cpu().data.numpy()]\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 423,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "5.289999999999999\n",
+      "         Mixed Linear Model Regression Results\n",
+      "========================================================\n",
+      "Model:            MixedLM Dependent Variable: Y         \n",
+      "No. Observations: 1000    Method:             ML        \n",
+      "No. Groups:       50      Scale:              0.6364    \n",
+      "Min. group size:  20      Log-Likelihood:     -1323.5319\n",
+      "Max. group size:  20      Converged:          Yes       \n",
+      "Mean group size:  20.0                                  \n",
+      "--------------------------------------------------------\n",
+      "             Coef.  Std.Err.    z    P>|z| [0.025 0.975]\n",
+      "--------------------------------------------------------\n",
+      "X1           -1.395    0.045 -30.859 0.000 -1.484 -1.306\n",
+      "X2            0.223    0.046   4.891 0.000  0.134  0.313\n",
+      "Group Var     5.869    1.517                            \n",
+      "========================================================\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "data = simulate(N=1000, G = 50, sd_re=2.3)\n",
+    "print(2.3**2)\n",
+    "dataset = pd.DataFrame(np.concatenate([data.X, data.Y, np.reshape(data.g, [data.g.shape[0],1])], axis = 1),\n",
+    "                       columns = [\"X1\", \"X2\", \"Y\", \"g\"]\n",
+    "                      )\n",
+    "md = smf.mixedlm(\"Y~0+X1+X2\", dataset, groups=dataset[\"g\"], )\n",
+    "mdf = md.fit(reml=False)\n",
+    "print(mdf.summary())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 428,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "tensor(1.2997, grad_fn=<StdBackward0>)\n",
+      "[inf, tensor([0.])]\n",
+      "tensor(175.1028, grad_fn=<LogdetBackward0>)\n",
+      "tensor(181.4552, grad_fn=<LogdetBackward0>)\n",
+      "tensor(186.3749, grad_fn=<LogdetBackward0>)\n",
+      "tensor(190.2208, grad_fn=<LogdetBackward0>)\n",
+      "tensor(193.2012, grad_fn=<LogdetBackward0>)\n",
+      "tensor(195.5276, grad_fn=<LogdetBackward0>)\n",
+      "tensor(197.3473, grad_fn=<LogdetBackward0>)\n",
+      "tensor(198.7539, grad_fn=<LogdetBackward0>)\n",
+      "tensor(199.8382, grad_fn=<LogdetBackward0>)\n",
+      "tensor(200.6704, grad_fn=<LogdetBackward0>)\n",
+      "tensor(201.2995, grad_fn=<LogdetBackward0>)\n",
+      "tensor(201.7719, grad_fn=<LogdetBackward0>)\n",
+      "tensor(202.1229, grad_fn=<LogdetBackward0>)\n",
+      "tensor(202.3760, grad_fn=<LogdetBackward0>)\n",
+      "tensor(202.5502, grad_fn=<LogdetBackward0>)\n",
+      "tensor(202.6577, grad_fn=<LogdetBackward0>)\n",
+      "tensor(202.7038, grad_fn=<LogdetBackward0>)\n",
+      "tensor(202.6926, grad_fn=<LogdetBackward0>)\n",
+      "tensor(202.6261, grad_fn=<LogdetBackward0>)\n",
+      "tensor(202.5042, grad_fn=<LogdetBackward0>)\n",
+      "tensor(202.3284, grad_fn=<LogdetBackward0>)\n",
+      "tensor(202.1013, grad_fn=<LogdetBackward0>)\n",
+      "tensor(201.8234, grad_fn=<LogdetBackward0>)\n",
+      "tensor(201.4984, grad_fn=<LogdetBackward0>)\n",
+      "tensor(201.1310, grad_fn=<LogdetBackward0>)\n",
+      "tensor(200.7227, grad_fn=<LogdetBackward0>)\n",
+      "tensor(200.2790, grad_fn=<LogdetBackward0>)\n",
+      "tensor(199.8050, grad_fn=<LogdetBackward0>)\n",
+      "tensor(199.3027, grad_fn=<LogdetBackward0>)\n",
+      "tensor(198.7776, grad_fn=<LogdetBackward0>)\n",
+      "tensor(198.2352, grad_fn=<LogdetBackward0>)\n",
+      "tensor(197.6769, grad_fn=<LogdetBackward0>)\n",
+      "tensor(197.1081, grad_fn=<LogdetBackward0>)\n",
+      "tensor(196.5340, grad_fn=<LogdetBackward0>)\n",
+      "tensor(195.9553, grad_fn=<LogdetBackward0>)\n",
+      "tensor(195.3772, grad_fn=<LogdetBackward0>)\n",
+      "tensor(194.8043, grad_fn=<LogdetBackward0>)\n",
+      "tensor(194.2368, grad_fn=<LogdetBackward0>)\n",
+      "tensor(193.6790, grad_fn=<LogdetBackward0>)\n",
+      "tensor(193.1351, grad_fn=<LogdetBackward0>)\n",
+      "tensor(192.6045, grad_fn=<LogdetBackward0>)\n",
+      "tensor(192.0908, grad_fn=<LogdetBackward0>)\n",
+      "tensor(191.5974, grad_fn=<LogdetBackward0>)\n",
+      "tensor(191.1226, grad_fn=<LogdetBackward0>)\n",
+      "tensor(190.6696, grad_fn=<LogdetBackward0>)\n",
+      "tensor(190.2406, grad_fn=<LogdetBackward0>)\n",
+      "tensor(189.8331, grad_fn=<LogdetBackward0>)\n",
+      "tensor(189.4492, grad_fn=<LogdetBackward0>)\n",
+      "tensor(189.0907, grad_fn=<LogdetBackward0>)\n",
+      "tensor(188.7537, grad_fn=<LogdetBackward0>)\n",
+      "tensor(188.4399, grad_fn=<LogdetBackward0>)\n",
+      "tensor(188.1501, grad_fn=<LogdetBackward0>)\n",
+      "tensor(187.8799, grad_fn=<LogdetBackward0>)\n",
+      "tensor(187.6303, grad_fn=<LogdetBackward0>)\n",
+      "tensor(187.4016, grad_fn=<LogdetBackward0>)\n",
+      "tensor(187.1890, grad_fn=<LogdetBackward0>)\n",
+      "tensor(186.9930, grad_fn=<LogdetBackward0>)\n",
+      "tensor(186.8141, grad_fn=<LogdetBackward0>)\n",
+      "tensor(186.6468, grad_fn=<LogdetBackward0>)\n",
+      "tensor(186.4917, grad_fn=<LogdetBackward0>)\n",
+      "[1.641299843788147, tensor(-0.0473, grad_fn=<DivBackward0>)]\n",
+      "tensor(186.3492, grad_fn=<LogdetBackward0>)\n",
+      "tensor(186.2140, grad_fn=<LogdetBackward0>)\n",
+      "tensor(186.0867, grad_fn=<LogdetBackward0>)\n",
+      "tensor(185.9677, grad_fn=<LogdetBackward0>)\n",
+      "tensor(185.8521, grad_fn=<LogdetBackward0>)\n",
+      "tensor(185.7405, grad_fn=<LogdetBackward0>)\n",
+      "tensor(185.6338, grad_fn=<LogdetBackward0>)\n",
+      "tensor(185.5270, grad_fn=<LogdetBackward0>)\n",
+      "tensor(185.4211, grad_fn=<LogdetBackward0>)\n",
+      "tensor(185.3172, grad_fn=<LogdetBackward0>)\n",
+      "tensor(185.2108, grad_fn=<LogdetBackward0>)\n",
+      "tensor(185.1030, grad_fn=<LogdetBackward0>)\n",
+      "tensor(184.9951, grad_fn=<LogdetBackward0>)\n",
+      "tensor(184.8828, grad_fn=<LogdetBackward0>)\n",
+      "tensor(184.7678, grad_fn=<LogdetBackward0>)\n",
+      "tensor(184.6514, grad_fn=<LogdetBackward0>)\n",
+      "tensor(184.5295, grad_fn=<LogdetBackward0>)\n",
+      "tensor(184.4041, grad_fn=<LogdetBackward0>)\n",
+      "tensor(184.2766, grad_fn=<LogdetBackward0>)\n",
+      "tensor(184.1432, grad_fn=<LogdetBackward0>)\n",
+      "tensor(184.0059, grad_fn=<LogdetBackward0>)\n",
+      "tensor(183.8665, grad_fn=<LogdetBackward0>)\n",
+      "tensor(183.7212, grad_fn=<LogdetBackward0>)\n",
+      "tensor(183.5721, grad_fn=<LogdetBackward0>)\n",
+      "tensor(183.4210, grad_fn=<LogdetBackward0>)\n",
+      "tensor(183.2644, grad_fn=<LogdetBackward0>)\n",
+      "tensor(183.1045, grad_fn=<LogdetBackward0>)\n",
+      "tensor(182.9430, grad_fn=<LogdetBackward0>)\n",
+      "tensor(182.7766, grad_fn=<LogdetBackward0>)\n",
+      "tensor(182.6073, grad_fn=<LogdetBackward0>)\n",
+      "tensor(182.4371, grad_fn=<LogdetBackward0>)\n",
+      "tensor(182.2626, grad_fn=<LogdetBackward0>)\n",
+      "tensor(182.0859, grad_fn=<LogdetBackward0>)\n",
+      "tensor(181.9089, grad_fn=<LogdetBackward0>)\n",
+      "tensor(181.7283, grad_fn=<LogdetBackward0>)\n",
+      "tensor(181.5461, grad_fn=<LogdetBackward0>)\n",
+      "tensor(181.3644, grad_fn=<LogdetBackward0>)\n",
+      "tensor(181.1796, grad_fn=<LogdetBackward0>)\n",
+      "tensor(180.9939, grad_fn=<LogdetBackward0>)\n",
+      "tensor(180.8093, grad_fn=<LogdetBackward0>)\n",
+      "tensor(180.6221, grad_fn=<LogdetBackward0>)\n",
+      "tensor(180.4346, grad_fn=<LogdetBackward0>)\n",
+      "tensor(180.2487, grad_fn=<LogdetBackward0>)\n",
+      "tensor(180.0608, grad_fn=<LogdetBackward0>)\n",
+      "tensor(179.8730, grad_fn=<LogdetBackward0>)\n",
+      "tensor(179.6872, grad_fn=<LogdetBackward0>)\n",
+      "tensor(179.4998, grad_fn=<LogdetBackward0>)\n",
+      "tensor(179.3128, grad_fn=<LogdetBackward0>)\n",
+      "tensor(179.1282, grad_fn=<LogdetBackward0>)\n",
+      "tensor(178.9422, grad_fn=<LogdetBackward0>)\n",
+      "tensor(178.7569, grad_fn=<LogdetBackward0>)\n",
+      "tensor(178.5742, grad_fn=<LogdetBackward0>)\n",
+      "tensor(178.3903, grad_fn=<LogdetBackward0>)\n",
+      "tensor(178.2073, grad_fn=<LogdetBackward0>)\n",
+      "tensor(178.0269, grad_fn=<LogdetBackward0>)\n",
+      "tensor(177.8454, grad_fn=<LogdetBackward0>)\n",
+      "tensor(177.6649, grad_fn=<LogdetBackward0>)\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "[1.9286086559295654,\n",
+       " 1.0407851934432983,\n",
+       " array([[-1.3586684 ],\n",
+       "        [ 0.22159573]], dtype=float32)]"
+      ]
+     },
+     "execution_count": 428,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fit_model(data, device=\"cpu:0\", det=True, batch_size=1000)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 379,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[ 0.17998728],\n",
+       "       [-1.49548422]])"
+      ]
+     },
+     "execution_count": 379,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "data.W\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "interpreter": {
+   "hash": "b2c8e16c7ccab309b6c7ff727be01f1772ce77360f8d9c2b6c3a69dbab5f4903"
+  },
+  "kernelspec": {
+   "display_name": "Python 3.7.1 64-bit (conda)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.1"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/README.Rmd b/README.Rmd
index 8b80d1d..e7b0bbf 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -106,8 +106,13 @@ plot(an)
 
 The anova shows the relative changes in the deviance of the groups and their intersections.
 
+We can also visualize the individual contributions to the species and site $$R^2$$:
 
 
+```{r,fig.height=7, fig.width=8, warning=FALSE}
+plot(an, internal=TRUE)
+```
+
 
 If it fails, check out the help of ?install_sjSDM, ?installation_help, and vignette("Dependencies", package = "sjSDM"). 
 
diff --git a/README.md b/README.md
index d76151b..7dec929 100644
--- a/README.md
+++ b/README.md
@@ -71,7 +71,7 @@ Simulate a community and fit model:
 library(sjSDM)
 ```
 
-    ## ── Attaching sjSDM ──────────────────────────────────────────────────── 1.0.0 ──
+    ## ── Attaching sjSDM ──────────────────────────────────────────────────── 1.0.1 ──
 
     ## ✓ torch <environment> 
     ## ✓ torch_optimizer  
@@ -89,73 +89,73 @@ model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1+X2+X3), spatial =
 summary(model)
 ```
 
-    ## LogLik:  -523.0463 
+    ## LogLik:  -517.1162 
     ## Regularization loss:  0 
     ## 
     ## Species-species correlation matrix: 
     ## 
     ##  sp1  1.0000                                 
-    ##  sp2 -0.3370  1.0000                             
-    ##  sp3 -0.1490 -0.3520  1.0000                         
-    ##  sp4 -0.1690 -0.3120  0.6950  1.0000                     
-    ##  sp5  0.6020 -0.3360 -0.0950 -0.1020  1.0000                 
-    ##  sp6 -0.2250  0.3830  0.1630  0.1670 -0.0940  1.0000             
-    ##  sp7  0.4590 -0.1390  0.1470  0.1230  0.4940  0.2370  1.0000         
-    ##  sp8  0.2910  0.1520 -0.4220 -0.4120  0.2240 -0.0590  0.1370  1.0000     
-    ##  sp9 -0.0990 -0.0580  0.0720  0.0620 -0.3110 -0.2330 -0.2330 -0.0850  1.0000 
-    ##  sp10     0.2280  0.3710 -0.6010 -0.5630  0.2410  0.0750  0.1310  0.4160 -0.2220  1.0000
+    ##  sp2 -0.3460  1.0000                             
+    ##  sp3 -0.1400 -0.3770  1.0000                         
+    ##  sp4 -0.1410 -0.3500  0.7070  1.0000                     
+    ##  sp5  0.6090 -0.3260 -0.0870 -0.0940  1.0000                 
+    ##  sp6 -0.2180  0.3640  0.1600  0.1490 -0.0550  1.0000             
+    ##  sp7  0.4700 -0.1420  0.1320  0.1190  0.5120  0.2390  1.0000         
+    ##  sp8  0.2530  0.1540 -0.4360 -0.4220  0.2340 -0.0540  0.1100  1.0000     
+    ##  sp9 -0.1060 -0.0420  0.0510  0.0610 -0.2900 -0.2650 -0.2250 -0.1180  1.0000 
+    ##  sp10     0.2070  0.3830 -0.6020 -0.5750  0.2350  0.0820  0.1270  0.4150 -0.2100  1.0000
     ## 
     ## 
     ## 
     ## Spatial: 
-    ##             sp1        sp2       sp3         sp4       sp5       sp6      sp7
-    ## X1:X2 0.2097643 -0.5288849 0.5181209 0.009623682 0.3563328 0.1980719 0.423929
-    ##             sp8       sp9     sp10
-    ## X1:X2 0.4698656 0.2615229 0.257226
+    ##             sp1        sp2       sp3         sp4       sp5       sp6       sp7
+    ## X1:X2 0.1831122 -0.5116199 0.5413426 -0.00264241 0.3427651 0.1870106 0.4059039
+    ##             sp8       sp9      sp10
+    ## X1:X2 0.4804094 0.2659622 0.2485598
     ## 
     ## 
     ## 
     ##                  Estimate  Std.Err Z value Pr(>|z|)    
-    ## sp1 (Intercept)  -0.03977  0.21206   -0.19  0.85124    
-    ## sp1 X1            0.82571  0.41443    1.99  0.04633 *  
-    ## sp1 X2           -1.62219  0.40169   -4.04  5.4e-05 ***
-    ## sp1 X3           -0.13454  0.35189   -0.38  0.70221    
-    ## sp2 (Intercept)  -0.01476  0.21101   -0.07  0.94422    
-    ## sp2 X1            0.89789  0.40524    2.22  0.02671 *  
-    ## sp2 X2            0.19296  0.38875    0.50  0.61965    
-    ## sp2 X3            0.46829  0.35138    1.33  0.18263    
-    ## sp3 (Intercept)  -0.36290  0.22196   -1.64  0.10204    
-    ## sp3 X1            1.01730  0.40704    2.50  0.01245 *  
-    ## sp3 X2           -0.26439  0.40040   -0.66  0.50905    
-    ## sp3 X3           -0.72271  0.38358   -1.88  0.05955 .  
-    ## sp4 (Intercept)  -0.05324  0.20073   -0.27  0.79082    
-    ## sp4 X1           -1.02942  0.38118   -2.70  0.00692 ** 
-    ## sp4 X2           -1.29268  0.37599   -3.44  0.00059 ***
-    ## sp4 X3           -0.32540  0.33982   -0.96  0.33828    
-    ## sp5 (Intercept)  -0.14809  0.20236   -0.73  0.46427    
-    ## sp5 X1            0.50222  0.39101    1.28  0.19899    
-    ## sp5 X2            0.44313  0.36177    1.22  0.22062    
-    ## sp5 X3           -0.43606  0.34987   -1.25  0.21263    
-    ## sp6 (Intercept)   0.16414  0.20103    0.82  0.41422    
-    ## sp6 X1            1.65042  0.41523    3.97  7.0e-05 ***
-    ## sp6 X2           -0.66862  0.38461   -1.74  0.08213 .  
-    ## sp6 X3            0.16584  0.32957    0.50  0.61482    
-    ## sp7 (Intercept)  -0.00151  0.20626   -0.01  0.99416    
-    ## sp7 X1           -0.28922  0.39346   -0.74  0.46230    
-    ## sp7 X2            0.28814  0.37885    0.76  0.44691    
-    ## sp7 X3           -1.06155  0.34752   -3.05  0.00225 ** 
-    ## sp8 (Intercept)   0.13099  0.15214    0.86  0.38925    
-    ## sp8 X1            0.24814  0.28800    0.86  0.38890    
-    ## sp8 X2            0.30580  0.28384    1.08  0.28131    
-    ## sp8 X3           -1.01204  0.26964   -3.75  0.00017 ***
-    ## sp9 (Intercept)   0.03149  0.17057    0.18  0.85352    
-    ## sp9 X1            1.07060  0.33145    3.23  0.00124 ** 
-    ## sp9 X2           -0.79677  0.32507   -2.45  0.01424 *  
-    ## sp9 X3            0.59462  0.28658    2.07  0.03800 *  
-    ## sp10 (Intercept) -0.06706  0.18177   -0.37  0.71219    
-    ## sp10 X1          -0.43999  0.33907   -1.30  0.19441    
-    ## sp10 X2          -0.94309  0.32855   -2.87  0.00410 ** 
-    ## sp10 X3          -0.35444  0.30660   -1.16  0.24766    
+    ## sp1 (Intercept)  -0.02294  0.20454   -0.11  0.91068    
+    ## sp1 X1            0.83191  0.40909    2.03  0.04199 *  
+    ## sp1 X2           -1.63252  0.37725   -4.33  1.5e-05 ***
+    ## sp1 X3           -0.12155  0.33693   -0.36  0.71828    
+    ## sp2 (Intercept)  -0.03995  0.22700   -0.18  0.86029    
+    ## sp2 X1            0.89410  0.41314    2.16  0.03045 *  
+    ## sp2 X2            0.15840  0.40862    0.39  0.69827    
+    ## sp2 X3            0.43637  0.39043    1.12  0.26371    
+    ## sp3 (Intercept)  -0.34723  0.22196   -1.56  0.11773    
+    ## sp3 X1            1.03418  0.40596    2.55  0.01085 *  
+    ## sp3 X2           -0.21915  0.40775   -0.54  0.59095    
+    ## sp3 X3           -0.68271  0.38546   -1.77  0.07654 .  
+    ## sp4 (Intercept)  -0.04001  0.20259   -0.20  0.84345    
+    ## sp4 X1           -1.02553  0.38775   -2.64  0.00817 ** 
+    ## sp4 X2           -1.27792  0.38960   -3.28  0.00104 ** 
+    ## sp4 X3           -0.27332  0.34129   -0.80  0.42323    
+    ## sp5 (Intercept)  -0.13282  0.19686   -0.67  0.49986    
+    ## sp5 X1            0.45864  0.37292    1.23  0.21875    
+    ## sp5 X2            0.42704  0.36548    1.17  0.24264    
+    ## sp5 X3           -0.40132  0.34380   -1.17  0.24308    
+    ## sp6 (Intercept)   0.17045  0.19915    0.86  0.39206    
+    ## sp6 X1            1.68681  0.40498    4.17  3.1e-05 ***
+    ## sp6 X2           -0.66825  0.38106   -1.75  0.07949 .  
+    ## sp6 X3            0.17839  0.33066    0.54  0.58954    
+    ## sp7 (Intercept)   0.00701  0.20640    0.03  0.97289    
+    ## sp7 X1           -0.28037  0.40207   -0.70  0.48560    
+    ## sp7 X2            0.29147  0.37519    0.78  0.43725    
+    ## sp7 X3           -1.06023  0.35214   -3.01  0.00261 ** 
+    ## sp8 (Intercept)   0.13752  0.14967    0.92  0.35816    
+    ## sp8 X1            0.26011  0.28608    0.91  0.36323    
+    ## sp8 X2            0.27259  0.27891    0.98  0.32840    
+    ## sp8 X3           -1.01502  0.26326   -3.86  0.00012 ***
+    ## sp9 (Intercept)   0.04048  0.16892    0.24  0.81061    
+    ## sp9 X1            1.10252  0.32760    3.37  0.00076 ***
+    ## sp9 X2           -0.81973  0.31488   -2.60  0.00923 ** 
+    ## sp9 X3            0.59836  0.28592    2.09  0.03637 *  
+    ## sp10 (Intercept) -0.09065  0.18045   -0.50  0.61542    
+    ## sp10 X1          -0.43825  0.33418   -1.31  0.18972    
+    ## sp10 X2          -0.94811  0.32707   -2.90  0.00375 ** 
+    ## sp10 X3          -0.37165  0.30311   -1.23  0.22015    
     ## ---
     ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 
@@ -163,78 +163,76 @@ summary(model)
 plot(model)
 ```
 
-    ## LogLik:  -523.0463 
+    ## LogLik:  -517.1162 
     ## Regularization loss:  0 
     ## 
     ## Species-species correlation matrix: 
     ## 
     ##  sp1  1.0000                                 
-    ##  sp2 -0.3370  1.0000                             
-    ##  sp3 -0.1490 -0.3520  1.0000                         
-    ##  sp4 -0.1690 -0.3120  0.6950  1.0000                     
-    ##  sp5  0.6020 -0.3360 -0.0950 -0.1020  1.0000                 
-    ##  sp6 -0.2250  0.3830  0.1630  0.1670 -0.0940  1.0000             
-    ##  sp7  0.4590 -0.1390  0.1470  0.1230  0.4940  0.2370  1.0000         
-    ##  sp8  0.2910  0.1520 -0.4220 -0.4120  0.2240 -0.0590  0.1370  1.0000     
-    ##  sp9 -0.0990 -0.0580  0.0720  0.0620 -0.3110 -0.2330 -0.2330 -0.0850  1.0000 
-    ##  sp10     0.2280  0.3710 -0.6010 -0.5630  0.2410  0.0750  0.1310  0.4160 -0.2220  1.0000
+    ##  sp2 -0.3460  1.0000                             
+    ##  sp3 -0.1400 -0.3770  1.0000                         
+    ##  sp4 -0.1410 -0.3500  0.7070  1.0000                     
+    ##  sp5  0.6090 -0.3260 -0.0870 -0.0940  1.0000                 
+    ##  sp6 -0.2180  0.3640  0.1600  0.1490 -0.0550  1.0000             
+    ##  sp7  0.4700 -0.1420  0.1320  0.1190  0.5120  0.2390  1.0000         
+    ##  sp8  0.2530  0.1540 -0.4360 -0.4220  0.2340 -0.0540  0.1100  1.0000     
+    ##  sp9 -0.1060 -0.0420  0.0510  0.0610 -0.2900 -0.2650 -0.2250 -0.1180  1.0000 
+    ##  sp10     0.2070  0.3830 -0.6020 -0.5750  0.2350  0.0820  0.1270  0.4150 -0.2100  1.0000
     ## 
     ## 
     ## 
     ## Spatial: 
-    ##             sp1        sp2       sp3         sp4       sp5       sp6      sp7
-    ## X1:X2 0.2097643 -0.5288849 0.5181209 0.009623682 0.3563328 0.1980719 0.423929
-    ##             sp8       sp9     sp10
-    ## X1:X2 0.4698656 0.2615229 0.257226
+    ##             sp1        sp2       sp3         sp4       sp5       sp6       sp7
+    ## X1:X2 0.1831122 -0.5116199 0.5413426 -0.00264241 0.3427651 0.1870106 0.4059039
+    ##             sp8       sp9      sp10
+    ## X1:X2 0.4804094 0.2659622 0.2485598
     ## 
     ## 
     ## 
     ##                  Estimate  Std.Err Z value Pr(>|z|)    
-    ## sp1 (Intercept)  -0.03977  0.21206   -0.19  0.85124    
-    ## sp1 X1            0.82571  0.41443    1.99  0.04633 *  
-    ## sp1 X2           -1.62219  0.40169   -4.04  5.4e-05 ***
-    ## sp1 X3           -0.13454  0.35189   -0.38  0.70221    
-    ## sp2 (Intercept)  -0.01476  0.21101   -0.07  0.94422    
-    ## sp2 X1            0.89789  0.40524    2.22  0.02671 *  
-    ## sp2 X2            0.19296  0.38875    0.50  0.61965    
-    ## sp2 X3            0.46829  0.35138    1.33  0.18263    
-    ## sp3 (Intercept)  -0.36290  0.22196   -1.64  0.10204    
-    ## sp3 X1            1.01730  0.40704    2.50  0.01245 *  
-    ## sp3 X2           -0.26439  0.40040   -0.66  0.50905    
-    ## sp3 X3           -0.72271  0.38358   -1.88  0.05955 .  
-    ## sp4 (Intercept)  -0.05324  0.20073   -0.27  0.79082    
-    ## sp4 X1           -1.02942  0.38118   -2.70  0.00692 ** 
-    ## sp4 X2           -1.29268  0.37599   -3.44  0.00059 ***
-    ## sp4 X3           -0.32540  0.33982   -0.96  0.33828    
-    ## sp5 (Intercept)  -0.14809  0.20236   -0.73  0.46427    
-    ## sp5 X1            0.50222  0.39101    1.28  0.19899    
-    ## sp5 X2            0.44313  0.36177    1.22  0.22062    
-    ## sp5 X3           -0.43606  0.34987   -1.25  0.21263    
-    ## sp6 (Intercept)   0.16414  0.20103    0.82  0.41422    
-    ## sp6 X1            1.65042  0.41523    3.97  7.0e-05 ***
-    ## sp6 X2           -0.66862  0.38461   -1.74  0.08213 .  
-    ## sp6 X3            0.16584  0.32957    0.50  0.61482    
-    ## sp7 (Intercept)  -0.00151  0.20626   -0.01  0.99416    
-    ## sp7 X1           -0.28922  0.39346   -0.74  0.46230    
-    ## sp7 X2            0.28814  0.37885    0.76  0.44691    
-    ## sp7 X3           -1.06155  0.34752   -3.05  0.00225 ** 
-    ## sp8 (Intercept)   0.13099  0.15214    0.86  0.38925    
-    ## sp8 X1            0.24814  0.28800    0.86  0.38890    
-    ## sp8 X2            0.30580  0.28384    1.08  0.28131    
-    ## sp8 X3           -1.01204  0.26964   -3.75  0.00017 ***
-    ## sp9 (Intercept)   0.03149  0.17057    0.18  0.85352    
-    ## sp9 X1            1.07060  0.33145    3.23  0.00124 ** 
-    ## sp9 X2           -0.79677  0.32507   -2.45  0.01424 *  
-    ## sp9 X3            0.59462  0.28658    2.07  0.03800 *  
-    ## sp10 (Intercept) -0.06706  0.18177   -0.37  0.71219    
-    ## sp10 X1          -0.43999  0.33907   -1.30  0.19441    
-    ## sp10 X2          -0.94309  0.32855   -2.87  0.00410 ** 
-    ## sp10 X3          -0.35444  0.30660   -1.16  0.24766    
+    ## sp1 (Intercept)  -0.02294  0.20454   -0.11  0.91068    
+    ## sp1 X1            0.83191  0.40909    2.03  0.04199 *  
+    ## sp1 X2           -1.63252  0.37725   -4.33  1.5e-05 ***
+    ## sp1 X3           -0.12155  0.33693   -0.36  0.71828    
+    ## sp2 (Intercept)  -0.03995  0.22700   -0.18  0.86029    
+    ## sp2 X1            0.89410  0.41314    2.16  0.03045 *  
+    ## sp2 X2            0.15840  0.40862    0.39  0.69827    
+    ## sp2 X3            0.43637  0.39043    1.12  0.26371    
+    ## sp3 (Intercept)  -0.34723  0.22196   -1.56  0.11773    
+    ## sp3 X1            1.03418  0.40596    2.55  0.01085 *  
+    ## sp3 X2           -0.21915  0.40775   -0.54  0.59095    
+    ## sp3 X3           -0.68271  0.38546   -1.77  0.07654 .  
+    ## sp4 (Intercept)  -0.04001  0.20259   -0.20  0.84345    
+    ## sp4 X1           -1.02553  0.38775   -2.64  0.00817 ** 
+    ## sp4 X2           -1.27792  0.38960   -3.28  0.00104 ** 
+    ## sp4 X3           -0.27332  0.34129   -0.80  0.42323    
+    ## sp5 (Intercept)  -0.13282  0.19686   -0.67  0.49986    
+    ## sp5 X1            0.45864  0.37292    1.23  0.21875    
+    ## sp5 X2            0.42704  0.36548    1.17  0.24264    
+    ## sp5 X3           -0.40132  0.34380   -1.17  0.24308    
+    ## sp6 (Intercept)   0.17045  0.19915    0.86  0.39206    
+    ## sp6 X1            1.68681  0.40498    4.17  3.1e-05 ***
+    ## sp6 X2           -0.66825  0.38106   -1.75  0.07949 .  
+    ## sp6 X3            0.17839  0.33066    0.54  0.58954    
+    ## sp7 (Intercept)   0.00701  0.20640    0.03  0.97289    
+    ## sp7 X1           -0.28037  0.40207   -0.70  0.48560    
+    ## sp7 X2            0.29147  0.37519    0.78  0.43725    
+    ## sp7 X3           -1.06023  0.35214   -3.01  0.00261 ** 
+    ## sp8 (Intercept)   0.13752  0.14967    0.92  0.35816    
+    ## sp8 X1            0.26011  0.28608    0.91  0.36323    
+    ## sp8 X2            0.27259  0.27891    0.98  0.32840    
+    ## sp8 X3           -1.01502  0.26326   -3.86  0.00012 ***
+    ## sp9 (Intercept)   0.04048  0.16892    0.24  0.81061    
+    ## sp9 X1            1.10252  0.32760    3.37  0.00076 ***
+    ## sp9 X2           -0.81973  0.31488   -2.60  0.00923 ** 
+    ## sp9 X3            0.59836  0.28592    2.09  0.03637 *  
+    ## sp10 (Intercept) -0.09065  0.18045   -0.50  0.61542    
+    ## sp10 X1          -0.43825  0.33418   -1.31  0.18972    
+    ## sp10 X2          -0.94811  0.32707   -2.90  0.00375 ** 
+    ## sp10 X3          -0.37165  0.30311   -1.23  0.22015    
     ## ---
     ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 
-![](README_files/figure-gfm/unnamed-chunk-5-1.png)<!-- -->
-
 We also support also other response families: Count data:
 
 ``` r
@@ -252,20 +250,27 @@ associations, and space) on the occurences:
 
 ``` r
 imp = importance(model)
+```
+
+    ## Warning: 'importance' is deprecated.
+    ## Use 'plot(anova(x, internal=TRUE))' instead.
+    ## See help("Deprecated")
+
+``` r
 print(imp)
 ```
 
     ##    sp       env      spatial     biotic
-    ## 1   1 0.8539057 3.194661e-04 0.14577481
-    ## 2   2 0.5897330 3.662053e-03 0.40660492
-    ## 3   3 0.7400899 3.535578e-03 0.25637448
-    ## 4   4 0.9069016 7.076025e-07 0.09309772
-    ## 5   5 0.5101848 2.644762e-03 0.48717043
-    ## 6   6 0.8644908 3.016936e-04 0.13520756
-    ## 7   7 0.5335647 1.650771e-03 0.46478456
-    ## 8   8 0.8930446 3.978094e-03 0.10297730
-    ## 9   9 0.7327270 6.191073e-04 0.26665390
-    ## 10 10 0.9761007 1.266072e-03 0.02263325
+    ## 1   1 0.8661759 2.439290e-04 0.13358012
+    ## 2   2 0.5512735 3.389213e-03 0.44533733
+    ## 3   3 0.7243540 3.892151e-03 0.27175385
+    ## 4   4 0.8974715 5.453921e-08 0.10252855
+    ## 5   5 0.4491234 2.478703e-03 0.54839790
+    ## 6   6 0.8595459 2.568011e-04 0.14019726
+    ## 7   7 0.5315694 1.515744e-03 0.46691483
+    ## 8   8 0.9476700 4.443077e-03 0.04788696
+    ## 9   9 0.7424834 6.183612e-04 0.25689825
+    ## 10 10 0.9604391 1.141932e-03 0.03841905
 
 ``` r
 plot(imp)
@@ -288,9 +293,9 @@ print(an)
     ## Terms added sequentially:
     ## 
     ##          Deviance Residual deviance R2 Nagelkerke R2 McFadden
-    ## Biotic  151.02063         521.48694       0.77914      0.1092
-    ## Abiotic 192.75159         328.73535       0.96786      0.2487
-    ## Spatial  -7.40771         336.14307       0.96539      0.2433
+    ## Biotic  153.69963         518.82466       0.78497      0.1112
+    ## Abiotic 180.74306         338.08160       0.96472      0.2419
+    ## Spatial  12.84565         325.23595       0.96897      0.2512
 
 ``` r
 plot(an)
@@ -301,6 +306,23 @@ plot(an)
 The anova shows the relative changes in the deviance of the groups and
 their intersections.
 
+We can also visualize the individual contributions to the species and
+site
+*R*<sup>2</sup>
+:
+
+``` r
+plot(an, internal=TRUE)
+```
+
+    ## Registered S3 methods overwritten by 'ggtern':
+    ##   method           from   
+    ##   grid.draw.ggplot ggplot2
+    ##   plot.ggplot      ggplot2
+    ##   print.ggplot     ggplot2
+
+![](README_files/figure-gfm/unnamed-chunk-9-1.png)<!-- -->
+
 If it fails, check out the help of ?install_sjSDM, ?installation_help,
 and vignette(“Dependencies”, package = “sjSDM”).
 
@@ -336,7 +358,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.358]Iter: 1/10  10%|#         | [00:00,  5.12it/s, loss=7.358]Iter: 1/10  10%|#         | [00:00,  5.12it/s, loss=7.331]Iter: 1/10  10%|#         | [00:00,  5.12it/s, loss=7.323]Iter: 1/10  10%|#         | [00:00,  5.12it/s, loss=7.324]Iter: 1/10  10%|#         | [00:00,  5.12it/s, loss=7.334]Iter: 1/10  10%|#         | [00:00,  5.12it/s, loss=7.318]Iter: 6/10  60%|######    | [00:00, 23.86it/s, loss=7.318]Iter: 6/10  60%|######    | [00:00, 23.86it/s, loss=7.326]Iter: 6/10  60%|######    | [00:00, 23.86it/s, loss=7.309]Iter: 6/10  60%|######    | [00:00, 23.86it/s, loss=7.304]Iter: 6/10  60%|######    | [00:00, 23.86it/s, loss=7.286]Iter: 10/10 100%|##########| [00:00, 26.03it/s, loss=7.286]
+    ## Iter: 0/10   0%|          | [00:00, ?it/s]Iter: 0/10   0%|          | [00:00, ?it/s, loss=7.233]Iter: 1/10  10%|#         | [00:00,  5.80it/s, loss=7.233]Iter: 1/10  10%|#         | [00:00,  5.80it/s, loss=7.237]Iter: 1/10  10%|#         | [00:00,  5.80it/s, loss=7.222]Iter: 1/10  10%|#         | [00:00,  5.80it/s, loss=7.204]Iter: 1/10  10%|#         | [00:00,  5.80it/s, loss=7.22] Iter: 1/10  10%|#         | [00:00,  5.80it/s, loss=7.2] Iter: 1/10  10%|#         | [00:00,  5.80it/s, loss=7.188]Iter: 7/10  70%|#######   | [00:00, 30.03it/s, loss=7.188]Iter: 7/10  70%|#######   | [00:00, 30.03it/s, loss=7.196]Iter: 7/10  70%|#######   | [00:00, 30.03it/s, loss=7.187]Iter: 7/10  70%|#######   | [00:00, 30.03it/s, loss=7.192]Iter: 10/10 100%|##########| [00:00, 30.66it/s, loss=7.192]
 
 Calculate Importance:
 
@@ -348,25 +370,25 @@ covX = fa.covariance( torch.tensor(Env).t() ).data.cpu().numpy()
 fa.importance(beta=Beta, covX=covX, sigma=Sigma)
 ```
 
-    ## {'env': array([[ 6.00627158e-04,  1.14630593e-03,  1.09053524e-02,
-    ##          1.05654253e-02,  1.45557038e-02],
-    ##        [ 9.56870615e-03,  4.21576016e-03,  1.06933614e-04,
-    ##          1.33172970e-03,  1.05470214e-02],
-    ##        [ 2.32625613e-03,  1.61409937e-02,  6.64134495e-05,
-    ##          9.74757495e-05,  2.18663458e-03],
-    ##        [ 8.52300692e-03,  1.07057637e-03,  1.46411667e-02,
-    ##          8.68435390e-03,  7.08195148e-03],
-    ##        [-7.17847797e-06, -1.23505345e-06,  1.45622957e-02,
-    ##          3.56136751e-03,  4.91820509e-03],
-    ##        [-1.31524590e-04,  1.22741107e-02,  1.22232996e-02,
-    ##          7.41305854e-03,  1.47087574e-02],
-    ##        [ 1.26086129e-02,  1.54245849e-04,  6.33575348e-03,
-    ##          3.37137189e-03,  9.59621649e-03],
-    ##        [ 8.97107311e-05,  5.07542631e-03,  9.64390172e-04,
-    ##          1.04126502e-02,  1.18228560e-02],
-    ##        [ 9.25059337e-03,  5.78574790e-03,  9.44872445e-05,
-    ##          8.52529984e-03,  3.14252358e-03],
-    ##        [ 1.40638428e-03,  6.01684162e-03,  4.81087307e-04,
-    ##          1.57201290e-02,  1.10209584e-02]], dtype=float32), 'biotic': array([0.9622266 , 0.9742299 , 0.9791822 , 0.95999897, 0.9769665 ,
-    ##        0.95351225, 0.96793383, 0.971635  , 0.97320133, 0.9653546 ],
+    ## {'env': array([[ 8.01924523e-03,  5.88351442e-03,  3.46450630e-04,
+    ##          2.41728313e-03,  2.93385610e-03],
+    ##        [ 2.51798570e-04,  1.30314089e-03,  4.08564508e-03,
+    ##          8.43194837e-04,  1.00540789e-02],
+    ##        [ 4.58702561e-04,  6.93971897e-03,  9.88384054e-05,
+    ##          1.51673518e-03,  4.04162938e-03],
+    ##        [-1.49338421e-05,  3.24389548e-05,  1.47342135e-03,
+    ##          3.60254053e-05,  7.97839928e-03],
+    ##        [ 7.79314898e-04,  1.94915745e-04,  6.61612488e-03,
+    ##          1.23276515e-02,  3.61960870e-03],
+    ##        [ 3.79100814e-03,  9.13846270e-06,  7.05215894e-03,
+    ##          8.15463625e-03,  2.77351518e-03],
+    ##        [ 3.37206497e-04,  9.80060431e-04,  6.51973824e-05,
+    ##          1.43132071e-04,  1.15610389e-02],
+    ##        [ 6.48319907e-03,  6.12059329e-03,  9.50422324e-03,
+    ##          8.78254510e-03,  9.26101953e-03],
+    ##        [ 8.99332296e-03,  1.01236743e-03,  3.25980218e-04,
+    ##          1.23339929e-02,  1.86911039e-03],
+    ##        [ 8.95618182e-03,  1.53849716e-03,  1.09001040e-03,
+    ##          6.02281233e-03,  7.11724581e-03]], dtype=float32), 'biotic': array([0.98039967, 0.98346215, 0.9869444 , 0.9904946 , 0.9764624 ,
+    ##        0.97821957, 0.9869134 , 0.9598484 , 0.97546524, 0.9752752 ],
     ##       dtype=float32)}
diff --git a/README_files/figure-gfm/unnamed-chunk-7-1.png b/README_files/figure-gfm/unnamed-chunk-7-1.png
index 185af38..9bb39f8 100644
Binary files a/README_files/figure-gfm/unnamed-chunk-7-1.png and b/README_files/figure-gfm/unnamed-chunk-7-1.png differ
diff --git a/README_files/figure-gfm/unnamed-chunk-8-1.png b/README_files/figure-gfm/unnamed-chunk-8-1.png
index 92e2d3c..348af80 100644
Binary files a/README_files/figure-gfm/unnamed-chunk-8-1.png and b/README_files/figure-gfm/unnamed-chunk-8-1.png differ
diff --git a/README_files/figure-gfm/unnamed-chunk-9-1.png b/README_files/figure-gfm/unnamed-chunk-9-1.png
index 41f2ac4..1a0e9d0 100644
Binary files a/README_files/figure-gfm/unnamed-chunk-9-1.png and b/README_files/figure-gfm/unnamed-chunk-9-1.png differ
diff --git a/sjSDM/DESCRIPTION b/sjSDM/DESCRIPTION
index 341338d..d57f103 100644
--- a/sjSDM/DESCRIPTION
+++ b/sjSDM/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sjSDM
 Type: Package
 Title: Scalable Joint Species Distribution Modeling
-Version: 1.0.0
-Date: 2022-01-05
+Version: 1.0.1
+Date: 2022-03-11
 Authors@R: 
     c(person(given = "Maximilian",
              family = "Pichler",
@@ -41,7 +41,8 @@ Imports:
     crayon,
     ggplot2,
     checkmate,
-    mathjaxr
+    mathjaxr,
+    ggtern
 Suggests:
     testthat,
     knitr,
diff --git a/sjSDM/NAMESPACE b/sjSDM/NAMESPACE
index 6e6aca6..fa78405 100644
--- a/sjSDM/NAMESPACE
+++ b/sjSDM/NAMESPACE
@@ -1,54 +1,54 @@
-# Generated by roxygen2: do not edit by hand
-
-S3method(anova,sjSDM)
-S3method(coef,sjSDM)
-S3method(getCov,sjSDM)
-S3method(getWeights,sjSDM)
-S3method(logLik,sjSDM)
-S3method(plot,sjSDM)
-S3method(plot,sjSDM_cv)
-S3method(plot,sjSDManova)
-S3method(plot,sjSDMimportance)
-S3method(predict,sjSDM)
-S3method(print,DNN)
-S3method(print,bioticStruct)
-S3method(print,linear)
-S3method(print,sjSDM)
-S3method(print,sjSDM_cv)
-S3method(print,sjSDManova)
-S3method(print,sjSDMimportance)
-S3method(setWeights,sjSDM)
-S3method(simulate,sjSDM)
-S3method(summary,sjSDM)
-S3method(summary,sjSDM_cv)
-S3method(update,sjSDM)
-export(AccSGD)
-export(AdaBound)
-export(Adamax)
-export(DNN)
-export(RMSprop)
-export(Rsquared)
-export(SGD)
-export(bioticStruct)
-export(generateSpatialEV)
-export(getCov)
-export(getSe)
-export(getWeights)
-export(importance)
-export(install_diagnostic)
-export(install_sjSDM)
-export(is_torch_available)
-export(linear)
-export(madgrad)
-export(setWeights)
-export(simulate_SDM)
-export(sjSDM)
-export(sjSDM.tune)
-export(sjSDMControl)
-export(sjSDM_cv)
-import(checkmate)
-import(ggplot2)
-import(graphics)
-import(mathjaxr)
-import(stats)
-importFrom(stats,simulate)
+# Generated by roxygen2: do not edit by hand
+
+S3method(anova,sjSDM)
+S3method(coef,sjSDM)
+S3method(getCov,sjSDM)
+S3method(getWeights,sjSDM)
+S3method(logLik,sjSDM)
+S3method(plot,sjSDM)
+S3method(plot,sjSDM_cv)
+S3method(plot,sjSDManova)
+S3method(plot,sjSDMimportance)
+S3method(predict,sjSDM)
+S3method(print,DNN)
+S3method(print,bioticStruct)
+S3method(print,linear)
+S3method(print,sjSDM)
+S3method(print,sjSDM_cv)
+S3method(print,sjSDManova)
+S3method(print,sjSDMimportance)
+S3method(setWeights,sjSDM)
+S3method(simulate,sjSDM)
+S3method(summary,sjSDM)
+S3method(summary,sjSDM_cv)
+S3method(update,sjSDM)
+export(AccSGD)
+export(AdaBound)
+export(Adamax)
+export(DNN)
+export(RMSprop)
+export(Rsquared)
+export(SGD)
+export(bioticStruct)
+export(generateSpatialEV)
+export(getCov)
+export(getSe)
+export(getWeights)
+export(importance)
+export(install_diagnostic)
+export(install_sjSDM)
+export(is_torch_available)
+export(linear)
+export(madgrad)
+export(setWeights)
+export(simulate_SDM)
+export(sjSDM)
+export(sjSDM.tune)
+export(sjSDMControl)
+export(sjSDM_cv)
+import(checkmate)
+import(ggplot2)
+import(graphics)
+import(mathjaxr)
+import(stats)
+importFrom(stats,simulate)
diff --git a/sjSDM/NEWS.md b/sjSDM/NEWS.md
index 8046acf..666f686 100644
--- a/sjSDM/NEWS.md
+++ b/sjSDM/NEWS.md
@@ -1,3 +1,27 @@
+# sjSDM 1.0.1
+
+## New Features
+
+* anova plots for internal meta-community structure (based on individual R-squared values)
+
+## Minor changes
+
+* first layer of DNN now always without an explicit bias (bias/intercept is passed by model/formula, if desired)
+* revised prediction function, improved stability
+* revised simulation function, samples now from a multivariate probit model
+
+## Bug fixes
+
+* unlisting of config objects in `sjSDM::sjSDM_cv` (thanks to Máté) (added unit tests)  #88
+* `sjSDM::Rsquared` bug for spatial models (thanks to Máté) #90
+* revised regularization behavior, l1 and l2 were not correctly imposed on DNN structure
+* revised and improved setWeights function
+* bugs in vignettes (thanks to Doug) #92
+* bugs in plot function for models with DNN objects
+
+
+
+
 # sjSDM 1.0.0
 
 ## Major changes
diff --git a/sjSDM/R/Rsquared.R b/sjSDM/R/Rsquared.R
index a0c97c5..e3f1196 100644
--- a/sjSDM/R/Rsquared.R
+++ b/sjSDM/R/Rsquared.R
@@ -28,11 +28,10 @@ Rsquared = function(model, method = c("Nagelkerke", "McFadden")) {
   
   if(is.null(model$Null)) {
     if(inherits(model, "spatial")) {
-      spatial_formula=stats::as.formula(~0)
+      model$Null = update(model, env_formula=~1, spatial_formula=~0)
     } else {
-      spatial_formula=NULL
+      model$Null = update(model, env_formula=~1)
     }
-    model$Null = update(model, env_formula=~1, spatial_formula=spatial_formula)
   }
   N0 = -logLik(model$Null)[[1]]
   N1 = -logLik(model)[[1]]
diff --git a/sjSDM/R/anova.R b/sjSDM/R/anova.R
index b306df9..d10f558 100644
--- a/sjSDM/R/anova.R
+++ b/sjSDM/R/anova.R
@@ -19,7 +19,10 @@
 #' \item{results}{Data frame of results.}
 #' \item{to_print}{Data frame, summarized results for type I anova.}
 #' \item{N}{Number of observations (sites).}
-#' \item{spatial}{Logical, spatial model or not}
+#' \item{spatial}{Logical, spatial model or not.}
+#' \item{species}{individual species R2s.}
+#' \item{sites}{individual site R2s.}
+#' \item{lls}{individual site by species negative-log-likelihood values.}
 #' 
 #' Implemented S3 methods are \code{\link{print.sjSDManova}} and \code{\link{plot.sjSDManova}}
 #'  
@@ -29,15 +32,22 @@
 
 anova.sjSDM = function(object, ...) {
   out = list()
-  full_m = logLik(object)[[1]]
+  individual = TRUE
+  
+  if(object$family$family$family == "gaussian") stop("gaussian not yet supported")
+  
+  full_m = logLik(object, individual=individual)[[1]]
+  object$settings$se = FALSE
   ### fit different models ###
   #e_form = stats::as.formula(paste0(as.character(object$settings$env$formula), collapse = ""))
   if(!inherits(object, "spatial")) {
     out$spatial = FALSE
-    null_m = logLik(update(object, env_formula = ~1, biotic=bioticStruct(diag = TRUE )))[[1]]
-    A_m =    logLik(update(object, env_formula = NULL, biotic=bioticStruct(diag = TRUE )))[[1]]
-    B_m =    logLik(update(object, env_formula = ~1, biotic=bioticStruct(diag = FALSE)))[[1]]
-    SAT_m =  logLik(update(object, env_formula = ~as.factor(1:nrow(object$data$X)), biotic=bioticStruct(diag = TRUE )))[[1]]
+    
+    
+    null_m = logLik(update(object, env_formula = ~1, biotic=bioticStruct(diag = TRUE )), individual=individual)[[1]]
+    A_m =    logLik(update(object, env_formula = NULL, biotic=bioticStruct(diag = TRUE )), individual=individual)[[1]]
+    B_m =    logLik(update(object, env_formula = ~1, biotic=bioticStruct(diag = FALSE)), individual=individual)[[1]]
+    SAT_m =  logLik(update(object, env_formula = ~as.factor(1:nrow(object$data$X)), biotic=bioticStruct(diag = TRUE )), individual=individual)[[1]]
     
     A =    A_m
     B =    B_m
@@ -48,18 +58,21 @@ anova.sjSDM = function(object, ...) {
     anova_rows = c("Null", "B", "Full")
     names(anova_rows) = c("Null", "Biotic", "Abiotic")
     results = data.frame(models = c("A", "B","AB", "Full", "Saturated", "Null"),
-                         ll = -c(A, B,AB, full, sat, null))
+                         ll = -c(sum(A), sum(B), sum(AB), sum(full), sum(sat), sum(null)))
+    results_ind = list(A = -A, B = -B, AB = - AB, Full = -full, Saturated = -sat, Null = - null)
+    results_predict = list(A = -A, B = -B, AB = - AB, Full = -full, Saturated = -sat, Null = - null)
+    
   } else {
     out$spatial = TRUE
     s_form = stats::as.formula(paste0(as.character(object$settings$spatial$formula), collapse = ""))
-    null_m = logLik(update(object, env_formula = ~1, spatial_formula= ~0, biotic=bioticStruct(diag = TRUE )))[[1]]
-    A_m =    logLik(update(object, env_formula = NULL, spatial_formula= ~0, biotic=bioticStruct(diag = TRUE )))[[1]]
-    B_m =    logLik(update(object, env_formula = ~1, spatial_formula= ~0, biotic=bioticStruct(diag = FALSE)))[[1]]
-    S_m =    logLik(update(object, env_formula = ~1, spatial_formula= NULL, biotic=bioticStruct(diag = TRUE)))[[1]]
-    AB_m =   logLik(update(object, env_formula = NULL, spatial_formula= ~0, biotic=bioticStruct(diag = FALSE )))[[1]]
-    AS_m =   logLik(update(object, env_formula = NULL, spatial_formula= NULL, biotic=bioticStruct(diag = FALSE )))[[1]]
-    BS_m =   logLik(update(object, env_formula = ~1, spatial_formula= NULL, biotic=bioticStruct(diag = TRUE )))[[1]]
-    SAT_m =  logLik(update(object, env_formula = ~as.factor(1:nrow(object$data$X)), spatial_formula= ~0, biotic=bioticStruct(diag = TRUE )))[[1]]
+    null_m = logLik(update(object, env_formula = ~1, spatial_formula= ~0, biotic=bioticStruct(diag = TRUE )), individual=individual)[[1]]
+    A_m =    logLik(update(object, env_formula = NULL, spatial_formula= ~0, biotic=bioticStruct(diag = TRUE )), individual=individual)[[1]]
+    B_m =    logLik(update(object, env_formula = ~1, spatial_formula= ~0, biotic=bioticStruct(diag = FALSE)), individual=individual)[[1]]
+    S_m =    logLik(update(object, env_formula = ~1, spatial_formula= NULL, biotic=bioticStruct(diag = TRUE)), individual=individual)[[1]]
+    AB_m =   logLik(update(object, env_formula = NULL, spatial_formula= ~0, biotic=bioticStruct(diag = FALSE )), individual=individual)[[1]]
+    AS_m =   logLik(update(object, env_formula = NULL, spatial_formula= NULL, biotic=bioticStruct(diag = FALSE )), individual=individual)[[1]]
+    BS_m =   logLik(update(object, env_formula = ~1, spatial_formula= NULL, biotic=bioticStruct(diag = TRUE )), individual=individual)[[1]]
+    SAT_m =  logLik(update(object, env_formula = ~as.factor(1:nrow(object$data$X)), spatial_formula= ~0, biotic=bioticStruct(diag = TRUE )), individual=individual)[[1]]
     
     A =    A_m
     B =    B_m
@@ -72,11 +85,12 @@ anova.sjSDM = function(object, ...) {
     null = null_m
     sat =  SAT_m
     
+    
     results = data.frame(models = c("A", "B","S","B+A","B+S","A+S","AB","AS", "BS", "ABS", "Full", "Saturated", "Null"),
-                         ll = -c(A, B,S,AB_m,BS_m, AS_m, AB, AS, BS, ABS, full, sat, null))
+                         ll = -c(sum(A), sum(B),sum(S),sum(AB_m),sum(BS_m), sum(AS_m), sum(AB), sum(AS), sum(BS), sum(ABS), sum(full), sum(sat), sum(null)))
     
+    results_ind = list("A"=-A, "B"=-B,"S"=-S,"B+A"=-AB_m,"B+S"=-BS_m,"A+S"=-AS_m,"AB"=-AB,"AS"=-AS, "BS"=-BS, "ABS"=-ABS, "Full"=-full, "Saturated"=-sat, "Null"=-null)
     
-    #1-exp(2/nrow(object$data$Y)*(-A))
     anova_rows = c("Null", "B", "B+A", "Full")
     names(anova_rows) = c("Null", "Biotic", "Abiotic", "Spatial")
   }
@@ -84,10 +98,21 @@ anova.sjSDM = function(object, ...) {
   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`
-  R2 = function(a, b) return(1-exp(2/(nrow(object$data$Y))*(-a+b)))
-  results$`R2 Nagelkerke` = R2(rep(-results$ll[which(results$models == "Null", arr.ind = TRUE)], length(results$ll)), - results$ll)
-  R2 = function(a, b) 1 - (b/a)
-  results$`R2 McFadden`= R2(rep(-results$ll[which(results$models == "Null", arr.ind = TRUE)], length(results$ll)), - results$ll)
+  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)
+  
+  # individual
+  
+  Residual_deviance_ind = lapply(results_ind, function(r) r - results_ind$Saturated)
+  Deviance_ind = lapply(Residual_deviance_ind, function(r) r - Residual_deviance_ind$Null)
+  R211 = function(a, b, n=1) return(1-exp(2/(n)*(-a+b)))   # divide by what?
+  R2_Nagelkerke_ind = lapply(results_ind, function(r) R211(-colSums(results_ind$Null), -colSums(r), n=nrow(object$data$Y)))
+  R2_Nagelkerke_sites = lapply(results_ind, function(r) R211(-rowSums(results_ind$Null), -rowSums(r), n=ncol(object$data$Y)))
+  R222 = function(a, b) 1 - (b/a)
+  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)))
   
   to_print = results
   rownames(to_print) = to_print$models
@@ -100,6 +125,13 @@ anova.sjSDM = function(object, ...) {
   out$results = results
   out$to_print = to_print
   out$N = nrow(object$data$Y)
+  out$species = list(Residual_deviance = Residual_deviance_ind,
+                     Deviance = Deviance_ind,
+                     R2_Nagelkerke = R2_Nagelkerke_ind,
+                     R2_McFadden = R2_McFadden_ind)
+  out$sites = list(R2_Nagelkerke = R2_Nagelkerke_sites,
+                   R2_McFadden = R2_McFadden_sites)
+  out$lls = list(results_ind)
   class(out) = "sjSDManova"
   return(invisible(out))
 }
@@ -129,60 +161,188 @@ print.sjSDManova = function(x, ...) {
 #' 
 #' @param x anova object from \code{\link{anova.sjSDM}}
 #' @param y unused argument
-#' @param type use of deviance or of Nagelkerke or McFadden R-squared
+#' @param type deviance, Nagelkerke or McFadden R-squared
+#' @param internal logical, plot internal or total structure
 #' @param cols colors for the groups
 #' @param alpha alpha for colors
+#' @param env_deviance environmental deviance
 #' @param ... Additional arguments to pass to \code{plot()}
 #' 
-#' @return The visualized matrix is silently returned
+#' 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 
+#' 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
-plot.sjSDManova = function(x, y, type = c("Deviance", "Nagelkerke", "McFadden"), cols = c("#7FC97F","#BEAED4","#FDC086"),alpha=0.15, ...) {
+plot.sjSDManova = function(x, 
+                           y, 
+                           type = c("Deviance", "Nagelkerke", "McFadden"), 
+                           internal = FALSE,
+                           cols = c("#7FC97F","#BEAED4","#FDC086"),
+                           alpha=0.15, 
+                           env_deviance = NULL,
+                           ...) {
   lineSeq = 0.3
   nseg = 100
   dr = 1.0
   type = match.arg(type)
+  out = list()
   
   oldpar = par(no.readonly = TRUE)
   on.exit(par(oldpar))
   
-  values = x$results
-  select_rows = 
-    if(x$spatial) { 
-      sapply(c("A", "B", "AB","S", "AS", "BS", "ABS"), function(i) which(values$models == i, arr.ind = TRUE))
-    } else {
-      sapply(c("A", "B", "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)
+  if(!x$spatial && internal) {
+    internal=FALSE
+    warning("'internal=TRUE' currently only supported for spatial models.")
+  }
   
-  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(!internal) {
+    
+    values = x$results
+    select_rows = 
+      if(x$spatial) { 
+        sapply(c("A", "B", "AB","S", "AS", "BS", "ABS"), function(i) which(values$models == i, arr.ind = TRUE))
+      } else {
+        sapply(c("A", "B", "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(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
+    if(type == "Deviance") {type = "R2_McFadden"
+    } else {
+      if(type == "McFadden") type = "R2_McFadden"
+      else type = "R2_Nagelkerke"
+    }
+    internals = list()
+    
+
+    df = data.frame(
+        env = ifelse(x$sites[[type]]$A<0, 0, x$sites[[type]]$A),
+        spa = ifelse(x$sites[[type]]$S<0, 0, x$sites[[type]]$S),
+        codist = ifelse(x$sites[[type]]$B<0, 0, x$sites[[type]]$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]]$A<0, 0, x$species[[type]]$A),
+        spa = ifelse(x$species[[type]]$S<0, 0, x$species[[type]]$S),
+        codist = ifelse(x$species[[type]]$B<0, 0, x$species[[type]]$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"
+      
+    
+    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 = "Co-Distribution",
+                          z = "S", 
+                          zarrow = "Spatial Autocorrelation") +
+            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
+      }
+    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
   }
-  return(invisible(values))
+  return(invisible(out))
 }
 
diff --git a/sjSDM/R/importance.R b/sjSDM/R/importance.R
index 155a735..0c93183 100644
--- a/sjSDM/R/importance.R
+++ b/sjSDM/R/importance.R
@@ -9,6 +9,7 @@
 #' @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).
 #' 
 #' @return
 #' 
diff --git a/sjSDM/R/methods.R b/sjSDM/R/methods.R
index 397a0d5..27e0b1d 100644
--- a/sjSDM/R/methods.R
+++ b/sjSDM/R/methods.R
@@ -50,7 +50,7 @@ getWeights.sjSDM= function(object) {
 #' 
 #' set layer weights and sigma in \code{\link{sjSDM}} with \code{\link{DNN}} object
 #' @param object object of class  \code{\link{sjSDM}} with \code{\link{DNN}} object
-#' @param weights list of layer weights and sigma, see \code{\link{getWeights}}
+#' @param weights list of layer weights:  \code{list(env=list(matrix(...)), spatial=list(matrix(...)), sigma=matrix(...))}, see \code{\link{getWeights}}
 #' 
 #' @return No return value, weights are changed in place. 
 #' 
@@ -63,6 +63,32 @@ setWeights = function(object, weights) UseMethod("setWeights")
 setWeights.sjSDM= function(object, weights = NULL) {
   if(is.null(weights)) weights = list(env = object$weights, spatial = object$spatial_weights)
   
-  object$model$set_env_weights(w = weights[[1]])
-  object$model$set_spatial_weights(w = weights[[2]])
+  e = weights[[1]]
+  if(!is.null(e)){
+    if(!inherits(e, "list")) e = list(e)
+    object$model$set_env_weights(w = e)
+  }
+  if(inherits(object, "spatial")) {
+    if(length(weights) > 1){
+      s = weights[[2]]
+      if(!is.null(s)){
+        if(!inherits(s, "list")) s = list(s)
+        object$model$set_spatial_weights(w = s)
+      }
+    }
+    if(length(weights) > 2) {
+      sig = weights[[3]]
+      if(inherits(sig, "list")) sig = unlist(sig)
+      object$model$set_sigma(w = sig)
+    }
+    return(invisible(NULL))
+  }
+  if(length(weights) > 1) {
+    sig = weights[[3]]
+    if(!is.null(sig)) {
+      if(inherits(sig, "list")) sig = unlist(sig)
+      object$model$set_sigma(w = sig)
+    }
+  }
+
 }
\ No newline at end of file
diff --git a/sjSDM/R/plot.R b/sjSDM/R/plot.R
index da4a60f..9fe1417 100644
--- a/sjSDM/R/plot.R
+++ b/sjSDM/R/plot.R
@@ -15,7 +15,8 @@
 #' @export
 #' 
 plot.sjSDM = function(x, ...) {
-  plotsjSDMcoef(x, ...)
+  if(inherits(x, "linear")) plotsjSDMcoef(x, ...)
+  if(inherits(x, "DNN")) plot.sjSDM.DNN(x, ...)
 }
 
 #' Internal coefficients plot
@@ -96,3 +97,19 @@ plotsjSDMcoef = function(object,wrap_col=NULL,group=NULL,col=NULL,slist=NULL) {
       geom_errorbar(aes(ymax = Estimate + Std.Err, ymin = Estimate - Std.Err), width = 0.3)+ scale_y_continuous(limits = c(miny-0.3,maxy+0.1))
   })
 }
+
+
+#' Training history
+#' 
+#' Plot training loss history
+#' 
+#' @param x a model fitted by \code{\link{sjSDM}} with \code{\link{DNN}} object
+#' @param ... passed to plot
+#' @example /inst/examples/plot.sjSDM-example.R
+#' 
+#' @return No return value, called for side effects.
+#' 
+#' @import graphics
+plot.sjSDM.DNN = function(x, ...) {
+  graphics::plot(x$history, ylab = "loss", xlab = "epoch", ...)
+}
diff --git a/sjSDM/R/simulation.R b/sjSDM/R/simulation.R
index 8ced5e8..125df9e 100644
--- a/sjSDM/R/simulation.R
+++ b/sjSDM/R/simulation.R
@@ -26,7 +26,7 @@
 #' \item{link}{Which link}
 #' \item{response_type}{Which response type}
 #' \item{response}{Species occurrence matrix}
-#' \item{correlation}{Species covariance matirx}
+#' \item{correlation}{Species covariance matrix}
 #' \item{species_weights}{Species-environment coefficients}
 #' \item{env_weights}{Environmental covariates}
 #' \item{corr_acc}{Method to calculate sign accurracy}
diff --git a/sjSDM/R/sjSDM.R b/sjSDM/R/sjSDM.R
index 06da425..d082242 100644
--- a/sjSDM/R/sjSDM.R
+++ b/sjSDM/R/sjSDM.R
@@ -201,6 +201,7 @@ sjSDM = function(Y = NULL,
   
   output = as.integer(ncol(Y))
   input = as.integer(ncol(env$X))
+  intercept = "(Intercept)" %in% colnames(env$X)
   
   out$get_model = function(){
     model = pkg.env$fa$Model_sjSDM( device = device, dtype = dtype)
@@ -218,7 +219,7 @@ sjSDM = function(Y = NULL,
     if(!is.null(env$dropout)) dropout_env = env$dropout
     else dropout_env = -99
     
-    model$add_env(input, output, hidden = hidden, activation = activation,bias = bias, l1 = env$l1, l2=env$l2, dropout=dropout_env)
+    model$add_env(input, output, hidden = hidden, activation = activation,bias = bias, l1 = env$l1, l2=env$l2, dropout=dropout_env, intercept=intercept)
     
     if(!is.null(spatial)) {
       
@@ -376,7 +377,7 @@ predict.sjSDM = function(object, newdata = NULL, SP = NULL, type = c("link", "ra
   } else {
     
     if(is.null(newdata)) {
-      return(force_r(object$model$predict(newdata = object$data$X, link=link)))
+      return(force_r(object$model$predict(newdata = object$data$X, link=link, ...)))
     } else {
       if(is.data.frame(newdata)) {
         newdata = stats::model.matrix(object$formula, newdata)
@@ -531,14 +532,15 @@ summary.sjSDM = function(object, ...) {
 
 #' Generates simulations from sjSDM model
 #'
-#' Simulate nsim responses from the fitted model
+#' Simulate nsim responses from the fitted model following a multivariate probit model. 
+#' So currently only supported for \code{family = stats::binomial("probit")}
 #'
 #' @param object a model fitted by \code{\link{sjSDM}}
 #' @param nsim number of simulations
 #' @param seed seed for random numer generator
 #' @param ... optional arguments for compatibility with the generic function, no functionality implemented
 #' 
-#' @return Array of simulated species occurrences.
+#' @return Array of simulated species occurrences of dimension order [nsim, sites, species]
 #' 
 #' @importFrom stats simulate
 #' @export
@@ -549,24 +551,29 @@ simulate.sjSDM = function(object, nsim = 1, seed = NULL, ...) {
     pkg.env$torch$cuda$manual_seed(seed)
     pkg.env$torch$manual_seed(seed)
   }
-  preds = abind::abind(lapply(1:nsim, function(i) predict.sjSDM(object)), along = 0L)
-  simulation = apply(preds, 2:3, function(p) stats::rbinom(nsim, 1L,p))
-  return(simulation)
+  sims = force_r(predict(object, simulate=TRUE, sampling = as.integer(nsim)))
+  return(apply(sims, 1:3, function(i) ifelse(i > 0, 1, 0)))
 }
 
 
 
-#' Extract Log-Likelihood from a fitted sjSDM model
+#' Extract negative-log-Likelihood from a fitted sjSDM model
 #'
 #' @param object a model fitted by \code{\link{sjSDM}}
-#' @param ... optional arguments for compatibility with the generic function, no functionality implemented
+#' @param individual returns internal ll structure, mostly for internal useage
+#' @param ... optional arguments passed to internal logLik function (only used if \code{individual=TRUE})
 #' 
-#' @return Numeric value
+#' @return Numeric value or numeric matrix if individual is true.
 #' 
 #' @importFrom stats simulate
 #' @export
-logLik.sjSDM <- function(object, ...){
-  return(object$logLik[[1]])
+logLik.sjSDM <- function(object, individual=FALSE,...){
+  if(!individual) return(object$logLik[[1]])
+  else {
+    object = checkModel(object)
+    if(!inherits(object, "spatial")) return(force_r(object$model$logLik(object$data$X, object$data$Y, individual = TRUE)))
+    else return(force_r(object$model$logLik(object$data$X, object$data$Y, object$spatial$X, individual = TRUE)))
+  }
 }
 
 
@@ -587,7 +594,7 @@ update.sjSDM = function(object, env_formula = NULL, spatial_formula = NULL, biot
   mf = match.call()
   if(!is.null(env_formula)){
     m = match("env_formula", names(mf))
-    if(class(mf[3]$env_formula) == "name") mf[3]$env_formula = eval(mf[3]$env_formula, envir = parent.env(environment()))
+    if(inherits(mf[m]$env_formula,"name")) mf[m]$env_formula = eval(mf[m]$env_formula, envir = parent.env(environment()))
     env_formula = stats::as.formula(mf[m]$env_formula)
   } else {
     env_formula = object$settings$env$formula
@@ -595,7 +602,7 @@ update.sjSDM = function(object, env_formula = NULL, spatial_formula = NULL, biot
   
   if(!is.null(spatial_formula)){
     m = match("spatial_formula", names(mf))
-    if(class(mf[4]$spatial_formula) == "name") mf[4]$spatial_formula = eval(mf[4]$spatial_formula, envir = parent.env(environment()))
+    if(inherits(mf[m]$spatial_formula, "name")) mf[m]$spatial_formula = eval(mf[m]$spatial_formula, envir = parent.env(environment()))
     spatial_formula = stats::as.formula(mf[m]$spatial_formula)
   } else {
     spatial_formula = object$settings$spatial$formula
diff --git a/sjSDM/R/sjSDM_configs.R b/sjSDM/R/sjSDM_configs.R
index 5753595..a89296d 100644
--- a/sjSDM/R/sjSDM_configs.R
+++ b/sjSDM/R/sjSDM_configs.R
@@ -32,7 +32,7 @@ linear = function(data = NULL, formula = NULL, lambda = 0.0, alpha = 0.5) {
     if(!is.null(formula)){
       mf = match.call()
       m = match("formula", names(mf))
-      if(class(mf[3]$formula) == "name") mf[3]$formula = eval(mf[3]$formula, envir = parent.env(environment()))
+      if(inherits(mf[3]$formula, "name")) mf[3]$formula = eval(mf[3]$formula, envir = parent.env(environment()))
       formula = stats::as.formula(mf[m]$formula)
       X = stats::model.matrix(formula, data)
     } else {
@@ -45,7 +45,7 @@ linear = function(data = NULL, formula = NULL, lambda = 0.0, alpha = 0.5) {
     if(!is.null(formula)) {
       mf = match.call()
       m = match("formula", names(mf))
-      if(class(mf[3]$formula) == "name") mf[3]$formula = eval(mf[3]$formula, envir = parent.env(environment()))
+      if(inherits(mf[3]$formula, "name")) mf[3]$formula = eval(mf[3]$formula, envir = parent.env(environment()))
       formula = stats::as.formula(mf[m]$formula)
       data = data.frame(data)
       X = stats::model.matrix(formula, data)
@@ -91,7 +91,7 @@ print.linear = function(x, ...) {
 #' @param formula formula object for predictors
 #' @param hidden hidden units in layers, length of hidden corresponds to number of layers
 #' @param activation activation functions, can be of length one, or a vector of activation functions for each layer. Currently supported: tanh, relu, leakyrelu, selu, or sigmoid
-#' @param bias whether use biases in the layers, can be of length one, or a vector (number of hidden layers + 1 (last layer)) of logicals for each layer.
+#' @param bias whether use biases in the layers, can be of length one, or a vector (number of hidden layers including (last layer) but not first layer (intercept in first layer is specified by formula)) of logicals for each layer.
 #' @param lambda lambda penalty, strength of regularization: \eqn{\lambda * (lasso + ridge)}
 #' @param alpha weighting between lasso and ridge: \eqn{(1 - \alpha) * |weights| + \alpha ||weights||^2}
 #' @param dropout probability of dropout rate 
@@ -105,7 +105,7 @@ print.linear = function(x, ...) {
 #' \item{l1_coef}{L1 regularization strength, can be -99 if \code{lambda = 0.0}}
 #' \item{l2_coef}{L2 regularization strength, can be -99 if \code{lambda = 0.0}}
 #' \item{hidden}{Integer vector of hidden neurons in the deep neural network. Length of vector corresponds to the number of hidden layers.}
-#' \item{activation}{Charactervector of activation functions.}
+#' \item{activation}{Character vector of activation functions.}
 #' \item{bias}{Logical vector whether to use bias or not in each hidden layer.}
 #' 
 #' Implemented S3 methods include \code{\link{print.DNN}}
@@ -130,7 +130,7 @@ DNN = function(data = NULL, formula = NULL, hidden = c(10L, 10L, 10L), activatio
     if(!is.null(formula)){
       mf = match.call()
       m = match("formula", names(mf))
-      if(class(mf[3]$formula) == "name") mf[3]$formula = eval(mf[3]$formula, envir = parent.env(environment()))
+      if(inherits(mf[3]$formula, "name")) mf[3]$formula = eval(mf[3]$formula, envir = parent.env(environment()))
       formula = stats::as.formula(mf[m]$formula)
       X = stats::model.matrix(formula, data)
     } else {
@@ -143,7 +143,7 @@ DNN = function(data = NULL, formula = NULL, hidden = c(10L, 10L, 10L), activatio
     if(!is.null(formula)) {
       mf = match.call()
       m = match("formula", names(mf))
-      if(class(mf[3]$formula) == "name") mf[3]$formula = eval(mf[3]$formula, envir = parent.env(environment()))
+      if(inherits(mf[3]$formula, "name")) mf[3]$formula = eval(mf[3]$formula, envir = parent.env(environment()))
       formula = stats::as.formula(mf[m]$formula)
       data = data.frame(data)
       X = stats::model.matrix(formula, data)
diff --git a/sjSDM/R/sjSDM_cv.R b/sjSDM/R/sjSDM_cv.R
index afc7091..5f4c0c2 100644
--- a/sjSDM/R/sjSDM_cv.R
+++ b/sjSDM/R/sjSDM_cv.R
@@ -281,7 +281,7 @@ sjSDM_cv = function(Y,
              summary = summary_results, 
              settings = list(tune_samples = tune_samples, CV = CV, tune = tune),
              data = list(Y = Y, env = env, biotic = biotic, spatial = spatial),
-             config = c(list(device=device, sampling=sampling), as.list(unlist(ellip))),
+             config = c(list(device=device, sampling=sampling), ellip),
              spatial = !is.null(spatial))
   class(out) = c("sjSDM_cv")
   return(out)
diff --git a/sjSDM/README.Rmd b/sjSDM/README.Rmd
index 665c834..12088e8 100644
--- a/sjSDM/README.Rmd
+++ b/sjSDM/README.Rmd
@@ -5,7 +5,7 @@ output: github_document
 
 
 ```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = TRUE, collapse = TRUE)
+knitr::opts_chunk$set(echo = TRUE, collapse = TRUE, fig.path = "man/figures/README-")
 ```
 
 
@@ -77,6 +77,14 @@ model2 <- update(model, env_formula = ~I(X1^2))
 summary(model2)
 ```
 
+Estimate and visualize meta-community structure:
+      
+```{r}
+anv <- anova(model)
+plot(anv, internal=TRUE)
+```
+
+
 Change linear part of model to a deep neural network:
 
 ```{r}
diff --git a/sjSDM/README.md b/sjSDM/README.md
index 6c5fc9c..731c4f7 100644
--- a/sjSDM/README.md
+++ b/sjSDM/README.md
@@ -60,7 +60,7 @@ Let’s first simulate a community data set:
 
 ``` r
 library(sjSDM)
-## ── Attaching sjSDM ──────────────────────────────────────────────────── 1.0.0 ──
+## ── Attaching sjSDM ──────────────────────────────────────────────────── 1.0.1 ──
 ## ✓ torch <environment> 
 ## ✓ torch_optimizer  
 ## ✓ pyro  
@@ -77,53 +77,53 @@ Estimate jSDM:
 ``` r
 model <- sjSDM(Y = Occ, env = linear(data = Env, formula = ~X1), spatial = linear(data = SP, formula = ~0+X1:X2), se = TRUE, family=binomial("probit"), sampling = 100L)
 summary(model)
-## LogLik:  -575.5497 
+## LogLik:  -569.5459 
 ## Regularization loss:  0 
 ## 
 ## Species-species correlation matrix: 
 ## 
 ##  sp1  1.0000                                 
-##  sp2 -0.3430  1.0000                             
-##  sp3 -0.0570 -0.3700  1.0000                         
-##  sp4  0.0770 -0.3130  0.6730  1.0000                     
-##  sp5  0.4870 -0.3630 -0.0450 -0.1260  1.0000                 
-##  sp6 -0.0970  0.3320  0.1880  0.1880 -0.0880  1.0000             
-##  sp7  0.4360 -0.2090  0.1800  0.1530  0.5090  0.2490  1.0000         
-##  sp8  0.2350  0.0560 -0.3290 -0.3380  0.3190 -0.0090  0.2010  1.0000     
-##  sp9  0.0400 -0.0330  0.0160  0.2050 -0.3410 -0.1920 -0.2980 -0.2220  1.0000 
-##  sp10     0.3090  0.3230 -0.5050 -0.3800  0.1930  0.1360  0.1890  0.3570 -0.0930  1.0000
+##  sp2 -0.3380  1.0000                             
+##  sp3 -0.0310 -0.3640  1.0000                         
+##  sp4  0.0890 -0.2770  0.6670  1.0000                     
+##  sp5  0.4590 -0.3710 -0.0140 -0.0960  1.0000                 
+##  sp6 -0.1190  0.3690  0.1730  0.1980 -0.0770  1.0000             
+##  sp7  0.4290 -0.1860  0.1700  0.1470  0.5150  0.2420  1.0000         
+##  sp8  0.2490  0.0510 -0.3270 -0.3430  0.3090 -0.0230  0.2280  1.0000     
+##  sp9  0.0450 -0.0290  0.0060  0.1660 -0.3590 -0.2000 -0.3070 -0.1840  1.0000 
+##  sp10     0.2910  0.3150 -0.5000 -0.3540  0.1800  0.1080  0.1860  0.3410 -0.0710  1.0000
 ## 
 ## 
 ## 
 ## Spatial: 
-##              sp1        sp2       sp3       sp4       sp5       sp6       sp7
-## X1:X2 0.08259914 -0.5019723 0.4352706 -0.169348 0.4040644 0.2000494 0.4183032
+##             sp1        sp2       sp3        sp4       sp5       sp6       sp7
+## X1:X2 0.1362036 -0.4616366 0.3906074 -0.1906729 0.4108464 0.1651098 0.4153814
 ##             sp8       sp9        sp10
-## X1:X2 0.3042413 0.2288286 -0.02296865
-## 
-## 
-## 
-##                  Estimate  Std.Err Z value Pr(>|z|)    
-## sp1 (Intercept)  -0.06842  0.20499   -0.33  0.73857    
-## sp1 X1            0.72483  0.39751    1.82  0.06824 .  
-## sp2 (Intercept)  -0.03794  0.20898   -0.18  0.85593    
-## sp2 X1            0.92104  0.40582    2.27  0.02323 *  
-## sp3 (Intercept)  -0.26596  0.21984   -1.21  0.22637    
-## sp3 X1            0.95088  0.41292    2.30  0.02129 *  
-## sp4 (Intercept)  -0.04988  0.20611   -0.24  0.80877    
-## sp4 X1           -1.15879  0.38965   -2.97  0.00294 ** 
-## sp5 (Intercept)  -0.15585  0.19463   -0.80  0.42326    
-## sp5 X1            0.50131  0.37329    1.34  0.17929    
-## sp6 (Intercept)   0.18499  0.19405    0.95  0.34043    
-## sp6 X1            1.55482  0.39826    3.90  9.5e-05 ***
-## sp7 (Intercept)   0.00426  0.17742    0.02  0.98084    
-## sp7 X1           -0.34412  0.34556   -1.00  0.31933    
-## sp8 (Intercept)   0.15873  0.13976    1.14  0.25607    
-## sp8 X1            0.14996  0.26153    0.57  0.56637    
-## sp9 (Intercept)   0.01891  0.17133    0.11  0.91211    
-## sp9 X1            1.09219  0.32656    3.34  0.00082 ***
-## sp10 (Intercept) -0.07258  0.17739   -0.41  0.68243    
-## sp10 X1          -0.49947  0.32436   -1.54  0.12359    
+## X1:X2 0.2773075 0.2295802 -0.02492988
+## 
+## 
+## 
+##                  Estimate Std.Err Z value Pr(>|z|)    
+## sp1 (Intercept)   -0.0405  0.1998   -0.20  0.83921    
+## sp1 X1             0.6836  0.3973    1.72  0.08527 .  
+## sp2 (Intercept)   -0.0101  0.2120   -0.05  0.96207    
+## sp2 X1             0.9244  0.4110    2.25  0.02450 *  
+## sp3 (Intercept)   -0.2737  0.2150   -1.27  0.20302    
+## sp3 X1             0.9865  0.4143    2.38  0.01726 *  
+## sp4 (Intercept)   -0.0560  0.2111   -0.27  0.79074    
+## sp4 X1            -1.1353  0.4182   -2.71  0.00663 ** 
+## sp5 (Intercept)   -0.1702  0.1978   -0.86  0.38964    
+## sp5 X1             0.4623  0.3792    1.22  0.22281    
+## sp6 (Intercept)    0.2006  0.1902    1.05  0.29152    
+## sp6 X1             1.5796  0.3863    4.09  4.3e-05 ***
+## sp7 (Intercept)    0.0180  0.1775    0.10  0.91941    
+## sp7 X1            -0.3499  0.3488   -1.00  0.31576    
+## sp8 (Intercept)    0.1588  0.1412    1.12  0.26069    
+## sp8 X1             0.1323  0.2627    0.50  0.61452    
+## sp9 (Intercept)    0.0264  0.1703    0.15  0.87699    
+## sp9 X1             1.1001  0.3269    3.36  0.00077 ***
+## sp10 (Intercept)  -0.0525  0.1735   -0.30  0.76199    
+## sp10 X1           -0.5166  0.3234   -1.60  0.11022    
 ## ---
 ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 ```
@@ -133,83 +133,97 @@ Update model (change main effects to quadratic effects):
 ``` r
 model2 <- update(model, env_formula = ~I(X1^2))
 summary(model2)
-## LogLik:  -608.1784 
+## LogLik:  -600.4297 
 ## Regularization loss:  0 
 ## 
 ## Species-species correlation matrix: 
 ## 
 ##  sp1  1.0000                                 
-##  sp2 -0.3030  1.0000                             
-##  sp3 -0.0130 -0.2560  1.0000                         
-##  sp4  0.0600 -0.3530  0.5650  1.0000                     
-##  sp5  0.4550 -0.3540 -0.0540 -0.0690  1.0000                 
-##  sp6 -0.0420  0.4110  0.2420  0.0420 -0.0790  1.0000             
-##  sp7  0.4220 -0.2130  0.1750  0.0780  0.5140  0.2090  1.0000         
-##  sp8  0.2440  0.0710 -0.3000 -0.3740  0.3210  0.0220  0.2360  1.0000     
-##  sp9  0.0920  0.0950  0.1040  0.1220 -0.3260  0.0200 -0.2400 -0.1540  1.0000 
-##  sp10     0.2760  0.2130 -0.5060 -0.3470  0.2340  0.0230  0.1450  0.3220 -0.1210  1.0000
+##  sp2 -0.2870  1.0000                             
+##  sp3 -0.0050 -0.2600  1.0000                         
+##  sp4  0.0600 -0.3630  0.5690  1.0000                     
+##  sp5  0.4710 -0.3350 -0.0270 -0.0520  1.0000                 
+##  sp6 -0.0310  0.4350  0.2190  0.0140 -0.0410  1.0000             
+##  sp7  0.4280 -0.2050  0.1700  0.0730  0.5290  0.2020  1.0000         
+##  sp8  0.2530  0.0580 -0.2890 -0.3710  0.3270  0.0220  0.2400  1.0000     
+##  sp9  0.0550  0.1120  0.0980  0.1210 -0.3220  0.0230 -0.2450 -0.1590  1.0000 
+##  sp10     0.2720  0.2180 -0.4960 -0.3520  0.2270  0.0360  0.1420  0.3030 -0.1140  1.0000
 ## 
 ## 
 ## 
 ## Spatial: 
-##             sp1        sp2       sp3        sp4       sp5       sp6      sp7
-## X1:X2 0.1553214 -0.4860287 0.4846192 -0.2180294 0.4662001 0.1784182 0.378258
-##             sp8       sp9        sp10
-## X1:X2 0.2950488 0.2057981 -0.03424472
+##             sp1       sp2       sp3        sp4       sp5       sp6       sp7
+## X1:X2 0.1946324 -0.462581 0.4739597 -0.1982521 0.4937007 0.1960442 0.4032983
+##             sp8       sp9         sp10
+## X1:X2 0.2994181 0.2302684 -0.006441596
 ## 
 ## 
 ## 
 ##                  Estimate Std.Err Z value Pr(>|z|)
-## sp1 (Intercept)   -0.1303  0.2849   -0.46     0.65
-## sp1 I(X1^2)        0.4434  0.6996    0.63     0.53
-## sp2 (Intercept)   -0.2555  0.2864   -0.89     0.37
-## sp2 I(X1^2)        0.7050  0.7215    0.98     0.33
-## sp3 (Intercept)   -0.2640  0.2952   -0.89     0.37
-## sp3 I(X1^2)        0.1099  0.6970    0.16     0.87
-## sp4 (Intercept)   -0.2272  0.2737   -0.83     0.41
-## sp4 I(X1^2)        0.7222  0.6841    1.06     0.29
-## sp5 (Intercept)   -0.1383  0.2487   -0.56     0.58
-## sp5 I(X1^2)        0.2325  0.6248    0.37     0.71
-## sp6 (Intercept)    0.1669  0.2690    0.62     0.53
-## sp6 I(X1^2)       -0.0286  0.6517   -0.04     0.97
-## sp7 (Intercept)    0.0521  0.2466    0.21     0.83
-## sp7 I(X1^2)       -0.0193  0.6236   -0.03     0.98
-## sp8 (Intercept)   -0.0328  0.2086   -0.16     0.88
-## sp8 I(X1^2)        0.7716  0.5249    1.47     0.14
-## sp9 (Intercept)    0.2076  0.2478    0.84     0.40
-## sp9 I(X1^2)       -0.6296  0.5992   -1.05     0.29
-## sp10 (Intercept)  -0.1308  0.2505   -0.52     0.60
-## sp10 I(X1^2)       0.2070  0.5999    0.35     0.73
+## sp1 (Intercept)   -0.1307  0.2899   -0.45     0.65
+## sp1 I(X1^2)        0.4196  0.7427    0.56     0.57
+## sp2 (Intercept)   -0.2218  0.3048   -0.73     0.47
+## sp2 I(X1^2)        0.6981  0.7684    0.91     0.36
+## sp3 (Intercept)   -0.3121  0.2963   -1.05     0.29
+## sp3 I(X1^2)        0.1102  0.7289    0.15     0.88
+## sp4 (Intercept)   -0.2547  0.2787   -0.91     0.36
+## sp4 I(X1^2)        0.7601  0.6738    1.13     0.26
+## sp5 (Intercept)   -0.1277  0.2489   -0.51     0.61
+## sp5 I(X1^2)        0.2279  0.6497    0.35     0.73
+## sp6 (Intercept)    0.1619  0.2606    0.62     0.53
+## sp6 I(X1^2)       -0.0363  0.6880   -0.05     0.96
+## sp7 (Intercept)    0.0430  0.2347    0.18     0.85
+## sp7 I(X1^2)       -0.0318  0.6161   -0.05     0.96
+## sp8 (Intercept)   -0.0196  0.2080   -0.09     0.93
+## sp8 I(X1^2)        0.7520  0.5211    1.44     0.15
+## sp9 (Intercept)    0.2057  0.2406    0.85     0.39
+## sp9 I(X1^2)       -0.6158  0.5886   -1.05     0.30
+## sp10 (Intercept)  -0.1004  0.2440   -0.41     0.68
+## sp10 I(X1^2)       0.2091  0.5959    0.35     0.73
+```
+
+Estimate and visualize meta-community structure:
+
+``` r
+anv <- anova(model)
+plot(anv, 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-7-1.png)<!-- -->
+
 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)
-## LogLik:  -564.1605 
+## LogLik:  -507.1421 
 ## Regularization loss:  0 
 ## 
 ## Species-species correlation matrix: 
 ## 
 ##  sp1  1.0000                                 
-##  sp2 -0.4150  1.0000                             
-##  sp3 -0.0480 -0.3220  1.0000                         
-##  sp4  0.0440 -0.3230  0.5590  1.0000                     
-##  sp5  0.5510 -0.2900 -0.0280 -0.0880  1.0000                 
-##  sp6 -0.1760  0.3670  0.2260  0.0560  0.0130  1.0000             
-##  sp7  0.5050 -0.1590  0.2130  0.0860  0.5020  0.3070  1.0000         
-##  sp8  0.2690  0.1100 -0.2800 -0.3930  0.2870  0.0620  0.2770  1.0000     
-##  sp9 -0.0630 -0.1070  0.0550  0.1380 -0.2600 -0.2420 -0.2260 -0.1470  1.0000 
-##  sp10     0.2720  0.2690 -0.5240 -0.3000  0.2340  0.0440  0.1590  0.2880 -0.1800  1.0000
+##  sp2 -0.4000  1.0000                             
+##  sp3 -0.1060 -0.3570  1.0000                         
+##  sp4 -0.0650 -0.3470  0.7530  1.0000                     
+##  sp5  0.5770 -0.3450 -0.1110 -0.0710  1.0000                 
+##  sp6 -0.2720  0.3630  0.2170  0.2060 -0.0600  1.0000             
+##  sp7  0.4220 -0.1520  0.1560  0.1970  0.4780  0.2430  1.0000         
+##  sp8  0.2150  0.1140 -0.4060 -0.4120  0.2190 -0.0540  0.0890  1.0000     
+##  sp9 -0.0390 -0.0630  0.0980  0.0970 -0.3440 -0.2930 -0.1830 -0.1220  1.0000 
+##  sp10     0.1270  0.3730 -0.6080 -0.5760  0.2280  0.0620  0.0820  0.3760 -0.2510  1.0000
 ## 
 ## 
 ## 
 ## Spatial: 
-##              sp1        sp2       sp3        sp4       sp5       sp6       sp7
-## X1:X2 0.09318528 -0.4935189 0.4328462 -0.1994821 0.4871828 0.1590104 0.3766337
-##             sp8       sp9         sp10
-## X1:X2 0.2689047 0.2027397 -0.002816148
+##             sp1        sp2       sp3        sp4       sp5       sp6      sp7
+## X1:X2 0.1791221 -0.5187001 0.5328111 -0.1147005 0.3768144 0.2153512 0.528778
+##             sp8       sp9      sp10
+## X1:X2 0.4690367 0.2996197 0.2277844
 ## 
 ## 
 ## 
diff --git a/sjSDM/cran-comments.md b/sjSDM/cran-comments.md
index b554eec..467d4bb 100644
--- a/sjSDM/cran-comments.md
+++ b/sjSDM/cran-comments.md
@@ -1,3 +1,45 @@
+## Version 1.0.1
+
+### Submission 1, 03/11/2022
+This is a minor update, fixing a few minor bugs and undesired behaviors reported 
+by users and adding a few new functionalities for the internal meta-community 
+structure
+
+Examples cannot be executed directly after the pkg installation
+because they need additional python dependencies which first have to be 
+installed by the `sjSDM::install_sjSDM()` function. Therefore we deiced to use
+the `\dontrun{}` flag
+
+
+### Successfull R CMD checks under
+
+* Locally: MacOS Monterey 12.2.1 (R x86_64 version)
+* Github actions: 
+  - MacOS Catalina 10.15 R-release 
+  - Ubuntu 20.04 R-release, R-oldrelease, and R-development
+  - Windows-latest R-release
+* Rhub:
+  - fedora 
+  - solairs
+* Win-builder R-release, R-development, and R-oldrelease
+
+Notes (from win-builder): 
+
+Possibly misspelled words in DESCRIPTION:
+  CPUs (21:262)
+  GPUs (21:271)
+  Hartig (21:320)
+  Pichler (21:310)
+  Scalable (3:8)
+  autocorrelation (21:488)
+  jSDMs (21:79, 21:565)
+  scalable (21:16)
+
+The spelling is intended.
+
+
+
+
 ## Version 1.0.0
 
 
@@ -64,7 +106,7 @@ Done. All plot functions have been updated accordingly.
 > Please ensure that you do not modify the global environment (e.g. by
 > using <<-) in your functions. This is not allowed by the CRAN policies.
 
-Done. We moved global variables into the pkg.env.
+Done. We haved moved global variables into the pkg.env.
 
 
 ### Successfull R CMD checks under
diff --git a/sjSDM/inst/examples/sjSDM-example.R b/sjSDM/inst/examples/sjSDM-example.R
index 2493d51..d2295ab 100644
--- a/sjSDM/inst/examples/sjSDM-example.R
+++ b/sjSDM/inst/examples/sjSDM-example.R
@@ -55,6 +55,8 @@ model = sjSDM(Y = com$response, env = linear(com$env_weights),
               iter = 50L) # increase iter for your own data 
 summary(model)
 predict(model, newdata = com$env_weights, SP = XY)
+R2 = Rsquared(model)
+print(R2)
 
 ## Using spatial eigenvectors as predictors to account 
 ## for spatial autocorrelation is a common approach:
@@ -65,7 +67,9 @@ model = sjSDM(Y = com$response, env = linear(com$env_weights),
 summary(model)
 predict(model, newdata = com$env_weights, SP = SPV)
 
-
+## Visualize internal meta-community structure
+an = anova(model)
+plot(an, internal=TRUE)
 
 ## non-linear(deep neural network) model
 model = sjSDM(Y = com$response, env = linear(com$env_weights), 
@@ -104,6 +108,9 @@ result = anova(model)
 print(result)
 plot(result)
 
+## visualize meta-community structure
+plot(result, internal=TRUE)
+
 
 # Deep neural network
 ## we can fit also a deep neural network instead of a linear model:
diff --git a/sjSDM/inst/python/sjSDM_py/model_sjSDM.py b/sjSDM/inst/python/sjSDM_py/model_sjSDM.py
index f6b6cae..d98b3ff 100644
--- a/sjSDM/inst/python/sjSDM_py/model_sjSDM.py
+++ b/sjSDM/inst/python/sjSDM_py/model_sjSDM.py
@@ -33,7 +33,11 @@ def l1_loss(tensor: torch.Tensor, l1: float):
         def l2_loss(tensor: torch.Tensor, l2: float):
             return tensor.pow(2.0).sum().mul(l2)
         
-        self.l1_l2 = [l1_loss, l2_loss]
+        @torch.jit.script
+        def l1_l2_loss(tensor: torch.Tensor, l1: float, l2: float):
+            return tensor.pow(2.0).sum().mul(l2)+tensor.abs().sum().mul(l1)
+        
+        self.l1_l2 = [l1_loss, l2_loss, l1_l2_loss]
 
     def __call__(self):
         pass
@@ -78,7 +82,8 @@ def add_env(self,
                 bias: List[bool] = [False], 
                 l1: float = -99, 
                 l2: float = -99, 
-                dropout: float = -99) -> None:
+                dropout: float = -99,
+                intercept=False) -> None:
         """Add environmental model
 
         Args:
@@ -96,11 +101,44 @@ def add_env(self,
         self.params.append(self.env.parameters())
         self.input_shape = input_shape
         self.output_shape = output_shape
-        for p in self.env.parameters():
-            if l1 > 0.0: 
-                self.losses.append(lambda: self.l1_l2[0](p, l1))
-            if l2 > 0.0:
-                self.losses.append(lambda: self.l1_l2[1](p, l2))
+
+        individual_losses = []
+        for index, p in enumerate(self.env.parameters()):
+            if index == 0:
+                if (l1 > 0.0) & (l2 > 0.0): 
+                    if intercept is False:
+                        individual_losses.append(lambda p: self.l1_l2[2](p, l1, l2))
+                    else:
+                        individual_losses.append(lambda p: self.l1_l2[2](p[:,1:], l1, l2))
+                    next
+                elif (l1 <= 0.0) & (l2 > 0.0):
+                    if intercept is False:
+                        individual_losses.append(lambda p: self.l1_l2[1](p, l2))
+                    else:
+                        individual_losses.append(lambda p: self.l1_l2[1](p[:,1:], l2))
+                elif (l1 > 0.0) & (l2 <= 0.0):
+                    if intercept is False:
+                        individual_losses.append(lambda p: self.l1_l2[0](p, l1))
+                    else:
+                        individual_losses.append(lambda p: self.l1_l2[0](p[:,1:], l1))
+            else:
+                if (l1 > 0.0) & (l2 > 0.0): 
+                    individual_losses.append(lambda p: self.l1_l2[2](p, l1, l2))
+                elif (l1 <= 0.0) & (l2 > 0.0):
+                    individual_losses.append(lambda p: self.l1_l2[1](p, l2))
+                elif (l1 > 0.0) & (l2 <= 0.0):
+                    individual_losses.append(lambda p: self.l1_l2[0](p, l1))
+        if len(individual_losses) > 0:
+            def loss():
+                tmp_loss = [
+                    individual_losses[index](p).reshape([1])
+                    for index, p in enumerate(self.env.parameters())
+                    if len(p.shape) > 1
+                ]
+                return torch.cat(tmp_loss).sum()
+            self.losses.append( loss )
+                    
+
 
     def add_spatial(self, 
                     input_shape: int, 
@@ -110,7 +148,8 @@ def add_spatial(self,
                     bias: List[bool] = [False], 
                     l1: float = -99, 
                     l2: float = -99, 
-                    dropout: float = -99) -> None:
+                    dropout: float = -99,
+                    intercept=False) -> None:
         """Add spatial model
 
         Args:
@@ -125,11 +164,25 @@ def add_spatial(self,
         """                    
         self.spatial = (self._build_NN(input_shape, output_shape, hidden, bias, activation, dropout))
         self.params.append(self.spatial.parameters())
+        
+        individual_losses = []
         for p in self.spatial.parameters():
-            if l1 > 0.0: 
-                self.losses.append(lambda: self.l1_l2[0](p, l1))
-            if l2 > 0.0:
-                self.losses.append(lambda: self.l1_l2[1](p, l2))
+            if (l1 > 0.0) & (l2 > 0.0): 
+                individual_losses.append(lambda p: self.l1_l2[2](p, l1, l2))
+            elif (l1 <= 0.0) & (l2 > 0.0):
+                individual_losses.append(lambda p: self.l1_l2[1](p, l2))
+            elif (l1 > 0.0) & (l2 <= 0.0):
+                individual_losses.append(lambda p: self.l1_l2[0](p, l1))
+        if len(individual_losses) > 0:
+            def loss():
+                tmp_loss = [
+                    individual_losses[index](p).reshape([1])
+                    for index, p in enumerate(self.spatial.parameters())
+                    if len(p.shape) > 1
+                ]
+                return torch.cat(tmp_loss).sum()
+            self.losses.append( loss )
+        
                 
     def _build_NN(self, 
                   input_shape: int, 
@@ -156,7 +209,9 @@ def _build_NN(self,
             activation = [activation[0] for _ in range(len(hidden))]
 
         if len(bias) == 1:
-            bias = [bias[0] for _ in range(len(hidden)+1)]  
+            bias = [bias[0] for _ in range(len(hidden))]
+            
+        bias.insert(0, False)
 
         if len(hidden) > 0:
             for i in range(len(hidden)):
@@ -234,7 +289,7 @@ def _get_DataLoader(self,
             else: 
                 data = torch.utils.data.TensorDataset(torch.tensor(X, dtype=self.dtype, device=torch.device('cpu')))            
 
-        DataLoader = torch.utils.data.DataLoader(data, batch_size=batch_size, shuffle=shuffle, num_workers=int(parallel), pin_memory=pin_memory, drop_last=drop_last)
+        DataLoader = torch.utils.data.DataLoader(data, batch_size=int(batch_size), shuffle=shuffle, num_workers=int(parallel), pin_memory=pin_memory, drop_last=drop_last)
         torch.cuda.empty_cache()
                 
         return DataLoader
@@ -464,7 +519,7 @@ def logLik(self,
             Tuple[float, float]: Tuple of log-Likelihood and regularization loss
         """               
         dataLoader = self._get_DataLoader(X = X, Y = Y, SP=SP, batch_size = batch_size, shuffle = False, parallel = parallel, drop_last = False)
-        loss_function = self._build_loss_function(train=train)
+        loss_function = self._build_loss_function(train=train, individual=individual)
         torch.cuda.empty_cache()
         any_losses = len(self.losses) > 0
 
@@ -512,7 +567,8 @@ def predict(self,
                 parallel: int = 0, 
                 sampling: int = 100, 
                 link: bool = True,
-                dropout: bool = False) -> np.ndarray:
+                dropout: bool = False,
+                simulate: bool = False) -> np.ndarray:
         """Predict with sjSDM
 
         Args:
@@ -524,13 +580,14 @@ def predict(self,
             sampling (int, optional): Number of MC-samples for each species. Defaults to 100.
             link (bool, optional): Linear or response scale. Defaults to True.
             dropout (bool, optional): Use dropout during predictions or not. Defaults to False.
+            simulate (bool, optional): Return simulated values on linear scale. Defaults to False.
 
         Returns:
             np.ndarray: Predictions
         """                
 
         dataLoader = self._get_DataLoader(X = newdata, Y = None, SP=SP, batch_size = batch_size, shuffle = False, parallel = parallel, drop_last = False)
-        loss_function = self._build_loss_function(train = False, raw=not link)
+        loss_function = self._build_loss_function(train = False, raw=not link, simulate=simulate)
 
         pred = []
         if self.device.type == 'cuda': # type: ignore
@@ -554,7 +611,7 @@ def predict(self,
                 sp = sp.to(self.device, non_blocking=True)
                 mu = self.env(x) + self.spatial(sp) # type: ignore
                 # loss = self._loss_function(mu, y, self.sigma, batch_size, sampling, df, alpha, device)
-                loss = loss_function(mu, self.sigma, x.shape[0], sampling, self.df, self.alpha, device)
+                loss = loss_function(mu, self.sigma, x.shape[0], sampling, self.df, self.alpha, device, self.dtype)
                 pred.append(loss)
             self.spatial.eval()
         else:
@@ -562,11 +619,14 @@ def predict(self,
                 x = x[0].to(self.device, non_blocking=True)
                 mu = self.env(x) # type: ignore
                 # loss = self._loss_function(mu, y, self.sigma, batch_size, sampling, df, alpha, device)
-                loss = loss_function(mu, self.sigma, x.shape[0], sampling, self.df, self.alpha, device)
+                loss = loss_function(mu, self.sigma, x.shape[0], sampling, self.df, self.alpha, device, self.dtype)
                 pred.append(loss)
 
         self.env.eval()
-        return torch.cat(pred, dim = 0).data.cpu().numpy()
+        cat_dim = 0
+        if simulate:
+            cat_dim = 1
+        return torch.cat(pred, dim = cat_dim).data.cpu().numpy()
 
     def se(self, 
            X: np.ndarray, 
@@ -744,16 +804,26 @@ def l2_ll(sigma: torch.Tensor, l2: float, diag: int):
             if l2 > 0.0:
                 self.losses.append( lambda: self.l1_l2[1](self.sigma, l2) )
 
-    def _build_loss_function(self, train: bool = True, raw: bool = False) -> Callable:
+    def _build_loss_function(self, train: bool = True, raw: bool = False, individual:bool = False, simulate: bool = False) -> Callable:
         """Build loss (likelihood) function
 
         Args:
             train (bool, optional): Train or evaulation mode Defaults to True.
             raw (bool, optional): Linear or response scale. Defaults to False.
+            individual(bool, optional): Return individual logLL values. Defaults to False.
+            simulate(bool, optional): Return simulated values. Defaults to False.
 
         Returns:
             Callable: loss function
-        """        
+        """
+        
+        if simulate:
+            #tmp(mu: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str):
+            @torch.jit.script
+            def tmp(mu: torch.Tensor,  sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):
+                noise = torch.randn(size = [sampling, batch_size, df], device=torch.device(device), dtype=dtype)
+                return torch.einsum("ijk, lk -> ijl", [noise, sigma]).add(mu)
+            return tmp
 
         if train:
             if self.link == "logit" or self.link == "probit":
@@ -788,7 +858,6 @@ def tmp(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int
             elif self.link == "normal":
                 def tmp(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):
                     return torch.distributions.MultivariateNormal(loc=mu, covariance_matrix=sigma.matmul(sigma.t()).add(torch.eye(sigma.shape[0]))).log_prob(Ys).neg()
-
         else:
             if self.link == "probit": 
                 link_func = lambda value: torch.distributions.Normal(0.0, 1.0).cdf(value)
@@ -804,13 +873,55 @@ def tmp(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int
             if raw:
                 link_func = lambda value: value
 
-            def tmp(mu: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str):
-                noise = torch.randn(size = [sampling, batch_size, df], device=torch.device(device))
+            def tmp(mu: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):
+                # noise = torch.randn(size = [sampling, batch_size, df], device=torch.device(device))
                 if self.link == "logit": 
-                    E = link_func(torch.einsum("ijk, lk -> ijl", [noise, sigma]).add(mu).mul(alpha)).mul(0.999999).add(0.0000005)
+                    E = link_func(mu.mul(alpha)).mul(0.999999).add(0.0000005)
                 else:
-                    E = link_func(torch.einsum("ijk, lk -> ijl", [noise, sigma]).add(mu)).mul(0.999999).add(0.0000005)
-                return E.mean(dim = 0)
+                    E = link_func(mu).mul(0.999999).add(0.0000005)
+                return E
+    
+    
+        if individual:
+            if self.link == "logit" or self.link == "probit":
+                @torch.jit.script
+                def tmp(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):
+                    noise = torch.randn(size = [sampling, batch_size, df], device=torch.device(device), dtype=dtype)
+                    E = torch.sigmoid(   torch.einsum("ijk, lk -> ijl", [noise, sigma]).add(mu).mul(alpha)   ).mul(0.999999).add(0.0000005)
+                    logprob = E.log().mul(Ys).add((1.0 - E).log().mul(1.0 - Ys)).neg().sum(dim = 2).neg()
+                    logprob2 = E.log().mul(Ys).add((1.0 - E).log().mul(1.0 - Ys))#.neg().sum(dim = 2).neg()
+                    Prop = logprob2.exp().mean(0).log().abs() # 
+                    Prop = Prop.multiply( (1.0/Prop.sum(dim=1)).reshape([-1, 1]).repeat_interleave(Ys.shape[1],1) )
+                    maxlogprob = logprob.max(dim = 0).values
+                    Eprob = logprob.sub(maxlogprob).exp().mean(dim = 0)
+                    return Eprob.log().neg().sub(maxlogprob).reshape([-1,1]).repeat_interleave(Ys.shape[1],1)*Prop
+            elif self.link == "linear":
+                @torch.jit.script
+                def tmp(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):
+                    noise = torch.randn(size = [sampling, batch_size, df], device=torch.device(device), dtype=dtype)
+                    E = torch.clamp(torch.einsum("ijk, lk -> ijl", [noise, sigma]).add(mu).mul(alpha), 0.0, 1.0).mul(0.999999).add(0.0000005)
+                    logprob = E.log().mul(Ys).add((1.0 - E).log().mul(1.0 - Ys)).neg().sum(dim = 2).neg()
+                    logprob2 = E.log().mul(Ys).add((1.0 - E).log().mul(1.0 - Ys))#.neg().sum(dim = 2).neg()
+                    Prop = logprob2.exp().mean(0).log().abs() # 
+                    Prop = Prop.multiply( (1.0/Prop.sum(dim=1)).reshape([-1, 1]).repeat_interleave(Ys.shape[1],1) )
+                    maxlogprob = logprob.max(dim = 0).values
+                    Eprob = logprob.sub(maxlogprob).exp().mean(dim = 0)
+                    return Eprob.log().neg().sub(maxlogprob).reshape([-1,1]).repeat_interleave(Ys.shape[1],1)*Prop
+            elif self.link == "count":
+                def tmp(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):
+                    noise = torch.randn(size = [sampling, batch_size, df], device=torch.device(device),dtype=dtype)
+                    E = torch.einsum("ijk, lk -> ijl", [noise, sigma]).add(mu).exp()
+                    logprob = torch.distributions.Poisson(rate=E).log_prob(Ys).sum(2)
+                    logprob2 = torch.distributions.Poisson(rate=E).log_prob(Ys)#.sum(2)
+                    Prop = logprob2.exp().mean(0).log().abs() # 
+                    Prop = Prop.multiply( (1.0/Prop.sum(dim=1)).reshape([-1, 1]).repeat_interleave(Ys.shape[1],1) )
+                    maxlogprob = logprob.max(dim = 0).values
+                    Eprob = logprob.sub(maxlogprob).exp().mean(dim = 0)
+                    return Eprob.log().neg().sub(maxlogprob).reshape([-1,1]).repeat_interleave(Ys.shape[1],1)*Prop
+
+            elif self.link == "normal":
+                def tmp(mu: torch.Tensor, Ys: torch.Tensor, sigma: torch.Tensor, batch_size: int, sampling: int, df: int, alpha: float, device: str, dtype: torch.dtype):
+                    return torch.distributions.MultivariateNormal(loc=mu, covariance_matrix=sigma.matmul(sigma.t()).add(torch.eye(sigma.shape[0]))).log_prob(Ys).neg()
 
         return tmp
 
diff --git a/sjSDM/man/DNN.Rd b/sjSDM/man/DNN.Rd
index 3aee925..a9602f1 100644
--- a/sjSDM/man/DNN.Rd
+++ b/sjSDM/man/DNN.Rd
@@ -24,7 +24,7 @@ DNN(
 
 \item{activation}{activation functions, can be of length one, or a vector of activation functions for each layer. Currently supported: tanh, relu, leakyrelu, selu, or sigmoid}
 
-\item{bias}{whether use biases in the layers, can be of length one, or a vector (number of hidden layers + 1 (last layer)) of logicals for each layer.}
+\item{bias}{whether use biases in the layers, can be of length one, or a vector (number of hidden layers including (last layer) but not first layer (intercept in first layer is specified by formula)) of logicals for each layer.}
 
 \item{lambda}{lambda penalty, strength of regularization: \eqn{\lambda * (lasso + ridge)}}
 
@@ -41,7 +41,7 @@ An S3 class of type 'DNN' including the following components:
 \item{l1_coef}{L1 regularization strength, can be -99 if \code{lambda = 0.0}}
 \item{l2_coef}{L2 regularization strength, can be -99 if \code{lambda = 0.0}}
 \item{hidden}{Integer vector of hidden neurons in the deep neural network. Length of vector corresponds to the number of hidden layers.}
-\item{activation}{Charactervector of activation functions.}
+\item{activation}{Character vector of activation functions.}
 \item{bias}{Logical vector whether to use bias or not in each hidden layer.}
 
 Implemented S3 methods include \code{\link{print.DNN}}
@@ -107,6 +107,8 @@ model = sjSDM(Y = com$response, env = linear(com$env_weights),
               iter = 50L) # increase iter for your own data 
 summary(model)
 predict(model, newdata = com$env_weights, SP = XY)
+R2 = Rsquared(model)
+print(R2)
 
 ## Using spatial eigenvectors as predictors to account 
 ## for spatial autocorrelation is a common approach:
@@ -117,7 +119,9 @@ model = sjSDM(Y = com$response, env = linear(com$env_weights),
 summary(model)
 predict(model, newdata = com$env_weights, SP = SPV)
 
-
+## Visualize internal meta-community structure
+an = anova(model)
+plot(an, internal=TRUE)
 
 ## non-linear(deep neural network) model
 model = sjSDM(Y = com$response, env = linear(com$env_weights), 
@@ -156,6 +160,9 @@ result = anova(model)
 print(result)
 plot(result)
 
+## visualize meta-community structure
+plot(result, internal=TRUE)
+
 
 # Deep neural network
 ## we can fit also a deep neural network instead of a linear model:
diff --git a/sjSDM/man/anova.sjSDM.Rd b/sjSDM/man/anova.sjSDM.Rd
index 8e8ea3a..7e0f040 100644
--- a/sjSDM/man/anova.sjSDM.Rd
+++ b/sjSDM/man/anova.sjSDM.Rd
@@ -17,7 +17,10 @@ An S3 class of type 'sjSDManova' including the following components:
 \item{results}{Data frame of results.}
 \item{to_print}{Data frame, summarized results for type I anova.}
 \item{N}{Number of observations (sites).}
-\item{spatial}{Logical, spatial model or not}
+\item{spatial}{Logical, spatial model or not.}
+\item{species}{individual species R2s.}
+\item{sites}{individual site R2s.}
+\item{lls}{individual site by species negative-log-likelihood values.}
 
 Implemented S3 methods are \code{\link{print.sjSDManova}} and \code{\link{plot.sjSDManova}}
 }
diff --git a/sjSDM/man/bioticStruct.Rd b/sjSDM/man/bioticStruct.Rd
index 633dddc..35ab290 100644
--- a/sjSDM/man/bioticStruct.Rd
+++ b/sjSDM/man/bioticStruct.Rd
@@ -102,6 +102,8 @@ model = sjSDM(Y = com$response, env = linear(com$env_weights),
               iter = 50L) # increase iter for your own data 
 summary(model)
 predict(model, newdata = com$env_weights, SP = XY)
+R2 = Rsquared(model)
+print(R2)
 
 ## Using spatial eigenvectors as predictors to account 
 ## for spatial autocorrelation is a common approach:
@@ -112,7 +114,9 @@ model = sjSDM(Y = com$response, env = linear(com$env_weights),
 summary(model)
 predict(model, newdata = com$env_weights, SP = SPV)
 
-
+## Visualize internal meta-community structure
+an = anova(model)
+plot(an, internal=TRUE)
 
 ## non-linear(deep neural network) model
 model = sjSDM(Y = com$response, env = linear(com$env_weights), 
@@ -151,6 +155,9 @@ result = anova(model)
 print(result)
 plot(result)
 
+## visualize meta-community structure
+plot(result, internal=TRUE)
+
 
 # Deep neural network
 ## we can fit also a deep neural network instead of a linear model:
diff --git a/sjSDM/man/figures/README-unnamed-chunk-7-1.png b/sjSDM/man/figures/README-unnamed-chunk-7-1.png
new file mode 100644
index 0000000..9ee2982
Binary files /dev/null and b/sjSDM/man/figures/README-unnamed-chunk-7-1.png differ
diff --git a/sjSDM/man/importance.Rd b/sjSDM/man/importance.Rd
index 4f6a7f4..9a6edfe 100644
--- a/sjSDM/man/importance.Rd
+++ b/sjSDM/man/importance.Rd
@@ -27,6 +27,7 @@ Computes standardized variance components with respect to abiotic, biotic, and s
 }
 \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).
 }
 \examples{
 \dontrun{
diff --git a/sjSDM/man/linear.Rd b/sjSDM/man/linear.Rd
index aef5f47..a2dc898 100644
--- a/sjSDM/man/linear.Rd
+++ b/sjSDM/man/linear.Rd
@@ -87,6 +87,8 @@ model = sjSDM(Y = com$response, env = linear(com$env_weights),
               iter = 50L) # increase iter for your own data 
 summary(model)
 predict(model, newdata = com$env_weights, SP = XY)
+R2 = Rsquared(model)
+print(R2)
 
 ## Using spatial eigenvectors as predictors to account 
 ## for spatial autocorrelation is a common approach:
@@ -97,7 +99,9 @@ model = sjSDM(Y = com$response, env = linear(com$env_weights),
 summary(model)
 predict(model, newdata = com$env_weights, SP = SPV)
 
-
+## Visualize internal meta-community structure
+an = anova(model)
+plot(an, internal=TRUE)
 
 ## non-linear(deep neural network) model
 model = sjSDM(Y = com$response, env = linear(com$env_weights), 
@@ -136,6 +140,9 @@ result = anova(model)
 print(result)
 plot(result)
 
+## visualize meta-community structure
+plot(result, internal=TRUE)
+
 
 # Deep neural network
 ## we can fit also a deep neural network instead of a linear model:
diff --git a/sjSDM/man/logLik.sjSDM.Rd b/sjSDM/man/logLik.sjSDM.Rd
index 9c1405a..9970a3b 100644
--- a/sjSDM/man/logLik.sjSDM.Rd
+++ b/sjSDM/man/logLik.sjSDM.Rd
@@ -2,18 +2,20 @@
 % Please edit documentation in R/sjSDM.R
 \name{logLik.sjSDM}
 \alias{logLik.sjSDM}
-\title{Extract Log-Likelihood from a fitted sjSDM model}
+\title{Extract negative-log-Likelihood from a fitted sjSDM model}
 \usage{
-\method{logLik}{sjSDM}(object, ...)
+\method{logLik}{sjSDM}(object, individual = FALSE, ...)
 }
 \arguments{
 \item{object}{a model fitted by \code{\link{sjSDM}}}
 
-\item{...}{optional arguments for compatibility with the generic function, no functionality implemented}
+\item{individual}{returns internal ll structure, mostly for internal useage}
+
+\item{...}{optional arguments passed to internal logLik function (only used if \code{individual=TRUE})}
 }
 \value{
-Numeric value
+Numeric value or numeric matrix if individual is true.
 }
 \description{
-Extract Log-Likelihood from a fitted sjSDM model
+Extract negative-log-Likelihood from a fitted sjSDM model
 }
diff --git a/sjSDM/man/plot.sjSDM.DNN.Rd b/sjSDM/man/plot.sjSDM.DNN.Rd
new file mode 100644
index 0000000..e4400a6
--- /dev/null
+++ b/sjSDM/man/plot.sjSDM.DNN.Rd
@@ -0,0 +1,36 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot.R
+\name{plot.sjSDM.DNN}
+\alias{plot.sjSDM.DNN}
+\title{Training history}
+\usage{
+\method{plot}{sjSDM.DNN}(x, ...)
+}
+\arguments{
+\item{x}{a model fitted by \code{\link{sjSDM}} with \code{\link{DNN}} object}
+
+\item{...}{passed to plot}
+}
+\value{
+No return value, called for side effects.
+}
+\description{
+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) 
+
+#create a group dataframe for plot
+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.sjSDManova.Rd b/sjSDM/man/plot.sjSDManova.Rd
index d668aa1..83192dd 100644
--- a/sjSDM/man/plot.sjSDManova.Rd
+++ b/sjSDM/man/plot.sjSDManova.Rd
@@ -8,8 +8,10 @@
   x,
   y,
   type = c("Deviance", "Nagelkerke", "McFadden"),
+  internal = FALSE,
   cols = c("#7FC97F", "#BEAED4", "#FDC086"),
   alpha = 0.15,
+  env_deviance = NULL,
   ...
 )
 }
@@ -18,17 +20,34 @@
 
 \item{y}{unused argument}
 
-\item{type}{use of deviance or of Nagelkerke or McFadden R-squared}
+\item{type}{deviance, Nagelkerke or McFadden R-squared}
+
+\item{internal}{logical, plot internal or total structure}
 
 \item{cols}{colors for the groups}
 
 \item{alpha}{alpha for colors}
 
-\item{...}{Additional arguments to pass to \code{plot()}}
+\item{env_deviance}{environmental deviance}
+
+\item{...}{Additional arguments to pass to \code{plot()}
+
+The \code{internal = TRUE} plot was heavily inspired by Leibold et al., 2022}
 }
 \value{
-The visualized matrix is silently returned
+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{
 Plot anova results
 }
+\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/setWeights.Rd b/sjSDM/man/setWeights.Rd
index bdb85a4..8f534ef 100644
--- a/sjSDM/man/setWeights.Rd
+++ b/sjSDM/man/setWeights.Rd
@@ -12,7 +12,7 @@ setWeights(object, weights)
 \arguments{
 \item{object}{object of class  \code{\link{sjSDM}} with \code{\link{DNN}} object}
 
-\item{weights}{list of layer weights and sigma, see \code{\link{getWeights}}}
+\item{weights}{list of layer weights:  \code{list(env=list(matrix(...)), spatial=list(matrix(...)), sigma=matrix(...))}, see \code{\link{getWeights}}}
 }
 \value{
 No return value, weights are changed in place.
diff --git a/sjSDM/man/simulate.sjSDM.Rd b/sjSDM/man/simulate.sjSDM.Rd
index e5e63c0..7a01ec2 100644
--- a/sjSDM/man/simulate.sjSDM.Rd
+++ b/sjSDM/man/simulate.sjSDM.Rd
@@ -16,8 +16,9 @@
 \item{...}{optional arguments for compatibility with the generic function, no functionality implemented}
 }
 \value{
-Array of simulated species occurrences.
+Array of simulated species occurrences of dimension order [nsim, sites, species]
 }
 \description{
-Simulate nsim responses from the fitted model
+Simulate nsim responses from the fitted model following a multivariate probit model. 
+So currently only supported for \code{family = stats::binomial("probit")}
 }
diff --git a/sjSDM/man/simulate_SDM.Rd b/sjSDM/man/simulate_SDM.Rd
index 0030675..c8de427 100644
--- a/sjSDM/man/simulate_SDM.Rd
+++ b/sjSDM/man/simulate_SDM.Rd
@@ -50,7 +50,7 @@ List of simulation results:
 \item{link}{Which link}
 \item{response_type}{Which response type}
 \item{response}{Species occurrence matrix}
-\item{correlation}{Species covariance matirx}
+\item{correlation}{Species covariance matrix}
 \item{species_weights}{Species-environment coefficients}
 \item{env_weights}{Environmental covariates}
 \item{corr_acc}{Method to calculate sign accurracy}
diff --git a/sjSDM/man/sjSDM.Rd b/sjSDM/man/sjSDM.Rd
index 0565dcb..16be929 100644
--- a/sjSDM/man/sjSDM.Rd
+++ b/sjSDM/man/sjSDM.Rd
@@ -230,6 +230,8 @@ model = sjSDM(Y = com$response, env = linear(com$env_weights),
               iter = 50L) # increase iter for your own data 
 summary(model)
 predict(model, newdata = com$env_weights, SP = XY)
+R2 = Rsquared(model)
+print(R2)
 
 ## Using spatial eigenvectors as predictors to account 
 ## for spatial autocorrelation is a common approach:
@@ -240,7 +242,9 @@ model = sjSDM(Y = com$response, env = linear(com$env_weights),
 summary(model)
 predict(model, newdata = com$env_weights, SP = SPV)
 
-
+## Visualize internal meta-community structure
+an = anova(model)
+plot(an, internal=TRUE)
 
 ## non-linear(deep neural network) model
 model = sjSDM(Y = com$response, env = linear(com$env_weights), 
@@ -279,6 +283,9 @@ result = anova(model)
 print(result)
 plot(result)
 
+## visualize meta-community structure
+plot(result, internal=TRUE)
+
 
 # Deep neural network
 ## we can fit also a deep neural network instead of a linear model:
diff --git a/sjSDM/vignettes/Dependencies.Rmd b/sjSDM/vignettes/Dependencies.Rmd
index d954cd7..b7021f1 100644
--- a/sjSDM/vignettes/Dependencies.Rmd
+++ b/sjSDM/vignettes/Dependencies.Rmd
@@ -16,12 +16,12 @@ editor_options:
   chunk_output_type: console
 ---
 
-```{r, echo = F, message = F}
+```{r, echo = F, message = F, purl=FALSE}
 set.seed(123)
 ```
 
 
-```{r global_options, include=FALSE}
+```{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)
 ```
 
@@ -44,7 +44,7 @@ Sometimes the automatic 'miniconda' installation (via
 Download and install the latest ['conda' version](https://www.anaconda.com/products/individual)
 
 Afterwards run:
-```{r,eval=FALSE}
+```{r,eval=FALSE, purl=FALSE}
 install_sjSDM(version = c("gpu")) # or "cpu" if you do not have a proper gpu device
 ```
 
@@ -63,7 +63,7 @@ Open the command window (cmd.exe - hit windows key + r and write cmd)
 
 Run in cmd.exe:
 
-```{bash, eval=FALSE}
+```{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
@@ -81,7 +81,7 @@ Restart R, try to run the example, and if this doesn't work:
 ## Linux - automatic installation:
  
 Run in R: 
-```{r,eval=FALSE}
+```{r,eval=FALSE, purl=FALSE}
 install_sjSDM(version = c("gpu")) # or "cpu" if you do not have a proper gpu device
 ```
 
@@ -99,7 +99,7 @@ We strongly advise to use a 'conda' environment but a virtual environment should
 Download and install the latest ['conda' version](https://www.anaconda.com/products/individual)
 
 Open your terminal and run:
-```{bash, eval=FALSE}
+```{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
@@ -117,7 +117,7 @@ Restart R try to run the example, if this doesn't work:
 ## MacOS - automatic installation:
  
 Run in R: 
-```{r,eval=FALSE}
+```{r,eval=FALSE, purl=FALSE}
 install_sjSDM()
 ```
 
@@ -135,7 +135,7 @@ We strongly advise to use a 'conda' environment but a virtual environment should
 Download and install the latest conda [conda version](https://www.anaconda.com/products/individual)
 
 Open your terminal and run:
-```{bash, eval=FALSE}
+```{bash, eval=FALSE, purl=FALSE}
 conda create --name r-sjsdm python=3.7
 conda activate r-sjsdm
 python -m pip install torch torchvision torchaudio 
@@ -163,6 +163,6 @@ community <- simulate_SDM(sites = 100, species = 10, env = 5)
 Env <- community$env_weights
 Occ <- community$response
 
-model <- sjSDM(X = Env, Y = Occ, formula = ~0+X1*X2 + X3 + X4)
+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/sjSDM.R b/sjSDM/vignettes/sjSDM.R
index 45e807c..092b186 100644
--- a/sjSDM/vignettes/sjSDM.R
+++ b/sjSDM/vignettes/sjSDM.R
@@ -7,15 +7,6 @@ knitr::opts_chunk$set(fig.width=7, fig.height=4.5, fig.align='center', warning=F
 ## ----eval=FALSE---------------------------------------------------------------
 #  library(sjSDM)
 
-## ---- eval = F----------------------------------------------------------------
-#  install_sjSDM()
-
-## ---- eval = F----------------------------------------------------------------
-#  install_sjSDM(method = "gpu")
-
-## ----eval = F-----------------------------------------------------------------
-#  vignette("Dependencies", package = "sjSDM")
-
 ## ----eval=FALSE---------------------------------------------------------------
 #  citation("sjSDM")
 
@@ -135,7 +126,7 @@ knitr::opts_chunk$set(fig.width=7, fig.height=4.5, fig.align='center', warning=F
 #  summary(model)
 
 ## ----eval=FALSE---------------------------------------------------------------
-#  zero_one = function(x) function(x) (x-min(x))/(max(x) -min(x))
+#  zero_one = function(x)  (x-min(x))/(max(x) -min(x))
 #  Y_pa = zero_one(log(Y+0.001))
 
 ## ----eval=FALSE---------------------------------------------------------------
diff --git a/sjSDM/vignettes/sjSDM.Rmd b/sjSDM/vignettes/sjSDM.Rmd
index 21521cd..51e7024 100644
--- a/sjSDM/vignettes/sjSDM.Rmd
+++ b/sjSDM/vignettes/sjSDM.Rmd
@@ -1,12 +1,13 @@
 ---
-title: "sjSDM: a scalable joint Species Distribution Model"
+title: "Getting started with sjSDM: a scalable joint Species Distribution Model"
 author: "Maximilian Pichler & Florian Hartig, Theoretical Ecology, University of Regensburg"
 date: "`r Sys.Date()`"
 output: 
+
   rmarkdown::html_vignette:
     toc: true
 vignette: >
-  %\VignetteIndexEntry{sjSDM: a scalable joint Species Distribution Model}
+  %\VignetteIndexEntry{Getting started with sjSDM: a scalable joint Species Distribution Model}
   \usepackage[utf8]{inputenc}
   %\VignetteEngine{knitr::rmarkdown}
 abstract: "sjSDM is a scalable joint species distribution model based on Monte-Carlo intergration of the joint likelihood. \n \n \n"
@@ -20,7 +21,7 @@ set.seed(123)
 
 
 ```{r global_options, include=FALSE}
-knitr::opts_chunk$set(fig.width=7, fig.height=4.5, fig.align='center', warning=FALSE, message=FALSE, cache = F)
+knitr::opts_chunk$set(fig.width=7, fig.height=4.5, fig.align='center', warning=FALSE, message=FALSE, cache = FALSE)
 ```
 
 
@@ -38,19 +39,19 @@ library(sjSDM)
 sjSDM depends on the Anaconda distribution of python and pytorch. You will get a warning if you don't have python installed. sjSDM depends on Anaconda which will be installed automatically if it is not available.
 
 To install the CPU version, run:
-```{r, eval = F}
+```{r, eval = F, purl=FALSE}
 install_sjSDM()
 ```
 
 For the GPU version (NVIDIA GPUs are only supported), run:
-```{r, eval = F}
+```{r, eval = F, purl=FALSE}
 install_sjSDM(method = "gpu")
 ```
 
 
 More details on and trouble-shooting for installing the necessary dependencies is explained in the vignette dependencies, call 
 
-```{r,eval = F}
+```{r,eval = F, purl=FALSE}
 vignette("Dependencies", package = "sjSDM")
 ```
 
@@ -159,8 +160,17 @@ an = anova(model)
 plot(an)
 ```
 
+## Internal meta-community structure
+We can also visualize the internal meta-community structure (only for models with space):
+```{r}
+plot(an, internal = TRUE)
+```
+
+
+
 ## Importance
-We can also partition the effects of the different components
+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)
@@ -271,19 +281,4 @@ summary(model)
 ```{r,eval=FALSE}
 model = sjSDM(log(Y+0.01), env = linear(X, ~.), se = TRUE, iter = 50L, family = gaussian("identity"))
 summary(model)
-```
-
-
-## Workaround to use the default binomial with probit or logit link:
-The log-normalization is one way to transform species counts into quasi-probabilities in the range of [0,1]:
-```{r,eval=FALSE}
-zero_one = function(x) function(x) (x-min(x))/(max(x) -min(x))
-Y_pa = zero_one(log(Y+0.001))
-```
-
-We can use the quasi probabilities the same way as presence-absence data:
-```{r,eval=FALSE}
-
-model = sjSDM(Y_pa, env = linear(X, ~.), se = TRUE, iter = 50L, family = binomial("logit"))
-summary(model)
 ```
\ No newline at end of file


Test patch:
diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
new file mode 100644
index 0000000..3f127f6
--- /dev/null
+++ b/.github/workflows/tests.yaml
@@ -0,0 +1,113 @@
+# For help debugging build failures open an issue on the RStudio community with the 'github-actions' tag.
+# https://community.rstudio.com/new-topic?category=Package%20development&tags=github-actions
+on:
+  schedule:
+  - cron: '0 0 * * *'
+  push:
+    branches:
+      - main
+      - master
+  pull_request:
+    branches:
+      - main
+      - master
+
+name: Tests
+
+jobs:
+  Tests:
+    runs-on: ${{ matrix.config.os }}
+
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: windows-latest, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+
+    env:
+      R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
+      RSPM: ${{ matrix.config.rspm }}
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      KMP_DUPLICATE_LIB_OK: TRUE
+    steps:
+      - uses: actions/checkout@v2
+
+      - uses: r-lib/actions/setup-r@v1
+        with:
+          r-version: ${{ matrix.config.r }}
+
+      - uses: r-lib/actions/setup-pandoc@v1
+      
+      - uses: r-lib/actions/setup-tinytex@master
+      
+      - name: Install makeindex
+        run: tlmgr install makeindex
+
+      - name: Query dependencies
+        run: |
+          setwd("sjSDM")
+          install.packages("pak", repos = "https://r-lib.github.io/p/pak/dev/")
+          saveRDS(pak::pkg_deps("local::.", dependencies = TRUE), "../.github/r-depends.rds")
+        shell: Rscript {0}
+
+      - name: Cache R packages
+        if: runner.os != 'Windows'
+        uses: actions/cache@v2
+        with:
+          path: ${{ env.R_LIBS_USER }}
+          key: ${{ matrix.config.os }}-${{ steps.install-r.outputs.installed-r-version }}-1-${{ hashFiles('.github/r-depends.rds') }}
+          restore-keys: ${{ matrix.config.os }}-${{ steps.install-r.outputs.installed-r-version }}-1-
+
+      - name: Install system dependencies
+        if: runner.os == 'Linux'
+        run: |
+          setwd("sjSDM")
+          pak::local_system_requirements(execute = TRUE)
+        shell: Rscript {0}
+          
+      - name: Install miniconda r-sjsdm env # taken from gluonts R package
+        uses: conda-incubator/setup-miniconda@v2
+        with:
+          activate-environment: r-sjsdm
+          python-version: 3.7.1
+          environment-file: _environment.yml
+
+      - name: Install dependencies
+        run: |
+          setwd("sjSDM")
+          pak::local_install_dev_deps(upgrade = TRUE)
+          pak::pkg_install("devtools")
+        shell: Rscript {0}
+        
+      - name: Install Pkg
+        run: |
+          cd sjSDM
+          Rscript -e "pak::local_install()"
+          Rscript -e "pak::pkg_install('rstudio/reticulate')"
+          Rscript -e "sjSDM::install_sjSDM(version='cpu')"
+          Rscript -e "reticulate::py_install('pytest', envname = 'r-sjsdm', pip=TRUE)"
+          
+      - name: tests
+        env:
+          _R_CHECK_CRAN_INCOMING_REMOTE_: false
+          CI: false
+        run: |
+          setwd("sjSDM")
+          options(crayon.enabled = TRUE)
+          devtools::test()
+        shell: Rscript {0}
+        
+      - name: testthat output
+        if: always()
+        run: find check -name 'testthat.Rout*' -exec cat '{}' \; || true
+        shell: bash
+        
+        
+      - name: Upload check results
+        if: failure()
+        uses: actions/upload-artifact@main
+        with:
+          name: ${{ matrix.config.os }}-r${{ matrix.config.r }}-results
+          path: check
diff --git a/Code/utils/test-8-vignette.R b/Code/utils/test-8-vignette.R
new file mode 100644
index 0000000..18f7ec5
--- /dev/null
+++ b/Code/utils/test-8-vignette.R
@@ -0,0 +1,66 @@
+source("utils.R")
+
+
+paths = c("../../vignettes", "../vignettes")
+if((!any(sapply(paths, dir.exists)))) testthat::skip("No vignettes found...")
+path = paths[sapply(paths, dir.exists)]
+rmds = c(paste0(path, "/Dependencies.Rmd"),
+         paste0(path, "/sjSDM.Rmd")
+)
+
+trim_name = function(name) {
+  name = stringr::str_replace_all( name , "'", " ")
+  name = stringr::str_replace_all( name , '"', " ")
+  name = stringr::str_remove_all(name, "-")
+  return(name)
+}
+
+create_test = function(input_file) {
+  name = strsplit(input_file, ".", fixed = TRUE)
+  name = name[[1]][length(name[[1]])-1]
+  name = strsplit(name, "/", fixed = TRUE)[[1]][3]
+  output = paste0(name, ".R")
+  knitr::opts_chunk$set(comment = NA)
+  knitr::purl(input_file, output = output, documentation = 1, quiet = TRUE)
+  knitr::opts_chunk$set(comment = "##")
+  input = file(output)
+  code = readLines(input)
+  close(input)
+  file.remove(output)
+  result = ''
+  open = FALSE
+  end = c( "list2env(as.list(environment()), envir = .GlobalEnv)", "}, NA)})" )
+  
+  # taken from https://stackoverflow.com/questions/6451152/how-to-catch-integer0
+  is.empty <- function(x, mode = NULL){
+    if (is.null(mode)) mode <- class(x)
+    identical(vector(mode, 1), c(x, vector(class(x), 1)))
+  }
+  
+  for(i in 1:length(code)){
+    start = grep("## ----", code[i])
+    if(is.empty(start)) start=FALSE
+    if(start == 1) start = TRUE
+    if(start){
+      if(open){
+        result = c(result, end)
+        result = c(result, paste0("testthat::test_that('",name,"__", trim_name(code[i]),
+                                  "', {testthat::expect_error( {"))
+      }else{
+        result = c(result, paste0("testthat::test_that('",name,"__", trim_name(code[i]),
+                                  "', {testthat::expect_error( {"))
+        open = TRUE
+      }
+    }
+    result = c(result, code[i])
+  }
+  
+  result = c(result, end) #, "rm(list=ls())", "gc()")
+  test_file = paste0("test-", name, ".R")
+  writeLines(result, con = test_file)
+  return(test_file)
+}
+
+vignette_tests = sapply(rmds, create_test)
+sapply(vignette_tests, source)
+file.remove(vignette_tests)
diff --git a/sjSDM/tests/testthat/test-1-anova.R b/sjSDM/tests/testthat/test-1-anova.R
index 02a5df8..939e572 100644
--- a/sjSDM/tests/testthat/test-1-anova.R
+++ b/sjSDM/tests/testthat/test-1-anova.R
@@ -33,8 +33,8 @@ test_model = function(occ = NULL, env, spatial=NULL, biotic = bioticStruct(),
   # iter, batch_size, se, link
   Funcs = list(
     list(2, 2, FALSE, binomial("logit")),
-    list(2, 23, FALSE, poisson()),
-    list(2, 40, FALSE, gaussian())
+    list(2, 23, FALSE, poisson())
+    #list(2, 40, FALSE, gaussian())
   )
   testthat::test_that("sjSDM anova Func", {
     testthat::skip_on_cran()
diff --git a/sjSDM/tests/testthat/test-2-sjSDM_cv.R b/sjSDM/tests/testthat/test-2-sjSDM_cv.R
index 2a675ed..a58dea3 100644
--- a/sjSDM/tests/testthat/test-2-sjSDM_cv.R
+++ b/sjSDM/tests/testthat/test-2-sjSDM_cv.R
@@ -49,4 +49,11 @@ testthat::test_that("sjSDM_cv", {
   #testthat::expect_error(suppressWarnings(plot(model)), NA)
   testthat::expect_error(summary(model), NA)  
   
+  
+  ### sjSDM.tune ###
+  testthat::expect_error(sjSDM.tune(sjSDM_cv(Y1, X1,spatial = linear(SP), iter = 1L, CV = 2L, tune_steps = 3L, n_cores = 2L, device=device, sampling=10L)), NA)
+  testthat::expect_error(sjSDM.tune(sjSDM_cv(Y1, X1,spatial = linear(SP), iter = 1L, family = binomial(),CV = 2L, tune_steps = 3L, n_cores = 2L, device=device, sampling=10L)), NA)
+  testthat::expect_error(sjSDM.tune(sjSDM_cv(Y1, X1,spatial = linear(SP), iter = 1L, family = binomial(), control = sjSDMControl(), CV = 2L, tune_steps = 3L, n_cores = 2L, device=device, sampling=10L)), NA)
+  testthat::expect_error(sjSDM.tune(sjSDM_cv(Y1, X1,spatial = linear(SP), CV = 2L, tune_steps = 3L, n_cores = 2L)), NA)
+  
 })
diff --git a/sjSDM/tests/testthat/test-4-sjSDM.R b/sjSDM/tests/testthat/test-4-sjSDM.R
index 7c06c0d..a0699e3 100644
--- a/sjSDM/tests/testthat/test-4-sjSDM.R
+++ b/sjSDM/tests/testthat/test-4-sjSDM.R
@@ -21,6 +21,9 @@ test_model = function(occ = NULL, env, spatial=NULL, biotic = bioticStruct(),
     testthat::expect_error({ .k = testthat::capture_output(coef(model)) }, NA)
     testthat::expect_error({.k = testthat::capture_output(print(model))}, NA)
     testthat::expect_error({ .k = testthat::capture_output(summary(model)) }, NA)
+    testthat::expect_error({ .k = testthat::capture_output(Rsquared(model)) }, NA)
+    testthat::expect_error({ .k = testthat::capture_output(update(model, env_formula=~1)) }, NA)
+    if(inherits(model, "spatial")) testthat::expect_error({ .k = testthat::capture_output(update(model, env_formula=~1, spatial_formula =~0)) }, NA)
     testthat::expect_false(any(is.na(model$history)))
     testthat::expect_error({ .k= testthat::capture_output(predict(model, batch_size=step_size)) }, NA)
     if(inherits(env, "matrix"))testthat::expect_error({ .k= testthat::capture_output(predict(model, newdata = env, batch_size=step_size)) }, NA)
@@ -256,4 +259,57 @@ testthat::test_that("sjSDM reload model", {
 })
 
 
+testthat::test_that("Changing weights", {
+  testthat::skip_on_cran()
+  testthat::skip_on_ci()
+  skip_if_no_torch()
+  sjSDM:::check_module()
+  XY = matrix(rnorm(100*2), 100L, 2L)
+  com = simulate_SDM(env = 3L, species = 5L, sites = 100L)
+  model = sjSDM(Y = com$response,env = com$env_weights, spatial = linear(XY), iter = 2L)
+  
+  setWeights(model, list(matrix(1.0, 5, 4)))
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$env_weights[[0]])), 1)
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$env_weights[[0]])), c(5, 4))
+  setWeights(model, list(matrix(2.0, 5, 4)))
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$env_weights[[0]])), 2)
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$env_weights[[0]])), c(5, 4))
+  
+  setWeights(model, list(matrix(3.0, 5, 4), matrix(4.0, 5, 3)))
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$env_weights[[0]])), 3)
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$spatial_weights[[0]])), 4)
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$env_weights[[0]])), c(5, 4))
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$spatial_weights[[0]])), c(5, 3))
+  
+  setWeights(model, list(matrix(5.0, 5, 4), list(matrix(6.0, 5, 3))))
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$env_weights[[0]])), 5)
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$spatial_weights[[0]])), 6)
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$env_weights[[0]])), c(5, 4))
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$spatial_weights[[0]])), c(5, 3))
+  
+  setWeights(model, list(matrix(5.0, 5, 4), list(matrix(6.0, 5, 3)), matrix(1.0, 5, 5)))
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$env_weights[[0]])), 5)
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$spatial_weights[[0]])), 6)
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$sigma$data$cpu()$numpy() )), 1)
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$env_weights[[0]])), c(5, 4))
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$spatial_weights[[0]])), c(5, 3))
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$sigma$data$cpu()$numpy() )), c(5, 5))
+  
+  setWeights(model, list(matrix(5.0, 5, 4), NULL, matrix(1.0, 5, 5)))
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$env_weights[[0]])), 5)
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$spatial_weights[[0]])), 6)
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$sigma$data$cpu()$numpy() )), 1)
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$env_weights[[0]])), c(5, 4))
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$spatial_weights[[0]])), c(5, 3))
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$sigma$data$cpu()$numpy() )), c(5, 5))
+  
+  setWeights(model, list(NULL, NULL, matrix(1.0, 5, 5)))
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$env_weights[[0]])), 5)
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$spatial_weights[[0]])), 6)
+  testthat::expect_equal(mean(reticulate::py_to_r(model$model$sigma$data$cpu()$numpy() )), 1)
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$env_weights[[0]])), c(5, 4))
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$spatial_weights[[0]])), c(5, 3))
+  testthat::expect_equal(dim(reticulate::py_to_r(model$model$sigma$data$cpu()$numpy() )), c(5, 5))
+  
+})
 
diff --git a/sjSDM/tests/testthat/utils.R b/sjSDM/tests/testthat/utils.R
index 95b628b..0e66ef3 100644
--- a/sjSDM/tests/testthat/utils.R
+++ b/sjSDM/tests/testthat/utils.R
@@ -9,6 +9,6 @@ force_r = function(x) {
 }
 
 is_gpu_available = function() {
-  if( force_r(pkg.env$torch$cuda$is_available()) ) return("gpu")
+  if( force_r(sjSDM:::pkg.env$torch$cuda$is_available()) ) return("gpu")
   else return("cpu")
 }
\ No newline at end of file