optimization - Error fitting a simulation model to data with optim in R -

good afternoon: have 30 years of population data on guanacos, , matrix simulation model; want estimate optimal values of 12 of model's parameters in r using optim, minimizing residual sums of squares (sum(obs-pred)^2). created function model. model working fine, if invoke function fixed parameter values exected results. when call optim sewt of initial parameters , function, following error message: "the argument a13 absent, no value omission"; also: message of lost warning: in if (ssq == 0) { : condition has length > 1 , first element used called from: top level (freely translated spanish).

please find below r code in 3 sections: (1) first function "guafit" declaration, (2) standalone run of model invoking "guafit", , (3) call of "optim" starting parameter values.

thank you, jorge rabinovich

# section (1): # clean     rm(list=ls())  ##################################################################### ####################### function "guafit"  ########################## #####################################################################      guafit <- function(ssq,a13,a21,a32,a33,ad1,ad2,ad3,ad4,bd1,bd2,bd3,bd4) {  # create field population (a 30-years time series) # ====================================================     tfield <- c(12334,10670,19078,11219,11771,12323,13027,14094,14604,17775,20774,16410,17626,21445,21111,20777,28978,27809,28935,38841,38363,32273,43128,58597,52456,33125,61334,60488,44773,56973)  # assign values vector of initial population # =====================================================     g <- matrix(c(0,0,0),ncol=3, nrow=1)     g[1]= 1542     g[2]= 740     g[3]= 3885  # load matrices initial values 30 time units (years) # =========================================================================     if (ssq == 0) {     a<-array(0,dim=c(3,3,30))     for(i in 1:29) {     a[1,3,i]= a13     a[2,1,i]= a21     a[3,2,i]= a32     a[3,3,i]= a33             }         } # initialize  variables # ==========================     tmod<-array(1,dim=c(1,30)); tmod <- as.numeric(tmod)     densprop<-array(1,dim=c(1,30)); densprop <- as.numeric(densprop)     fdeltafe<-array(1,dim=c(1,30)); fdeltafe <- as.numeric(fdeltafe)     fdeltasc<-array(1,dim=c(1,30)); fdeltasc <- as.numeric(fdeltasc)     fdeltasj<-array(1,dim=c(1,30)); fdeltasj <- as.numeric(fdeltasj)     fdeltasa<-array(1,dim=c(1,30)); fdeltasa <- as.numeric(fdeltasa)  # n0 initial population vector # multiplied 2 represewnt both sexes # =============================================== # transfer guanacos (g) vector 3 age classes     n0 <- g     tmod[1] <- (n0[1]+n0[2]+n0[3]) * 2  # declaration of initial simulation conditions # ================================================ # ng number of female individuals per age class (dim 3x30) # tmod total (both sexes) population (sum of 3 age classes * 2)     ng <- matrix( 0, nrow= 3, ncol=30)     ng[,1] <- n0  # assume constant carrying capacity (k= 60000 individuals)     carcap= 60000  # start simulation 30 years     for(i in 1:29) { #set density-dependent coefficients      densprop[i] <- tmod[i] / carcap  # calculate density-dependent factors     fdeltafe[i]= 1/(1+exp((densprop[i]-ad1)*bd1))     fdeltasc[i]= 1/(1+exp((densprop[i]-ad2)*bd2))     fdeltasj[i]= 1/(1+exp((densprop[i]-ad3)*bd3))     fdeltasa[i]= 1/(1+exp((densprop[i]-ad4)*bd4))  # apply density-dependent factors each coefficient in own age class      a[1,3,i]= a[1,3,i] * fdeltafe[i]     a[2,1,i]= a[2,1,i] * fdeltasc[i]     a[3,2,i]= a[3,2,i] * fdeltasj[i]     a[3,3,i]= a[3,3,i] * fdeltasa[i]  # project total population matrix operation      ng[,i+1] <- a[,,i]%*%ng[,i]     tmod[i+1] <- (ng[1,i+1] + ng[2,i+1] + ng[3,i+1]) * 2  # end of 30-years simulation loop             }  # calculate residual sum of squares (ssq)                    ssq = sum((tmod - tfield)^2)                return(list(outm=tmod, outc=tfield, elssq=ssq, matrices=a,     losgua=g, losguaxe=ng))  # end of function guafit         } #################################################################################  # section (2):  # initialize conditions , parameters before calling function guafit      ssq <- 0  # initialize 8 density-dependent coefficients (2 coefficients per age class) # ==============================================================================     ad1= 1.195680167     ad2= 1.127219245     ad3= 1.113739384     ad4= 1.320456815     bd1= 10.21559509     bd2= 9.80201883     bd3= 9.760834107     bd4= 10.59390027  # assign initial values transition matrix coefficients # ============================================================     a21= 0.6     a32=0.8     a33=0.9     a13=0.37  # initialization of conditions finished, can  call function guafit # test, call function guafit once check results     myresults <- guafit(ssq,a13,a21,a32,a33,ad1, ad2, ad3, ad4, bd1, bd2, bd3, bd4) # save results of interest of function guafit new variables      restmod <- myresults$outm     tfield <- myresults$outc     ssqvalue <- myresults$elssq     lasmatrices <- myresults$matrices     reslosgua <- myresults$losgua     reslosguaxe <- myresults$losguaxe     ssqvalue # graph results     axisx <- c(1982:2011)     plot(axisx,tfield)     lines(axisx,restmod)  #################################################################################  # section (3):  # try optim function   # first creating initial parameter values variable pass argument      startparam <- c(ssq, a13,a21,a32,a33,ad1, ad2, ad3, ad4, bd1, bd2, bd3, bd4)     optim(startparam, guafit) # , got error message mentioned above.  # tried calling as:     optim(par=startparam, fn=guafit)  # , got same error message  # tried calling optim passing values directly list of values:     startparam <- c(ssq=0, a13=0.37, a21=0.6, a32=0.8, a33=0.9, ad1=1.1, ad2=1.1, ad3=1.1, ad4=1.1, bd1=10, bd2=10, bd3=10, bd4=10)     optim(startparam, guafit)     optim(par=startparam, fn=guafit)   # , got same error message 

