niche_neutral_GLMM_tutorial.Rmd

---
title: "Tutorial on niche-neutral GLMMs to understand community assembly"
author: "Sara Mortara, Alexandre Adalardo & Paulo Prado"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: 
 html_document:
    toc: true
    toc_float: true

---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

### Presentation

The general idea of our framework is to translate niche and neutral dynamics on community assembly into fixed and random effects of generalized mixed models. 
Within our framework one can simultaneously test the importance of niche solely, neutrality solely and combinations of niche and neutral effects on community structure (abundance or occurrence). 

We use the model selection approach to test multiple hypothesis and we use conditional and marginal $R^{2}$ values to quantify the relative importance of niche and neutrality on community structure. 

Main criticism of our framework rely on: (1) Are fixed and random effects actually capturing niche and neutral dynamics? (2) Are random effects only capturing neutral dynamics or over-dispersion in data and uninformed species traits are inflating random effects? 

In order to show how we implemented our framework and how model components are translated into niche and neutral dynamics we perform simulations as described below. We do not intend to scan all parameter space of type of communities and traits, neither make a comprehensive investigation of sampling effects. Our goal with this tutorial is to address main criticisms of our framework and show how promising our framework is to quantify the relative importance of niche and neutral dynamics.           

## A test of our modelling framework

Here, we present a working example of our modelling framework. We use simulated communities based on niche and neutral dynamics in order to show how fixed and random effects capture different ecological processes. Also, we use strong and week traits in order to show how fixed and random effects capture neutral dynamics when traits are strong. Therefore, we simulate stochastic and deterministic meta-communities and use strong or weak traits in model selection. 

We simulate meta-communities with the same data structure as our abundance data of ferns in three mountain chains in southern Brazil. Then, we make Poisson and negative binomial samples of communities and use samples from the meta-communities in our model selection framework.

In order to access main criticisms of our framework, we look into four scenarios:

* Deterministic community with right traits from Poisson sample
* Deterministic community with wrong traits from Poisson sample
* Stochastic community with right traits from Poisson sample
* Stochastic community with right traits from Negative Binomial sample


### 1. Building simulated communities

#### 1.1. Getting data for simulations

The simulations are based on several packages that should be installed in the user's computer before starting the tutorial. The above packages should be installed as regular CRAN packages before loaded. To check if they are already installed run the above code and if not installed, run the code below.

```{r,  message=FALSE}
packages <- c("devtools", "ade4", "vegan", "lme4", "bbmle", "optimx", "piecewiseSEM", "sads", "ggplot2", "dplyr", "MASS", "knitr")
instpack <- packages[!packages %in% installed.packages()]
if(length(instpack)>0)
{
    install.packages(packages[!packages %in% installed.packages()])
}
```

One of the packages needed in this tutorial, MCSim (MetaCommunity Simulation) is not available at CRAN and should be installed directly from the development repository github. Follow the above code to install it. Information about this package can be found at <a href="http://rstudio-pubs-static.s3.amazonaws.com/159425_80725873417e42fdb13821c10a198281.html"> MCSim </a>

```{r,  message=FALSE}
if(! "MCSim" %in% installed.packages())
{
    devtools::install_github("sokole/MCSim")
}


```

Attach all packages needed to running this tutorial.


```{r packages, results='hide', message=FALSE}
library(devtools)
library(MCSim)
library(ade4)
library(vegan)
library(lme4)
library(bbmle)
library(optimx)
library(piecewiseSEM)
library(sads)
library(ggplot2)
library(dplyr)
library(MASS)
library(knitr)
```

Here we specify the random number generator seed:

```{r seed}
set.seed(42)
```

To evaluate our methods, we generate the meta-community data. Part of the code used here is adapted from <b> <a href="http://rstudio-pubs-static.s3.amazonaws.com/159425_80725873417e42fdb13821c10a198281.html">  Sokol's MCSim tutorial </a> </b>. We create data from 30 sites across three mountain chains. We simulate abundance values for 20 different species.

For both deterministic and stochastic scenarios, we simulate 10 different meta-communities until time-step 100 and then take mean values from the 100^ht^ last time-step.

