
Hi, there

I am a graduate student new to coding in S who is hitting a bit of a wall 
at present using an "optim" function.  I am running into some troubles, and 
was hoping someone might be able to recognize where I am going wrong.

As background:  I have constructed a loop that carries out a 365-day 
calculation for a mass-balance model. Basically, the model depends on 2 
variables (p, ACT) to solve the 365-day iteration such that the modelled 
estimates of Wt, Ct match the observed estimates.  Each iteration in the 
model depends on the value in the previous iteration. I have this model 
functioning in excel, using "solver" to find values for p, ACT.

The 365-day loop translated into R operates well on it's own- I've used the 
values for p, ACT from my excel model to ensure that I get the same answers 
back, and I do.  (For comparison,  I've provided my output from this file 
in the first block of output below).

However, I have difficulties generating anything when I try to place this 
365-day iteration inside a function that will find the values of p, ACT for 
me using 'optim' (see second output list, below).   I get errors in the 
iteration loop that I don't get when I run the iteration loop on it's 
own.  In fact, the error I get is when I don't specify any starting values 
for p, ACT.  Therefore, this implies that my specification for p, ACT 
starting values in the 'optim' function (c(1,1), eqn, etc....)  isn't 
working.  Any idea why?

I am new to writing code in R so I am wondering if there is a simple fix 
that you guys might be able to see in my code that I am not aware of.  The 
reason I want this in R as opposed to excel is because I need to automate 
this procedure.  Once I get the optimization function working, I will be 
generating distributions from the means of a few of my input variables, and 
getting the program to calculate p, ACT values 1000 times so I can generate 
means and standard errors around my output data, and rbinding this to my 
input variables and values of p, ACT for each iteration.  This last bit is 
something I am somewhat familiar with, so if I can get past the 
optimization stage, I should be good to go.

Any help anyone could offer me would be greatly appreciated- I hope to hear 
from either of you when it is convenient.  The output for the files I have 
is listed below.  Sorry for listing all of my output- I've tried listing 
simplified examples before regarding this question, but have received 
requests to see all the code I am using.

Sincerely,

Mike Rennie

OUTPUT 1- successful program, which successfully goes through 365 days of 
calculations (coding for optim functions are turned off);


R : Copyright 2001, The R Development Core Team
Version 1.4.0  (2001-12-19)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.

