Comparing quantities between mizer and stock-assessment models

[michaelspence]

Is there a problem with how yield is calculated? I think it could be inconsistent with how fishing mortality works when the model is not in a stationary state.

In the McKendrick-von Foerster equation \mu is the mortality rate (per time) and is the sum fishing mortality and all other mortality. These are rate parameters and not proportions but when the yield is calculated using the getYield() function, the fishing mortality is used similar to a harvest rate. To demonstrate look at the McKendrick-von Foerster equation when growth is zero. This would be
dN/dt = -\mu N
(I could be wrong but the partial derivative on t becomes a full derivative as it is no longer a function of w). The solution to this model is
N_{t+1} = exp(-\mu)N_{t}
and the catch is calculated by Baranov equation
C = (F/\mu) * (1 - exp(-\mu))N_t,
where F is the fishing mortality.

However, in size-spectrum models, including mizer (with getYield()), catch (or yield) is calculated by
C=\int{} F*N(w)*wdw
or
C=\int{} F * B(w)dw
(equation 5.7 in Ken's recent book). This treats F more like a harvest ratio.

I could be wrong, but in size-spectrum models, including mizer, the dynamics have fishing mortality as a rate, whereas the calculation of yield is done using a harvest ratio. This means that the numbers are not consistent with one another.

I have created an illustrative example attached. mizer_yield_issue.txt

I created the community setting the selectivity to 1 for all sizes, all interactions to 0, kappa to 1e20 and background mortality to 0. This means that all species at all sizes experience the same fishing mortality, which is the only mortality that they experience, hence any fish that die are caught. kappa is set large so that all fish are fully fed. First I run the simulation for 50 years with no effort, saving the final years densities sim@n and sim@n_pp.

I turned off growth (but this still works if its left on) by running params@psi[,] <- 1 and then run forward 1 year with effort = 1. I found getYield(sim1) == getBiomass(sim1) but getBiomass(sim1) != 0, which what would be expected as the whole biomass has just been fished out and there's no growth. plot(sim1) shows that there are still large individuals in the system.

I then run forward one year without recruitment and calculated the fishing rate by comparing the numbers before num1 and after num2 and the total catch catch with new_biomass, the biomass at the current time.

In mizer for a single time step, e.g. dt=1, the proportion catch/new_biomass==effort, that is numbers of dead (in weight) divided by the new biomass of the stock is equal to the F. However the fishing mortality applied to the fish is not equal to effort. For example, I set dt=1 and effort=1, then catch/new_biomass==1, so catch would return the same value as getYield(), but the mortality rate was 0.693. When dt=0.001, catch/new_biomass==1.716924, very different to the getYield() function but the mortality rate was 0.9995003, almost 1.

newF is smaller that r$mort[1] * dt, although as dt gets smaller and is run more times newF approaches r$mort[1].

A quick fix would be to run getYield() on every time-step, this would require t_save=1/dt and a slight change to getYield, by including a t_save arguement:
getYield <- function (sim,t_save) {
assert_that(is(sim, "MizerSim"))
biomass <- sweep(sim@n, 3, sim@params@w * sim@params@dw,
"*")
f <- getFMort(sim, drop = FALSE) * t_save
yield <- apply(f * biomass, c(1, 2), sum)
return(yield)
}
Alternatively I wonder if there is a Baranov equation equivalent equation for the McKendrick-von Foerster equation?

NB. dt is changed on line 44. Also, I'm not sure this is an issue if the community are in an equilibrium state.

I could be talking complete nonsense here like so sorry if I am.

[gustavdelius]

You are right that what getYield() calculates is the rate at which biomass is fished. The biomass fished during a particular time period is the integral of this rate over that time period. That is not made very clear in the documentation. That should be improved.

I am not sure I understand what you mean by "Is there a problem with how yield is calculated? I think it could be inconsistent with how fishing mortality works when the model is not in a stationary state". Is the issue that the term "Yield" is usually used for the biomass fished during a year, not for the instantaneous rate? You are right that the two would agree numerically only when the rate is constant throughout the year, i.e., in a steady state. However they would still have different dimensions, one being biomass per year and the other being just biomass. Is there a better term than "yield" for the rate at which biomass is fished?

There is of course nothing like a Baranov equation in a proper dynamic model where the growth and mortality rates change over time in response to changes in the abundances. The equations need to be integrated numerically.

[Kenhasteandersen]

Hi Michael, Thanks for your careful analysis. I agree with Gustav: the point here is that Mizer calculated what we could call the “instantaneous yield rate”, that is the rate of biomass extraction. That is the correct way to calculate yield, in my view. Age-based fisheries models operates with annual time steps. They therefore need to correct for the exponential decline in biomass due to fishing within that year. That’s where Baranov’s catch equation comes in. We have small time steps and don’t need that correction. Are you concerned about the calculation when the model is not in a steady state (ie one an oscillating solution)? In that case one should indeed average the yield over some time period, for example over the oscillating period to get a better average estimate. Hope this settles your concern, else please write back. Ken

…

--------------------------------------------------------------------------------------------------------------- Ken H. Andersen, http://ken.haste.dk<http://ken.haste.dk/>, twitter: @69kno Professor in theoretical marine ecology, head of section and phd school, and deputy director of Centre for Ocean Life http://www.oceanlifecentre.dk<http://www.oceanlifecentre.dk/> Try my plankton simulator: http://oceanlife.dtuaqua.dk/Plankton/ New book: Fish Ecology, Evolution, and Exploitation - A New Theoretical Synthesis https://press.princeton.edu/titles/13516.html ***@***.*** On 16 Sep 2021, at 23.19, Gustav W Delius ***@***.******@***.***>> wrote: You are right that what getYield() calculates is the rate at which biomass is fished. The biomass fished during a particular time period is the integral of this rate over that time period. That is not made very clear in the documentation. That should be improved. Is the issue that the term "Yield" is usually used for the biomass fished during a year, not for the instantaneous rate? You are right that the two would agree numerically when the rate at which biomass is fished is constant throughout the year, i.e., in a steady state. However they would still have different dimensions, one being biomass per year and the other being just biomass. There is of course nothing like a Baranov equation in a proper dynamic model where the growth and mortality rates change over time in response to changes in the abundances. The equations need to be integrated numerically. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#233 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ADFHLAMXMPOUGSFCKZUX4ZLUCJNNXANCNFSM5EE74C6A>. Triage notifications on the go with GitHub Mobile for iOS<https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android<https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[gustavdelius]

@michaelspence , I have updated the help page for getYield(), see https://sizespectrum.org/mizer/reference/getYield.html. Please feel free to suggest further changes. Or perhaps you want to discuss the issue in a blog post on the mizer blog? The fact that in a continuous-time and continuous -size model one has to work with rates and densities and obtain more conventional quantities as integrals I find quite challenging to communicate to non-mathematical scientists. So your help is very much appreciated.

[michaelspence]

Thanks for the quick reply. I see now thanks!

[michaelspence]

The attached shows that the catch is the biomass at the end of the update times the fishing mortality applied in that update.
getYield.pdf (sorry, this document maybe all over the place as I wrote it when I was trying to understand stuff above)
However the getYield calculates the yield as the biomass at the beginning of the time step times the fishing mortality during that time step. Of course, this disappears as the size of the time step decreases but would it be more accurate if we use the effort from the previous time or set effort=effort[i,] rather than effort=effort[i-1,] in the project_simple call in the project function.

This could have issues when using stock-assessment inputs to drive fishing mortality, which would make sudden jumps at the beginning of the years (something that is not realistic I know). For example if I want to run the model from with fishing mortality rates from 1984 to 2020 I will run the simulation to 2021 but the fishing mortality applied in that final time-step will be that of 2020 but sim@effort["2021",] wouldn't have a value, nor would it be important.

In the example in the getYield function, the actual catch between 2010 and 2011 is sum(getYield(sim)[2:11, "Herring"] / 10), i.e. the 2011 catch is the amount caught between times 2010.9 and 2011, however this was calculated with 2011 fishing mortality, which in this example is the same as 2010 fishing mortality. The difference between the two is very small log(sum(getYield(sim)[1:10, "Herring"] / 10)) - log(sum(getYield(sim)[2:11, "Herring"] / 10)) [1] 0.004297044 but this may not be the case if there were a big difference between fishing mortality between 2010 and 2011.

[gustavdelius]

@michaelspence . Thank you very much for thinking about this. If I understand your latest post correctly, there are three issues discussed in there.

First of all there is the issue of discretisation error coming from using discrete time steps. You are proposing in the linked pdf file to try to adjust for that by assuming an exponential decay of the biomass within a timestep. However, in practice, when there is growth into and out of each size class as well as mortality, the biomass in each size class can either increase or decrease during the time step and thus there is no recipe for automatically reducing the discretisation error, other than decreasing the time step. In practice, the biomass does not change much within a time step of 0.1 years, so I think the discretisation error is not very important. I am actually much more worried about the discretisation errors coming from discretising in the size direction, and there I have been thinking of making corrections using power-law assumptions about the size-dependence of various quantities, very much along the lines you were thinking along. Worth revisiting at some point in the future.

The second issue is about how to interpret the time index in the value returned by getYield(). In the example you mention from the help page for getYield(). The entry labelled "2010" is actually the rate (in grams/year) at which fish biomass is caught during the time step that starts at 2010 and in this example extends to 2010.1 (So roughly the period from 1 January to 5 Feburary 2020). I do not understand why you say that: "however this was calculated with 2011 fishing mortality". Can you explain?

You also raise the issue that the entries in the effort array for the final time is not actually used in a simulation, because the simulation stops at that time. That is correct, but not really a problem.

[michaelspence]

Sorry it has been a while since I thought about this.

I think what I meant was that the catch should be calculated at the end of the time step rather than the beginning of the time-step. In the simplified model, the one with out growth, the catch between time $t$ and time $t+dt$ is $B(t+dt)Fdt$ , where $F$ is the fishing mortality between $t$ and $t+dt$. In an example where $F$ is the same annually, which I have used from the assessments, then the time period 2010-2010.1 has 2010 fishing mortality. In the numerical approximation, then I think it would be better to calculate the yield in this period by by using the 2010.1 biomass multiplied by the F in 2010. Instead mizer does this the other way around. It uses the 2010 biomass (the biomass at the beginning of the time step) and multiplies by the $F$ that hasn't yet happened.

I understand that this is just numerical errors, which will disappear when dt is small, but it is more accurate, and I would say more intuitive to calculate it at the end of the time step. For example, in the example getBiomass(sim)[1, "Herring"] does not depend on the $F$ used in the first time step. For high $F$, specifically $F &gt; 1/dt$, then more biomass is caught in that time step than is in the system.

So going back to the example herring yield is calculated by getYield(sim)[1:10, "Herring"] / 10) but I think it should be getYield(sim)[2:11, "Herring"] / 10), with getYield(sim) calculating the yield rate in that time-step rather than the last. This would mean that getYield was
function (sim) {
assert_that(is(sim, "MizerSim"))
biomass <- sweep(sim@n[-1,,], 3, sim@params@w * sim@params@dw,
"*")
f <- getFMort(sim, drop = FALSE)[-(dim(sim@n)[1]),,]
yield <- apply(f * biomass, c(1, 2), sum)
return(yield)
}