```{r generate data}

## Number of sites
Nsites <- 30
## Number of regions
Nregions <- 3
## Number of species
Nspp <- 20

## Sites attributes ##
## Here are the xy-coordinates of sites: 3 regions with 10 sites each
## Distances between regions is an order of amgnitude distance bewteen sites within regions
sites <- expand.grid(x = seq(0,120, length=Nregions), y=seq(1,5, length=Nsites/Nregions))
sites$x <- jitter(sites$x)
sites$y <- jitter(sites$y)
## Each set of 10 points at same x coordinate is labelled as a region
sites$region <- rep(letters[1:Nregions], Nsites/Nregions)
## Enviromental variable: 10 states, that repeat at each site (e.g. altitude)
sites$env <- rep(seq(1,5, length=Nsites/Nregions), each=Nregions)

## Calculate niches for species: optimal values along the enviromental variable
sp.opt <- runif(Nspp, min = 1, max = 5)

```

Here, we generate the sites x species matrix with community initial conditions to run the simulations.

```{r initial community}
## Initial condition ##
## Initial condition: matrix of sites x spp

## Random values of species abundances 
m0b <- matrix(rlnorm(Nspp*Nsites), Nsites, Nspp)
## Round values of species abundance to represent discrete values of number of individuals
m0b <- round(m0b)
## Splitting species in 3 fractions that are exclusive of each region
R.ind <- sample(letters[1:Nregions], Nspp, replace=TRUE)
for(i in letters[1:Nregions])
  m0b[sites$region==i,R.ind!=i] <- 0
## Calculating Relative abundances
m0b <- sweep(m0b, 1, apply(m0b,1,sum), FUN="/")
```

#### 1.2 Deterministic community

##### 1.2.1 Building the landscape

```{r landscape, results="hide", message=FALSE}
## Following Sokol's tutorial we arbitrarily chose JM = 1e6
JM <- 1e6

## We set m=0.5 to allow some dispersal limitation. Still a half of the deaths are replaced by locals

simulation_landscape_det <- MCSim::fn.make.landscape(
  site.coords = sites[,1:2],
  Ef = sites$env,
  m = 0.5, 
  JM = JM)
## It seems that in R 3.5 it is necessary to convert the dist mat in this object form dataframe to a matrix
simulation_landscape_det$dist.mat <-  as.matrix(simulation_landscape_det$dist.mat)
```
##### 1.2.2 Simulating abundances of deterministic communities

```{r simu det, results="hide", message=FALSE}

## Data frames to store simulations
id <- data.frame(site=rep(1:Nsites, Nspp), sites[,3:4], spp=rep(paste("sp",1:Nspp, sep="."), each=Nsites))
## To store simulation results
det.resu <- matrix(NA, nrow=nrow(id), ncol=10)

for(i in 1:10){
    simu.det <- MCSim::fn.metaSIM( # simulation of deterministic community w correctly observed traits
                           landscape = simulation_landscape_det,
                           ##output.dir.path = 'FERN_SIM_RESULTS_det',
                           ##scenario.ID = 'fern_det',  
                           ##sim.ID = 'det',
                           trait.Ef = sp.opt,
                           trait.Ef.sd = 0.5,
                           J.t0 = m0b,
                           ##gamma.abund = gama.init/(sum(gama.init)),
                           n.timestep = 100, # increased to 100 time steps; initial conditions seems to persist after t=30
                           W.r = 0,
                           nu = 0,
                           speciation.limit = 0,
                           save.sim = FALSE
                       )
    det <- subset(simu.det$J.long, timestep=="100")[,4]
    det.resu[,i] <- det 
}

det <- data.frame(id, count=rowMeans(det.resu))
```


#### 1.3 Stochastic community

##### 1.3.1 Simulating the landscape

  
```{r simu sto, results="hide", message=FALSE}

simulation_landscape_sto <- MCSim::fn.make.landscape(
  site.coords = sites[,1:2],
  Ef = sites$env,
  m = 0.5, # reducing immigration parameter, allowing limited dispersal
  JM = JM)
## It seems that in R 3.5 it is necessary to convert the dist mat in this object form dataframe to a matrix
simulation_landscape_sto$dist.mat <-  as.matrix(simulation_landscape_sto$dist.mat)
```

