[R-sig-eco] Fitting Lotka-Volterra Competition Model !

Alexandre F. Souza alexsouza.cb.ufrn.br at gmail.com
Sun Jan 17 18:39:30 CET 2016


Dear friends,

I am organizing an undergraduate course on population biology in Natal,
northeastern Brazil. I am willing to use the functions and code posted by
Noam Ross aimed at finding the parameters of the lotka-volterra competition
equation below (
https://github.com/noamross/Working-Code/blob/master/OldStuff/logisticFit.R
).

# Data-fitting to a continuous-time model
# here we'll fit the paramecium data to a logistic model, then you'll fit
the competition paramters

# make sure the library is loaded
library(deSolve) # loads the ODE numerical integration package

# the approach: minimize square deviations
# run the logistic model under a given r and K, find the squared deviations
from the data
# use optim to find the least squared deviation

# so we need the logistic model ready to run with lsoda:
logInt = function(t, n, parms) {
  with(as.list(parms), { # extract parameters from parms vector
    dn = r*n*(K-n)/K # logistic dn/dt
    return(list(dn)) # return dn/dt as a list
  })
}

# plus a function to find the square deviations
logIntMin = function(parms) {
  out = as.data.frame(lsoda(n0, times, logInt, parms)) # run ode
  mse = mean((out$n-nTrue)^2) # calculate mean squared error between
simulated and original data
  return(mse) # return mean squared error
}

# parameters and data
nTrue =
as.matrix(c(2,10,17,29,39,63,185,258,267,392,510,570,650,575,650,550,480,520,520,500))
# actual data
n0 = c(n=nTrue[1]) # initial population size - us first data point, make
sure to label
tf = length(nTrue) # time
times = 1:tf # vector of times to run over
rguess = 1 # guess for growth rate
Kguess = 500 # guess for carrying capacity
parms0 = c(r=rguess, K=Kguess)

# find fit and run simulations
optimOut = optim(parms0, logIntMin) # give optim initial guess and function
parms = optimOut$par # extract parameters
nSim = as.data.frame(lsoda(n0,times,logInt,parms)) # run ode to get
simulated population

# plot results - actual and simulated data
plot(times, nTrue, type="p", xlab="Time", ylab="Abundance")
lines(nSim$time, nSim$n)

# task: fitting paramecium data to Lotka-Volterra competition model
# you get r's and K's - this is in the "parameciumrk.txt" file - laoding
this under,
# the default for read.table will lead to column names V1, V2, etc. so the
evnetual labels,
# for the r's and K's will be r.V1, r.V2, etc.; same for n's
# you'll have to find the alpha's to the data in the "paramecium2.txt" file
# note that this means you have to reconstruct the parameters vector within
optim
# start with parmsrk=r and K, then do parms = c(parmsrk, a=a) to add a's
within minimizing,
# function

parameciumrk.txt

0.7816  0.6283
559.6860 202.4931

I run it and read it carefully. Yet, I could not progress to finding the
alpha (a) coefficients for the two-species version you mention at the end
of the code, and it is just these coefficients that I need to include in
the class.

If yes, my questions are:

1 - To build this extension I should work with each species at a time? You
say I should load the data and r an K parameters of the two species, but
the function logInt itself at lines 12-7 works for one species at a time,
right?

2 - This function by the way needs to be modified in order to include the
a, right? Is the modification just to alter the math part itself at line 14
from

dn = r*n*(K-n)/K # logistic dn/dt
to

dn = r*n1*(K-n1-a*n2)/K
However this second equation invokes the second dataset, and I could not
figure it out how to include both datasets inthe code, my efforts to adapt
the remaining of the code to this two-speces situation failed.

The point is that it is not difficult to create a class on the theory of
the Lotka-Volterra competition model, but students generally feel that it
is quite abstract and detached from real-world data and tend not to
consider it a concrete and applicable piece of knowledge. I think this code
makes this useful and concrete, and I would very much like to base the
whole class on its application to concrete time series.

Thank you very much in advance for any assistance,

Sincerely,

Alexandre

-- 
Dr. Alexandre F. Souza
Professor Adjunto III
Universidade Federal do Rio Grande do Norte
CB, Departamento de Ecologia
Campus Universitário - Lagoa Nova
59072-970 - Natal, RN - Brasil
lattes: lattes.cnpq.br/7844758818522706
http://www.docente.ufrn.br/alexsouza

	[[alternative HTML version deleted]]



More information about the R-sig-ecology mailing list