About a fortnight ago I posed a question regarding fitting predator-prey
models. Thanks to Nick Ellis who posted an answer, which was what I
required. Also to Dave Fournier who suggested using his software "AD
Model Builder" which is "hundreds of times faster than splus" for my problem.
Nick's answer follows, in case you missed it ....
------- Forwarded message follows -------
From: "Nick Ellis" <nick.ellis@marine.csiro.au>
Date sent: Tue, 10 Apr 2001 12:54:25 +1000
The problem is non-linear optimization for which there exist tools nls(),
nlmin() and nlminb(). I had a go at simulating some data and estimating with
nlminb with a guess at the variance model (normal in N and P, sigma=1/time).
I attach my code below. Hope it's helpful.
Nick Ellis
LotkaVolterra <- function(n, a, b, cc, r, K, N0, P0) {
N <- P <- numeric(n)
N[1] <- N0 + r*N0*(1-N0/K) - a*P0*N0
P[1] <- P0 + b*P0*(a*N0-cc)
for (i in 2:n) {
N[i] <- N[i-1] + r*N[i-1]*(1-N[i-1]/K) - a*P[i-1]*N[i-1]
P[i] <- N[i-1] + b*P[i-1]*(a*N[i-1]-cc)
}
cbind(N=N,P=P)
}
# check it works
plot(LotkaVolterra(20,.1,.1,.1,1,10,1,1),type='l')
# generate some data with normal errors with sigma=1/time
dat <- LotkaVolterra(20,.1,.1,.1,1,10,1,1)+
rnorm(40,0,matrix(rep(1/(1:20),2),ncol=2))
# nls didn't work for me
nls(dat ~ LotkaVolterra(20,a,b,cc,r,10,1,1),
start=list(a=.1,b=.1,cc=.1,r=1), trace=T, control=list(minscale=0.0001,
maxiter=200))
# try different approach with nlminb costraining all parameters to be
positive
# need auxiliary functions to package parms into single vector
fLVaux <- function(x,dat) { # evaluates LotkaVolterra()
n <- nrow(dat)
LotkaVolterra(n,x[1],x[2],x[3],x[4],x[5],1,1)
}
fLV <- function(x, dat) { # weighted sum of squares
n <- nrow(dat)
sum((dat - LotkaVolterra(n, x[1], x[2], x[3], x[4], x[5], 1, 1))^2 *
(1:n)^2)
}
# do optimization, get predicted values and plot
res <- nlminb(c(.1,.1,.1,1,10),fLV,lower=c(0,0,0,0,0),dat=dat)
pred <- fLVaux(res$parameters, dat)
plot(dat,type='l')
lines(pred,lty=3,col=2)
> -----Original Message-----
> I'm using S-PLUS 2000 on Windows 95, and wish to
> fit a predator/prey model to some mite data (Background: the prey is a 'pest'
> species, and the obligate predator is introduced each season in an attempt at
> control. The hope of the modeling is to be able to predict whether the
> predator will bring the prey under control and thus avoid the necessity of
> spraying). Initially, I would like to fit a basic Lotka-Volterra type model:
>
> delta N = rN(1 - N/K) - aPN
> &
> delta P = bP(aN - c)
>
> where delta N and delta P are the increments (from time t to t+1) in prey and
> predator numbers respectively; N and P are the numbers of prey and predators
> at time t; and r, K, a, b and c are coefficients to be estimated.
>
> Given a set of estimates of prey and predator numbers (for t in 0:n), my
> question is how to estimate the coefficients using S-Plus? I'd appreciate
> pointers to appropriate literature as well as more direct advice.
Peter Alspach
HortResearch Phone +64-3-528 0708
PO Box 220, Motueka FAX +64-3-528 7813
NEW ZEALAND Internet palspach@hortresearch.co.nz
____________________________________________________
The contents of this e-mail are privileged and/or confidential to the
named recipient and are not to be used by any other person and/or
organisation. If you have received this e-mail in error, please notify
the sender and delete all material pertaining to this e-mail.
____________________________________________________
|