##### 1.3.2 Simulating abundances of stochastic communities

```{r sto, results="hide", message=FALSE}
## A matrix for the results
sto.resu <- matrix(NA, nrow=nrow(id), ncol=10)

for(i in 1:10){
    simu.sto <- MCSim::fn.metaSIM(
                           landscape = simulation_landscape_sto,
                           #output.dir.path = 'FERN_SIM_RESULTS_sto_RT',
                           #scenario.ID = 'fern__sto_RT',  
                           #sim.ID = 'sto_RT',
                           trait.Ef = sp.opt,
                           trait.Ef.sd = 1000, # niche deviation changed for neutral dynamics
                           J.t0 = m0b,
                           ##gamma.abund = gama.init/(sum(gama.init)),
                           n.timestep = 100, 
                           W.r = 200, # Dispersal Kernel no longer flat
                           nu = 0,
                           speciation.limit = 0,
                           save.sim = FALSE
                       )
    sto <- subset(simu.sto$J.long, timestep=="100")[,4]
    sto.resu[,i] <- sto 
}

sto <- data.frame(id, count=rowMeans(sto.resu))


```

#### 1.4 Taking samples from simulated communities

Then, we take Poisson and negative binomial samples of different scenarios. 

Here, we are sampling 50% of simulated communities following Poisson and negative binomial distributions. Our objective here is to generate a scenario in which data is sampled from a Poisson distribution, representing the structure of species abundance data and from a negative binomial sample, representing the structure of overdispersed species abundance data. 

```{r sampling}
# Poisson samples
## Deterministic community
det.pois <- rpois(length(det$count), det$count*0.5)
## Stochastic community
sto.pois <- rpois(length(sto$count), sto$count*0.5)

# Negative binomial sample
## Stochastic community
sto.nb <- rnbinom(n= length(sto$count), 
                  size=1, #arbitrally chosen
                  #size=gamma.shape(glm(sto$count ~ 1, family=Gamma))$alpha, 
                  mu=sto$count*0.5)

## Looking into results from simulation and sample
par(mfrow=c(2,3))
plot(rad(sto$count), main="Stochastic simulation")
plot(rad(sto.pois), main="Stochastic Poisson sample")
plot(rad(sto.nb), main="Stochastic Negative Binomial sample")
plot(rad(det$count), main="Deterministic simulation")
plot(rad(det.pois), main="Deterministic Poisson sample")
par(mfrow=c(1,1))
```


### 2. Building hypothesis and models 

We aim to simultaneously test these four main hypotheses:

1. Solely niche dynamics affect species abundance

2. Solely neutral dynamics affect species abundance

3. Niche and neutral dynamics affect species abundance

4. Species abundances vary randomly across regions (null model)

In order to represent each hypothesis, we will build a set of models using GLMM with fixed and random effects as below. In general, niche dynamics are represented by fixed effects, whereas neutral dynamics are represented by random effects. Given that species traits are measured at species level, we included species as a random effect on all models, including niche dynamics models. Therefore, this random effect does not represent neutral dynamics, but a random intercept for species abundances.

Hypothesis | Fixed effects | Random Effects
-----------|---------------|----------------
Solely niche dynamics | trait and environment| species
Solely neutral dynamics | - | species, species within sites and species within regions 
Niche and neutral dynamics | trait and environment | species, species within sites and species within regions
Null model | - | species and region independently

### 3. Performing model selection

#### 3.1. Preparing data table for model selection

##### Preparing abundance, trait and environmental data

Here, we combine species abundance data sampled from deterministic and stochastic simulations with traits and environmental data. 

```{r selecting comm, results="hide", message=FALSE}

# Binding all data togheter
data <- data.frame(site=id[,"site"], spp=id[,"spp"], det=det.pois, sto=sto.pois, 
                   sto.nb=sto.nb) 
    
# Now we need species traits, gradient and spacial info

## Traits 
## A vector with wrong traits with correlation of less than 0.01 with the true traits
cor.t <- 1
while(cor.t>0.01){
    wrong.t <- runif(length(sp.opt), min(sp.opt), max(sp.opt))
    cor.t <- abs(cor(wrong.t, sp.opt))
}
## Trait data
trait.data <- data.frame(spp=unique(id$spp), trait= scale(sp.opt),
                    trait.wr=scale(wrong.t)) # creating vector w/ wrong traits
## Gradient
env.data <- data.frame(site=unique(id$site), grad=scale(sites$env), region=sites$region)

## Preparing data table for model selection
all.data <- merge(data, env.data, by=c("site"))
all.data <- merge(all.data, trait.data, by=c("spp"))
```

