This code produces only one of the scenarios in Figure 6.11. The parameter nu
needs to be altered to generate the other scenarios.
WARNING 1: this code takes a while to run (~20 minutes or more)
WARNING 2: More than half the code is post-simulation calculations in order to generate graphics. The end result (a relative abundance distribution) is needed for the code in Online Box 8
## specify parameters, initial conditions, and output vector
num.years <- 10000
J <- 10000
sp <- 1 # a counter to keep track of added species
COM <- vector(length = J); COM[] <- sp
nu <- 5e-04 # speciation rate
year <- 2
COM.out <- matrix(nrow = (num.years/100+1), ncol = J); COM.out[1,] <- COM
## run simulation
for (i in 1:(J*(num.years))) {
dead.indiv <- ceiling(J*runif(1))
## speciation (replace dead individual with new species)
if (runif(1) < nu) {
COM[dead.indiv] <- sp+1
sp <- sp + 1
} else {
## not speciation (replace dead individual from community)
COM[dead.indiv] <- COM[ceiling(J*runif(1))]
}
## record data every 100 years
if (i %% (J*100) == 0) {
COM.out[year,] <- COM
year <- year+ 1
} }
################################
# Calculations to allow plotting relative abundance distributions
## First calculate the number of species (SR) at each time
SR.out <- vector(length = dim(COM.out)[1])
for (t in 1:dim(COM.out)[1]) {
COM.hist <- hist(COM.out[t,],c(1:max(COM.out[t,])),plot=FALSE)
abund.vec <- COM.hist$counts[COM.hist$counts>0]
SR.out[t] <- sum(abund.vec>0)
}
## Calculate species abundances
abund.out <- matrix(nrow = (num.years/100+1), ncol = max(SR.out))
abund.out[] <- 0
for (t in 1:dim(COM.out)[1]) {
COM.hist <- hist(COM.out[t,],c(1:max(COM.out[t,])),plot=FALSE)
abund.vec <- COM.hist$counts[COM.hist$counts>0]
abund.out[t,1:SR.out[t]] <- abund.vec
}
abund.out <- abund.out/J
## The final rank-order of relative abundances
final.rad <- sort(abund.out[(num.years/100+1),1:SR.out[(num.years/100+1)]],decreasing=TRUE)
## Write the final relative abundance distribution to a file for use in Online Box 8 (island biogeography model)
write.csv(final.rad, file = "RAD1.csv")
## graph the results
plot(log(final.rad),type="l",
xlab="Species rank order", ylab="log(Relative abundance)")
num.years
, J
, and COM
are defined as in previous Boxes. nu
is the speciation rate. With probability nu
,
the dead individual will be replaced by a species new to the
metacommunity. Since we will have an undetermined number of species, sp
will keep track of species identities. The first species is species 1, and each new species will be numbered in sequence.num.years <- 100
J <- 100
sp <- 1 # a counter to keep track of added species
COM <- vector(length = J); COM[] <- sp
nu <- 1e-04 # speciation rate
ncol = J
because we will keep track of the identity of each individual in each
year of recorded data. In order to reduce the duration of the
simulation, we will record data only every 100 years, which is why nrow = (num.years/100+1)
COM.out <- matrix(nrow = (num.years/100+1), ncol = J); COM.out[1,] <- COM
COM
can be updated, we select an individual to die before figuring out if
it will be replaced by a speciation even or by a local reproducer. dead.indiv <- ceiling(J*runif(1))
nu
), we update the sp
counter, keeping track of how many species have entered the community. if (runif(1) < nu) {
COM[dead.indiv] <- sp+1
sp <- sp + 1
} else {
J*100
deaths) instead of every year (after J
deaths) as in previous Boxes). if (i %% (J*100) == 0) {
COM.out[year,] <- COM
year <- year+ 1
} }
dim(COM.OUT)[1]
returns the first dimension (number of rows) of COM.OUT
, which is the number of time steps at which we recorded data. The hist
(histogram) function returns, in the $counts
object, a vector of abundances of each species in the community, most
of which will be zero (because most species entering via speciation go
quickly extinct). We specify plot=FALSE
so that R doesn’t produce a histogram graph each time. abund.vec
will contain only the non-zero elements, and abund.vec>0
returns a vector with a ‘1’ for each species with non-zero abundance, the sum of which is species richness.SR.out <- vector(length = dim(COM.out)[1])
for (t in 1:dim(COM.out)[1]) {
COM.hist <- hist(COM.out[t,],c(1:max(COM.out[t,])),plot=FALSE)
abund.vec <- COM.hist$counts[COM.hist$counts>0]
SR.out[t] <- sum(abund.vec>0)
}
J
returns relative abundances that sum to one.abund.out <- matrix(nrow = (num.years/100+1), ncol = max(SR.out))
abund.out[] <- 0
for (t in 1:dim(COM.out)[1]) {
COM.hist <- hist(COM.out[t,],c(1:max(COM.out[t,])),plot=FALSE)
abund.vec <- COM.hist$counts[COM.hist$counts>0]
abund.out[t,1:SR.out[t]] <- abund.vec
}
abund.out <- abund.out/J
final.rad
is the rank-order relative abundance vector
at the final time step. This needs to be saved if we want to use it as
the metacommunity relative abundance distribution in Online Box 8. Note
that you will need to create/set the directory on your own (i.e.,
replaced c:/temp/RAD1.csv
with whatever directory and filename you want).final.rad <- sort(abund.out[(num.years/100+1),1:SR.out[(num.years/100+1)]],decreasing=TRUE)
## Write the final relative abundance distribution to a file for use in Online Box 8 (island biogeography model)
write.csv(final.rad, file = "RAD1.csv")