Evolutionary Biology Online Journal Club

Dr. Matthew Forister provides R code for today’s paper!

3 Comments

I e-mailed Dr. Matthew Forister about getting a copy of the R code he used for his paper that we are discussing today entitled: “Ant association facilitates the evolution of diet breadth in a lycaenid butterfly“. About 1 hour later I got the following response:

Hi Thomas,

Sure, I’m happy to provide the code! See attached. It’s been a while since I looked at this, but I just verified that it does work, and I added a bit of annotation at the top that should help. Please let me know if you or your group has any questions.

thanks for your interest,
Matt

I specifically asked if I could post the code on our blog so our members could run the code and test the simulations. I have posted the code below. You can find out more about Dr. Forister’s research on his faculty webpage, or his lab webpage.

####### Code for running simulations in Forister et al. 2011, "Ant association facilitates..."
## Below there are two functions: runAntsPlants and plotAntsPlants
## After these functions have been loaded, you can run the simulation with defaults thus:
## result <- runAntsPlants()
## Of course parameters can be changed by adding specifications in those parentheses.
## And then time-plots of the results can be generated with:
## plotAntsPlants(result)
##
## Note that there is one external package that needs to be loaded (Hmisc) which provides
## an error bar function for the graphing.

library(Hmisc)
runAntsPlants<-function(nRep=1000,nGen=100,elSurvMin=0.15,elSurvMax=0.65,lpSurvMin=0.1,lpSurvMax=0.55,
paSurvMin=0.35,paSurvMax=0.8,aeMin=10,aeMax=180,AFmin=1,AFmax=1,MFmin=1,MFmax=1,
Tend=1,freqMed=0.5,eggK=50000,startPropK=0.75,astragFlowerFactor=1.95,medQualFactor=0.053,
antFactor=2.32,flowerFactor=7.53,medFecFactorWF=0.67,medFecFactorNF=0.17){
cat(“ants and plants model running \n”)
## set up data objects
eggsMed<-array(dim=c(nRep,nGen))
eggsAst<-array(dim=c(nRep,nGen))
larvaMed<-array(dim=c(nRep,nGen))
larvaAst<-array(dim=c(nRep,nGen))
pupaMed<-array(dim=c(nRep,nGen))
pupaAst<-array(dim=c(nRep,nGen))
adultMed<-array(dim=c(nRep,nGen))
adultAst<-array(dim=c(nRep,nGen))

## row1=median, row2=ub, row3=lb
summaryEM<-array(dim=c(3,nGen))
summaryEA<-array(dim=c(3,nGen))
summaryLM<-array(dim=c(3,nGen))
summaryLA<-array(dim=c(3,nGen))
summaryPM<-array(dim=c(3,nGen))
summaryPA<-array(dim=c(3,nGen))
summaryAM<-array(dim=c(3,nGen))
summaryAA<-array(dim=c(3,nGen))

## vectors for lambda
curLambda<-numeric(nGen-1)
lambda<-numeric(nRep)

## begin simulation and cycle through reps and generations
for (i in 1:nRep){
cat(“rep”,i,”\n”)

## draw AF and MF for initialization
MF<-MFmin + rbeta(1,2,3) * (MFmax – MFmin)
AF<-AFmin + rbeta(1,2,3) * (AFmax – AFmin)

## figure out proportion of eggs on Medicago
regOut<-lm(c(1/astragFlowerFactor,1) ~ c(0,1))
prefMed<-regOut$coefficients[2] * AF + regOut$coefficients[1]
propMed<-freqMed * prefMed

## startPropK gives the ratio of the initial egg population to carrying capacity for eggs
eggsMed[i,1]<-eggK * startPropK * propMed
eggsAst[i,1]<-eggK * startPropK *(1 – propMed)

## begin looping through generations
for (j in 1:nGen){

## draw AF and MF for current generation
MF<-MFmin + rbeta(1,2,3) * (MFmax – MFmin)
AF<-AFmin + rbeta(1,2,3) * (AFmax – AFmin)

## figure out proportion of eggs that will be laid this gen on each host plant
prefMed<-regOut$coefficients[2] * AF + regOut$coefficients[1]
propMed<-freqMed * prefMed

## figure out relatvie larval survival on medicago (based on ants and flower, once per gen)
medSurvival<-medQualFactor * (antFactor * Tend + 1 – Tend) * (flowerFactor * MF + 1 – MF)

## if there is at least one egg go ahead
if (eggsMed[i,j]>=1 | eggsAst[i,j]>=1){

## survival of eggs to newly hatched larvae, host, flowers and ants don’t matter here
if (eggsMed[i,j]>=1) larvaMed[i,j]<-eggsMed[i,j] * (elSurvMin + rbeta(1,2,3) * (elSurvMax – elSurvMin))
else larvaMed[i,j:nGen]<-rep(0,nGen-j+1)
if (eggsAst[i,j]>=1) larvaAst[i,j]<-eggsAst[i,j] * (elSurvMin + rbeta(1,2,3) * (elSurvMax – elSurvMin))
else larvaAst[i,j:nGen]<-rep(0,nGen-j+1)

## survival of newly hatched larvae to pupation, depends on host, and if med, ants and flowers
if (larvaMed[i,j]>=1) pupaMed[i,j]<-medSurvival *
(larvaMed[i,j] * (lpSurvMin + rbeta(1,2,3) * (lpSurvMax – lpSurvMin)))
else pupaMed[i,j:nGen]<-rep(0,nGen-j+1)
if (larvaAst[i,j]>=1) pupaAst[i,j]<-larvaAst[i,j] * (lpSurvMin + rbeta(1,2,3) * (lpSurvMax – lpSurvMin))
else pupaAst[i,j:nGen]<-rep(0,nGen-j+1)

## survival of pupae to adults, host, flowers and andts don’t matter
if (pupaMed[i,j]>=1) adultMed[i,j]<-pupaMed[i,j] * (paSurvMin + rbeta(1,2,3) * (paSurvMax – paSurvMin))
else adultMed[i,j:nGen]<-rep(0,nGen-j+1)
if (pupaAst[i,j]>=1) adultAst[i,j]<-pupaAst[i,j] * (paSurvMin + rbeta(1,2,3) * (paSurvMax – paSurvMin))
else adultAst[i,j:nGen]<-rep(0,nGen-j+1)

## reproduction, requires at least two adults, and doesn’t occur in the last generation
if ((adultMed[i,j] + adultAst[i,j]) >= 2){
if (j < nGen){
astragEggs<-aeMin + rbeta(((adultAst[i,j])/2),2,3) * (aeMax – aeMin)
medWFeggs<-aeMin + rbeta(((adultMed[i,j])/2 * MF),2,3) * ((aeMax * medFecFactorWF) – aeMin)
medNFeggs<-aeMin + rbeta(((adultMed[i,j])/2 * (1 – MF)),2,3) * ((aeMax * medFecFactorNF) – aeMin)
totalEggs<-sum(c(astragEggs,medWFeggs,medNFeggs))
if (totalEggs > eggK) totalEggs<-eggK
eggsMed[i,j+1]<-totalEggs * propMed
eggsAst[i,j+1]<-totalEggs * (1 – propMed)
}
}
else {
eggsMed[i,j:nGen]<-rep(0,nGen-j+1)
eggsAst[i,j:nGen]<-rep(0,nGen-j+1)
larvaMed[i,j:nGen]<-rep(0,nGen-j+1)
larvaAst[i,j:nGen]<-rep(0,nGen-j+1)
pupaMed[i,j:nGen]<-rep(0,nGen-j+1)
pupaAst[i,j:nGen]<-rep(0,nGen-j+1)
adultMed[i,j:nGen]<-rep(0,nGen-j+1)
adultAst[i,j:nGen]<-rep(0,nGen-j+1)
}
}
else {
eggsMed[i,j:nGen]<-rep(0,nGen-j+1)
eggsAst[i,j:nGen]<-rep(0,nGen-j+1)
larvaMed[i,j:nGen]<-rep(0,nGen-j+1)
larvaAst[i,j:nGen]<-rep(0,nGen-j+1)
pupaMed[i,j:nGen]<-rep(0,nGen-j+1)
pupaAst[i,j:nGen]<-rep(0,nGen-j+1)
adultMed[i,j:nGen]<-rep(0,nGen-j+1)
adultAst[i,j:nGen]<-rep(0,nGen-j+1)
}
if (j != 1){
if ((adultMed[i,j-1] + adultAst[i,j-1]) > 0) curLambda[j-1]<-(adultMed[i,j] + adultAst[i,j])/(adultMed[i,j-1] + adultAst[i,j-1])
else curLambda[i]<-NA
}
## end of generation loop
}
lambda[i]<-exp(mean(log(curLambda),na.rm=TRUE))
## end of rep loop
}
## calculate median and quantiles for each generation, I can make this mean and sd if prefered
lb<-0.25
ub<-0.75
for (j in 1:nGen){
summaryEM[,j]<-quantile(eggsMed[,j],prob=c(lb,0.5,ub),na.rm=TRUE)
summaryEA[,j]<-quantile(eggsAst[,j],prob=c(lb,0.5,ub),na.rm=TRUE)
summaryLM[,j]<-quantile(larvaMed[,j],prob=c(lb,0.5,ub),na.rm=TRUE)
summaryLA[,j]<-quantile(larvaAst[,j],prob=c(lb,0.5,ub),na.rm=TRUE)
summaryPM[,j]<-quantile(pupaMed[,j],prob=c(lb,0.5,ub),na.rm=TRUE)
summaryPA[,j]<-quantile(pupaAst[,j],prob=c(lb,0.5,ub),na.rm=TRUE)
summaryAM[,j]<-quantile(adultMed[,j],prob=c(lb,0.5,ub),na.rm=TRUE)
summaryAA[,j]<-quantile(adultAst[,j],prob=c(lb,0.5,ub),na.rm=TRUE)
}
## determine proportion of replicates ending in extinction
propExtinct<-sum((adultMed[,nGen] + adultAst[,nGen]) == 0)/nRep
cat(“Proportion of replicates extinct: “,propExtinct,”\n”)
cat(“lamda: median=”, quantile(lambda,probs=0.5),” lb=”, quantile(lambda,probs=0.025),” ub=”,quantile(lambda,probs=0.975),”\n”)
## make list to be returned
resultsSummary<-list(summaryEM,summaryEA,summaryLM,summaryLA,summaryPM,summaryPA,summaryAM,summaryAA)
names(resultsSummary)<-c(“eggsMed”,”eggsAst”,”larvaMed”,”larvaAst”,”pupaMed”,”pupaAst”,”adultMed”,
“adultAst”)
out<-list(resultsSummary,eggsMed,eggsAst,larvaMed,larvaAst,pupaMed,pupaAst,adultMed,adultAst,
propExtinct,lambda)
names(out)<-c(“summary”,”eggsMed”,”eggsAst”,”larvaMed”,”larvaAst”,”pupaMed”,”pupaAst”,”adultMed”,
“adultAst”,”propExtinct”,”lambda”)
return(out)
## end of function
}
## Hmisc needed for errorbar plots, library location should be adjusted
library(Hmisc,lib.loc=”~/Library/R/2.8/library”)
## this plots the summary data, which is probably the most useful for now
plotAntsPlants<-function(results=NULL,pdf=FALSE,myFile=”./antsPlants.pdf”,myfile2=”./antsPlantsLam.pdf”){
## make errbar plots
if (pdf == TRUE) pdf(file=myFile,width=16,height=11)
else quartz()
par(mfrow=c(2,4))
par(pty=”s”)
errbar(1:dim(results$eggsMed)[2],results$summary$eggsMed[2,],results$summary$eggsMed[3,],results$summary$eggsMed[1,],
xlab=”Generation number”,ylab=”No. eggs on Medicago”)
axis(1)
axis(2)
box()
errbar(1:dim(results$eggsMed)[2],results$summary$eggsAst[2,],results$summary$eggsAst[3,],results$summary$eggsAst[1,],
xlab=”Generation number”,ylab=”No. eggs on Astragalus”)
axis(1)
axis(2)
box()
errbar(1:dim(results$eggsMed)[2],results$summary$larvaMed[2,],results$summary$larvaMed[3,],results$summary$larvaMed[1,],
xlab=”Generation number”,ylab=”No. larva on Medicago”)
axis(1)
axis(2)
box()
errbar(1:dim(results$eggsMed)[2],results$summary$larvaAst[2,],results$summary$larvaAst[3,],results$summary$larvaAst[1,],
xlab=”Generation number”,ylab=”No. larva on Astragalus”)
axis(1)
axis(2)
box()
errbar(1:dim(results$eggsMed)[2],results$summary$pupaMed[2,],results$summary$pupaMed[3,],results$summary$pupaMed[1,],
xlab=”Generation number”,ylab=”No. pupa on Medicago”)
axis(1)
axis(2)
box()
errbar(1:dim(results$eggsMed)[2],results$summary$pupaAst[2,],results$summary$pupaAst[3,],results$summary$pupaAst[1,],
xlab=”Generation number”,ylab=”No. pupa on Astragalus”)
axis(1)
axis(2)
box()
errbar(1:dim(results$eggsMed)[2],results$summary$adultMed[2,],results$summary$adultMed[3,],results$summary$adultMed[1,],
xlab=”Generation number”,ylab=”No. adults from Medicago”)
axis(1)
axis(2)
box()
errbar(1:dim(results$eggsMed)[2],results$summary$adultAst[2,],results$summary$adultAst[3,],results$summary$adultAst[1,],
xlab=”Generation number”,ylab=”No. adults from Astragalus”)
axis(1)
axis(2)
box()
#if (pdf == TRUE) dev.off()
## lambda plot
#if (pdf == TRUE) pdf(file=myFile2,width=16,height=11)
#else quartz()
#boxplot(results$lambda,ylab=”lambda”)
#if(pdf == TRUE) dev.off()

Advertisements

Author: hossiet

PhD student at Carleton University. I study predator-prey interactions. My current project examines the ecology and evolution of caterpillar eyespots.

3 thoughts on “Dr. Matthew Forister provides R code for today’s paper!

  1. This is awesome Tom, good job!

  2. Pingback: Forister et al 2011 Proc B – Paper summary « Evolutionary Biology Online Journal Club

  3. Pingback: Archive of Papers Discussed « Evolutionary Biology Online Journal Club

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s