Abundances curves of species along the environmental gradient are
as expected:

```{r checking abundance x grad plots}
## A sample of species, a panel for each region
## Niche
all.data %>%
    filter(as.integer(spp)<11) %>%
    ggplot(aes(grad, det)) + geom_line(aes(colour=spp)) +
    scale_y_log10()  +
    facet_wrap(~region) +
   labs(x="Environmental gradient", y="Species abundance", title="Deterministic community from Poisson sample")
all.data %>%
    filter(spp=="sp.10")%>%
    ggplot(aes(grad, det)) + geom_point() + facet_wrap(~region) +
    labs(x="Environmental gradient", y="Species abundance", title="Species sp.10 of deterministic community from Poisson sample")
## Neutral
all.data %>%
#    filter(as.integer(spp)<21) %>%
    ggplot(aes(grad, sto)) + geom_line(aes(colour=spp)) +
    scale_y_log10() +
    facet_wrap(~region) +
    labs(x="Environmental gradient", y="Species abundance", title="Stochastic community from Poisson sample")
all.data %>%
    filter(spp=="sp.10")%>%
    ggplot(aes(grad, sto)) + geom_point() + facet_wrap(~region) +
    labs(x="Environmental gradient", y="Species abundance", title="Species sp.10 of stochastic community from Poisson sample")
## Neutral from Negative Binomial sample 
all.data %>%
#    filter(as.integer(spp)<21) %>%
    ggplot(aes(grad, sto.nb)) + geom_line(aes(colour=spp)) +
    scale_y_log10() +
    facet_wrap(~region) +
    labs(x="Environmental gradient", y="Species abundance", title="Stochastic community from Negative binomial sample")
all.data %>%
    filter(spp=="sp.10")%>%
    ggplot(aes(grad, sto.nb)) + geom_point() + facet_wrap(~region) + 
    labs(x="Environmental gradient", y="Species abundance", title="Species sp.10 of stochastic community from Poisson sample")
```


#### 3.2 Performing model selection

First, we create simple functions to perform model selection. Although the logic is generic, functions below are specific to run this example. We create a function for each of our hypothesis. 

```{r packages model selection}

## 1. Hypothesis of solely niche dynamics
# Niche dynamics models in which abundances are a function of species ecological strategies (traits) or ecological strategies interacting with the environmental gradient (tradeoff model). 
m.niche <- function(ab, spp, trait, grad, family, ...){
      es <- glmer(ab ~ trait + (1|spp), family=family, 
               control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
      es.tradeoff <- glmer(ab ~ trait + I(grad^2) + trait:grad + (1|spp),
                         family=family, 
                    control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
      return(list(es=es, es.tradeoff=es.tradeoff))
}

## 2. Hypothesis of solely neutral dynamics 
# Neutral model in which variation in species abundances are restricted to sites and regions, representing limited dispersal at local and regional scales, respectively.
m.neu <- function(ab, spp, region, site, family, ...){
    drift <- glmer(ab ~ 1 + (1|spp:region) + (1|spp), family=family, 
                   control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
#    drift2 <- glmer(ab ~ 1 + (1|spp:region) + (1|spp:site) + (1|spp), family=family, 
                #   control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
    return(drift)
}

## 3. Hypothesis of niche and neutral dynamics
# Models in which niche and neutral dynamics models are combined = niche dynamic models + random term for limited dispersal
m.nineu <- function(ab, trait, grad, spp, region, family, ...){
    es.drift <- glmer(ab ~ trait + (1|spp:region) + (1|spp), family=family, 
                      nAGQ=1, control=glmerControl(optimizer="bobyqa", 
                                              optCtrl=list(maxfun=2e5)))
    es.tradeoff.drift <- glmer(ab ~ trait + I(grad^2) + trait:grad + 
                                  (1|spp:region) + (1|spp), family=family, 
                                control=glmerControl(optimizer="bobyqa", 
                                                     optCtrl=list(maxfun=2e5)))
 #   es.drift2 <- glmer(ab ~ trait + (1|spp:region) + (1|spp:site) 
  #                    + (1|spp), family=family, 
   #                   nAGQ=1, control=glmerControl(optimizer="bobyqa", 
    #                                          optCtrl=list(maxfun=5e5)))
    #es.tradeoff.drift2 <- glmer(ab ~ trait + I(grad^2) + trait:grad + 
     #                             (1|spp:region) + (1|spp:site) + 
      #                            (1|spp), family=family, 
       #                         control=glmerControl(optimizer="bobyqa", 
        #                                             optCtrl=list(maxfun=5e5)))
       return(list(es.drift=es.drift, es.tradeoff.drift=es.tradeoff.drift))#, 
         #          es.drift=es.drift2, es.tradeoff.drift2=es.tradeoff.drift)) 
}

## 4. Null hypothesis
# Null model in which species and sites are independent random intercepts
m.null <- function(ab, region, spp, family, ...){
      null <- glmer(ab ~ 1 + (1|region) + (1|spp), family=family,
                   control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=5e5)))
      return(null)
}
```

