Welcome to the vignette of GUILDS. GUILDS is an R package aimed at providing the user with easy to use functions that help dealing with the Unified Neutral Theory of Biodiversity and Biogeography. More specifically, the package bundles together the sampling functions for the Neutral Model including multiple guilds that differ in dispersal ability (Janzen et al. 2015), and the Etienne sampling formula for the standard neutral model. Furthermore, the package contains functions to generate artificial datasets, and to calculate the expected abundances.

To generate data with the standard neutral model, we can make use of
the function `generate.ESF`

, where ESF stands for the Etienne
Sampling Formula. generate.ESF makes use of the urn scheme described in
Etienne (2005), where the avid reader can also find more details on the
Etienne Sampling Formula. The ESF combines the universal biodiversity
number theta, and the universal dispersal number I. Typically, the
migration rate is a bit easier to interpret, so we will make use of that
first, then convert it to I, and then generate the data:

```
library(GUILDS)
set.seed(42)
= 100
theta <- 0.1
m <- 10000
J <- m * (J - 1) / (1 - m)
I <- generate.ESF(theta, I, J)
abund abund
```

```
## [1] 875 238 212 177 176 172 159 158 149 145 144 126 124 122 117 112 111 110
## [19] 110 109 108 108 106 101 99 95 93 92 86 80 80 79 78 77 76 75
## [37] 70 68 66 63 62 60 60 59 59 58 57 57 56 56 55 54 54 53
## [55] 52 51 51 50 48 47 46 46 45 44 44 43 43 43 43 42 42 41
## [73] 40 40 39 39 39 37 37 37 36 36 35 35 34 34 34 33 33 31
## [91] 31 30 30 30 30 30 29 29 28 28 28 27 26 26 26 26 26 26
## [109] 26 25 25 24 24 24 24 24 24 23 22 22 22 21 21 21 21 20
## [127] 20 20 20 20 20 20 20 19 19 19 19 18 18 18 18 17 17 17
## [145] 17 17 17 16 16 16 16 16 16 15 15 15 15 15 14 14 14 13
## [163] 13 13 12 12 12 12 12 12 12 11 11 11 11 11 11 10 10 10
## [181] 10 10 9 9 9 9 9 9 9 9 8 8 8 8 8 8 8 8
## [199] 8 8 8 8 8 7 7 7 7 7 7 7 7 7 7 7 7 6
## [217] 6 6 6 6 6 6 6 5 5 5 5 5 5 5 5 5 5 5
## [235] 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4
## [253] 4 4 4 4 4 4 3 3 3 3 3 3 3 3 3 3 3 3
## [271] 3 3 3 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2
## [289] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [307] 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [325] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [343] 1
```

The package has now generated a nice species abundance distribution.
Using the function ``preston_plot`

we can visualize our
abundance distribution

`preston_plot(abund)`

We can even do better, using the function `expected.SAD`

we can calculate the expected number of individuals per species.
Furthermore, we can easily add this to our plot:

```
<- expected.SAD(theta, m, J)
abund.expect preston_plot(abund, abund.expect)
```

We can clearly see that the simulated abundance distribution closely matches the expected distribution (the black line), although the stochastic nature of the simulation has caused some small deviations. Using the Etienne Sampling Formula we can now estimate theta and m, using maximum likelihood. Because we already know the real values (those are the values that we used to simulate the data), we can nicely crosscheck whether the Etienne Sampling Formula yields good results. We start the maximum likelihood at three different starting points, to compare convergence.

```
<- maxLikelihood.ESF(init_vals = c(theta, m), abund)
LL1 <- maxLikelihood.ESF(init_vals = c(100, 0.01), abund)
LL2 <- maxLikelihood.ESF(init_vals = c(10, 0.5), abund)
LL3 $par LL1
```

`## [1] 102.0395279 0.1159656`

`$par LL2`

`## [1] 101.9890695 0.1150816`

`$par LL3`

`## [1] 102.0662586 0.1158091`

Regardless of the starting point, we recover the same parameter estimates that are close to the real values theta = 100 and m = 0.1 .

The guilds model, as described in Janzen et al. (2015), asks a slightly different question from the standard neutral model. How would species abundances look if there would be two guilds of species, where the guilds would only differ in dispersal ability? Imagine for instance a forest, where some tree species disperse using wind-aided seed dispersal, and other tree species disperse using animal-aided seed dispersal. The model assumes that the two guilds share a local community, and the metacommunity, and share the same speciation dynamics (e.g. have the same theta). Dispersal is different however, and in contrast to the standard neutral model, where we would infer migration rates, here we asses the dispersal ability: the guilds model postulates that migration (m) is the product of the dispersal ability of the species (alpha) and the species’ frequency in the metacommunity (p).

In the paper, we show that in order for parameter estimation to work reliably, we have to condition our sampling formula’s on guild size. The package also contains sampling formula’s without guild size, in order for the user to reproduce our train of thought, but here we will focus solely on situations where we have conditioned our sampling formula’s on guild size, and hence we assume that the user knows the number of individuals per guild.

