Open andrechalom opened 8 years ago
Definitely the most important and highest priority enhancement. We need to allow that at least birth and/or death rates vary according to a Wienner process (random walk) with a parameter for the variance (environmental variability). I see three issues here:
1 - How to calculate the time ellapsed between steps when rates vary. The exponential distribution does not apply in this case. Code in Time-Varying_rates.R
is a first step to implement Gibson & Bruck's recipe to the pdf of the time distribution, given a reaization of the random walk of rates along time. To be checked, see comments there.
2 - From this pdf we could then draw a time for the event for each individual and then proceed with the Gillespie algorithm.
A sidenote: both birth and death rates are bounded to zero. Probably a random walk of the log of the rates is simple solution (this means a log-normal random-walk).
A zero-order interpolation of the environmental variability would be to use the bdm() process considering fixed birth and death rates, and after the time step is completed, the birth and death rates can be updated to reflect the environmental variability.
We can use a zero-order interpolation of the environmental variability and use Marion & Renshaw's equation 2.5 in order to validate the interpolation error. We can also use the expected derivative of the environmental factor to compare with the inverse of the sum of the Gillespie weights. If the ISGW is one or two orders of magnitude larger than the derivative of the environmental fluctuations, the zero-order interpolation errors can be regarded as negligible.
The environmental stochasticity should be given as a data.frame with time x birth/death rates which will then be interpolated to give instantaneous birth/death rates.
[EDIT: reimplemented this as K fluctuations, which affect the death rate. The "first" and "third" issue below have been solved]
Implemented a first version of the zero-order interpolation for the birth rates. The stochasticity is clearly affecting the simulation results as can be seen in the plots below (The code for this examples can be found on the vignette, just run some times the first chunck until you find a nice stochastic variation). This is the birth rate multiplier:
And this spike in birth rate and subsequent decline can be seen here:
There are three main concerns for me now. First, I'm currently concerned that this structure might be too hamstrung if we decide to incorporate other kinds of environmental stochasticity. So I'd like to take some time to delineate which kinds of stochasticity can be used in the foreseeable future of this package, to try to get them all in the same framework.
Second, there's the question of assessing the cost of this approximation. Are the errors (in respect to Marion's eq 2.5, for instance) sufficiently small?
Third, there are some performance problems in the current code. But I think it's really best to focus on the first issue before attempting any kind of optimization here.
A note on the theory behind this issue: Marion's 2000 paper makes a distinction between i) Constant Environment; ii) Deterministic Variation and iii) Continuous Stochastic Variation. Equation 2.5, mentioned above, is derived under the Deterministic Variation section. Is it still valid for CSV?
[EDIT: implemented]
Different species may be affected by the stochastic component in different ways. Maybe there could be a vector that associates how much the environmental stochasticity affects each specific species - so a[25] = 1
means that species 25 is 100% positevely affected by the stochastic factor. a[35] = 0
means that species 35 is totally unaffected, and a[99] = -0.1
means that species 99 suffers a small detrimental effect.
Thinking about the relation between the characteristic time for the environmental stochasticity and the mean time for a simulated event to occur, our worst case scenario is when the community is empty, in which case the simulated event is expected to occur in 1 / sum (migration rates)
; maybe we should warn the users that the stochastic functions are better approximation for climax communities only.
Can we write a diagnostic function for the expected error of the zero-order interpolation of the stochastic proccess?
Different sensitivities to environmental stochasticity implemented. The effect of the environmental sensitivity can be seen on the following plot, which traces the population evolution over time for four species, ranging from fully affected by the environmental fluctuations to one that is reversely affected (environmental strengths = c(1, 0.5, 0, -1)
)
According to Chesson's theory of coexistence only non-linear relationships of species competitive ability to environent can contribute to coexistence. So, can we allow the effect of environmental fluctuations on k be nonlinear? The simplest way that occurred to me is to allow a linear relationship at log-log scale, that implies in a power-law at arithmetic sclae. Linear relationhip will be then a special case.
Can we incorporate environmental stochasticity in our models? How is the Gillespie processed changed? See Gibson & Bruck, 2000.