Now, we will apply all the models to different sets of data. 

##### Deterministic community w/ right traits, from poisson sample

```{r models det rt}
# Solely niche dynamics
niche.det.rt <- m.niche(ab=all.data$det,
                        spp=all.data$spp, 
                        trait=all.data$trait, 
                        grad=all.data$grad,
                       family = "poisson") 
# Solely neutral dynamics
neu.det <- m.neu(ab=all.data$det, 
                      spp=all.data$spp, 
                      region=all.data$region, 
                      site=all.data$grad,
                      "poisson")

# Niche and neutral dynamics
nineu.det.rt <- m.nineu(ab=all.data$det, 
                        spp=all.data$spp, 
                        region=all.data$region,
                        trait=all.data$trait, 
                        grad=all.data$grad, 
                        site=all.data$grad,
                        "poisson")
# Null hypothesis
null.det <- m.null(ab=all.data$det, 
                      spp=all.data$spp, 
                      region=all.data$region, 
                      "poisson")
# BIC for each model
mod.det.rt.BIC <- BIC(null.det, 
                      neu.det[[1]], neu.det[[2]], 
                   niche.det.rt[[1]], niche.det.rt[[2]],
                   nineu.det.rt[[1]], nineu.det.rt[[2]], 
                   nineu.det.rt[[3]], nineu.det.rt[[4]]) 
mod.det.rt.BIC <- mod.det.rt.BIC[order(mod.det.rt.BIC[,2]),]
mod.det.rt.BIC
# R2 for the best model explaining data
sem.model.fits(niche.det.rt[[2]])

```

##### Deterministic community w/ wrong traits, from poisson sample

```{r models det wt}
# Solely niche dynamics
niche.det.wt <- m.niche(ab=all.data$det,
                        spp=all.data$spp, 
                        trait=all.data$trait.wr, 
                        grad=all.data$grad, family="poisson") 
# Niche and neutral dynamics
nineu.det.wt <- m.nineu(ab=all.data$det, 
                        spp=all.data$spp, 
                        region=all.data$region,
                        trait=all.data$trait.wr, 
                        grad=all.data$grad, 
                        site=all.data$grad, family="poisson")
# BIC for each model
mod.det.wt.BIC <- BIC(null.det, neu.det[[1]], neu.det[[2]], 
                   niche.det.wt[[1]], niche.det.wt[[2]],
                   nineu.det.wt[[1]], nineu.det.wt[[2]],
                   nineu.det.wt[[3]], nineu.det.wt[[4]]) 
mod.det.wt.BIC <- mod.det.wt.BIC[order(mod.det.wt.BIC[,2]),]
mod.det.wt.BIC
# R2 for the best model explaining data
sem.model.fits(niche.det.wt[[2]])
```

##### Stochastic community w/ right traits, from poisson sample


