VirtualCom2014 Files

Manual for running examples with VirtualCom in R

First steps

Download the VirtualCom.zip file and unpack it into your preferred working directory.

Open an R consol. And then load the required add-on packages.

pkgs <- c('geiger', 'picante', 'apTreeshape','OUwie','FD', 'Hmisc')
for(pkg in pkgs) library(pkg, character.only = TRUE)

Load all functions and datasets of VirtualCom.
my.path <- "~/myWorkingDirectory"
for(fun in list.files(my.path)) source(file.path(my.path, fun))

Each simulation of VirtualCom starts with the generation of a species pool through a trait evolution process on a simulated phylogeny (function: create.pool). In a second step and over several years, native communities assemble from the species pool of natives according to the chosen assembly rules. Then invasions can be added by allowing invasive species from the species pool to enter communities for some additional years (function: tamaure). These two steps can be either performed independently, or a simulation experiment can be run over all steps at once (function: simulation.experiment). One can get a first idea on how these functions work by running simple examples.

Simple examples for the function create.pool:

# create.pool with trait evolution based on Brownian motion
pool_BM <- create.pool(n.species.pool=500, n.invader.pool=1, evol.model="BM", min.phyl.signal=NA, evol.model.param=0.1, nrep=1)
hist(pool_BM$func$niche_evol)

# create.pool with trait evolution based on an OU-model with 3 different optima
pool_OU <- create.pool(n.species.pool=500, n.invader.pool=1, evol.model="OUwie.sim", min.phyl.signal=NA, evol.model.param=0.1, nrep=1, OUwie.sigma.sq=rep(0.00001,3), OUwie.theta=c(1, 50, 100), OUwie.alpha = rep(10,3))
hist(pool_OU$func$niche_evol)

# plot the simulated phylogenetic trees
plot(pool_BM$phy)
plot(pool_OU$phy)

# calculate the resulting phylogenetic signals
phylosignal(pool_BM$func$niche_evol, pool_BM$phy, checkdata=FALSE)
phylosignal(pool_OU$func$niche_evol, pool_BM$phy, checkdata=FALSE)



Simple examples for the function tamaure:

# extract the simulated traits from the species pool simulated with create.pool 
niche.optima <- pool_BM$func$niche_evol
names(niche.optima) <- pool_BM$func$SpeciesID
# run community assembly with only an environmental filter, with only a symmetric competition filter or with only a hierarchical competition filter and plot results
env.filter <- tamaure(niche.breadth=2, niche.optima=niche.optima, env=50, beta.env=1, beta.comp=0, beta.abun=10, years=10, K=100, plot=TRUE)
sym.comp.filter <- tamaure(niche.breadth=2, niche.optima=niche.optima, env=50, beta.env=0, beta.comp=10, beta.abun=10, years=10, K=100, plot=TRUE, competition="symmetric", intra.sp.com=1)
hier.comp.filter <- tamaure(niche.breadth=2, niche.optima=niche.optima, env=50, beta.env=0, beta.comp=0.5, beta.abun=10, years=10, K=100, plot=TRUE, competition="hierarchical", intra.sp.com=1)

If the argument plot=TRUE in tamaure then community composition is plotted, together with the strengths of the different filters and the evolution of the functional mean pairwise distances (mpd) over time. Species IDs on the x-axes are sorted by trait values. These plots help to understand how assembly processes influence final community composition. For example, when only environmental filtering drives community assembly, then species are very similar (clustering pattern); when only symmetric competition drives community assembly, then species are very dissimilar (over-dispersion pattern); when only hierarchic competition drives community assembly, then species tend to have high trait values and are very similar (clustering pattern).

To answer a scientific question, one often wants to run repetitions of scenarios and compare different combinations of parameters (Please note that these simulations can take a few minutes). To run such a virtual experiment we implemented an additional function: 
Simple examples for the function simulation.experiment:
# look at the parameter table created with the file make.parameter.file.R
print(param.virtual.exp)

# run the first line of this parameter table, i.e. the first scenario
single.run <- simulation.experiment(param.virtual.exp[1,]) # run first line
str(single.run)

