# [R] Excel can do what R can't?????

Michael Rennie mrennie at utm.utoronto.ca
Thu Jul 17 02:51:23 CEST 2003

```Hmmm.

I tried entering 'Hgtmod = Hgt' at the end of my 'optim' function, but that
didn't help me any- still getting poor optimizations.  Perhaps this isn't
working to set a condition, as I was hoping it to.  I think that if I can set
the condition Hgmod = Hgt, then it should be able to find reasonable solutions
to this problem set, since that seems to be the trick in my Excel 'solver'
function.

Also, I'm a little hesitant to simplify this too much in terms of reducing the
model, becuase I need this thing to work over 365 iterations, as it does in
excel.  I've at least cleaned up some of the commenting so it should be a bit
more straightforward.  I tried making the arrays more clear with tabs, but lost
them upon pasting them into this file.

I've included below the code I am currently using with my temp.dat file
attached (it's also below so someone can copy and paste it into a text file
named 'temp.dat')- that's all anyone should need to play around with this if
they are feeling so inclined by cutting and pasting into R.

Thanks again,

Mike

#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

#Bring in temp file

temper <- scan("temp.dat", na.strings = ".", list(Day=0, jday=0, Temp=0))

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 <- 1 #  0.558626306252032
ACT <- 2 #  1.66764519286918

q<-c(p,ACT)

#introduce function to solve
f <- function (q, Hgtmod)
{

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

C[i]<- q[1]*CA*(W[i]^CB)*((comp[i,4])^Xc)*(exp(Xc*(1-(comp[i,4]))))*Pc

ASMR[i]<- q[2]*RA*(W[i]^RB)*((comp[i,5])^Xa)*(exp(Xa*(1-(comp[i,5]))))

SMR[i]<- (ASMR[i]/q[2])

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)

}

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"))

Wtmod<- bioday [length(W),2]
Wtmod

Hgtmod<- bioday [length(Hg),14]
Hgtmod

q

f <- 1000000000*((((Wt-Wtmod)^2)/Wt) + (((Hgt-Hgtmod)^2)/Hgt)) ; f
}

optim(q, f, method = "L-BFGS-B",
lower = c(0.2, 2), upper=c(2, 3),
Hgtmod = Hgt)

#-----------------------------

Temp.dat:

1       153     9.4
2       154     9.6
3       155     9.8
4       156     10
5       157     10.2
6       158     10.4
7       159     10.6
8       160     10.8
9       161     11
10      162     11.2
11      163     11.4
12      164     11.6
13      165     11.8
14      166     12
15      167     12.3
16      168     12.5
17      169     12.7
18      170     12.9
19      171     13.1
20      172     13.4
21      173     13.6
22      174     13.8
23      175     14
24      176     14.2
25      177     14.5
26      178     14.7
27      179     14.9
28      180     15.1
29      181     15.4
30      182     15.6
31      183     15.8
32      184     16
33      185     16.2
34      186     16.5
35      187     16.7
36      188     16.9
37      189     17.1
38      190     17.3
39      191     17.5
40      192     17.7
41      193     17.9
42      194     18.1
43      195     18.3
44      196     18.5
45      197     18.7
46      198     18.9
47      199     19
48      200     19.2
49      201     19.4
50      202     19.5
51      203     19.7
52      204     19.9
53      205     20
54      206     20.2
55      207     20.3
56      208     20.4
57      209     20.5
58      210     20.7
59      211     20.8
60      212     20.9
61      213     21
62      214     21.1
63      215     21.2
64      216     21.3
65      217     21.3
66      218     21.4
67      219     21.5
68      220     21.5
69      221     21.6
70      222     21.6
71      223     21.6
72      224     21.7
73      225     21.7
74      226     21.7
75      227     21.7
76      228     21.7
77      229     21.7
78      230     21.7
79      231     21.6
80      232     21.6
81      233     21.6
82      234     21.5
83      235     21.5
84      236     21.4
85      237     21.3
86      238     21.3
87      239     21.2
88      240     21.1
89      241     21
90      242     20.9
91      243     20.8
92      244     20.7
93      245     20.5
94      246     20.4
95      247     20.3
96      248     20.2
97      249     20
98      250     19.9
99      251     19.7
100     252     19.5
101     253     19.4
102     254     19.2
103     255     19
104     256     18.9
105     257     18.7
106     258     18.5
107     259     18.3
108     260     18.1
109     261     17.9
110     262     17.7
111     263     17.5
112     264     17.3
113     265     17.1
114     266     16.9
115     267     16.7
116     268     16.5
117     269     16.2
118     270     16
119     271     15.8
120     272     15.6
121     273     15.4
122     274     15.1
123     275     14.9
124     276     14.7
125     277     14.5
126     278     14.2
127     279     14
128     280     13.8
129     281     13.6
130     282     13.4
131     283     13.1
132     284     12.9
133     285     12.7
134     286     12.5
135     287     12.3
136     288     12
137     289     11.8
138     290     11.6
139     291     11.4
140     292     11.2
141     293     11
142     294     10.8
143     295     10.6
144     296     10.4
145     297     10.2
146     298     10
147     299     9.8
148     300     9.6
149     301     9.4
150     302     9.3
151     303     9.1
152     304     8.9
153     305     8.7
154     306     8.6
155     307     8.4
156     308     8.2
157     309     8.1
158     310     7.9
159     311     7.8
160     312     7.6
161     313     7.5
162     314     7.3
163     315     7.2
164     316     7
165     317     6.9
166     318     6.8
167     319     6.7
168     320     6.5
169     321     6.4
170     322     6.3
171     323     6.2
172     324     6.1
173     325     6
174     326     5.8
175     327     5.7
176     328     5.6
177     329     5.5
178     330     5.5
179     331     5.4
180     332     5.3
181     333     5.2
182     334     5.1
183     335     5
184     336     5
185     337     4.9
186     338     4.8
187     339     4.7
188     340     4.7
189     341     4.6
190     342     4.5
191     343     4.5
192     344     4.4
193     345     4.4
194     346     4.3
195     347     4.3
196     348     4.2
197     349     4.2
198     350     4.1
199     351     4.1
200     352     4
201     353     4
202     354     4
203     355     3.9
204     356     3.9
205     357     3.8
206     358     3.8
207     359     3.8
208     360     3.8
209     361     3.7
210     362     3.7
211     363     3.7
212     364     3.6
213     365     3.6
214     366     3.6
215     1       3.2
216     2       3.2
217     3       3.2
218     4       3.2
219     5       3.2
220     6       3.2
221     7       3.2
222     8       3.2
223     9       3.2
224     10      3.2
225     11      3.2
226     12      3.2
227     13      3.2
228     14      3.2
229     15      3.2
230     16      3.2
231     17      3.2
232     18      3.2
233     19      3.2
234     20      3.2
235     21      3.2
236     22      3.2
237     23      3.2
238     24      3.2
239     25      3.2
240     26      3.2
241     27      3.2
242     28      3.2
243     29      3.2
244     30      3.2
245     31      3.2
246     32      3.2
247     33      3.2
248     34      3.2
249     35      3.2
250     36      3.2
251     37      3.2
252     38      3.2
253     39      3.2
254     40      3.2
255     41      3.2
256     42      3.2
257     43      3.2
258     44      3.2
259     45      3.2
260     46      3.2
261     47      3.2
262     48      3.2
263     49      3.2
264     50      3.2
265     51      3.2
266     52      3.2
267     53      3.2
268     54      3.3
269     55      3.3
270     56      3.3
271     57      3.3
272     58      3.3
273     59      3.3
274     60      3.3
275     61      3.3
276     62      3.3
277     63      3.3
278     64      3.3
279     65      3.3
280     66      3.3
281     67      3.3
282     68      3.3
283     69      3.3
284     70      3.3
285     71      3.4
286     72      3.4
287     73      3.4
288     74      3.4
289     75      3.4
290     76      3.4
291     77      3.4
292     78      3.4
293     79      3.5
294     80      3.5
295     81      3.5
296     82      3.5
297     83      3.5
298     84      3.5
299     85      3.6
300     86      3.6
301     87      3.6
302     88      3.6
303     89      3.6
304     90      3.7
305     91      3.7
306     92      3.7
307     93      3.8
308     94      3.8
309     95      3.8
310     96      3.8
311     97      3.9
312     98      3.9
313     99      4
314     100     4
315     101     4
316     102     4.1
317     103     4.1
318     104     4.2
319     105     4.2
320     106     4.3
321     107     4.3
322     108     4.4
323     109     4.4
324     110     4.5
325     111     4.5
326     112     4.6
327     113     4.7
328     114     4.7
329     115     4.8
330     116     4.9
331     117     5
332     118     5
333     119     5.1
334     120     5.2
335     121     5.3
336     122     5.4
337     123     5.5
338     124     5.5
339     125     5.6
340     126     5.7
341     127     5.8
342     128     6
343     129     6.1
344     130     6.2
345     131     6.3
346     132     6.4
347     133     6.5
348     134     6.7
349     135     6.8
350     136     6.9
351     137     7
352     138     7.2
353     139     7.3
354     140     7.5
355     141     7.6
356     142     7.8
357     143     7.9
358     144     8.1
359     145     8.2
360     146     8.4
361     147     8.6
362     148     8.7
363     149     8.9
364     150     9.1
365     151     9.3
366     152     9.3

--
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
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: temp.dat
Url: https://stat.ethz.ch/pipermail/r-help/attachments/20030716/1285b704/temp.pl
```