R is a collaborative project with many contributors.
Type `contributors()' for more information.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

 > ####################################
 > #    perch.R                       #
 > # Hewett and Johnson bioenergetics #
 > # model combined with              #
 > # Trudel MMBM to estimate          #
 > # Consumption in perch in R code   #
 > # Execute with                     #
 > # R --vanilla < perch.R > perch.out#
 > ####################################
 >
 > #USER INPUT BELOW
 >
 > #Weight at time 0
 > Wo<- 9.2
 >
 > #Hg concentration at time 0 (ugHg/g wet weight)
 > Hgo<- 0.08
 >
 > #Weight at time t
 > Wt<- 32.2
 >
 > #Hg concentration at time t (ugHg/g wet weight)
 > Hgt<- 0.110
 >
 > #Prey methylmercury concentration (as constant)
 > Hgp<- 0.033
 >
 > #Prey caloric value (as constant)
 > Pc<- 800
 >
 > #Energy density of fish (as constant, calories)
 > Ef <- 1000
 >
 > #Maturity status, 0=immature, 1=mature
 > Mat<- 0
 >
 > #Sex, 1=male, 2=female
 > Sex<- 1
 >
 > #USER INPUT ABOVE
 >
 > #Bioenergetics parameters for perch
 > CA <- 0.25
 > CB <- 0.73  #same as 1+(-0.27)- convert g/g/d to g/d * Pc to get cal/d
 > CQ <- 2.3
 > CTO <- 23
 > CTM <- 28
 > Zc<- (log(CQ))*(CTM-CTO)
 > Yc<- (log(CQ))*(CTM-CTO+2)
 > Xc<- ((Zc^2)*(1+(1+40/Yc)^0.5)^2)/400
 >
 > RA <- 34.992  #0.0108*3240 cal/g 02, converting weight of 02 to cal
 > RB <- 0.8   #same as 1+(-0.2) see above...
 > RQ <- 2.1
 > RTO <- 28
 > RTM <- 33
 > Za <- (log(RQ))*(RTM-RTO)
 > Ya<- (log(RQ))*(RTM-RTO+2)
 > Xa<- ((Za^2)*(1+(1+40/Ya)^0.5)^2)/400
 >
 > S <- 0.172
 >
 > FA <- 0.158
 > FB <- -0.222
 > FG <- 0.631
 >
 > UA<- 0.0253
 > UB<- 0.58
 > UG<- -0.299
 >
 > #Mass balance model parameters
 > EA <- 0.002938
 > EB <- -0.2
 > EQ <- 0.066
 > a <- 0.8
 >
 > #Specifying sex-specific parameters
 >
 > GSI<- NULL
 >
 > if (Sex==1) GSI<-0.05 else
+ if (Sex==2) GSI<-0.17
 >
 > # Define margin of error functions
 > #merror <- function(phat,M,alpha) # (1-alpha)*100% merror for a proportion
 > #    {
 > #    z <- qnorm(1-alpha/2)
 > #    merror <- z * sqrt(phat*(1-phat)/M)  # M is (Monte Carlo) sample size
 > #    merror
 > #    }
 >
 > #Bring in temp file
 >
 > temper <- scan("temp.dat", na.strings = ".", list(Day=0, jday=0, Temp=0))
Read 366 records
 >
 > Day<-temper$Day ; jday<-temper$jday ; Temp<-temper$Temp ;
 >
 > temp<- cbind (Day, jday, Temp)
 > #Day = number of days modelled, jday=julian day, Temp = daily avg. temp.
 > #temp [,2]
 >
 > Vc<-(CTM-(temp[,3]))/(CTM-CTO)
 > Vr<-(RTM-(temp[,3]))/(RTM-RTO)
 >
 > comp<- cbind (Day, jday, Temp, Vc, Vr)
 >
 > #comp
 >
 > bio<-matrix(NA, ncol=13, nrow=length(Day))
 > W<-NULL
 > C<-NULL
 > ASMR<-NULL
 > SMR<-NULL
 > A<-NULL
 > F<-NULL
 > U<-NULL
 > SDA<-NULL
 > Gr<-NULL
 > Hg<-NULL
 > Ed<-NULL
 > GHg<-NULL
 > K<-NULL
 > Expegk<-NULL
 > EGK<-NULL
 > p<-NULL
 > ACT<-NULL
 >
 >
 > p <- 0.558626306252032
 > ACT <- 1.66764519286918
 >
 > q<-c(p,ACT)
 >
 > #introduce function to solve
 > #f <- function (q)
 > #{
 >
 > M<- length(Day) #number of days iterated
 >
 > for (i in 1:M)
+ {
+
+ #Bioenergetics model
+
+ if (Day[i]==1) W[i] <- Wo else
+ if (jday[i]==121 && Mat==1) W[i] <- (W[i-1]-(W[i-1]*GSI*1.2)) else
+ W[i] <- (W[i-1]+(Gr[i-1]/Ef))
+ #W
+
+ #W<-Wo
+
+ C[i]<- p*CA*(W[i]^CB)*((comp[i,4])^Xc)*(exp(Xc*(1-(comp[i,4]))))*Pc
+
+ ASMR[i]<- ACT*RA*(W[i]^RB)*((comp[i,5])^Xa)*(exp(Xa*(1-(comp[i,5]))))
+
+ SMR[i]<- (ASMR[i]/ACT)
+
+ A[i]<- (ASMR[i]-SMR[i])
+
+ F[i]<- (FA*((comp[i,3])^FB)*(exp(FG*p))*C[i])
+
+ U[i]<- (UA*((comp[i,3])^UB)*(exp(UG*p))*(C[i]-F[i]))
+
+ SDA[i]<- (S*(C[i]-F[i]))
+
+ Gr[i]<- (C[i]-(ASMR[i]+F[i]+U[i]+SDA[i]))
+
+ #Trudel MMBM
+
+ if (Day[i]==1) Hg[i] <- Hgo else Hg[i] <- 
a*Hgp*(C[i-1]/Pc/W[i-1])/EGK[i-1]*(1-Expegk[i-1])+(Hg[i-1]*Expegk[i-1])
+
+ Ed[i]<- EA*(W[i]^EB)*(exp(EQ*(comp[i,3])))
+
+ GHg[i] <- Gr[i]/Ef/W[i]
+
+ if (Sex==1) 
K[i]<-(((0.1681*(10^(1.3324+(0.000453*Hg[i])))/1000)/Hg[i])*GSI)/M else
+ if (Sex==2) 
K[i]<-(((0.1500*(10^(0.8840+(0.000903*Hg[i])))/1000)/Hg[i])*GSI)/M
+ # = dw/ww conversion * gonad ~ body conc'n function(ng/g) / convert to ug/g
+ # then express as Q times GSI gives K / M gives daily K
+
+ EGK[i] <- (Ed[i] + GHg[i] + (K[i]*Mat))
+
+ Expegk[i] <- exp(-1*EGK[i])
+
+ bio<- cbind(W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)
+
+ }
 >
 > #warnings()
 >
 > dimnames (bio) <-list(NULL, c("W", "C", "ASMR", "SMR", "A", "F", "U", 
"SDA", "Gr", "Ed", "GHg", "EGK", "Hg"))
 >
 > bioday<-cbind(jday, W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)
 >
 > dimnames (bioday) <-list(NULL, c("jday", "W", "C", "ASMR", "SMR", "A", 
"F", "U", "SDA", "Gr", "Ed", "GHg", "EGK", "Hg"))
 >
 > #bioday
 >
 > Wtmod<- bioday [length(W),2]
 > Wtmod
[1] 32.20906
 >
 > Hgtmod<- bioday [length(Hg),14]
 > Hgtmod
[1] 0.1099862
 >
 > f <- (((Wt-Wtmod)^2 + (Hgt-Hgtmod)^2)/2)^2 ; f
[1] 1.682521e-09
 >
 > q
[1] 0.5586263 1.6676452
 >
 > #warnings()
 >
 > write.table (bioday, file = "perch.csv", append = FALSE, sep=",", na = 
NA, col.names = TRUE)
 >
 > #}
 >
 > #nlm(f,c(1,1))
 >
 > #optim(c(1,1), f, method = "L-BFGS-B",
 > #     lower = c(0, 0), upper=c(2, 10))
 >
 > #warnings()
 >
 > bioday
        jday         W        C     ASMR       SMR         A        F         U
   [1,]  153  9.200000 233.6647 107.5640  64.50050  43.06345 31.93755 15.840142



OUTPUT 2- program that doesn't work, and gets screwed up in the daily 
iterations- cqan't recognize starting conditions for q, even though it 
worked fine before I placed it within the 'optim' function;

R : Copyright 2001, The R Development Core Team
Version 1.4.0  (2001-12-19)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.

R is a collaborative project with many contributors.
Type `contributors()' for more information.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

 > ####################################
 > #    perch.R                       #
 > # Hewett and Johnson bioenergetics #
 > # model combined with              #
 > # Trudel MMBM to estimate          #
 > # Consumption in perch in R code   #
 > # Execute with                     #
 > # R --vanilla < perch.R > perch.out#
 > ####################################
 >
 > #USER INPUT BELOW
 >
 > #Weight at time 0
 > Wo<- 9.2
 >
 > #Hg concentration at time 0 (ugHg/g wet weight)
 > Hgo<- 0.08
 >
 > #Weight at time t
 > Wt<- 32.2
 >
 > #Hg concentration at time t (ugHg/g wet weight)
 > Hgt<- 0.110
 >
 > #Prey methylmercury concentration (as constant)
 > Hgp<- 0.033
 >
 > #Prey caloric value (as constant)
 > Pc<- 800
 >
 > #Energy density of fish (as constant, calories)
 > Ef <- 1000
 >
 > #Maturity status, 0=immature, 1=mature
 > Mat<- 0
 >
 > #Sex, 1=male, 2=female
 > Sex<- 1
 >
 > #USER INPUT ABOVE
 >
 > #Bioenergetics parameters for perch
 > CA <- 0.25
 > CB <- 0.73  #same as 1+(-0.27)- convert g/g/d to g/d * Pc to get cal/d
 > CQ <- 2.3
 > CTO <- 23
 > CTM <- 28
 > Zc<- (log(CQ))*(CTM-CTO)
 > Yc<- (log(CQ))*(CTM-CTO+2)
 > Xc<- ((Zc^2)*(1+(1+40/Yc)^0.5)^2)/400
 >
 > RA <- 34.992  #0.0108*3240 cal/g 02, converting weight of 02 to cal
 > RB <- 0.8   #same as 1+(-0.2) see above...
 > RQ <- 2.1
 > RTO <- 28
 > RTM <- 33
 > Za <- (log(RQ))*(RTM-RTO)
 > Ya<- (log(RQ))*(RTM-RTO+2)
 > Xa<- ((Za^2)*(1+(1+40/Ya)^0.5)^2)/400
 >
 > S <- 0.172
 >
 > FA <- 0.158
 > FB <- -0.222
 > FG <- 0.631
 >
 > UA<- 0.0253
 > UB<- 0.58
 > UG<- -0.299
 >
 > #Mass balance model parameters
 > EA <- 0.002938
 > EB <- -0.2
 > EQ <- 0.066
 > a <- 0.8
 >
 > #Specifying sex-specific parameters
 >
 > GSI<- NULL
 >
 > if (Sex==1) GSI<-0.05 else