i got code running, not sure whether returns correct parameter estimates. tried 2 different approaches.

with first approach did following:

1. placed data , initializations outside function.
2. removed constants optim statement
3. assumed ad1, ad2, ad3 , ad4 constants , not parameters estimated. guess assumed ad1, ad2, ad3 , ad4 constant offsets.

this first approach returned estimates of bd1, bd2, bd3 , bd4 closely matched initial values.

with second approach did following:

1. placed data , initializations outside function.
2. removed constants optim statement
3. assumed ad1, ad2, ad3 , ad4 parameters estimated.

this second approach returned estimates of 8 parameters: ad1, ad2, ad3, ad4, bd1, bd2, bd3 , bd4 not think estimates correct. note 2 columns in hessian 0. se not estimated second approach , estimates of bd1, bd2, bd3 , bd4 not close initial values.

the r code both approaches below. check code see whether either approach doing want. if treating ad1, ad2, ad3, , ad4 parameters critical perhaps model code have revised.

here r code , output first approach:

set.seed(2345)  guafit <- function(param) {            bd1 <- param[1]           bd2 <- param[2]           bd3 <- param[3]           bd4 <- param[4]  # start simulation 30 years     for(i in 1:29) {  # set density-dependent coefficients     densprop[i] <- tmod[i] / carcap  # calculate density-dependent factors     fdeltafe[i]= 1/(1+exp((densprop[i]-ad1)*bd1))     fdeltasc[i]= 1/(1+exp((densprop[i]-ad2)*bd2))     fdeltasj[i]= 1/(1+exp((densprop[i]-ad3)*bd3))     fdeltasa[i]= 1/(1+exp((densprop[i]-ad4)*bd4))  # apply density-dependent factors each coefficient in own age class     a[1,3,i]= a[1,3,i] * fdeltafe[i]     a[2,1,i]= a[2,1,i] * fdeltasc[i]     a[3,2,i]= a[3,2,i] * fdeltasj[i]     a[3,3,i]= a[3,3,i] * fdeltasa[i]  # project total population matrix operation     ng[,i+1]  <- a[,,i]%*%ng[,i]     tmod[i+1] <- (ng[1,i+1] + ng[2,i+1] + ng[3,i+1]) * 2  # end of 30-years simulation loop             }  # calculate residual sum of squares (ssq)                    ssq = sum((tmod - tfield)^2)  # end of function guafit         }  # create field population (a 30-years time series) # ====================================================     tfield <- c(12334,10670,19078,11219,11771,12323,13027,14094,14604,17775,                 20774,16410,17626,21445,21111,20777,28978,27809,28935,38841,                 38363,32273,43128,58597,52456,33125,61334,60488,44773,56973)  # initialize conditions , parameters before calling function guafit     ssq <- 0  # assign initial values transition matrix coefficients # ============================================================     a21 = 0.6     a32 = 0.8     a33 = 0.9     a13 = 0.37  # assign values vector of initial population # =====================================================     g <- matrix(c(0,0,0),ncol=3, nrow=1)     g[1] = 1542     g[2] = 740     g[3] = 3885  # load matrices initial values 30 time units (years) # =========================================================================      <- array(0,dim=c(3,3,30))      for(i in 1:29) {     a[1,3,i]= a13     a[2,1,i]= a21     a[3,2,i]= a32     a[3,3,i]= a33             }  # initialize  variables # ==========================     tmod<-array(1,dim=c(1,30)); tmod <- as.numeric(tmod)      densprop <- array(1,dim=c(1,30))      densprop <- as.numeric(densprop)      fdeltafe <- array(1,dim=c(1,30))      fdeltafe <- as.numeric(fdeltafe)      fdeltasc <- array(1,dim=c(1,30))     fdeltasc <- as.numeric(fdeltasc)      fdeltasj <- array(1,dim=c(1,30))     fdeltasj <- as.numeric(fdeltasj)      fdeltasa <- array(1,dim=c(1,30))     fdeltasa <- as.numeric(fdeltasa)      ng <- matrix( 0, nrow= 3, ncol=30)  # assume constant carrying capacity (k= 60000 individuals)     carcap= 60000  # n0 initial population vector # multiplied 2 represewnt both sexes # =============================================== # transfer guanacos (g) vector 3 age classes      n0 <- g      tmod[1] <- (n0[1]+n0[2]+n0[3]) * 2  # declaration of initial simulation conditions # ================================================ # ng number of female individuals per age class (dim 3x30) # tmod total (both sexes) population (sum of 3 age classes * 2)      ng[,1] <- n0   ad1 <- 1.195680167 ad2 <- 1.127219245 ad3 <- 1.113739384 ad4 <- 1.320456815   fit <- optim ( c(10.21559509,                    9.80201883,                    9.760834107,                   10.59390027), fn=guafit, method='bfgs', hessian=true)  fit 

here output first approach:

$par [1] 11.315899 11.886347 11.912239  9.675885  $value [1] 1177381086  $counts function gradient       150       17   $convergence [1] 0  $message null  $hessian          [,1]      [,2]      [,3]      [,4] [1,] 417100.8  306371.2  326483.2  941186.8 [2,] 306371.2  516636.4  370602.5 1061540.5 [3,] 326483.2  370602.5  577929.7 1135695.9 [4,] 941186.8 1061540.5 1135695.9 3720506.4 

here 4 parameter estimates , standard errors:

           mle          se [1,] 11.315899 0.002439888 [2,] 11.886347 0.002240669 [3,] 11.912239 0.002153876 [4,]  9.675885 0.001042035 

here r code second approaches, estimates 8 parameters. se not estimated second approach:

set.seed(2345)  guafit <- function(param) {            ad1 <- param[1]           ad2 <- param[2]           ad3 <- param[3]           ad4 <- param[4]           bd1 <- param[5]           bd2 <- param[6]           bd3 <- param[7]           bd4 <- param[8]  # start simulation 30 years     for(i in 1:29) {  # set density-dependent coefficients     densprop[i] <- tmod[i] / carcap  # calculate density-dependent factors     fdeltafe[i]= 1/(1+exp((densprop[i]-ad1)*bd1))     fdeltasc[i]= 1/(1+exp((densprop[i]-ad2)*bd2))     fdeltasj[i]= 1/(1+exp((densprop[i]-ad3)*bd3))     fdeltasa[i]= 1/(1+exp((densprop[i]-ad4)*bd4))  # apply density-dependent factors each coefficient in own age class     a[1,3,i]= a[1,3,i] * fdeltafe[i]     a[2,1,i]= a[2,1,i] * fdeltasc[i]     a[3,2,i]= a[3,2,i] * fdeltasj[i]     a[3,3,i]= a[3,3,i] * fdeltasa[i]  # project total population matrix operation      ng[,i+1]  <- a[,,i]%*%ng[,i]     tmod[i+1] <- (ng[1,i+1] + ng[2,i+1] + ng[3,i+1]) * 2  # end of 30-years simulation loop             }  # calculate residual sum of squares (ssq)                    ssq = sum((tmod - tfield)^2)  # end of function guafit         }  # create field population (a 30-years time series) # ====================================================     tfield <- c(12334,10670,19078,11219,11771,12323,13027,14094,14604,17775,                 20774,16410,17626,21445,21111,20777,28978,27809,28935,38841,                 38363,32273,43128,58597,52456,33125,61334,60488,44773,56973)  # initialize conditions , parameters before calling function guafit     ssq <- 0  # assign initial values transition matrix coefficients # ============================================================     a21 = 0.6     a32 = 0.8     a33 = 0.9     a13 = 0.37  # assign values vector of initial population # =====================================================     g <- matrix(c(0,0,0),ncol=3, nrow=1)     g[1] = 1542     g[2] = 740     g[3] = 3885  # load matrices initial values 30 time units (years) # =========================================================================      <- array(0,dim=c(3,3,30))      for(i in 1:29) {     a[1,3,i]= a13     a[2,1,i]= a21     a[3,2,i]= a32     a[3,3,i]= a33             }  # initialize  variables # ==========================     tmod<-array(1,dim=c(1,30)); tmod <- as.numeric(tmod)      densprop <- array(1,dim=c(1,30))      densprop <- as.numeric(densprop)      fdeltafe <- array(1,dim=c(1,30))      fdeltafe <- as.numeric(fdeltafe)      fdeltasc <- array(1,dim=c(1,30))     fdeltasc <- as.numeric(fdeltasc)      fdeltasj <- array(1,dim=c(1,30))     fdeltasj <- as.numeric(fdeltasj)      fdeltasa <- array(1,dim=c(1,30))     fdeltasa <- as.numeric(fdeltasa)      ng <- matrix( 0, nrow= 3, ncol=30)  # assume constant carrying capacity (k= 60000 individuals)     carcap= 60000  # n0 initial population vector # multiplied 2 represewnt both sexes # =============================================== # transfer guanacos (g) vector 3 age classes      n0 <- g      tmod[1] <- (n0[1]+n0[2]+n0[3]) * 2  # declaration of initial simulation conditions # ================================================ # ng number of female individuals per age class (dim 3x30) # tmod total (both sexes) population (sum of 3 age classes * 2)      ng[,1] <- n0   fit <- optim ( c( 1.195680167,                    1.127219245,                    1.113739384,                    1.320456815,                   10.21559509,                    9.80201883,                    9.760834107,                   10.59390027), fn=guafit, method='bfgs', hessian=true)   fit 

here output second approach:

$par [1]   0.9191288   1.0266562   1.1754382   0.9973042 248.5229404  75.2993413 380.2445394 460.2091730  $value [1] 1031918725  $counts function gradient       705       76   $convergence [1] 0  $message null  $hessian              [,1]          [,2] [,3]        [,4]          [,5]          [,6] [,7]        [,8] [1,] 7.688038e+09  4.293539e+07    0  2.89082527  6.248865e+04  4.860839e+04    0  0.02980232 [2,] 4.293539e+07  2.761336e+06    0  0.29802322 -3.247231e+03  1.245818e+04    0 -0.02980232 [3,] 0.000000e+00  0.000000e+00    0  0.00000000  0.000000e+00  0.000000e+00    0  0.00000000 [4,] 2.890825e+00  2.980232e-01    0 55.25350571 -2.980232e-02  0.000000e+00    0  0.02980232 [5,] 6.248865e+04 -3.247231e+03    0 -0.02980232  8.401275e+01 -3.635883e+00    0 -0.02980232 [6,] 4.860839e+04  1.245818e+04    0  0.00000000 -3.635883e+00  3.185868e+01    0 -0.02980232 [7,] 0.000000e+00  0.000000e+00    0  0.00000000  0.000000e+00  0.000000e+00    0  0.00000000 [8,] 2.980232e-02 -2.980232e-02    0  0.02980232 -2.980232e-02 -2.980232e-02    0  0.02980232