(This article was first published on

**me nugget**, and kindly contributed to R-bloggers)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 on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...