Lotaka-Volterra_Example.R

## Markus Gesmann, March 2012

library(simecol)
library(latticeExtra)

HaresLynxObservations <- "Year	 Hares.x.1000	 Lynx.x.1000
1900	30	4
1901	47.2	6.1
1902	70.2	9.8
1903	77.4	35.2
1904	36.3	59.4
1905	20.6	41.7
1906	18.1	19
1907	21.4	13
1908	22	8.3
1909	25.4	9.1
1910	27.1	7.4
1911	40.3	8
1912	57	12.3
1913	76.6	19.5
1914	52.3	45.7
1915	19.5	51.1
1916	11.2	29.7
1917	7.6	15.8
1918	14.6	9.7
1919	16.2	10.1
1920	24.7	8.6"

## Data sourced from
## http://www-rohan.sdsu.edu/~jmahaffy/courses/f00/math122/labs/labj/q3v1.htm
## Originally published in E. P. Odum (1953), Fundamentals of Ecology,
## Philadelphia, W. B. Saunders.
## Natural oscillations have been observed between the populations of
## predators and their prey. The predator-prey model was developed
## independently in the 1920s by Alfred J. Lotka and Vito
## Volterra. This problem examines the Lotka-Volterra model with data
## on lynx and snowshoe hares from North Canada. Between 1845 and
## 1935, the Hudson Bay company kept good records on the numbers of
## animals trapped in Canada

hl.obs <- read.table(text=HaresLynxObservations, header=TRUE)
names(hl.obs) <- c("Year", "Hares", "Lynx") ## Population in 000's
hl.obs$Type <- "Observations"

Plot <- xyplot(Hares + Lynx ~ Year, data=hl.obs, t="l", groups=Type,
               auto.key=list(points=FALSE, lines=TRUE),
               scales=list(alternating=1), 
               ylab="Population in 000's", xlab="Year")  
asTheEconomist(Plot)


## Load Lotka-Volterra model
data(lv)
## Modify model
main(lv) <- function (time, init, parms) 
{
  x <- init
  p <- parms
  dH <- p["a1"] * x[1] - p["a2"] * x[1] * x[2]
  dL <- -p["b1"] * x[2] + p["b2"] * x[1] * x[2]
  list(c(dH, dL))
}
init(lv) <-c(Hares=30, Lynx=4)
times(lv) <- c(from=1900, to=1920, by=0.1)
parms(lv) <- c(a1=0.5, a2=0.02, b1=0.9, b2=0.03)

output <- sim(lv)
plot(output)

## Fit ode mode

my.years <- hl.obs$Year
my.obs <- hl.obs[, c("Hares", "Lynx")]
myfit <- fitOdeModel(lv, whichpar=c("a1", "a2", "b1","b2"),
                     obstime=my.years, yobs=my.obs,
                     method="BFGS",
                     lower=c(a1=0, a2=0, b1=0, b2=0),
                     upper=c(a1=1, a2=1, b1=1, b2=1))

parms(lv) <-  myfit$par

sim.run <- sim(lv)
plot(sim.run)
sim.data <- out(sim.run)
names(sim.data)[1] <- "Year"
sim.data$Type <- "Model simulation"

allData <- rbind(hl.obs, sim.data)

## Compare real data with model simulation
Plot <- xyplot(Hares + Lynx ~ Year, data=allData, t="l", groups=Type,
               auto.key=list(points=FALSE, lines=TRUE),
               scales=list(alternating=1), 
               ylab="Population in 000's", xlab="Year")  
asTheEconomist(Plot)

phase.Plot <- xyplot(Hares ~ Lynx,  data=allData, t="l", groups=Type,
                     auto.key=list(points=FALSE, lines=TRUE),
                     scales=list(alternating=1), 
                     ylab="Hares vs. Lynx")  
asTheEconomist(phase.Plot)