NB. When doing this I actually edited the project code rather that the getYield code so that sim@effort would be correct rather than just dropping the last term from getFMort.

[gustavdelius]

I agree with you that only reporting rates at discrete points in time fails to give much information of what the rates do in between those times. There is nothing one can do about that other than decreasing the time steps. It is important that the rates reported by mizer are actually the correct instantaneous rates. The yield rate for example at a particular time is the product of the biomass at that time and the fishing mortality rate at that time. If I understand correctly, you are proposing to give the product of the biomass at that time and the fishing mortality at the next time point. That however is not the instantaneous yield rate at any point in time.

I also don't see the benefit of using the biomass at the next time point. It seems that you feel that the biomass at the end of a time step is a better representation of the biomass during the time step than the biomass at the start of a time step. But why would that be the case?

If you want a yield rate that is closer to the average yield rate over the period between two neighbouring time points, you could consider taking the average of the yield rates at those two time points. But really, if you have substantial changes in yield rate between time points, it would be better to increase the number of time points.

I also somewhat regret that the function for getting the yield rate is called getYield() instead of getYieldRate(). That may lead people to misinterpret the result from getYield() as the yield obtained between one time point and the next rather than the instantaneous yield rate at one time point. However other functions for getting rates, like getMort() for example, also don't include "Rate" in their names.

