Required data
Although in typical use, the user will have his own data available, for clarity and portability, we will generate artificial data, and use that to perform STEPCAM. Using the function generate.Artificial.Data we generate a dataset, with 3 traits, 10 plots and 10 species:
library(STEPCAM)
set.seed(42)
simul_data <- generate.Artificial.Data(n_species = 10, n_traits = 3, n_communities = 10,
occurence_distribution = 0.5, average_richness = 0.8,
sd_richness = 0.25, mechanism_random = FALSE)For the sake of this vignette we want to display the ability of STEPCAM to estimate habitat filtering. Habitat Filtering causes species with too dissimilar trait values to be selected against in a community, or in a plot. Because generate.Artificial.Data generates random trait data, for our example it is more convenient to enter specific trait values that will aid in detecting Habitat Filtering:
data_species <- simul_data$traits
data_species$trait1 <- c(runif(8,0,1), 5, 10)
data_species$trait2 <- c(runif(8,0,1), -10, 30)
data_species$trait3 <- c(runif(8,0,1), -20, 40)The three traits selected here do not represent a specific trait, but we can clearly see that the two last entries (of species9 and species10) have very deviant values from the values of the rest of the species:
data_species## species trait1 trait2 trait3
## 1 species1 0.42137842 0.14957900 0.7987603
## 2 species2 0.06763682 0.49927288 0.1103351
## 3 species3 0.56143796 0.94056488 0.5397983
## 4 species4 0.07072189 0.33423133 0.5712339
## 5 species5 0.21139194 0.18843433 0.6189515
## 6 species6 0.54962041 0.26971618 0.7148537
## 7 species7 0.48198145 0.53074408 0.1233006
## 8 species8 0.15946985 0.02145023 0.3110496
## 9 species9 5.00000000 -10.00000000 -20.0000000
## 10 species10 10.00000000 30.00000000 40.0000000To complete the construction of our artificial dataset, we have to set the abundance data to something a bit more meaningful. Our focus here is on plot 1, and for plot one we want to “trick” the dataset into showing strong Habitat Filtering. In order to do so we include all species, except species 9 and 10:
data_abundances <- simul_data$abundances
for (i in 1:8) {
data_abundances[1,i] <- 1
}
data_abundances[1,9] <- 0
data_abundances[1,10] <- 0
data_abundances## species1 species2 species3 species4 species5 species6 species7 species8
## 1 1 1 1 1 1 1 1 1
## 2 0 0 1 0 1 1 1 1
## 3 1 1 1 1 1 1 1 1
## 4 1 0 0 0 0 1 1 1
## 5 0 0 0 1 1 1 1 0
## 6 1 0 1 1 1 1 1 1
## 7 1 0 0 1 1 1 1 0
## 8 1 1 1 1 1 1 1 1
## 9 1 0 1 1 0 1 1 1
## 10 1 1 1 1 1 1 1 1
## species9 species10
## 1 0 0
## 2 1 1
## 3 1 1
## 4 1 1
## 5 0 1
## 6 1 1
## 7 0 1
## 8 1 1
## 9 0 1
## 10 1 1We now have all the basic ingredients to start a STEPCAM analysis. What STEPCAM basically does, is pick a random combination of parameters (DA, HF, LS), and simulate the stepwise removal (especially for a high relative contribution of DA, the outcome can be very stochastic). After stepwise removal, the resulting species pool is used to estimate four trait-distribution summary statistics, which in turn are compared with the values obtained for the original data. If they are similar enough, the parameter combination is accepted, and the algorithm moves on. As you will see, the algorithm goes through a range of iterations. Each iteration, the threshold of similarity is lowered (e.g. the difference between the simulated and observed data has to be smaller). As the threshold is lowered, the acceptance rate often drops, and good results are best obtained for extremely low acceptance rates, if computation time permits.
Starting STEPCAM
output <- STEPCAM_ABC(data_abundances, data_species,
numParticles = 100, n_traits = 3, plot_number = 1,
stopRate = 0.05, stop_at_iteration = 12, continue_from_file = FALSE)##
## Generating Particles for iteration 1
## 0--------25--------50--------75--------100
## ***************************************** 0.59 0.74 0.67 accept rate = 1
##
## Generating Particles for iteration 2
## 0--------25--------50--------75--------100
## ***************************************** 0.56 0.75 0.69 accept rate = 1
##
## Generating Particles for iteration 3
## 0--------25--------50--------75--------100
## ***************************************** 0.63 0.9 0.47 accept rate = 1
##
## Generating Particles for iteration 4
## 0--------25--------50--------75--------100
## ***************************************** 0.46 1.07 0.47 accept rate = 0.7769231
##
## Generating Particles for iteration 5
## 0--------25--------50--------75--------100
## ***************************************** 0.09 1.91 0 accept rate = 0.3136646
##
## Generating Particles for iteration 6
## 0--------25--------50--------75--------100
## ***************************************** 0.07 1.93 0 accept rate = 0.6778523
##
## Generating Particles for iteration 7
## 0--------25--------50--------75--------100
## ***************************************** 0.08 1.92 0 accept rate = 0.6558442
##
## Generating Particles for iteration 8
## 0--------25--------50--------75--------100
## ***************************************** 0.1 1.9 0 accept rate = 0.6392405
##
## Generating Particles for iteration 9
## 0--------25--------50--------75--------100
## ***************************************** 0.07 1.93 0 accept rate = 0.6516129
##
## Generating Particles for iteration 10
## 0--------25--------50--------75--------100
## ***************************************** 0.02 1.98 0 accept rate = 0.7112676
##
## Generating Particles for iteration 11
## 0--------25--------50--------75--------100
## ***************************************** 0.06 1.94 0 accept rate = 0.7372263
##
## Generating Particles for iteration 12
## 0--------25--------50--------75--------100
## * 50.5 50.5 50.5 accept rate = 1The results above show that the simulation ran for 9 iterations, and at the 10th iteration the acceptance rate of newly proposed parameter values dropped below 1 in 20 (5%). Right next to the progressbars, the simulation provides us with some intermediate output: the mean estimates for DS, HF and LS, and the acceptance rate. Furthermore, we see that while the acceptance rate is 1, the relative contribution of our three parameters of interest is more or less equal to the prior. What is the prior then? The prior in this case is the non-informative starting point that STEPCAM takes as parameter combination in the first iteration. In our case we are removing two species in a stepwise fashion (remember that we set the presence of two species to 0?), and hence the initial expected frequency of the three processes is 2*(1/3,1/3,1/3), which is about (0.67,0.67,0.67). From iteration 7 onwards we see that the acceptance rate starts to drop, and that the relative contribution of the three processes starts to deviate from an equal distribution. By iteration 11 we have a very nice result of c(0.02,1.98,0), which is very much in line with our expectation of (0,2,0) (DS,HF,LS)! STEPCAM has thus correctly inferred the relative contribution of the three processes.
We can plot our results in a ternary plot:
TernPlot(output)
Because we have only 2 species to remove, only a limited number of combinations of (DS, HF, LS) are possible, hence the low number of points in the ternaryplot. For a more exciting example of the plotting abilities of the STEPCAM package, please see the relevant vignette.