#SIMULATION OF DECREASING TEMPORAL RESOLUTION IN TIME SERIES OF CONSECUTIVE SAMPLES
#ABUNDANCES OF DEAD INDIVIDUALS ARE TRACKED IN SIMULATIONS
#THREE MODELS OF LOCAL COMMUNITY DYNAMICS-NEUTRAL, NEGATIVE CONSPECIFIC DENSITY DEPENDENCE AND TRADE-OFF MODEL (WITHOUT DENSITY DEPENDENCE)
#ONE NEUTRAL LOG-SERIES-LIKE METACOMMUNITY
#EACH MODEL ENDS WITH A TWO-DIMENSIONAL ARRAY "MATRIXOFDEADS" OF DEAD INDIVIDUALS (SPECIES IDENTITY TRACKED BY THEIR NUMBER)
#THAT DIED OVER N YEARS OF TIME SERIES (ROWS) IN N SIMULATIONS (COLUMNS)
#THIS ARRAY IS THEN USED FOR AVERAGING PROCEDURE WHERE INDIVIDUALS ARE POOLED, WITH TWO OPTIONS - DEAD INDIVIDUALS ARE EITHER
#SAMPLED FROM UNIFORM (PER CAPITA PRESERVATION RATE = 1) OR EXPONENTIAL AGE-FREQUENCY DISTRIBUTIONS (PRESERVATION RATE < 1)
#SAMPLE SIZE OF AVERAGED ASSEMBLAGES IS FIXED WITH INCREASING TIME AVERAGING
#SAMPLING RATE IS DEFINED BY THE SAMPLE SIZE OF AVERAGED ASSEMBLAGES RELATIVE TO THE TOTAL NUMBER OF INDIVIDUALS DIED OVER
#THE DURATION OF TIME AVERAGING
#SOURCE: TOMASOVYCH A. AND KIDWELL, S.M. 2010. THE EFFECTS OF TEMPORAL RESOLUTION ON SPECIES TURNOVER AND ON TESTING METACOMMUNITY MODELS. AMERICAN NATURALIST.
#LIBRARIES
library(vegan)
library(untb)
#SOME ARBITRARY INPUTS
timetoequilibrium<-200 #TIME TO ATTAIN STEADY STATE - "disequilibrium" years do not go into output
abundances<-c(28435,15050,10991,8227,6368,5180,4250,3438,2872,2402,1961,1702,
1407,1212,962,802,658,566,488,436,360,305,271,229,199,171,148,134,116,95,80,69,
60,50,43,37,33,29,25,22,19,16,13,11,9,8,7,6,5,4,3,3,2,2,2,2,1,1,1,1,1) #SPECIES ABUNDANCES IN METACOMMUNITY
REPLACE<-30 #NUMBER OF DIED/BORN INDIVIDUALS PER TIME STEP (DEPENDS ON LIFESPAN)
m<-0.1 #PER CAPITA IMMIGRATION PROBABILITY FROM METACOMMUNITY PER TIME STEP
Nlive<-50 #NUMBER OF INDIVIDUALS IN NON-AVERAGED LOCAL COMMUNITY-SAMPLE SIZE
Ndead<-100 #NUMBER OF INDIVIDUALS IN AVERAGED LOCAL COMMUNITY-SAMPLE SIZE
simulations<-50 #NUMBER OF SIMULATIONS
TAvalues<-c(5,10,50,100) #SELECTED TIME AVERAGING VALUES
timewindow<-300 #DURATION OF TIME SERIES IN TIME STEPS
events<-9 #NUMBER OF CONSECUTIVE SAMPLES IN ONE TIME SERIES
time<-timewindow+timetoequilibrium+max(TAvalues) #TOTAL DURATION OF TIME STEPS IN SIMULATIONS-MUST EXCEED timewindow+timetoequilibrium+max(TAvalues)
a<-0.4 #STRENGTH OF DENSITY DEPENDENCE FOR DENSITY-DEPENDENT MODEL
lambda<-0.1 #PER CAPITA DESTRUCTION RATE FOR SAMPLING FROM EXPONENTIAL AFD
#############################################################################################################################################
#INITIAL PARAMETERS
REPLACE<-REPLACE
m<-m
timetoequilibrium<-timetoequilibrium
abundances<-abundances
Nlive<-Nlive
Ndead<-Ndead
simulations<-simulations
TAvalues<-TAvalues
time<-time
timewindow<-timewindow
events<-events
numberofspecies<-length(abundances)
#BUILD INDIVIDUALS IN METACOMMUNITY USING ABUNDANCE VECTOR (OBSERVED OR PRODUCED BY MODELING)
metacommunity<-rep(0,0);
for (i in 1:length(abundances)) {
v<-rep(i, times=abundances[i])
metacommunity<-append(metacommunity, v, after = length(metacommunity))
}
local<-rep(NA, Nlive) #vector of species in local living (non-averaged) community
matrixofdeads<-array(0,dim=c((REPLACE*time)-(timetoequilibrium*REPLACE), simulations))
#ARRAY WITH VECTORS CONTAINING INDIVIDUALS THAT DIED DURING N TIME STEPS
#(MINUS N OF INITIAL DISEQUILIBRIUM YEARS)
#IS DETERMINED BY THE NUMBER OF INDIVIDUALS DIED DURING ONE STEP
#TIMES NUMBER OF TIME STEPS
deadinputs<-array(0,dim=c(REPLACE, time,simulations)) #inputs of dead species for N time steps
liveinputs<-array(0, dim=c(Nlive, time)) #inputs of live species for N time steps
##########################################################################################################################################
#NEUTRAL MODEL
#LOCAL COMMUNITY DRIFT - DYING INDIVIDUALS CAN CONTRIBUTE TO OFFSPRING IN THE NEXT GENERATION - MORAN MODEL
for (j in 1:simulations) { #repeat local dynamics for N simulations
for (i in 1:Nlive) {
local[i]<-sample(metacommunity, 1, replace=T)
} #initialize live assemblage from metacommunity
simulationdeadvector<-rep(0,0) #vector for died individuals for one simulation
for (s in 1:time) { #loop for N time steps
yeardeadvector<-rep(0,0) #vector for died individuals for one time step
for (t in 1:REPLACE) { #simulate local dynamics within one year, for the number of replacements
random<-runif(1,0,1)
#DENSITY-INDEPENDENT MORTALITY
deadind<-sample(Nlive, 1, replace=F) #density-independent mortality
deads<-local[deadind] #identify species identity of dead individuals
#LOOP FOR MIGRATION
if (random<m)
{local[deadind]<-sample(metacommunity, 1, replace=T)
} #immigration from metacommunity
else {local[deadind]<-sample(local, 1, replace=T)
} #recruitment from local live assemblage
yeardeadvector<-append(yeardeadvector, deads, after = length(yeardeadvector)) #build vector of dead inds
}
liveinputs[,s]<-local #yearly inputs of live species for N years
deadinputs[1:(length(yeardeadvector)),s,j]<-yeardeadvector #yearly inputs of dead species for N years
#VECTOR OF ALL INDIVIDUALS THAT DIED DURING N TIME STEPS
simulationdeadvector<-append(simulationdeadvector, yeardeadvector, after = length(simulationdeadvector))
}
matrixofdeads[1:length(simulationdeadvector[((timetoequilibrium*REPLACE)+1):(time*REPLACE)]),j]<-simulationdeadvector[((timetoequilibrium*REPLACE)+1):(time*REPLACE)]
}
##############################################################################################################################################
#TRADE-OFF MODEL
deads<-rep(NA,REPLACE) #vector of individuals died during one time step
deadind<-rep(NA, 1);dispind<-rep(NA, 1);drawnind<-rep(NA,1);drawnind2<-rep(NA, 1); #bookkeeping identity of dying and immigrating/born individuals
mortalities<-rep(NA, numberofspecies);disp<-rep(NA, numberofspecies);
mortalitymatrix<-matrix(0,nrow=numberofspecies, ncol=simulations)
for (j in 1:simulations) { #simulate local dynamics
for (i in 1:numberofspecies) {
mortalities[i]<-runif(1,0,1) #assign predefined species-specific mortalities to individual species according to tradeoff
disp[i]<-mortalities[i] #species-specific dispersal rates are inversely related to mortalities
}
for (i in 1:Nlive) {
local[i]<-sample(metacommunity, 1, replace=T) #make initial local live sample from metacommunity
}
mortalitymatrix[,j]<-mortalities
simulationdeadvector<-rep(0,0)
for (t in 1:time) {
for (k in 1:REPLACE) {
repeat {
deadind<-sample(Nlive, 1, replace=F) #randomly draw deads from sample with n individuals(T,F)
drawnind<-local[deadind]
random<-runif(1,0,1)
if (random<mortalities[drawnind]) #if mortality of individual[species] is true
{deads[k]<-drawnind}
if (random<mortalities[drawnind])
{break}
}
if (random<m) {
repeat {
dispind<-sample(length(metacommunity), 1,replace=F) #randomly draw deads from sample with n individuals
drawnind2<-metacommunity[dispind]
if (random<disp[drawnind2]) #if mortality of individual[species] is true
{local[deadind]<-drawnind2}
if (random<disp[drawnind2])
{break}
}
}
else {repeat
{
dispind<-sample(length(local), 1, replace=F) #randomly draw deads from sample with n individuals
drawnind2<-local[dispind]
if (random<disp[drawnind2]) #if mortality of individual[species] is true
{local[deadind]<-drawnind2}
if (random<disp[drawnind2])
{break}
}
}
}
deads<-deads[!is.na(deads)]
deadinputs[1:(round(REPLACE)),t,j]<-deads
simulationdeadvector<-append(simulationdeadvector, deadinputs[,t,j], after = length(simulationdeadvector))
}
matrixofdeads[1:length(simulationdeadvector[((timetoequilibrium*REPLACE)+1):(time*REPLACE)]),j]<-simulationdeadvector[((timetoequilibrium*REPLACE)+1):(time*REPLACE)]
}
###########################################################################################################################################
#NEGATIVE INTRASPECIFIC DENSITY DEPENDENCE
for (j in 1:simulations) { #repeat local dynamics for N simulations
for (i in 1:Nlive) {
local[i]<-sample(metacommunity, 1, replace=T)
} #initialize live assemblage from metacommunity
simulationdeadvector<-rep(0,0) #vector for died individuals for one simulation
for (s in 1:time) { #loop for N time steps (years)
yeardeadvector<-rep(0,0) #vector for died individuals for one time step
repeat { #simulate local dynamics within one time step, for the number of replacements
random<-runif(1,0,1)
#DENSITY-DEPENDENT MORTALITY
deadind<-sample(Nlive, 1, replace=F) #initially random density-independent mortality
deads<-local[deadind] #identify species identity of dead individuals
proportion<-species.table(local)[deads]/Nlive
a<-a #strength of density dependence - MUST BE SPECIFIED
w<-1-a*proportion #survivorship-probability of to escape density-dep. mortality
migrandom<-runif(1,0,1)
if (w < random) { #is density-dependent survivorship for the given individual small?
yeardeadvector<-append(yeardeadvector, deads, after = length(yeardeadvector)) #build vector of dead inds
#LOOP FOR MIGRATION
if (migrandom<m)
{local[deadind]<-sample(metacommunity, 1, replace=T)
} #immigration from metacommunity
else {local[deadind]<-sample(local, 1, replace=T)
} #recruitment from local live assemblage
}
if (length(yeardeadvector) == REPLACE) {break}
}
liveinputs[,s]<-local #yearly inputs of live species for N years
deadinputs[1:(length(yeardeadvector)),s]<-yeardeadvector #yearly inputs of dead species for N years
#VECTOR OF ALL INDIVIDUALS THAT DIED DURING N YEARS
simulationdeadvector<-append(simulationdeadvector, yeardeadvector, after = length(simulationdeadvector))
}
matrixofdeads[1:length(simulationdeadvector[((timetoequilibrium*REPLACE)+1):(time*REPLACE)]),j]<-simulationdeadvector[((timetoequilibrium*REPLACE)+1):(time*REPLACE)]
}
#############################################################################################################################################
#SIMULATING TIME AVERAGING
uniform=TRUE #SAMPLING FROM UNIFORM AGE-FREQUENCY DISTRIBUTION OR NOT?
frame<-array(NA, dim=c(length(TAvalues),9)) #OUTPUT TABLE FOR PLOTTING - NINE VARIABLES VERSUS TIME AVERAGING
timeseries<-seq(length=events,from=1,to=timewindow) #selecting time steps that correspond to N equally-spaced events
for (x in 1:length(TAvalues)) {
TA<-TAvalues[x] #max TA is given by (REPLACE*TA)+(timewindow*REPLACE)>>time*REPLACE
TADEADINPUT<-(REPLACE*TA) #number of individuals died during time averaging
#IF TIME AVERAGING DOES NOT ALLOW ACCUMULATION OF NDEAD INDIVIDUALS, REDUCE THE NDEAD TO THE MAXIMUM POSSIBLE N
if (Ndead > TADEADINPUT) {N <- TADEADINPUT}
if (Ndead == TADEADINPUT) {N <- Ndead}
if (Ndead < TADEADINPUT) {N <- Ndead}
#SAMPLE PRESERVED DEAD INDIVIDUALS UNIFORMLY FROM THE COMPLETE POOL OF DEAD INDIVIDUALS
preserveddeadlist<-array(NA, dim=c(N, timewindow,simulations))
absmatrix<-array(0,dim=c(numberofspecies, timewindow,simulations))
percentmatrix<-array(0,dim=c(numberofspecies, timewindow,simulations))
for (i in 1:simulations) {
for (j in 1:timewindow) {
deadlist<-matrixofdeads[(j*REPLACE):(((j*REPLACE)+TADEADINPUT)-1),i] #take dead individuals during one interval after attaining equilibrium
#deadlist - list of died individuals in order they died, within the range of TA
#SAMPLING FROM UNIFORM AGE-FREQUENCY DISTRIBUTION
if (uniform==TRUE) {preserveddeadlist[,j,i]<-sample(deadlist,N,replace=F,prob=NULL)} #individuals of preserved species in matrix of local samples one column - one time step
#SAMPLING FROM EXPONENTIAL AGE-FREQUENCY DISTRIBUTION, DESTRUCTION RATE LAMBDA MUST BE SPECIFIED
else {
ND<-rep(NA,time)
lambda<-lambda
for (xx in time:1) {ND[xx]=1*exp(-lambda*xx)}
ND<-ND[1:TA] #SELECT ALL COHORTS WITHIN THE GIVEN RANGE OF TA
P<-ND/sum(ND) #PROPORTION OF PRESERVED INDIVIDUALS IN EACH COHORT
#RELATIVE TO ALL PRESERVED SHELLS
Pvector<-rep(P,each=REPLACE)
preserveddeadlist[,j,i]<-sample(deadlist,N,replace=F,prob=Pvector)
}
}
}
#SPECIES ABUNDANCES IN DEATH ASSEMBLAGES
for (i in 1:simulations) {
for (j in 1:timewindow) {
absmatrix[1:length(species.table(preserveddeadlist[1:N,j,i])),j,i]<-species.table(preserveddeadlist[1:N,j,i])
}
}
for (i in 1:simulations) {
for (j in 1:timewindow) {
percentmatrix[,j,i]<-abstimematrix[,j,i]/sum(abstimematrix[,j,i])
}
}
absolute<-absmatrix[,seq(length=events,from=timewindow,to=1),]
relative<-percentmatrix[,seq(length=events,from=timewindow,to=1),]
medianbc<-rep(NA, simulations); medianhorn<-rep(NA, simulations); medianjaccard<-rep(NA, simulations);
for (i in 1:simulations) {
medianbc[i]<-median(1-vegdist(t(relative[,,i]), method="bray"))
medianhorn[i]<-median(1-vegdist(t(relative[,,i]), method="horn"))
medianjaccard[i]<-median(1-vegdist(t(relative[,,i]), method="jaccard"))
}
frame[x,]<-as.vector(c(mean(medianbc), quantile(medianbc, c(0.025)), quantile(medianbc, c(0.975)),
mean(medianhorn), quantile(medianhorn, c(0.025)), quantile(medianhorn, c(0.975)),
mean(medianjaccard), quantile(medianjaccard, c(0.025)), quantile(medianjaccard, c(0.975))))
}
####################################################################################################################################
#PLOTS FOR ONE OF THE THREE MODELS
op <- par(mfrow=c(2, 2))
par(pty="s")
plot(TAvalues,frame[,1], type="n",pch=21,bg="white",log="x",
xlab="Time-averaging (years)", ylab="Bray-Curtis similarity", lwd=2, ylim=c(0,1))
lines(TAvalues,frame[,1], type="b",lty=2,pch=21,bg="white",lwd=2)
lines(TAvalues,frame[,2], type="l",lty=2,pch=21,bg="white",lwd=2)
lines(TAvalues,frame[,3], type="l",lty=2,pch=21,bg="white",lwd=2)
plot(TAvalues,frame[,1], type="n",pch=21,bg="white",log="x",
xlab="Time-averaging (years)", ylab="Horn-Morisita similarity", lwd=2, ylim=c(0,1))
lines(TAvalues,frame[,4], type="b",lty=2,pch=21,bg="white",lwd=2)
lines(TAvalues,frame[,5], type="l",lty=2,pch=21,bg="white",lwd=2)
lines(TAvalues,frame[,6], type="l",lty=2,pch=21,bg="white",lwd=2)
plot(TAvalues,frame[,1], type="n",pch=21,bg="white",log="x",
xlab="Time-averaging (years)", ylab="Jaccard similarity", lwd=2, ylim=c(0,1))
lines(TAvalues,frame[,7], type="b",lty=2,pch=21,bg="white",lwd=2)
lines(TAvalues,frame[,8], type="l",lty=2,pch=21,bg="white",lwd=2)
lines(TAvalues,frame[,9], type="l",lty=2,pch=21,bg="white",lwd=2)