# create a wrapper function to run the full experiment 
wrapper <- function(a, my.path){
    pkgs <- c('geiger', 'picante', 'apTreeshape','OUwie','FD', 'Hmisc')
    for(pkg in pkgs) library(pkg, character.only = TRUE)
    for (fun in list.files(my.path)) source(file.path(my.path, fun))
	   print(a)
	   return(try(simulation.experiment(param.virtual.exp[a,])))	
}  
   
# running the experiment either with lapply (this can take a few minutes)
output <- lapply(1:4, wrapper, my.path=my.path) 

# or running the experiment with sfLapply using snowfall to allow for parallel computing (this can take a few minutes)
library(snowfall)
sfInit(parallel=TRUE, cpus=8)
output <- sfLapply(1:4, wrapper, my.path=my.path)
sfStop()

# Analysis of results:
# 1. extract results from list
result.table <- get.results(output=output, myvar="obs", invader="FALSE")

# 2. define a variable that contains the information of the underlying assembly process   
result.table$process <- ifelse(result.table$beta.env==0 & result.table$beta.comp==0, "Random", ifelse(result.table$beta.env!=0 & result.table$beta.comp==0, "Env", ifelse(result.table$beta.env==0 & result.table$beta.comp!=0, "Comp", "Both") ))

# 3. plot the funtional diversity (mpd) of communities in dependence on assembly processes
library(ggplot2)
ggplot(data=result.table, aes(x=process, y=FD_ab_mpd)) + geom_boxplot(width=0.8) +  xlab("Assembly rules") +  ylab("Functional diversity")


Example 1: Effects of trait evolution on diversity patterns
In this example we asked whether phylogenetic diversity patterns (mpdphy), simulated from species pools with high phylogenetic signal or moderate phylogenetic signal, equally reveal true community assembly rules. To answer this question we compared mpdphy patterns when the phylogeny was or was not transformed via a delta-model (resulting in high vs. moderate phylogenetic signal) under 3 assembly scenarios: random, environmental filtering or competition. The parameters for these different scenarios can be loaded from the package. Each row of this data frame contains all the parameter values required for one simulation scenario.

# look at the parameter table created with the file make.param1.R
print(param1)

To run all the parameter combinations we define a function called wrapper_ex1.

wrapper_ex1 <- function(a, my.path){
    pkgs <- c('geiger', 'picante', 'apTreeshape','OUwie','FD', 'Hmisc')
    for(pkg in pkgs) library(pkg, character.only = TRUE)    
    for (fun in list.files(my.path)) source(file.path(my.path, fun))
    return(try(simulation.experiment(param1[a,])))	
}  

The simulation experiment can then easily be run with sfLapply{snowfall}.

library(snowfall)
sfInit(parallel=TRUE, cpus=8)
output1 <- sfLapply(1:8, wrapper_ex1, my.path)
sfStop()

Results can then be collected with get.results and be plotted to analyze the results of the experiment.

### Analysis of results:
# 1. extract results from list
result.table1 <- get.results(output=output1, myvar="obs", invader="FALSE")
# define a variable that contains the information of whether or not trees were transformed   

# 2. define variables that contain the information of the tree transformation and the underlying assembly process   
result.table1$evol.model.param <- ifelse(result.table1$evol.model.param==1, "no", "yes")
result.table1$process <- ifelse(result.table1$beta.env==0 & result.table1$beta.comp==0, "Random", ifelse(result.table1$beta.env!=0 & result.table1$beta.comp==0, "Env", ifelse(result.table1$beta.env==0 & result.table1$beta.comp!=0, "Comp", "Both") ))

# 3. remove scenarios with interacting competition and environmental filters   
result.table1 <- result.table1[result.table1$process != "Both",]
head(result.table1)

# 4. calculate phylogenetic signal for the two different species pools   
(Ks <- aggregate(result.table1$Blom_K, list(result.table1$evol.model.param), mean))

# 5. plot the results
library(ggplot2)
ggplot(data=result.table1, aes(x=process, y=FD_ab_mpd))+ geom_boxplot(width=0.8) + xlab("Assembly rules") + ylab("Functional diversity")
	
