[This article was first published on me nugget, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
I had a reader ask me recently to help understand how to modify the code of an individual-based model (IBM) that I posted a while back. It was my first attempt at an IBM in R, and I realized that I have made some significant changes to the way that I code such models nowadays. Most of the changes are structural, but seem to help a lot in clearly documenting the model and its underlying processes.
Basically, I follow a structure that a friend taught me (a very experienced modeller, specializing in IBMs). Granted R isn’t the best language for such models but, depending on your computational needs, it can also be quite easy to implement if you already have experience in basic R programming. The idea is to code important processes of the IBM as functions, which accept an object of individuals as the main argument. The functions can be as simple or elaborate as needed, but the key is that when you finally set up your simulation, you only need to call these functions again. This results in easily legible code, where you are less likely to get lost in the model and can concentrate more on the processes to be performed during each iteration (e.g. growth, reproduction, mortality):
inds <- grow.inds(inds)
inds <- reproduce.inds(inds)
inds <- kill.inds(inds)
Since the previous post, I have gone towards using data.frame objects to store individuals and their attributes since I find it easier to apply functions for extracting summary statistics (e.g. histogram of age distribution in the final population):
The example shown here is a modification of the genetic drift example that I showed before – only this time I have included 4 color phenotypes. The model runs until either one phenotype dominates the population or a maximum number of iterations is reached. Reproduction allows for a individuals to have more than one offspring per time step, but skewed towards zero. Death is modeled by a constant instantaneous mortality rate. I have used such a setup with more complicated models of fish genetics and found the performance to be quite fast, even with population sizes of 300,000 individuals. The key is to maintain as much vectorization in the functions as possible.
To reproduce the example:
######################################
### Population dynamics parameters ###
######################################
Colors <- c("blue", "green", "orange", "red") # population phenotypes
tmax <- 2000 # maximum time duration of simulation
n.initial <- 50*length(Colors) # inital population size
CC <- 500 # carrying capacity of offspring
M <- 0.5 # natural mortality
#########################
### Process functions ###
#########################
# Creation
make.inds <- function(id=NaN, color=NaN, age=NaN){
inds <- data.frame(id = id, color = color, age = 0) # give then starting attributes (id, color, age)
inds
}
# Growth
grow.inds <- function(inds){
inds$age <- inds$age + 1 # advance age by one
inds
}
# Reproduction (breeding, inheritance)
reproduce.inds <- function(inds){
n.eggs <- rlnorm(nrow(inds), mean=0.01, sd=0.7) # how many offspring per adult
n.eggs <- round(n.eggs * (1 - sum(n.eggs)/CC)) # mediate by carrying capacity
if(sum(n.eggs) > 0){
motherRow <- rep(seq(nrow(inds)), times=n.eggs) # mother's row
offspring <- make.inds(id=seq(max(inds$id)+1, length.out=sum(n.eggs))) # make offspring
offspring$color <- inds$color[motherRow] # inherit color
inds <- rbind(inds, offspring) # add new offspring to adults
}
inds
}
# Mortality
kill.inds <- function(inds){
pDeath <- 1 - exp(-M) # prob. of death
kill <- which(runif(nrow(inds)) < pDeath) # kill these ones
inds <- inds[-kill,]
inds
}
###################
### Simulation ####
###################
# Initial population
inds <- make.inds(id=1:n.initial, color=as.factor(array(Colors, dim=n.initial)))
head(inds, 10)
# Object for storing results (population size in time)
N <- array(NaN, dim=c(tmax, length(Colors)))
colnames(N) <- Colors
N[1,] <- summary(inds$color)
# Run
set.seed(1) # optional starting seed
for(t in 2:tmax){
# population processes
inds <- grow.inds(inds)
inds <- reproduce.inds(inds)
inds <- kill.inds(inds)
# store results
N[t,] <- summary(inds$color)
print(paste("time =", t, "; n = ", sum(N[t,])))
# break when one phenotype dominates
if(sum(N[t,] > 0) == 1) {break}
}
N <- N[1:t,] # remove excess rows of results object
##############
#### Plots ###
##############
# Population sizes
png("pop_size.png", width=5, height=5, res=300, units="in", type="cairo")
op <- par(mar=c(5,5,1,1))
ylim <- 1.2 * range(N)
Time <- seq(nrow(N))
plot(Time, N[,1], t="n", lwd=2, col=4, ylim=ylim, ylab="Population size")
for(i in seq(Colors)){
lines(Time, N[,i], col=Colors[i], lwd=2)
}
legend("topleft", legend=paste(Colors, "pop."), lwd=2, col=Colors, bty="n")
par(op)
dev.off()
# Age distribution of final population
png("age_dist.png", width=5, height=5, res=300, units="in", type="cairo")
op <- par(mar=c(5,5,2,1))
hist(inds$age, freq=FALSE, col=8, xlab="Age [years]", main="Age distribution")
par(op)
dev.off()To leave a comment for the author, please follow the link and comment on their blog: me nugget.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.