[michaelspence]

If I understand correctly, you are proposing to give the product of the biomass at that time and the fishing mortality at the next time point. That however is not the instantaneous yield rate at any point in time.

What you are describing is what currently happens, although I disagree, this is the rate at some point. I propose that the instantaneous yield rate is the product of the biomass at the time and the Fishing mortality of the previous time. I agree this is not the rate at any point but it is for an infinitely small time before hand.

Suppose F(s) to be the fishing mortality at time s with F(s)=F0 for t<=s<t+dt, F(t+dt)=F1, also let B(s) be the biomass time s. The yield rate at time s is therefore F(s)B(s). So the yield rate at time t is B(t)F0 and at t+dt B(t+dt)F1, however for some small value e>0 and e<dt, then the yield rate is B(t+dt-e)F0. If we take the limit as e goes to 0 this becomes B(t+dt)F0, as B(t+dt-e)=B(t+dt), hence that rate happens (could be wrong here but it very nearly happens) at some point.

If I understand correctly, you are proposing to give the product of the biomass at that time and the fishing mortality at the next time point. That however is not the instantaneous yield rate at any point in time.

In "getYield.pdf" (see my post on 22nd September), in the discretized model with no growth and no natural mortality (equation 2), the total biomass lost, which is entirely down to fishing, and is therefore the catch, is
B(t) - B(t+dt) = FB(t+dt)dt,
i.e. the biomass at the end of the time step times the fishing mortality in that time step (see the last equation on the first page). In this model, multiplying the fishing mortality by the biomass at the end of the time step gives exactly the right yield. Although there is growth and other mortality, mizer actually implements this discretized model, therefore I think this is a better representation of the rate of fishing that occurred during that time window.

I also somewhat regret that the function for getting the yield rate is called getYield() instead of getYieldRate(). That may lead people to misinterpret the result from getYield() as the yield obtained between one time point and the next rather than the instantaneous yield rate at one time point. However other functions for getting rates, like getMort() for example, also don't include "Rate" in their names.

Would it be worth including a getCatch() function?