David Zelený

Script for creating simulated community data, structured along two virtual ecological gradients

simul2.r
```## Simulation of ecological data structured by two ecological gradients
## (based on the script of Jason Fridley (Fridley et al. 2007, Appendix S2), modified by David Zelený ([email protected])

# Basic setting:
totS<-300	#total species in simulation
Np<-500		#number of sample plots
gr1.length <- 5000 # length of the first gradient (arbitrary units)
gr2.length <- 2000 # length of the second gradient (arbitrary units)

# in case you want to add random species into the matrix, change following:
no.of.random.sp <- 0  # no of randomly added species to matrix
max.fr.random <- 0  # maximum frequency of randomly added species
# (random species are added with randomly generated frequency <0; max.fr.random> into the matrix

#Creates coenoclines
#Creates survey designs (random, uniform, or biased exponentially along gradient)
#Creates plot-spp matrix based on survey design
#Uses matrix in GS algorithm
#Outputs graphs
#Loops through "random" and "exponential (biased)" survey designs

#Choose niche distribution type: "normal", "skewed", or "random"
niche.type<-"random"

#This is beta function for generating niches
curve <- function(Ao,m,r,a,g,x)	{
(Ao*((((x-m)/r)+(a/(a+g)))^a)*((1-(((x-m)/r)+(a/(a+g))))^g))/(((a/(a+g))^a)*((1-(a/(a+g)))^g))
}
x1 <- seq(1, gr1.length, by = 1) # gradient 1
x2 <- seq(1, gr2.length, by = 1) # gradient 2
# species values for random niches
if(niche.type=="random") {
S1<-totS					#number of species
Ao1<-rlnorm(S1,2,1)			#amplitude vector (lognormal distribution)
m1<-sample(seq(5:max(x1)-5),S1, replace = T)		#location of optima
a1<-(runif(S1,min=.1,max=4))		#shape parameter (alpha)
g1<-(runif(S1,min=.1,max=4))		#shape parameter (gamma)

S2<-totS					#number of species
Ao2<-rlnorm(S2,2,1)			#amplitude vector (lognormal distribution)
m2<-sample(seq(5:max(x2)-5),S2, replace = T)		#location of optima
a2<-(runif(S2,min=.1,max=4))		#shape parameter (alpha)
g2<-(runif(S2,min=.1,max=4))		#shape parameter (gamma)
}

# species values for normal niches
if(niche.type=="normal") {
S1<-totS					#number of species
Ao1<-rlnorm(S1,2,1)			#amplitude vector (lognormal distribution)
m1<-sample(seq(5:max(x1)-5),S1, replace = T)		#location of optima
a1<-rep(1.99,S1)				#shape parameter (alpha)
g1<-rep(1.99,S1)				#shape parameter (gamma)

S2<-totS					#number of species
Ao2<-rlnorm(S2,2,1)			#amplitude vector (lognormal distribution)
m2<-sample(seq(5:max(x2)-5),S2, replace = T)		#location of optima
a2<-rep(1.99,S2)				#shape parameter (alpha)
g2<-rep(1.99,S2)				#shape parameter (gamma)
}

# species values for skewed niches
if(niche.type=="skewed") {
S1<-totS					#number of species
Ao1<-rlnorm(S1,2,1)			#amplitude vector (lognormal distribution)
m1<-sample(seq(5:max(x1)-5),S1, replace = T)		#location of optima
#produce skews in either direction, randomly
a.1<-rep(1.99,S1)
g.1<-rep(.25,S1)
samp<-sample(c(1:S1),S1/2,replace=FALSE)
a.1[samp]<-.25
g.1[samp]<-1.99
a1<-a.1				#shape parameter (alpha)
g1<-g.1				#shape parameter (gamma)

S2<-totS					#number of species
Ao2<-rlnorm(S2,2,1)			#amplitude vector (lognormal distribution)
m2<-sample(seq(5:max(x2)-5),S2, replace = T)		#location of optima
#produce skews in either direction, randomly
a.1<-rep(1.99,S2)
g.1<-rep(.25,S2)
samp<-sample(c(1:S2),S2/2,replace=FALSE)
a.1[samp]<-.25
g.1[samp]<-1.99
a2<-a.1				#shape parameter (alpha)
g2<-g.1				#shape parameter (gamma)
}

#------Summary plotting:
#par(mfrow=c(2,2))

#	for(L in 2:S) {
#		lines(x,A[,L])
#	}
#lines(x,A[,51],lwd=3,col=2)
#lines(x,A[,52],lwd=3,col=3)

#plot(x,richness)

#plot(c(1:S),rev(sort(ranges)))

##SAMPLING FROM GENERATED CURVES
#---------------------------------

#sppmat <- data.frame(c(1:S),paste("S",c(1:S),sep=""))

#eq.sample.x <- trunc(seq(min(x),max(x)-1,length=Np))		#vector of sample locations along gradient

rand.sample.x1 <- trunc(sample(c(2:max(x1))-1,Np, replace = T))
rand.sample.x2 <- trunc(sample(c(2:max(x2))-1,Np, replace = T))

#Biased samples at one end of gradient
#based on exponential function; more samples at low end of gradient
#	exp.x<-sort(rexp(Np,rate=20))
#	exp.sample.x<-trunc(	exp.x/(max(exp.x))*(max(x)-1))
#remove duplicate Np x values
#		while( length(unique(exp.sample.x))!=Np ) {
#		exp.sample.x[duplicated(exp.sample.x)]<-exp.sample.x[duplicated(exp.sample.x)]+1
#	}
#	exp.sample.x<-rev(max(x)-exp.sample.x)	#switch to other side of gradient

A1 <- matrix(0,nrow=length (rand.sample.x1),ncol=S1) #response abundances
for(L in 1:S1){
A1[,L] <- curve(Ao1[L],m1[L],r1[L],a1[L],g1[L], rand.sample.x1)
}
A2 <- matrix(0,nrow=length (rand.sample.x2),ncol=S2) #response abundances
for(L in 1:S2){
A2[,L] <- curve(Ao2[L],m2[L],r2[L],a2[L],g2[L], rand.sample.x2)
}
A12 <- A1*A2

#Summary stats for each species
p.mat<-A12
p.mat.NA<-is.na(p.mat)
p.mat[p.mat.NA]<-0
a.mat<-p.mat			#abundance matrix with zeros
p.mat[p.mat>0]<-1		#presence-absence matrix

#ranges<-colSums(p.mat)		#total ranges (niche breadth) for each species

#Distribution of # inds per local community along gradient
#from beta distribution (see above)
density.curve.1<-curve(100,max(x1)/2,max(x1),.25,.25,x1)
density.curve.2<-curve(100,max(x2)/2,max(x2),.25,.25,x2)

#Number of inds in each plot for each above sampling schemes
#draws.eq <- round(rnorm(Np,mean=density.curve[eq.sample.x],sd=1))
draws.rand.1 <- round(rnorm(Np,mean=density.curve.1[rand.sample.x1],sd=1))
draws.rand.2 <- round(rnorm(Np,mean=density.curve.2[rand.sample.x2],sd=1))
draws.rand <- round ((draws.rand.1 + draws.rand.2)/2)
#draws.exp <- round(rnorm(Np,mean=density.curve[exp.sample.x],sd=1))

#output data frames
#samp.out.eq <- matrix(0,nrow=Np,ncol=S)
samp.out.rand <- matrix(0,nrow=Np,ncol=totS)
#samp.out.exp <- matrix(0,nrow=Np,ncol=S)

#Sampling for equal-sample-interval
#for(i in 1:Np) {
#	samp.prob<-a.mat[eq.sample.x[i],]		#probabilities of sampling each species in given location (based on rel abundance)
#	tab.samp <- table(sample(c(1:S),size=draws.eq[i],prob=samp.prob,replace=T))	#tabulated vector of spp identities after choosing "draws" no. of individuals
#	samp.out.eq[i,][as.numeric(names(tab.samp))] <- tab.samp
#}
#eq.pa <- samp.out.eq
#eq.pa[eq.pa>0] <- 1		#presence-absence version
#converted to 2-column list:
#spp.vec<-NULL
#plot.vec<-NULL
#for(i in 1:Np) {
#	vec.true<-as.logical(eq.pa[i,])
#	plot.vec<-c(plot.vec,rep(i,length=sum(eq.pa[i,])))
#	spp.vec<-c(spp.vec,c(1:S)[vec.true])
#}
#out.eq <- data.frame(plot.vec,spp.vec)

#Sampling for random-sample-interval
for(i in 1:Np) {
samp.prob<-a.mat[i,]		#probabilities of sampling each species in given location (based on rel abundance)
tab.samp <- table(sample(c(1:totS),size=draws.rand[i],prob=samp.prob,replace=T))	#tabulated vector of spp identities after choosing "draws" no. of individuals
samp.out.rand[i,][as.numeric(names(tab.samp))] <- tab.samp
}

rand.pa <- samp.out.rand
colnames (rand.pa) <- paste ('spec_', 1:dim (rand.pa)[2], sep = '')
rand.pa[rand.pa>0] <- 1		#presence-absence version

if (no.of.random.sp > 0)
{
random.part <- matrix (0, ncol = no.of.random.sp, nrow = Np, dimnames = list (NULL, paste ('rand_', 1:no.of.random.sp, sep = '')))
for (i in seq (1, no.of.random.sp))  random.part[sample (c (1:Np), round (runif (1)*max.fr.random)),i] <- 1
simul.table <- cbind (rand.pa, random.part)
} else simul.table <- rand.pa

######################################################

## RESULTS
# m = vector of species optimas
# rand.pa = sample-species matrix
# rand.sample.x = position of sample locations
# simul.table = simulated data (if random species were added, they are appended in the end of the table)```