estimate r from initial R

[sbfnk]

Description

This PR closes #920. For now it's using the existing estimate for initial number of infections. As an initial step I'm interested in how doing this affects performance and/or recovery of estimates.

Initial submission checklist

• My PR is based on a package issue and I have explicitly linked it.
• I have tested my changes locally (using devtools::test() and devtools::check()).
• I have added or updated unit tests where necessary.
• I have updated the documentation if required and rebuilt docs if yes (using devtools::document()).
• I have followed the established coding standards (and checked using lintr::lint_package()).
• I have added a news item linked to this PR.

After the initial Pull Request

• I have reviewed Checks for this PR and addressed any issues as far as I am able.

[seabbs]

it seems like there are some issues when accumulating to weekly

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if eada335 is merged into main:

• ❗🐌default: 22.4s -> 31.8s [+27.84%, +56.7%]
• ❗🐌no_delays: 24.4s -> 35s [+28.06%, +58.29%]
• ❗🐌random_walk: 9.08s -> 12.7s [+34.19%, +45.36%]
• ❗🐌stationary: 12.9s -> 19.2s [+27.9%, +70.83%]
• ❗🐌uncertain: 36.7s -> 47.2s [+16.65%, +40.78%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[sbfnk]

This is now ready for initial review but it might require a bit more testing (and doc updates / news item etc) before releasing it in the wild. A few observations so far

• I couldn't get the package to compile when using the algebraic newton solver either locally or on GHA. I could get it to compile with the Powell solver but this led to occasional failures. It all seems to work well (and possibly faster) hand-rolled Newton solver. For now these two methods are both available as stan functions.
• Initialisation of the number of infections remains a critical issue. If doing this independently of estimated growth we get a strong correlation with the initial R as predicted by @SamuelBrand1 - for now what's estimated is roughly the initial number of infections related to the first week of data, which is then scaled based on the growth rate and length of the seeding period; doing it this way reduces the correlation and improves fits, but means we can no longer interpret the corresponding parameter as being the initial number of infections. If going ahead with this we'll have to adjust the simulation functions.
• touchstone is failing now for unknown reasons, which is annoying because it would be good to measure the performance impact of the latest version.

[seabbs]

touchstone is failing now for unknown reasons, which is annoying because it would be good to measure the performance impact of the latest version.

I see this in epinowcast as well. Precisely what is happening is not clear to me.

[seabbs]

t a high level I think this looks good and it seems reasonable to rescale initial infections I think. A few comments and questions about output changes. I have just site read so far so need to have a play with the models.

[sbfnk]

touchstone is failing now for unknown reasons, which is annoying because it would be good to measure the performance impact of the latest version.

I see this in epinowcast as well. Precisely what is happening is not clear to me.

I now see the same error locally on a linux box. Pretty sure it's unrelated to touchstone or package code but why exactly it happens is a mystery to me, too (as the action goes from success to failure with nearly identical logs, i.e. same package versions, same runner version etc.).

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 75628ff is merged into main:

• ❗🐌default: 20.7s -> 23.7s [+5.37%, +22.74%]
• ✔️no_delays: 24.1s -> 24.6s [-6.11%, +9.54%]
• ❗🐌random_walk: 9.03s -> 9.8s [+0.92%, +16.17%]
• ✔️stationary: 12.5s -> 11.7s [-21.33%, +8.78%]
• ✔️uncertain: 34.1s -> 36.2s [-5.76%, +18.11%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 0449171 is merged into main:

• ✔️default: 22.2s -> 22.6s [-6.78%, +10.93%]
• ✔️no_delays: 23.5s -> 24.2s [-8.72%, +14.6%]
• ✔️random_walk: 9.34s -> 14.7s [-65.73%, +180.3%]
• ✔️stationary: 12.7s -> 13.6s [-20.85%, +35.9%]
• ✔️uncertain: 34.3s -> 37.2s [-4.16%, +21.23%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[sbfnk]

I think this is good to go - potentially of interest to @SamuelBrand1 @seabbs I had to introduce a truncation step to avoid numeric overflows when the first value of R came very close to zero.

EpiNow2/inst/stan/functions/rt.stan

Line 103 in dba0c2e

real r = fmax((R - 1) / (R * mean_gt), -1);

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if dba0c2e is merged into main:

• ❗🐌default: 19.5s -> 24s [+14.6%, +31.83%]
• ❗🐌no_delays: 23.9s -> 27.3s [+0.54%, +27.41%]
• ❗🐌random_walk: 9.1s -> 9.45s [+0.03%, +7.76%]
• ✔️stationary: 12.7s -> 11.8s [-21.57%, +8.88%]
• ✔️uncertain: 33.6s -> 36.9s [-25.35%, +44.61%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if dba0c2e is merged into main:

• ✔️default: 23.4s -> 22.5s [-9.79%, +2.42%]
• ✔️no_delays: 32.3s -> 26.4s [-76.64%, +40.15%]
• ✔️random_walk: 9.75s -> 9.6s [-7.67%, +4.61%]
• 🚀stationary: 13.9s -> 12.9s [-11.27%, -3.26%]
• 🚀uncertain: 36.3s -> 32.8s [-16.55%, -2.49%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[seabbs]

I think this is good to go - potentially of interest to @SamuelBrand1 @seabbs I had to introduce a truncation step to avoid numeric overflows when the first value of R came very close to zero.

Just looking now. Did you compare this approach to padding the initial guess as in theory fmax should be a bit of a no no?

[sbfnk]

I think this is good to go - potentially of interest to @SamuelBrand1 @seabbs I had to introduce a truncation step to avoid numeric overflows when the first value of R came very close to zero.

Just looking now. Did you compare this approach to padding the initial guess as in theory fmax should be a bit of a no no?

Padding how? You mean add a small number? I looked at this one but it would be expected to be slower because you now always start in the "wrong" place -- what would you expect the issue to be with fmax? It's not applied to anything being sampled.

real r = (R - 1) / (R * mean_gt + 1);

[seabbs]

Yes that definitely is the tradeoff so I am not clear it is worth it. It is being applied to something being sampled though, the first R index.

[sbfnk]

The last two comments tell you the speed difference between the two versions.

I also remember seeing more fitting failures (initial conditions leading to NaN) with that version but could check again. Do you have any other suggestion?

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 6a36058 is merged into main:

• ✔️default: 21.9s -> 25.2s [-8.24%, +37.99%]
• ✔️no_delays: 24.7s -> 24.7s [-11.71%, +11.55%]
• ✔️random_walk: 9.62s -> 9.46s [-8.06%, +4.69%]
• ✔️stationary: 12.9s -> 12.1s [-20.13%, +7.61%]
• ✔️uncertain: 36.4s -> 37.3s [-8.88%, +13.51%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[sbfnk]

The last two comments tell you the speed difference between the two versions.

Overall seems indistinguishable from stochasticity at the used sample size (n = 5)

I also remember seeing more fitting failures (initial conditions leading to NaN) with that version but could check again.

Don't see these upon checking again - perhaps got muddled up with other rstan issues.

I'm just updating the snapshots to see if checks pass then, in which case I'd be happy to go with this version.

[seabbs]

Do you have any other suggestion?

Aside from the padding you are trying no not really.

I'm just updating the snapshots to see if checks pass then, in which case I'd be happy to go with this version.

I really don't feel strongly especially if there is no real difference.

[seabbs]

I really like this. I've checked back against the Julia code and agree it aligns. Running locally all seems basically the same as the version on main which is good.

In terms of the padding vs max I don't feel very strongly and agree there is a bit of a tradeoff

The only thing that is missing which could come in a follow up PR is unit tests of the R to r functions in a similar vein to the julia code and the original issue with the test approach

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 35740e5 is merged into main:

• ✔️default: 21.1s -> 22.8s [-5.97%, +21.92%]
• ✔️no_delays: 23.8s -> 26.4s [-6.86%, +29.12%]
• ✔️random_walk: 9.81s -> 26.8s [-221.48%, +566.99%]
• ✔️stationary: 13s -> 19.9s [-15.9%, +121.69%]
• ✔️uncertain: 33.3s -> 34.1s [-8.37%, +12.96%]
  Further explanation regarding interpretation and methodology can be found in the documentation.