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]]