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]] _______________________________________________ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology