Setting up a mizer model that agrees with observations used to be difficult. That is not a surprise, because we have seen how all the species influence each other, as well as the resource, and how the reproduction rates of all species depend on the size spectra of all species and vice versa. So if you make changes at one end of the model to make it agree with some observation, things change at the other end to mess things up again.
There are three dynamic processes in action in a mizer model at the same time: size-spectrum dynamics, reproduction dynamics and resource dynamics. These are fully interacting with each other, with feedback loops. So for example the resource spectrum depends on the consumption by fish, which depends on the fish size spectra, which depend on the fish growth rates, which depends on the resource spectrum. Similarly the reproduction rate depends on the number of mature fish and on their energy income, which depends among other things on the rate at which new individuals are recruited, which depends among other things on the reproduction rate. And all of these processes depend on the model parameters that we are supposed to choose in a way that reproduces observed behaviour. It seems hopeless!
The way we have arrived at a simple process for the creation of a viable mizer model is to decouple the tuning of the size spectrum dynamics from the tuning of the reproduction dynamics and resource dynamics. So, as we have done last week, initially we turn off reproduction dynamics and resource dynamics. We set the constant reproduction rate to a level that produces the observed biomasses and we set the constant resource spectrum according to observations or, in the absence of observations, we set it to a Sheldon power law. We then use that size spectrum dynamics on its own quickly settles down to its steady state, so that we can interactively tune parameters to get the steady state to agree with observations.
Once we are happy with the steady state of the model, we turn the reproduction and resource dynamics back on, but with parameter choices that do not modify the steady state of the coupled system. We then have to tune the remaining parameters of the reproduction dynamics and resource dynamics to achieve the correct sensitivity of the system to perturbations away from its steady state. But by separating this tuning of the dynamics from the tuning of the steady state, the whole process becomes much more manageable.
We will concentrate on building models with the correct steady state this week and then tune the behaviour away from steady state next week.
In this tutorial we will take the species parameters that we assembled in the previous tutorial for the Curonian lagoon and use the newMultispeciesParams() function to build a mizer model with them. We will let mizer choose most of the defaults and then adjust a few things so that the model has a steady state that has the observed species biomasses and growth curves. We will not do any fine-tuning of the steady state because the first big improvement that we will make in the next tutorial will be to add multiple resource spectra, as is appropriate for this shallow coastal ecosystem.
Warning: cannot remove prior installation of package 'rlang'
Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying C:
\Users\Gustav\AppData\Local\R\win-library\4.2\00LOCK\rlang\libs\x64\rlang.dll
to C:\Users\Gustav\AppData\Local\R\win-library\4.2\rlang\libs\x64\rlang.dll:
Permission denied
Warning: restored 'rlang'
package 'htmltools' successfully unpacked and MD5 sums checked
Warning: cannot remove prior installation of package 'htmltools'
Warning in file.copy(savedcopy, lib, recursive = TRUE):
problem copying C:\Users\Gustav\AppData\Local\R\win-
library\4.2\00LOCK\htmltools\libs\x64\htmltools.dll to C:
\Users\Gustav\AppData\Local\R\win-library\4.2\htmltools\libs\x64\htmltools.dll:
Permission denied
Warning: restored 'htmltools'
package 'Rcpp' successfully unpacked and MD5 sums checked
package 'fontawesome' successfully unpacked and MD5 sums checked
package 'pillar' successfully unpacked and MD5 sums checked
package 'tibble' successfully unpacked and MD5 sums checked
package 'farver' successfully unpacked and MD5 sums checked
package 'stringi' successfully unpacked and MD5 sums checked
Warning: cannot remove prior installation of package 'stringi'
Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying C:
\Users\Gustav\AppData\Local\R\win-library\4.2\00LOCK\stringi\libs\icudt69l.dat
to C:\Users\Gustav\AppData\Local\R\win-library\4.2\stringi\libs\icudt69l.dat:
Invalid argument
Warning in file.copy(savedcopy, lib, recursive = TRUE):
problem copying C:\Users\Gustav\AppData\Local\R\win-
library\4.2\00LOCK\stringi\libs\x64\stringi.dll to C:
\Users\Gustav\AppData\Local\R\win-library\4.2\stringi\libs\x64\stringi.dll:
Permission denied
Warning: restored 'stringi'
package 'shiny' successfully unpacked and MD5 sums checked
package 'deSolve' successfully unpacked and MD5 sums checked
package 'rintrojs' successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\Gustav\AppData\Local\Temp\Rtmp6fMkfA\downloaded_packages
* checking for file 'C:\Users\Gustav\AppData\Local\Temp\Rtmp6fMkfA\remotes33c865fe11d4\sizespectrum-mizerExperimental-f2d9e2f/DESCRIPTION' ... OK
* preparing 'mizerExperimental':
* checking DESCRIPTION meta-information ... OK
* checking for LF line-endings in source and make files and shell scripts
* checking for empty or unneeded directories
* building 'mizerExperimental_2.3.1.tar.gz'
We collected the species parameters and the interaction matrix for our model in the previous tutorial and saved them as .csv files. So here we only need to read in those files.
We use rownames = 1 to let read.csv know that the first column in the spreadsheet, which contains the predator names, should be used as the names of the rows of the interaction matrix.
If you want to check that the data was read in correctly, you can click on curonian_species_params and curonian_interaction in the “Environment” tab. That will open the data frames in your editor window.
We will now set up a multi-species mizer model using the function newMultispeciesParams(). Besides the species parameters and the interaction matrix, the other information that flows into a multi-species model are the resource parameters, the gear parameters and the fishing effort.
We let mizer choose defaults for the resource parameters and gear parameters. By default, the resource carrying capacity will be set to a power law N_R(w) = \kappa w^{-\lambda} with \lambda = 2.05, as we are already familiar with from week 1.
The default fishing gear selects all individuals above maturity size of all species with a catchability of 1. So if we set the fishing effort to 0.3 that imposes a fishing mortality of 0.3 per year on those fish.
Last week we discussed that natural choice for the allometric exponents n (growth exponent) and p (metabolic exponent) is to take them both equal to 3/4. That is the choice we will use.
With this information we call the function newMultispeciesParams() which returns a MizerParams object that we save in the variable cur_model (lazy shorthand for “Curonian model”):
cur_model<-newMultispeciesParams(species_params =curonian_species_params, interaction =curonian_interaction, initial_effort =0.3, lambda =2.05, n =3/4, p =3/4)
No h provided for some species, so using f0 and k_vb to calculate it.
No ks column so calculating from critical feeding level.
Using z0 = z0pre * w_inf ^ z0exp for missing z0 values.
As you see from the messages it has printed, the newMultispeciesParams() function has made choices for some model parameters based on the information we supplied.
So for example it chose the coefficient h of the maximum intake rate for each species to produce growth curves that are roughly in agreement with the von Bertalannfy growth parameters k_vb that we supplied and with a larval feeding level f0 (which we did not supply but which mizer chose as f0 = 0.6).
Mizer chose the coefficient ks of the metabolic rate for each species so that a proportion fc of the maximum intake rate would be enough to cover the metabolic expenses. We did not specifiy this critical feeding level fc so mizer chose fc = 0.2.
Finally mizer has set a constant background mortality z0 that scales with the species’ asymptotic size. Let’s look at the result:
There is a lot wrong here. The species spectra lack the characteristic bulge at adult sizes. Also the species spectra do not line up nicely with the abundance of the resource. But most importantly, these spectra are not close to their steady state values.
We will now project to the steady state, which will finally give us realistic species spectra. To do this we use the special function steady() which implements our trick of keeping the reproduction rate and the resource spectrum constant.
cur_model<-steady(cur_model)
Convergence was achieved in 15 years.
Warning in setBevertonHolt(params, reproduction_level = old_reproduction_level):
The following species require an unrealistic reproductive efficiency greater
than 1: roach
We can ignore the warning from the setBevertonHolt() function about unrealistic reproductive efficiencies. Those warnings are an artefact of how the reproduction level is set by default by the matchBiomasses() function. We could fix those defaults, but we are not yet concerned with the reproduction dynamics so we don’t have to do that and just ignore the warnings.
Now let us look at the spectra in the steady state:
The steady() function is not guaranteed to find the steady state. It may be that the steady state is actually unstable. In that case the system evolves towards an oscillating state instead. This is luckily rare for realistic parameters, but may well happen while you are still trying to find the correct parameters. Let us illustrate this in a modified model. We decrease the width of the feeding kernel for all species to sigma = 1 and increase the slope of the resource spectrum to lambda = 2.2. We know from mathematical studies that both of these changes have a destabilising effect.
modified_species_params<-curonian_species_paramsmodified_species_params$sigma<-1modified_model<-newMultispeciesParams(species_params =modified_species_params, interaction =curonian_interaction, initial_effort =0.3, lambda =2.2, n =3/4, p =3/4)
No h provided for some species, so using f0 and k_vb to calculate it.
No ks column so calculating from critical feeding level.
Using z0 = z0pre * w_inf ^ z0exp for missing z0 values.
Now let us see what happens when we try to run this system to steady state:
modified_model<-steady(modified_model)
Simulation run did not converge after 99 years. Value returned by the distance function was: 0.759089044351571
sim<-steady(modified_model, return_sim =TRUE)
Simulation run did not converge after 99 years. Value returned by the distance function was: 1.7392491172154
Mizer is agnostic of whether you want to measure biomass per square meter, or per square kilometer, or for the entire area of the fishery or whatever. So initially it had set things up on an arbitrarily chosen scale. We can see this if we compare the biomasses of the species in the model with the observed biomasses from your species parameter file with the plotBiomassVsSpecies() function:
This shows for each species the model biomass (open circle) and the observed biomass (filled square) on a logarithmic y-axis. The line connecting model value and observed value is for visual purposes only. We see that model values and observed values are many orders of magnitude apart.
Using your supplied biomass observations, mizer can now calibrate its scale so that the total biomass in the model coincides with the total observed biomass, summed over all species.
Of course for the individual species the model biomasses will still disagree with the observed biomasses, with some being too high while others are too low. Just the total summed over all species agrees between model and observation.
To fix the discrepancy between the model biomasses and the observed biomasses we simply need to rescale the abundance of each species by the appropriate factor. The matchBiomasses() function does this for us.
Now the circles and squares lie on top of each other.
There are similar functions matchNumbers() and matchYields() that you would use in case either total numbers of individuals or fisheries yields are known instead of total biomasses.
Step 5: Rinse and repeat
Rescaling the spectra of the individual species has not created another steady state. All species now experience a new prey distribution and a new predator distribution and so their growth and death rates have changed, which requires us to run the dynamics again to find the new steady state. So we essentially go back to step 2 and project to steady state:
cur_model<-steady(cur_model)
Simulation run did not converge after 99 years. Value returned by the distance function was: 0.182533427975559
The function tell us that running the dynamics for 99 years was not enough to converge sufficiently close to the steady state. In other words, the state was still changing after 99 years. There are now two possibilities: 1) the dynamics will never converge to a steady state but instead continue to oscillate for ever or 2) we simply need to be a bit more patient. Let’s hope for the second option and call steady() again.
cur_model<-steady(cur_model)
Simulation run did not converge after 99 years. Value returned by the distance function was: 0.0280180720197587
We were lucky and only had to wait for another 13.5 years. We’ll discuss what to do when we are not lucky in the When things don’t go smoothly section.
But while the system settled down to the new steady state, the species biomasses have changed, so they now no longer agree with observation.
We appear to be in a bind: If we match the biomasses we are no longer at steady state, if we run to steady state we no longer match the biomasses. But notice that the discrepancies are not as big as previously. So we don’t give up but simply keep iterating.
After repeating the cycle of matching biomasses and running to steady state 10 more times we get this picture:
Warning in setBevertonHolt(params): For the following species `erepro` has been
increased to the smallest possible value: erepro[predator_fish] = 3.11e-06
Convergence was achieved in 58.5 years.
Warning in setBevertonHolt(params): For the following species `erepro` has
been increased to the smallest possible value: erepro[ruffe] = 0.00582;
erepro[perch] = 0.000903; erepro[pikeperch] = 3.52e-05; erepro[burbot] = 4e-05;
erepro[predator_fish] = 8.99e-07
Convergence was achieved in 30 years.
Warning in setBevertonHolt(params): For the following species `erepro` has been
increased to the smallest possible value: erepro[ruffe] = 0.00128; erepro[perch]
= 0.000189; erepro[pikeperch] = 1.09e-05; erepro[burbot] = 1.2e-05;
erepro[predator_fish] = 4.88e-07
Convergence was achieved in 22.5 years.
Warning in setBevertonHolt(params): For the following species `erepro` has
been increased to the smallest possible value: erepro[ruffe] = 0.000507;
erepro[perch] = 9.73e-05; erepro[pikeperch] = 7.46e-06; erepro[burbot] =
8.38e-06; erepro[predator_fish] = 4.36e-07
Convergence was achieved in 16.5 years.
Warning in setBevertonHolt(params): For the following species `erepro` has
been increased to the smallest possible value: erepro[ruffe] = 0.000375;
erepro[perch] = 8.15e-05; erepro[pikeperch] = 7.33e-06; erepro[burbot] =
8.36e-06; erepro[predator_fish] = 4.86e-07
Convergence was achieved in 12 years.
Warning in setBevertonHolt(params): For the following species `erepro` has been
increased to the smallest possible value: erepro[ruffe] = 0.00036; erepro[perch]
= 8.18e-05; erepro[pikeperch] = 7.68e-06; erepro[burbot] = 8.8e-06;
erepro[predator_fish] = 5.28e-07
Convergence was achieved in 9 years.
Convergence was achieved in 7.5 years.
Convergence was achieved in 6 years.
Convergence was achieved in 1.5 years.
Convergence was achieved in 1.5 years.
Now use your species parameter data frame that you assembled in the exercise in the previous tutorial. Go through the 5 steps that we went through above to build your own mizer model based on your species parameters.
Use the worksheet for this tutorial to submit your result.
