Comparison of EpiEstim::wallinga_teunis and R0::est.R0.TD

[tamas-ferenci]

The package R0 has a function called est.R0.TD which - I believe - does the same as EpiEstim's wallinga_teunis.

Indeed, results are almost identical if you consider the following example:

data("Germany.1918", package = "R0")
res1 <- R0::est.R0.TD(Germany.1918,generation.time("gamma",c(3, 1.5)),begin=1,end=126)
res2 <- EpiEstim::wallinga_teunis(Germany.1918, method="parametric_si",
                       config = list(t_start = 1:126, t_end = 1:126,
                                     mean_si = 3, std_si = 1.5,
                                     n_sim = 100))
ggplot(res2$R[1:124,],aes(x=t_start,y=`Mean(R)`)) + geom_line() + 
geom_line(data=data.frame(x=2:125,y=res1$R[2:125]),aes(x=x,y=y),color="red")

However, the story is different is you use only the initial part of the epidemic:

res1 <- R0::est.R0.TD(Germany.1918[1:30],generation.time("gamma",c(3, 1.5)),
begin=1,end=30)
res2 <- EpiEstim::wallinga_teunis(Germany.1918[1:30], method="parametric_si",
                                  config = list(t_start = 1:30, t_end = 1:30,
                                                mean_si = 3, std_si = 1.5,
                                                n_sim = 100))
ggplot(res2$R[1:28,],aes(x=t_start,y=`Mean(R)`)) + geom_line() + 
geom_line(data=data.frame(x=2:29,y=res1$R[2:29]),aes(x=x,y=y),color="red")

The result from EpiEstim starts to go to zero as we are approaching the end of the data, while R0 keeps a realistic trajectory.

I tried both functions on many epidemic curves, but the overall impression was the same: if we are in the middle of an epidemic, EpiEstim's results will go to zero at the end of the observation period, while R0 provides feasible results even at the end of the data.

I don't know the mechanics of both functions in detail, but I was wondering if this is an error in EpiEstim or there is some other explanation for this phenomenon.

[jstockwin]

Hi @tamas-ferenci, thanks for submitting the issue. We've also recently seen problems with this function in #92, which I naively try to fix in #93 but unfortunately I am not really familiar with the methods involved here. You could try my fix, but I think I'd be surprised if this changes things for you.

@zkamvar or @robin-thompson are either of you familiar with this method and would be able to review the code? Otherwise, do you know somebody who could?

[tamas-ferenci]

Thanks @jstockwin , I tried it, but indeed to no avail. As otherwise it is producing the same results as R0, and the difference only appears when we are in the middle of an epidemic, and even it that case only at the end, I'm pretty sure it is not a bug in the strict sense of the word, rather some conceptual difference (perhaps R0 has a method to estimate realistic values even when we are running out of future observations...?).

[annecori]

Thanks for reporting this and sorry for the slow reply. I think what happens is that EpiEstim ignores the incidence after day 30 while R0 does not. I will try to fix this asap.

[annecori]

@tamas-ferenci, so actually I don't think what we are doing is wrong, as when using only the truncated incidence as input as you do in. your example, the estimated effective reproduction number should indeed go to zero toward the present (this is a known issue of right censoring in the WT method). There have been some methods developed to fix this right censoring issue but we implemented the most simple version of WT here. If you use the whole curve you get a result much more similar to the R0 package, see my example below. Therefore I think the R0 package must somehow make assumptions about future incidence to fix the right censoring issue, which we don't do. Please let me know if you are happy with my explanation and I will close this issue as I actually don't think this is a bug in EpiEstim. Thanks!

Please note atm you need to run the below code from branch wallinga-fix.

data("Germany.1918", package = "R0")
res1 <- R0::est.R0.TD(Germany.1918[1:30],generation.time("gamma",c(3, 1.5)),
                      begin=1,end=30)
res2 <- EpiEstim::wallinga_teunis(Germany.1918[1:30], method="parametric_si",
                                  config = list(t_start = 1:30, t_end = 1:30,
                                                mean_si = 3, std_si = 1.5,
                                                n_sim = 100))
res3 <- EpiEstim::wallinga_teunis(Germany.1918, method="parametric_si",
                                  config = list(t_start = 1:30, t_end = 1:30,
                                                mean_si = 3, std_si = 1.5,
                                                n_sim = 100))
plot(res2$R$t_start, res2$R$`Mean(R)`, type = "l", xlab = "Time", ylab = "Mean R")
lines(2:29, res1$R[2:29], col="red")
lines(2:29, res3$R$`Mean(R)`[1:28], col="blue")
legend("topright", c("R0", "EpiEstim on early curve", "EpiEstim on whole curve"),
       lty = 1,
       col = c("red", "black", "blue"))

[tamas-ferenci]

Thanks @annecori for the answer, and sorry for the late reply! My best guess is that the solution is the following sentence from the R0's manual: "Correction for estimation in real time is implemented as in Cauchemez et al, AJE (2006)." Unfortunately they omitted the bibliographic details of this reference, so I can only speculate that this is the article they referred to. (The AJE paper is likely this, but methodologically it simply refers back to this paper in EID.) Interestingly, it involves a complicated process (at least at first glance), while the code of est.R0.TD only seems to contain minor differences, but this is perhaps a question rather for the authors of R0.

I agree that this is definitely not a bug, so at this point I think it's really up to you to either close this issue, or change it from "bug" to "enhancement", if you may decide at some point in the future to implement this extension of WT by Cauchemez to account for right censoring (which would be a nice addition in my opinion!).

Thanks again for your feedback!