ggplot(data=result.table1, aes(x=process, y=PD_ab_mpd)) +	geom_boxplot(aes(fill= evol.model.param), width=0.8) + xlab("Assembly rules") + ylab("Phylogenetic diversity") + theme(legend.justification=c(1,0), legend.position=c(0.42,0)) + scale_fill_discrete(name = "Tree transformed?", labels=c(paste("no; K~", round(Ks[1,2],2),sep=""), paste("yes; K~", round(Ks[2,2],2),sep=""))) 


Example 2: Different indices for revealing invasion processes
In this example, we compared the performance of different indices (Thuiller et al. 2010)at detecting invasion processes. To do so, we assembled native communities either under an environmental-filtering scenario or under a competition assembly scenario, and then invaded these communities following the same assembly processes. The parameters for the different scenarios (environmental filtering vs. competition) can be loaded from the package. Each row of this data frame contains all the parameter values required for one simulation scenario.

# look at the parameter table created with the file make.param2.R
print(param2)

To run all the parameter combinations we define a function called wrapper.

wrapper_ex2 <- function(a, my.path){
    pkgs <- c('geiger', 'picante', 'apTreeshape','OUwie','FD', 'Hmisc')
    for(pkg in pkgs) library(pkg, character.only = TRUE)
    for (fun in list.files(my.path)) source(file.path(my.path, fun))
    return(try(simulation.experiment(param2[a,])))	
}  

The simulation experiment can then easily be run with sfLapply{snowfall}.

library(snowfall)
sfInit(parallel=TRUE, cpus=9)
output2 <- sfLapply(1:9, wrapper_ex2, my.path)
sfStop()

Results can then be collected with get.results and be plotted to analyze the results of the experiment.

### Analysis of results:
# 1. extract results from list
result.comm <- get.results(output2, myvar="rankNULL", invader="FALSE")
# 2. define a variable that contains the information of the underlying assembly process   
result.comm$process <- ifelse(result.comm$beta.env==0 & result.comm$beta.comp==0, "Random", ifelse(result.comm$beta.env!=0 & result.comm$beta.comp==0, "Env", ifelse(result.comm$beta.env==0 & result.comm$beta.comp!=0, "Comp", "Both") ))
# 3. keep only results for a pure competition or environmental filter   
result.comm <- result.comm[(result.comm$process == "Comp") | (result.comm$process == "Env"),]

# 4. get invader-based indices
result.inv <- get.results(output2, myvar="rankNULL", invader="TRUE")
result.inv <- result.inv[! is.na(result.inv$fun_MDNC),]		
# 5. define a variable that contains the information of the underlying assembly process   
result.inv$process <- ifelse(result.inv$beta.env==0 & result.inv$beta.comp==0, "Random", ifelse(result.inv$beta.env!=0 & result.inv$beta.comp==0, "Env", ifelse(result.inv$beta.env==0 & result.inv$beta.comp!=0, "Comp", "Both") ))
# 6. keep only results for a pure competition or environmental filter   
result.inv <- result.inv[(result.inv$process == "Comp") | (result.inv$process == "Env"),]

# 7. combine results for community-based and invader-based indices
library(reshape2)
long.comm <- melt(result.comm[, c("process","FD_pa_mpd", "FD_pa_mntd")])
long.inv <- melt(result.inv[, c("process","fun_MDNC", "fun_MDNN")])
long <- cbind(rbind(long.comm, long.inv), Indices=c(rep("comm", nrow(long.comm)), rep("inv", nrow(long.inv))))
long$variable <- ifelse(long$variable == "FD_pa_mpd", "mpd", ifelse(long$variable == "FD_pa_mntd", "mntd", ifelse(long$variable == "fun_MDNC", "mdnc", "mdnn")))

# 8. plot the results to analyze the outcome of the experiment
library(ggplot2)
ggplot(long[long$process=="Env",], aes(x=variable, y=value)) + geom_boxplot(aes(fill= Indices),width=0.8) + xlab("Indices") + ylab("Rank") + geom_hline(aes(yintercept=0.05), linetype="dashed") + guides(fill=FALSE) 

ggplot(long[long$process=="Comp",], aes(x=variable, y=value)) + geom_boxplot(aes(fill= Indices), width=0.8) + xlab("Indices") +  ylab("Rank") + theme(legend.justification=c(1,0), legend.position=c(0.2,0.1)) + geom_hline(aes(yintercept=c(0.95,0.05)), linetype="dashed") 



1