Feature/core eco basic ricker

[clairecwinfrey]

Creates a new script, coreEcoBasicRicker.R, for use in the population ecology unit for Core Ecology to accompany discussion of Melbourne & Hastings, 2008
Script overview:

1. deterministic Ricker models from original 12_basic_ricker_fit (with extra commentary and questions)
2. Explanation and hands-on coding with Poisson and Bernoulli/binomial distributions
3. Plot and explore Suppl. Fig. 1, specifically Poisson, Poisson-binomial, and negative binomial-binomial gamma models. Code from 05_plotprodfunc.R with needed functions from source/ explicitly added in. Questions to walk students through models and how they reflect biology

[brettmelbourne]

Add to gitignore:

*.rproj
.gitignore

and delete gitignore from this branch

[brettmelbourne]

Delete this file from commit (suggested adding to gitignore)

[brettmelbourne]

It's worth showing how transforming the model allows this linear fit. Also, it's a rough fit because of not handling the various data issues (-Inf etc)
N_{t+1} = R N_t e^(-alpha N_t)
N_{t+1} / N_t = R e^(-alpha N_t)
log(N_{t+1} / N_t) = log(R e^(-alpha N_t))
r = log(R) - alpha N_t
r = r_0 - alpha N_t

[brettmelbourne]

Try also increasing R to see stable limit cycles and chaos. The code in bifurcation.R will make some plots that show where the bifurcations occur. This is one of the interesting dynamical properties that emerge from these discrete generation models.

[brettmelbourne]

By setting the seed the same, this 100 will be same as the first 100 in previous 1000. Perhaps that's what you wanted to show? Otherwise, suggest not resetting the seed.

[brettmelbourne]

I think here a good thing to do might be to walk through the lifecycle of a single individual.

R = B(1-m), where B is birth rate and m is mortality probability

1. It gives birth

births <- rpois(n=1, lambda=B)

2. Now each of these individuals survives or dies. For a single individual, it has to survive both cannibalism (a density-dependent process) and density-independent causes of death. Let's look at one individual:

First, does it survive density-dependent cannibalism?

p_sdd <- exp(-1 * alpha * Nt)
survive <- rbinom(n=1, size=1, prob = p_sdd  )

Now, if it does survive cannibalism, does it also survive density-independent causes of death?

p_sdi <- 1 - m
if ( survive )  {
    survive <- rbinom(n=1, size=1, prob = p_sdi  ) 
}

A simple emergent stochastic process arises from these two different binomial stochastic processes. The emergent process is also binomial. Combining them together gives the following stochastic process for the survival of a single individual:

p_sdd <- exp(-1 * alpha * Nt)
p_sdi <- 1 - m
survive <- rbinom(n=1, size=1, prob = p_sdd * p_sdi  ) 

Now summing up the survival of all the individuals that were born:

survivors <- 0
if ( births > 0 )  {
  for ( i in 1:births)  {
    survive <- rbinom(n=1, size=1, prob = p_sdd * p_sdi  ) 
    if ( survive )  {
      survivors <- survivors + 1
    }
  }
}

So, how many made it to become the next generation , N_{t+1}:

print(survivors)

For the whole population of N_t, we add up the contribution from all the individuals. The function following collapses these steps into two emergent overall stochastic distributions and is a vectorized version of this, and is further vectorized to accept multiple populations at once (i.e. Nt can be a vector).

[brettmelbourne]

One neat thing about the above, is that the emergent stochastic process for that original individual, even though a sequence of different stochastic processes (one Poisson and two binomial events), is simply Poisson:
Ntp1_i <- rpois(n=1, lambda = R * exp(-1 * alpha * Nt))

[brettmelbourne]

It could be worth replacing length(Nt) with n=1, as that's all we're doing here. The length(Nt) allows multiple populations to be done at once but it might be confusing at first.

[brettmelbourne]

What might be cool here instead is to add the steps back in (since the rpois is combining the two steps from the previous function, so it's less obvious what is happening biologically):

females <- rbinom(1, size=Nt, prob=p)
births <- rpois(1, lambda=females * (1/p) * R))  #why do we multiply R * 1/p?
survivors <- rbinom(1, size=births, prob= exp(-1 * alpha * Nt)

[brettmelbourne]

Then modify the function to add back in the length(Nt) argument, which allows multiple populations to be simulated at once in a vectorized way. Then simulate a bunch of populations with the same Nt

Ntp1 <- RickerStBS(Nt=rep(10, 10000),  R,  alpha)
hist(Ntp1)
var(Ntp1)

Make a histogram of N_{t+1}. This is the distribution of the stochastic process. What is the variance in N_{t+1}?

[brettmelbourne]

So this is the math solution to the variance. It's kinda neat ('cause math ;) but perhaps not obvious directly. We could add in the simulation of the variance suggested above and show that it matches the math.

[brettmelbourne]

vars <- rep(NA,101)
i = 1
for ( Nt in seq(0, 1000, 10) ) {
  Ntp1 <- RickerStBS(rep(Nt, 10000),  R,  alpha)
  vars[i] <- var(Ntp1)
  i <- i + 1
}
plot(seq(0, 1000, 10) , sqrt(vars), xlab="Nt", ylab="Standard deviation in N_{t+1}")

[brettmelbourne]

It would be good to simulate the dynamics over time now. Code is simple and is in plottimeseries.R. Do a few runs to see the variability. Compare to the deterministic model. Compare Poisson to say a model with environmental variation.

[brettmelbourne]

Some people may be familiar with the negative binomial. This is the same K. This is where a negative binomial comes from: it's a compound Poisson-gamma distribution. K can be thought of as the dispersion parameter, much like the variance of the Normal is the dispersion parameter of the Normal. The only twist is that magnitude is reversed: small K is large dispersion, infinite K is zero dispersion and returns the Poisson model. Biologically, the way to think about this is: small K = large variation in birth rates, large K = small variation in birth rates, infinite K = no variation in birth rates among individuals or environments (reduces to Poisson variation, i.e. pure demographic stochasticity).

[brettmelbourne]

For me, the neat discovery was that these variances simply add up (and had simple equations for the components). The reason the Ricker function is there is because it gives the mean N_{t+1} of the stochastic process. In the Poisson, the variance is equal to the mean of the stochastic process. I forget the derivation of the others (I remember discovering it and being surprised, actually after the paper was published, so the results in the paper were actually by simulation) but it's the expected variance for that component.

[brettmelbourne]

This is great! I've made a few suggestions.

[brettmelbourne]

Thanks for this, it looks wonderful! I made a few suggestions.