[R] optim() not finding optimal values

[Derek Ogle]

I am trying to use optim() to minimize a sum-of-squared deviations function based upon four parameters.  The basic function is defined as ...

SPsse <- function(par,B,CPE,SSE.only=TRUE)  {
  n <- length(B)                             # get number of years of data
  B0 <- par["B0"]                            # isolate B0 parameter
  K <- par["K"]                              # isolate K parameter
  q <- par["q"]                              # isolate q parameter
  r <- par["r"]                              # isolate r parameter
  predB <- numeric(n)
  predB[1] <- B0
  for (i in 2:n) predB[i] <- predB[i-1]+r*predB[i-1]*(1-predB[i-1]/K)-B[i-1]
  predCPE <- q*predB
  sse <- sum((CPE-predCPE)^2)
  if (SSE.only) sse
    else list(sse=sse,predB=predB,predCPE=predCPE)
}

My call to optim() looks like this

# the data
d <- data.frame(catch= c(90000,113300,155860,181128,198584,198395,139040,109969,71896,59314,62300,65343,76990,88606,118016,108250,108674), cpe=c(109.1,112.4,110.5,99.1,84.5,95.7,74.1,70.2,63.1,66.4,60.5,89.9,117.0,93.0,116.6,90.0,105.1))

pars <- c(800000,1000000,0.0001,0.17)                   # put all parameters into one vector
names(pars) <- c("B0","K","q","r")                      # name the parameters
( SPoptim <- optim(pars,SPsse,B=d$catch,CPE=d$cpe) )    # run optim()


This produces parameter estimates, however, that are not at the minimum value of the SPsse function.  For example, these parameter estimates produce a smaller SPsse,

parsbox <- c(732506,1160771,0.0001484,0.4049)
names(parsbox) <- c("B0","K","q","r")        
( res2 <- SPsse(parsbox,d$catch,d$cpe,SSE.only=FALSE) )

Setting the starting values near the parameters shown in parsbox even resulted in a movement away from (to a larger SSE) those parameter values.

( SPoptim2 <- optim(parsbox,SPsse,B=d$catch,CPE=d$cpe) )    # run optim()


This "issue" most likely has to do with my lack of understanding of optimization routines but I'm thinking that it may have to do with the optimization method used, tolerance levels in the optim algorithm, or the shape of the surface being minimized.

Ultimately I was hoping to provide an alternative method to fisheries biologists who use Excel's solver routine.

If anyone can offer any help or insight into my problem here I would be greatly appreciative.  Thank you in advance.