1
$\begingroup$

I'm playing with the JHU nCOV data and looking to calculate the doubling rate in my region (Western Australia) - I can get it down to an integer value via a kind of brute force (halve the current case value, use excel maxifs to look at the most recent date where the case value is =< take the difference)

Is there a better way?

$\endgroup$

2 Answers 2

1
$\begingroup$

In the early periods of an infection, we can model the number of people sick as

$$ y = \beta_0 \exp(\beta_1 t) $$

We are looking for a time, $\Delta t$, so that

$$ 2 y = \beta_0 \exp(\beta_1 t + \beta_1\Delta t)$$

This means that

$$ 2 = \exp(\beta_1\Delta t)$$

or

$$ \dfrac{\log(2)}{\beta_1} = \Delta t$$

So, the time required for the exponential to double is $\log(2) / \beta_1$. We need an estimate of $\beta_1$, $\hat{\beta}_1$. If you have infection data, the easiest way to do this is to do a linear regression on the log scale where $\hat{\beta}_1$ will be the coefficient of time. Generally speaking, this is not integer valued. My own estimates from the United States are approximately 2.6 days (though we shouldn't take this estimate very seriously). Since this isn't integer valued, you have to round up to 3 days. Why round up? After 2 days, the number of infections hasn't quite doubled. After 3 days, the number of infections is a little more than double. So, in order for the infection to double, we have to wait $3 = \lceil 2.6\rceil$ days.

$\endgroup$
8
  • $\begingroup$ Why do any rounding at all? $\endgroup$
    – whuber
    Mar 31, 2020 at 15:14
  • $\begingroup$ @whuber No need to, but reports of infection are not in continuous time. Data are probably updated daily, so although 2.6 is a fine estimate of the doubling time, every.3 days we will see the counts double in the data. $\endgroup$ Mar 31, 2020 at 15:37
  • 1
    $\begingroup$ The difference between 2.6 and 3 is large enough to make one cautious about rounding. For instance, after a month of 30 days the prediction in one case is an increase by a factor of $2^{30/2.6} = 3000$ and in the other case only a factor of $2^{30/3}=1000.$ $\endgroup$
    – whuber
    Mar 31, 2020 at 16:07
  • $\begingroup$ In the short term, I see no problem rounding the doubling time up to the nearest integer, especially since data is reported daily and not continuously. Extrapolating 30 days would make any exponential model fit to the data look bad. There are no uncertainty estimates here either, which to me is a bigger concern than me rounding my estimate. $\endgroup$ Mar 31, 2020 at 16:25
  • 2
    $\begingroup$ It is usually (far) better not to round overly much, especially when the value you report is likely to be used in further calculations or analyses that might be sensitive to the rounding error, as in my previous example. $\endgroup$
    – whuber
    Mar 31, 2020 at 17:57
0
$\begingroup$

It might be interesting to express the doubling rate as a function of time. It is definitely not a constant.

For this purpose, I plotted the curve on a logarithmic scale. The slope indicates how fast the number doubles.

Below I use the inverse of the slope/rate to determine the number of days that are necessary to double the number of cases

$$\tau_{\text{doubling time}} \approx \frac{\Delta t}{ log_2(N(t+\Delta t))-log_2(N(t))}$$

Where $\tau_{\text{doubling time}}$ is the doubling time-scale, $N(t)$ is the number of cases at time $t$, $\Delta t$ is a small-time difference between which we compare the change.

Instead of using a difference with $\Delta t$ we could use a derivative. The $\Delta t$ is not very precise and does not exactly express the instantaneous growth rate, and instead, it is the growth over the last $\Delta t$ period of time. But to use differentiation you would need to fit a curve to the data, and if this curve is not a good model (ie. biased) then this will introduce error in a different way.

plot

# John Hopkins
dat <- read.csv("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv")

# extract country data
y <- as.numeric(dat[dat$Province.State == "Western Australia", -c(1:4)])
d <- 3
tr <- max(which(y<1)):(length(y)-d)  ##trimming data to only the numbers >1

## plot of cases on logarithmic scale
plot(-1,-1, ylim = c(1,512), xlim = c(0,40), 
     ylab = "cases", xlab = "date", log = "y",xaxt = "n", yaxt = "n", 
     main = "plot of cases on logarithmic scale")

for (i in 0:10) {
  lines(c(-10,100),c(1,1)*2^i, col = "gray", lty=2)
}
axis(2,2^c(0:10), las = 2)
axis(1,at = c(1)+c(0:5)*7, labels = c("29 feb","7 Mar", "14 Mar", "21 mar", "28 Mar", "4 Apr") )

points(y[tr+d], pch = 21, col = 1,  bg =  floor(log(y[tr+d],2)))

# plot the number of days neccesary to increase the logarithm
plot(-1,-1, ylim = c(1,20), xlim = c(0,40), 
     ylab = "cases", xlab = "date", log = "y", yaxt = "n", xaxt = "n",
     main= "number of days neccesary to increase the logarithm \n inverse slope 3/log(cases[t]) - log(cases[t-3])")
axis(2,2^c(0:10), las = 2)
axis(1,at = c(1)+c(0:5)*7, labels = c("29 feb","7 Mar", "14 Mar", "21 mar", "28 Mar", "4 Apr") )
for (i in 1:20) {
  lines(c(-10,100),c(1,1)*i, col = "gray", lty=2)
}

lines( d/(log(y[tr+d],2)-log(y[tr],2)))
points( d/(log(y[tr+d],2)-log(y[tr],2)), pch = 21, col = 1,  bg =  floor(log(y[tr+d],2)) )

You can see this also in the many plots that go around in the media. Below is an example for Australian provinces. Here the slopes for various doubling rates are plotted explicitly.

multiple

## plot of multiple cases on logarithmic scale
plot(-1,-1, ylim = c(8,2^12), xlim = c(0,40), 
     ylab = "cases", xlab = "days since first 8 cases", log = "y",xaxt = "n", yaxt = "n", 
     main = "plot of multiple territories")

axis(1, 1:31, 0:30)
for (i in 0:12) {
  lines(c(-10,100),c(1,1)*2^i, col = "gray", lty=3)
}
axis(2,2^c(3:12), las = 2)
#axis(1,at = c(1)+c(0:5)*7, labels = c("29 feb","7 Mar", "14 Mar", "21 mar", "28 Mar", "4 Apr") )


lines(c(1,2*14), 8*c(1,2^14), col = "dark gray", lty = 2)
lines(c(1,3*14), 8*c(1,2^14), col = "dark gray", lty = 2)
lines(c(1,4*14), 8*c(1,2^14), col = "dark gray", lty = 2)
lines(c(1,7*14), 8*c(1,2^14), col = "dark gray", lty = 2)

text(log(4096/8,2)*2,4096, "doubling every 2 days", pos = 2, srt = 57, col = "dark gray")
text(log(4096/8,2)*3-0.25,4096, "doubling every 3 days", pos = 2, srt = 45, col = "dark gray")
text(log(4096/8,2)*4-0.5,4096, "doubling every 4 days", pos = 2, srt = 38, col = "dark gray")
text(40,8*2^(40/7)*1.1, "doubling every week", pos = 2, srt = 22, col = "dark gray")

for (province in c("Victoria","Tasmania","South Australia",
                   "Queensland","Northern Territory",
                   "New South Wales", "Australian Capital Territory","Western Australia")){
  # extract country data
  y <- as.numeric(dat[dat$Province.State == province, -c(1:4)])
  tr <- (max(which(y<8))+1):(length(y))  ##trimming data to only the numbers >1

  col = rgb(0.4,0.4,0.4)
  if (province == "Western Australia") {
    col = "red"
  }
  lines(y[tr], pch = 21, col = col,  bg =  floor(log(y[tr+d],2)), lwd = 2)
  if (province == "Western Australia") {
    points(y[tr], pch = 21, col = 1,  bg =  "white")
    text(length(y[tr]), tail(y[tr],1), "Western Australia", col =2, pos = 4)
  }

}

This type of data analysis is just exploratory data analysis. Beware of overinterpreting the differences. There are many potential causes for the differences and the data does not allow us to make strong conclusions about the underlying theories that might explain the differences.

Note that it is possible to make models with the same epidemiological parameters that will turn out different epidemiological curves due to different seeding, distance from the epicenter, random variations, etc.

$\endgroup$
1
  • $\begingroup$ this is a great answer!.. just later than the other answer.. IMO both are valid $\endgroup$
    – baradhili
    Apr 10, 2020 at 3:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.