+ if (Sex==2) GSI<-0.17
 >
 > # Define margin of error functions
 > #merror <- function(phat,M,alpha) # (1-alpha)*100% merror for a proportion
 > #    {
 > #    z <- qnorm(1-alpha/2)
 > #    merror <- z * sqrt(phat*(1-phat)/M)  # M is (Monte Carlo) sample size
 > #    merror
 > #    }
 >
 > #Bring in temp file
 >
 > temper <- scan("temp.dat", na.strings = ".", list(Day=0, jday=0, Temp=0))
Read 366 records
 >
 > Day<-temper$Day ; jday<-temper$jday ; Temp<-temper$Temp ;
 >
 > temp<- cbind (Day, jday, Temp)
 > #Day = number of days modelled, jday=julian day, Temp = daily avg. temp.
 > #temp [,2]
 >
 > Vc<-(CTM-(temp[,3]))/(CTM-CTO)
 > Vr<-(RTM-(temp[,3]))/(RTM-RTO)
 >
 > comp<- cbind (Day, jday, Temp, Vc, Vr)
 >
 > #comp
 >
 > bio<-matrix(NA, ncol=13, nrow=length(Day))
 > W<-NULL
 > C<-NULL
 > ASMR<-NULL
 > SMR<-NULL
 > A<-NULL
 > F<-NULL
 > U<-NULL
 > SDA<-NULL
 > Gr<-NULL
 > Hg<-NULL
 > Ed<-NULL
 > GHg<-NULL
 > K<-NULL
 > Expegk<-NULL
 > EGK<-NULL
 > p<-NULL
 > ACT<-NULL
 >
 >
 > #p <- 0.558626306252032
 > #ACT <- 1.66764519286918
 >
 > q<-c(p,ACT)
 >
 > #introduce function to solve
 > f <- function (q)