```{r models sto rt}
# Solely niche dynamics
niche.sto.rt <- m.niche(ab=all.data$sto,
                        spp=all.data$spp, 
                        grad=all.data$grad, 
                        trait=all.data$trait, 
                        family="poisson")
# Solely neutral dynamics
neu.sto.rt <- m.neu(ab=all.data$sto, 
                    spp=all.data$spp, 
                    region=all.data$region,
                    site=all.data$grad,
                    family="poisson")
# Niche and neutral dynamics
nineu.sto.rt <- m.nineu(ab=all.data$sto, 
                        spp=all.data$spp, 
                        region=all.data$region,
                        grad=all.data$grad, 
                        trait=all.data$trait,
                        site=all.data$grad,
                        family="poisson") 
# Null hypothesis
null.sto.rt <- m.null(ab=all.data$sto, 
                      spp=all.data$spp, 
                      region=all.data$region, 
                      family="poisson")                        
# BIC for each model
mod.sto.rt.BIC <- BIC(null.sto.rt, neu.sto.rt[[1]], neu.sto.rt[[2,]], 
                   niche.sto.rt[[1]], niche.sto.rt[[2]],
                   nineu.sto.rt[[1]], nineu.sto.rt[[2]], 
                   nineu.sto.rt[[3]], nineu.sto.rt[[4]])
mod.sto.rt.BIC <- mod.sto.rt.BIC[order(mod.sto.rt.BIC[,2]),]
mod.sto.rt.BIC
# R2 for the best model explaining data
sem.model.fits(neu.sto.rt)
```

##### Stochastic community w/ right traits, from negative binomial

```{r models sto nb rt}
# Solely niche dynamics
niche.stonb.rt <- m.niche(ab=all.data$sto.nb,
                        spp=all.data$spp, 
                        grad=all.data$grad, 
                        trait=all.data$trait, 
                        family="poisson") 
# Solely neutral dynamics
neu.stonb.rt <- m.neu(ab=all.data$sto.nb, 
                    spp=all.data$spp, 
                    region=all.data$region, 
                    family="poisson")
# Niche and neutral dynamics
nineu.stonb.rt <- m.nineu(ab=all.data$sto.nb, 
                        spp=all.data$spp, 
                        region=all.data$region,
                        grad=all.data$grad, 
                        trait=all.data$trait, 
                        family="poisson") 
# Null hypothesis
null.stonb.rt <- m.null(ab=all.data$sto.nb,
                      spp=all.data$spp, 
                      region=all.data$region, 
                      family="poisson")                        
# BIC for each model
mod.stonb.rt.BIC <- BIC(null.stonb.rt, neu.stonb.rt, 
                   niche.stonb.rt[[1]], niche.stonb.rt[[2]],
                   nineu.stonb.rt[[1]], nineu.stonb.rt[[2]])
mod.stonb.rt.BIC <- mod.stonb.rt.BIC[order(mod.stonb.rt.BIC[,2]),]
mod.sto.rt.BIC
# R2 for the best model explaining data
sem.model.fits(neu.stonb.rt)
```

### 4. Addressing main criticisms of our framework

We now will compare results from model selection of the different scenarios we created. We consider the model with lowest BIC (Bayesian Information Criterion) as the model that best explains the community data. We will depict all best models in terms of their adjusted $R^{2}$. Adjusted $R^{2}$, as they measure the relative importance of fixed and random effects in the model, are used here a proxy of correspondence of community processes to terms in the model. Therefore, conditional $R^{2}$ represents the influence of niche dynamics whereas marginal $R^{2}$ is partitioned into each random term and can represent either niche (in the case of the term (1|spp)) or, more commonly, neutral dynamics (for other terms).  

