estimate growth rate r using Poisson glm instead of log-linear model
Opened this issue · 4 comments
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")
}
similar issue, different package:
epiforecasts/EpiNow#74
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.
yes, if there is overdispersion quasipoisson or negative binomial regression would be best.
Related issue I just posted: epiforecasts/EpiNow#75