+ {
+
+ M<- length(Day) #number of days iterated
+
+ for (i in 1:M)
+ {
+
+ #Bioenergetics model
+
+ if (Day[i]==1) W[i] <- Wo else
+ if (jday[i]==121 && Mat==1) W[i] <- (W[i-1]-(W[i-1]*GSI*1.2)) else
+ W[i] <- (W[i-1]+(Gr[i-1]/Ef))
+ #W
+
+ #W<-Wo
+
+ C[i]<- p*CA*(W[i]^CB)*((comp[i,4])^Xc)*(exp(Xc*(1-(comp[i,4]))))*Pc
+
+ ASMR[i]<- ACT*RA*(W[i]^RB)*((comp[i,5])^Xa)*(exp(Xa*(1-(comp[i,5]))))
+
+ SMR[i]<- (ASMR[i]/ACT)
+
+ A[i]<- (ASMR[i]-SMR[i])
+
+ F[i]<- (FA*((comp[i,3])^FB)*(exp(FG*p))*C[i])
+
+ U[i]<- (UA*((comp[i,3])^UB)*(exp(UG*p))*(C[i]-F[i]))
+
+ SDA[i]<- (S*(C[i]-F[i]))
+
+ Gr[i]<- (C[i]-(ASMR[i]+F[i]+U[i]+SDA[i]))
+
+ #Trudel MMBM
+
+ if (Day[i]==1) Hg[i] <- Hgo else Hg[i] <- 
a*Hgp*(C[i-1]/Pc/W[i-1])/EGK[i-1]*(1-Expegk[i-1])+(Hg[i-1]*Expegk[i-1])
+
+ Ed[i]<- EA*(W[i]^EB)*(exp(EQ*(comp[i,3])))
+
+ GHg[i] <- Gr[i]/Ef/W[i]
+
+ if (Sex==1) 
K[i]<-(((0.1681*(10^(1.3324+(0.000453*Hg[i])))/1000)/Hg[i])*GSI)/M else
+ if (Sex==2) 
K[i]<-(((0.1500*(10^(0.8840+(0.000903*Hg[i])))/1000)/Hg[i])*GSI)/M
+ # = dw/ww conversion * gonad ~ body conc'n function(ng/g) / convert to ug/g
+ # then express as Q times GSI gives K / M gives daily K
+
+ EGK[i] <- (Ed[i] + GHg[i] + (K[i]*Mat))
+
+ Expegk[i] <- exp(-1*EGK[i])
+
+ bio<- cbind(W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)
+
+ }
+
+ #warnings()
+
+ dimnames (bio) <-list(NULL, c("W", "C", "ASMR", "SMR", "A", "F", "U", 
"SDA", "Gr", "Ed", "GHg", "EGK", "Hg"))
+
+ bioday<-cbind(jday, W, C, ASMR, SMR, A, F, U, SDA, Gr, Ed, GHg, EGK, Hg)
+
+ dimnames (bioday) <-list(NULL, c("jday", "W", "C", "ASMR", "SMR", "A", 
"F", "U", "SDA", "Gr", "Ed", "GHg", "EGK", "Hg"))
+
+ #bioday
+
+ Wtmod<- bioday [length(W),2]
+ Wtmod
+
+ Hgtmod<- bioday [length(Hg),14]
+ Hgtmod
+
+ f <- (((Wt-Wtmod)^2 + (Hgt-Hgtmod)^2)/2)^2 ; f
+
+ q
+
+ #warnings()
+
+ write.table (bioday, file = "perch.csv", append = FALSE, sep=",", na = 
NA, col.names = TRUE)
+
+ }
 >
 > #nlm(f,c(1,1))
 >
 > optim(c(1,1), f, method = "L-BFGS-B",
+       lower = c(0, 0), upper=c(2, 10))
Error in "[<-"(*tmp*, i, value = (W[i - 1] + (Gr[i - 1]/Ef))) :
         nothing to replace with
Execution halted


Michael Rennie
M.Sc. Candidate
University of Toronto at Mississauga
3359 Mississauga Rd. N.
Mississauga, ON  L5L 1C6
Ph: 905-828-5452  Fax: 905-828-3792
	[[alternative HTML version deleted]]