```{r r-squared table}
# Function to calculate R-squared for each random effect based on code from Melina Leite and Nakagawa
# Works only for Poisson GLMMs
r2.table <- function(model){
# Function to calculate the null model, model with all random terms 
best.null <- function(model) {
parens <- function(x) paste0("(",x,")")
onlyBars <- function(form) reformulate(sapply(findbars(form),
                                              function(x)  parens(deparse(x))),
                                              response=".")
onlyBars(formula(model))
best.null <- update(model,onlyBars(formula(model)))
return(best.null)
}
# Calculates null model
m0 <- best.null(model)
# Variance for fixed effects
VarF <- var(as.vector(fixef(model) %*% t(model@pp$X)))
# Denominator for R2GLMM formula works for Poisson distribution only
deno <- (VarF + sum(unlist(VarCorr(model))) +  
              log(1 + 1/exp(as.numeric(fixef(m0)))))
# R2GLMM(m) - marginal R2GLMM 
r2f <- VarF/deno
# R2GLMM(c) - conditional R2GLMM for full model
r2t <- (VarF + sum(unlist(VarCorr(model))))/deno
# R2 random effects only
r2rand <- r2t-r2f
## R2 Residuals
r2res <- 1-r2t
## Partitioning R2 GLMM for each random effect
r2rand.part <- unlist(VarCorr(model))/deno
r2.tab <- t(as.data.frame(c(conditional = r2t,
      fixed = r2f,
      random = r2rand,
      r2rand.part)))
row.names(r2.tab) <- "model"
return(as.data.frame(r2.tab))
}

```

#### 4.1. Are fixed and random effects actually capturing niche and neutral dynamics?

In order to address this question, we compare results of the selected model for deterministic and stochastic communities with right traits, both sampled from Poisson distribution.


```{r test 4.1}
#Calculating R2
r2.det.rt <- r2.table(niche.det.rt[[2]])
r2.sto.rt <- r2.table(neu.sto.rt)
#Creating table
table1 <- bind_rows(r2.det.rt, r2.sto.rt)
row.names(table1) <- c("Deterministic community w/ right traits", 
                   "Stochastic community w/ right traits")
kable(table1)
```

For communities built by known niche dynamics conditional $R^{2}$ value is composed basically by fixed effects, representing niche dynamics. For communities built by known neutral dynamics, conditional $R^{2}$ is built by random effects, specially the effect represented by the term (1|region:spp), which represents of limited dispersal of species within regions.

#### 4.2. Are random effects only capturing neutral dynamics or overdipersion in data and uninformed species traits are inflating random effects? 

First, we will address the question of over-dispersion in data. In  order to examine over-dispersion issues, we compare the results of the selected model from stochastic communities with right traits sampled from Poisson and negative binomial distributions.   
```{r test 4.2 overdisp}
#Calculating R2
# r2.sto.rt
r2.stonb.rt <- r2.table(neu.stonb.rt)
#Creating table
table2 <- bind_rows(r2.sto.rt, r2.stonb.rt)
row.names(table2) <- c("Stochastic community w/ right traits from Poisson", 
                       "Stochastic community w/ right traits from Negative Binomial")
kable(table2)
```


For stochastic communities sampled from both Poisson and Negative Binomial distributions, conditional $R^{2}$ are basically the same, indicating that over-dispersion in data is being well captured by Poisson GLMM's with no prejudice on the interpretation of which processes are affecting species abundances. For stochastic communities, the random term (1|spp) should be particularly small.  

Second, in order to examine how random effects reflect uninformed traits, we compare the results of the selected model from deterministic communities with right and wrong traits. 
```{r test 4.2 right-wrong}
#Calculating R2
# r2.det.rt
r2.det.wt <- r2.table(niche.det.wt[[2]])
#Creating table
table3 <- bind_rows(r2.det.rt, r2.det.wt)
row.names(table3) <- c("Deterministic community w/ right traits", 
                       "Deterministic community w/ wrong traits")
kable(table3)
```

Only looking into model selection results, one cannot examine difference between deterministic communities built by uninformative and informative species traits. However, by examining $R^{2}$ values, one can detect if traits used in the model are actually influencing species abundances. Uninformative traits inflate conditional $R^{2}$ because of the importance of the term (1|spp) and minimizes the influence of the environment. When traits are informative and niche dynamics are preponderant, $R^{2}$ for fixed effects should be particularly high. 

With this tutorial we show how to implement our modeling framework and how robust our proposal is to reflect niche and neutral dynamics. We recognize that further investigation of how fixed and random effects respond to variation on community dynamics is needed to evaluate the robustness of the proposal. However, with this tutorial, we already show that main criticism of our work does not apply.