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]] ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help