estimate growth rate r using Poisson glm instead of log-linear model

[ffinger]

The way the exponential growth rate is estimated here currently (in fit()) is to use a log-linear model (linear model similar to lm(log(incidence) ~ date). For the log to work you need to remove the 0s or replace them with a small positive value.

A better way to estimate the growth rate would be to use a Poisson glm of the incident cases (glm(incidence ~ date, family = poisson) or similar). In addition to the more appropriate error model for count data this can handle 0s in the data natively.

Here are some references:
https://besjournals.onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2010.00021.x
https://bmcmedinformdecismak.biomedcentral.com/articles/10.1186/1472-6947-12-147

And here is how this is done in the R0 package, where you can chose between poisson and log-linear models, with a default to poisson (https://github.com/cran/R0/blob/master/R/est.R0.EG.R):

  # Method 2 == Poisson regression
  else if (reg.met == "poisson") {
    #tmp <- glm(incid ~ t.glm, family=poisson(), data=epid)
    tmp <- glm(incid ~ t, family=poisson(), data=epid)
    Rsquared = (tmp$null.deviance-tmp$deviance)/(tmp$null.deviance)
    r <- coefficients(tmp)[2]
    confint = confint(tmp)[2,]
    pred= predict(tmp,type="response")
  }

[ffinger]

similar issue, different package:
epiforecasts/EpiNow#74

related:
#93
#28
#94

[thibautjombart]

Yes, definitely agreed! This is what I am currently using - quasi-Poisson actually. I was not aware of EpiNow .. there is some duplication there, and I am not sure how to handle this going forward. incidence should be about data handling and plots really, and reforms we discussed before were aiming at moving fitting to a separate package.

[ffinger]

yes, if there is overdispersion quasipoisson or negative binomial regression would be best.

[thibautjombart]

Related issue I just posted: epiforecasts/EpiNow#75