exponential growth rate calculation

[davidearn]

Irena and were just talking about how to estimate r for the vaxified model. The existing rExp function fails if compartment sizes are zero (e.g., some compartments in the vaxprotect layer). To fix that, we added a check so that if r_last == r_nextlast then we return 1 for the ratio rather than 0. Then the log returns 1, which is what we want (not NaN). My question now is why is the estimated r the mean of the growth rates of all the growth rates for each epidemiological class? The fastest growth will dominate this mean if it is much larger than the rest, but what is the logic for using the mean here? Irena also wonders if it would be better to condense over vaccination layers before estimating the growth rate (this would also almost certainly avoid any empty compartments).

[dushoff]

Presumably the logic for using the mean is that this is meant to be used on a linearized model where the differences between growth rates are meant to be numerical errors. I would nonetheless argue that it's more robust to take the natural weighted average -- which is exactly the same as what you would get by adding first and dividing once, afaics

On Wed, Jun 2, 2021 at 12:25 PM David Earn @.***> wrote:

Irena and were just talking about how to estimate r for the vaxified model. The existing rExp function fails if compartment sizes are zero (e.g., some compartments in the vaxprotect layer). To fix that, we added a check so that if r_last == r_nextlast then we return 1 for the ratio rather than 0. Then the log returns 1, which is what we want (not NaN). My question now is why is the estimated r the mean of the growth rates of all the growth rates for each epidemiological class? The fastest growth will dominate this mean if it is much larger than the rest, but what is the logic for using the mean here? Irena also wonders if it would be better to condense over vaccination layers before estimating the growth rate (this would also almost certainly avoid any empty compartments).

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/bbolker/McMasterPandemic/issues/51, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAW5E6IW7SYLKKXQXUQVCB3TQZLO3ANCNFSM457CNNCQ .

-- Jonathan Dushoff (https://tinyurl.com/jd-pronouns) McMaster University Department of Biology https://mac-theobio.github.io/ https://twitter.com/jd_mathbio http://jd-mathbio.blogspot.com/

[bbolker]

Exactly.

On 6/2/21 1:14 PM, Jonathan Dushoff wrote:

Presumably the logic for using the mean is that this is meant to be used on a linearized model where the differences between growth rates are meant to be numerical errors. I would nonetheless argue that it's more robust to take the natural weighted average -- which is exactly the same as what you would get by adding first and dividing once, afaics

On Wed, Jun 2, 2021 at 12:25 PM David Earn @.***> wrote:

Irena and were just talking about how to estimate r for the vaxified model. The existing rExp function fails if compartment sizes are zero (e.g., some compartments in the vaxprotect layer). To fix that, we added a check so that if r_last == r_nextlast then we return 1 for the ratio rather than 0. Then the log returns 1, which is what we want (not NaN). My question now is why is the estimated r the mean of the growth rates of all the growth rates for each epidemiological class? The fastest growth will dominate this mean if it is much larger than the rest, but what is the logic for using the mean here? Irena also wonders if it would be better to condense over vaccination layers before estimating the growth rate (this would also almost certainly avoid any empty compartments).

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/bbolker/McMasterPandemic/issues/51, or unsubscribe

https://github.com/notifications/unsubscribe-auth/AAW5E6IW7SYLKKXQXUQVCB3TQZLO3ANCNFSM457CNNCQ .

-- Jonathan Dushoff (https://tinyurl.com/jd-pronouns) McMaster University Department of Biology https://mac-theobio.github.io/ https://twitter.com/jd_mathbio http://jd-mathbio.blogspot.com/

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/bbolker/McMasterPandemic/issues/51#issuecomment-853231998, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAATIRWUVCT6OFE3J3TPKYLTQZRHTANCNFSM457CNNCQ.

[papsti]

OK, so this makes me think that computing the growth rate as the mean over the estimated rates for the vaxified compartments is not correct, because some compartments have growth rate zero while others don't (and the difference is larger than just numerical error). Should I condense the results over vax layers first, before taking the ratio?

[davidearn]

While chatting with Ben I think I got my head around why this is OK. My current thinking is:

• by discarding susceptible outflow terms in the ODEs, we turn the system into a linear system.

• providing the linear system is non-degenerate (irreducible Jacobian matrix), the linear system converges as t-->infty to exponential growth along the dominant eigenvector (which exists because of nondegeneracy).

• provided we integrate long enough (e.g, "many" GIs), we will be close to the asymptotic behaviour in which each component of the state vector grows with exactly the rate we are trying to estimate, hence we want to average out any differences between components, and if we integrate for sufficiently long then it should make any difference how we "average out"

• all of this works iff the initial state vector does not contain any zero components

It seems to me that the solution is not the hack Irena and I implemented earlier today, but to replace zero components in the initial state with something small but positive. Then the question becomes how small should this be in order to reach the asymptotic growth behaviour by the end of the integration?

Alternatively, condensing should be fine provided we end up with a strictly positive state vector.

David

On Wed, 2 Jun 2021, Irena Papst wrote:

OK, so this makes me think that computing the growth rate as the mean over the estimated rates for the vaxified compartments is not correct, because some compartments have growth rate zero while others don't (and the difference is larger than just numerical error). Should I condense the results over vax layers first, before taking the ratio?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, orunsubscribe.[ADMXZJ37OHPDA4LUM2SAQ3DTQ2Q3ZA5CNFSM457CNNC2YY3PNVWWK3TUL52HS4 DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOGLO6RQY.gif]

[papsti]

Ah, OK. The default for vaxifying a state is that everyone starts in the unvaccinated layer, so just distributing people evenly among the vaccination strata in rExp's initial state should do the trick. I could also just integrate for longer than 100 steps, but that's a less efficient solution.

I'll implement the first fix I suggested and revert the hack we implemented this morning.

[dushoff]

It's still probably better to condense. Maybe also worth thinking about some sort of check for whether you've integrated long enough (small difference between r_halfway and r_final, for example)?

On Wed, Jun 2, 2021 at 6:48 PM Irena Papst @.***> wrote:

Ah, OK. The default for vaxifying a state is that everyone starts in the unvaccinated layer, so just distributing people evenly among the vaccination strata in rExp's initial state should do the trick. I could also just integrate for longer than 100 steps, but that's a less efficient solution.

I'll implement the first fix I suggested and revert the hack we implemented this morning.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/bbolker/McMasterPandemic/issues/51#issuecomment-853431944, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAW5E6PBA4QOTBNMEDAVGODTQ2YLLANCNFSM457CNNCQ .

-- Jonathan Dushoff (https://tinyurl.com/jd-pronouns) McMaster University Department of Biology https://mac-theobio.github.io/ https://twitter.com/jd_mathbio http://jd-mathbio.blogspot.com/

[papsti]

Yeah, actually I realized that after trying the uniform distribution fix and realizing that nothing ever flows into the hospital compartments of the vaxprotect layer when we turn off severe infections for vaccinated individuals (the default). Even if we allow a small number of severe infections for vaccinated individuals, this flow will be slow and I'm not sure 100 steps will be enough to integrate.

I'm adding the condense step now. I'll add a warning if the variance of the log(r_last/r_nextlast) is above some threshold.

[bbolker]

I've seen some commits on your branch about this; based on our conversation yesterday (4 June) and your subsequent steps, don't forget to comment here/explain the solutions you ended up with/close the issue.