So, without verder ado, let’s simulate some data and then apply the relevant sampling formula’s. We start by simulating data using the “D0” model, which assumes no differences in dispersal ability between the guilds (e.g. the guilds model then reduces to the standard neutral model)

```
set.seed(666)
<- 200
theta <- 0.1
alpha_x <- 10000
J <- 8000
J_x <- J - J_x
J_y
<- generate.Guilds.Cond(theta, alpha_x, alpha_x, J_x, J_y) abund
```

Again, we can calculate the expected distributions, and visualize our distributions. This time however, a small part of the calculation of the expected distribution can not be treated analytically, and has to be done via the means of simulation, so a larger number of replicates means a more accurate expected distribution. Here we will use a low number of replicates for computational reasons.

```
<- expected.SAD.Guilds.Conditional(theta, alpha_x, alpha_x,
abund.expected n_replicates = 5)
J_x, J_y,
par(mfrow = c(1, 2))
preston_plot(abund$guildX, abund.expected$guildX, main = "Guild X")
preston_plot(abund$guildY, abund.expected$guildY, main = "Guild Y")
```

Before, we set the size of Guild X to 8000 individuals, and the size of Guild Y to 2000 individuals. From the Preston plots we notice that indeed Guild X contains more species, but that apart from that, there are not many striking differences (we assumed equal dispersal of course!) Now that we have our data, we can apply the guilds’ sampling formula, and see whether it can accruately infer parameter values.

```
<- maxLikelihood.Guilds.Conditional(init_vals = c(theta, alpha_x),
LL model = "D0",
sadx = abund$guildX,
sady = abund$guildY,
verbose = FALSE)
$par LL
```

`## [1] 180.88797516 0.09923591`

We see that the maximum likelihood algorithm accurately infers the value of alpha, and gets relatively close to inferring the right value of theta. If we would want to infer the rate of migration, we could, in this case, also apply the standard neutral model (after all, there were no differences in dispersal, rendering the guilds construction obsolete):

```
<- maxLikelihood.ESF(init_vals = c(theta, alpha_x),
LL1 abund = c(abund$guildX, abund$guildY),
verbose = FALSE)
$par LL1
```

`## [1] 246.0613192 0.2132893`

Using the ESF, we find a slightly higher estimate for theta (but notice that the mean of the two estimates becomes really close to the value we put in), and we find a higher value for m.

So, what happens if we assume differences in dispersal? Let’s, for fun, assume that there are two guilds (guild X and Y), where guild X is much bigger, but where species from guild X have a much lower dispersal ability:

```
set.seed(666 + 42)
<- 200
theta <- 0.01
alpha_x <- 0.1
alpha_y <- 1000
J <- 800
J_x <- J - J_x
J_y
<- generate.Guilds.Cond(theta, alpha_x, alpha_y, J_x, J_y) abund
```

This generates the following dataset:

```
<- expected.SAD.Guilds.Conditional(theta, alpha_x, alpha_y,
abund.expected n_replicates = 5)
J_x, J_y,
par(mfrow = c(1, 2))
preston_plot(abund$guildX, abund.expected$guildX, main = "Guild X")
preston_plot(abund$guildY, abund.expected$guildY, main = "Guild Y")
```

Clearly, there are now strong differences in the community distributions: the larger guild has a much broader distribution, with much less singletons than the smaller guild. Again, we can use maximum likelihood to infer the parameters from the data:

```
<- maxLikelihood.Guilds.Conditional(init_vals = c(theta, alpha_x, alpha_y),
ML1 model = "D1",
sadx = abund$guildX, sady = abund$guildY,
verbose = FALSE)
```

Maximum likelihood estimation shows that we can accurately infer the underlying parameters. Furthermore, we can test whether the model assuming differences in dispersal between guilds better explains the data, than the model without. We do this by also fitting the model without differences to the data, and compare the estimated likelihoods:

```
<- maxLikelihood.Guilds.Conditional(init_vals = c(theta, alpha_x),
ML2 model = "D0",
sadx = abund$guildX, sady = abund$guildY,
verbose = FALSE)
```

To do an honest comparison, we have to correct for the extra parameter (alpha_y) in the “D1” model, we do this by calculating the AIC value:

```
= 2 * 3 - 2 * -ML1$value
AIC_D1 = 2 * 2 - 2 * -ML2$value
AIC_D0 AIC_D1
```

`## [1] 153.3964`

` AIC_D0`

`## [1] 194.3017`

The AIC of the D1 model is lower, indicating a much better fit - as expected.

This concludes the vignette. Go out there, find abundance data and check whether the data has guild structure! If not, use the ESF, but if you can come up with a potential guild structure, apply the conditional guilds sampling formula! Don’t forget to apply both the model with, and without differences in dispersal ability, in order to verify that the guilds really differ in dispersal ability!