reconhub / earlyR

Estimation of infectiousness in the early stage of an outbreak, with short-term predictions.
http://www.repidemicsconsortium.org/earlyR
Other
9 stars 3 forks source link

Potential issue in log likelihood calculation #7

Open Gulfa opened 5 years ago

Gulfa commented 5 years ago

I was trying to replicate the results of the ML estimate for R in this package without being successful. I have narrowed down the difference between my implementation and the one in this package to line 175 in get_R.R.

sum(stats::dpois(x[-1], lambda = R * x_lambdas, log = TRUE))

I might be misunderstand something, but I think that the length of x[-1] and x_lambdas are different in this line of code and the dpois function then recycles values from x[-1] which is shorter. From my understanding x_lambdas has one entry too many and the final element should be removed. I tried this in a fork and then the estimated value for R matched my python implementation. In the example in https://www.reconlearn.org/post/practical-ebola-reconstruction.html, this changes the estimated value of R from 1.286 to 1.336.

Happy to submit a PR from my fork if that would be helpful and I am not completely mistaken about this.

annecori commented 5 years ago

@thibautjombart I agree with the above issue spotted by @Gulfa - maybe worth having a look? One way of solving this I think is to not do the subsetting [-1] on line 170 but instead do it inside line 176 with sum(stats::dpois(x[-1], lambda = R * x_lambdas[-1], log = TRUE))

thibautjombart commented 5 years ago

I think this is now fixed by the recent PR: https://github.com/reconhub/earlyR/blob/master/R/get_R.R#L177-L206

zkamvar commented 5 years ago

I don't think this is exactly fixed. The first lambda is now NA

thibautjombart commented 5 years ago

@annecori will confirm, but I think the first lambda is NA: the force of infection prior to the first case is unknown, and the first non-negative incidence actually does not enter the likelihood calculation.

annecori commented 5 years ago

yes that's correct

zkamvar commented 5 years ago

@Gulfa, can you check the most recent github version of earlyR to check if it's correct?