---
# System prepended metadata

title: Coronavirus situation and projections in New Zealand
tags: [coronavirus, ' covid19', ' prediction']

---

---
title: "Coronavirus situation and projections in New Zealand"
author: "Arindam Basu"
date: '2020-04-14'
output: word_document
tags: coronavirus, covid19, prediction
---

## Coronavirus predictions in NZ 

New Zealand is one of the few countries that seem to have got a control over coronavirus outbreak. Here we present the incidence of Coronavirus in New Zealand and some predictions as to how and when we can expect to get to a point of 0 incident cases in NZ. The codes are presented along with the conclusions so that you can revise them. Please let me know what you think. 

```{r}
## first load the libraries
library(tidyverse)
library(lubridate)
library(incidence)
library(projections)
library(epitrix)
library(distcrete)

## then obtain the data from JHU CSSE archives on cases

fulldata <- read_csv("https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv")

fulldata %>% head(3) # gives us the headers and some data
```


```{r}
## now select NZ and process NZ data and convert it to an incidence object

## create a New Zealand data object first

nz_data <- fulldata %>%
   rename(province = `Province/State`,
          country = `Country/Region`) %>%
   filter(country == "New Zealand") %>%
   select(-Lat, -Long) %>%
   gather(dates, count, -province, -country) %>%
   mutate(date1 = mdy(dates),
          inc_cases = c(0, diff(count, lag = 1)),
          days_since = as.numeric(date1 - min(date1))) %>%
   select(date1, inc_cases, days_since)

nz_data %>% head()
```


```{r}
nz_inc <- fulldata %>%
   rename(province = `Province/State`,
          country = `Country/Region`) %>%
   filter(country == "New Zealand") %>%
   select(-Lat, -Long) %>%
   gather(dates, count, -province, -country) %>%
   mutate(date1 = mdy(dates),
          inc_cases = c(0, diff(count, lag = 1))) %>%
   select(date1, inc_cases) %>%
   uncount(inc_cases) %>%
   pull(date1) %>%
   incidence()
   
# nz_inc # will provide detailed information about the incidence object
```


```{r}
# plot incidence
plot(nz_inc, border = "white") +
   ggsave("incidence_plot.jpg")

```

## The first plot

![incidence plot](https://i.imgur.com/FQSTFQs.jpg)


This shows that the number of cases increased steadily till the first week of April and then progressively dropped. This suggests that we had reached a peak in the number of cases in New Zealand and then there has  been a progressive drop in the number of new cases we continued to find despite significantly increasing the number of tests. Intuitively, this would mean that we would find more new cases even with somewhat relaxed testing criteria that the Ministry of Health used for testing (instead of tight testing criteria for overseas travel or contact with someone with history of overseas travel and physical symptoms, now the  new criteria included anyone with suspicious symptoms and signs). We have fitted an optimised model to find a peak date and patterns of fitted lines on the existing data. The `find_peak()` function from the incidence package returns that the peak was reached around 4th April, 2020. 

Next, we fitted an optimised Weibull regression model for the data using the `incidence` package and ran the regression model against time with no other coviariates. 

```{r}
peak_day <- find_peak(nz_inc) # find the day of the peak

reg_results <- fit_optim_split(nz_inc) # fit an optimised regression model

fit_info <- reg_results$fit # returns the results of the fit 

# fit_info # will provide details of the model

```

## Results of the regression  model

In the model, we regressed log-incidence over time, splitting the dates into two halves: one before 4th April, 2020, and the ones after 4th April, 2020. Each of these graphs were exponential graphs, and had different growth rates. We obtained the following statistics from the regression model:

| Parameter | Estimate | 95% confidence interval |
|-----------|----------|-------------------------|
| Daily Growth Rate before 4th April | 0.16 | 0.13, 0.10 |
| Daily growth rate after 4th April | - 0.15 | - 0.22, -0.08 |
| Doubling time of infection before 4/4/2020 | 4.31 days | 3.64, 5.28 |
| Doubling (down) time of infection after 4/4/2020 | 4.35 | 3.02, 7.75 |

This gives us several interesting insights into the process of the growth and shrinkage of this epidemic. Prior to 4th April, the number of cases grew rapidly at the rate of doubling every nearly five days; likewise, after the peak period was reached, when the epidemic started on a downhill climb, it did so equally rapidly as evidenced over the next few days, new cases halving roughly every five days. We can graphically examine these trends.

![Fit information](https://i.imgur.com/OduCrRH.jpg)


As can be seen in this graph, the trend of the cases grew rapidly till 2nd April, then there was a plateau, and then the number of cases steadily decreased starting 3rd April. In the absence of any intervention, the number of cases would have gone upward. New Zealand started "lock down" or level 4 of the intervention to "eliminate" the infection in and around 25th March, 2020. This meant that the effects to reduce the new number of cases even with increased surveillance and relaxed criteria took about 10 days to take effect. Note that this roughly or on average the time it takes an infection to manifest from first infection. This suggests that despite limitation in the number of tests being conducted, the tests were on track to identify and estimate the true number of cases. 



```{r}
plot(nz_inc, border = "white",
     fit = reg_results$fit) +
  ggsave("fit_info.jpg")
```
```{r}

## Find the average number of cases and their standard deviation

mu <- mean(nz_inc$counts, na.rm = T) # average number of cases
sigma <- sd(nz_inc$counts, na.rm = T) # standard deviation in the number of cases
cv <- mu / sigma # coefficient of variation


params <- gamma_mucv2shapescale(mu, cv)

serial_interval <- distcrete("gamma",
                             shape = params$shape,
                             scale = params$scale,
                             w = 0,
                             interval = 1)
```

Here are the details of these parameter estimates as of now:

| Parameter | Value |
|-----------|-------|
| Average incident cases | 29.8|
| standard deviation | 33.1 |
| Coefficient of variation | 0.90 |
|

While the numbers continue to decline, to get to a 0 incidence rate, it will be weeks before that can happen. 


```{r}
pred_pre_lockdown <- project(nz_inc[1:30], R = 2.4, 
                             si = serial_interval,
                             model = "poisson",
                             n_days = 20,
                             R_fix_within = T,
                             n_sim = 5)

pred_pre <- as.data.frame(pred_pre_lockdown)
plot(pred_pre_lockdown) +
  ggsave("prelockdown.jpg")

```

If we consider the trend lines based on our estimation of the possible rate of rise in the number of new cases assuming the trend of the growth in the number of cases were to continue as before, then we see from our plotting that by April 16th onwards, we would see around 70 new cases everyday. But we do not see that, instead we have trickled down to around 18 cases in a day (some days may still be high, this is on an average). 

![Prelockdown](https://i.imgur.com/Mc3uRfP.jpg)


It is reasonable to believe that when a country goes for suppression of infection, then the aim is to bring down the effective reproduction number to below 1.0. In theory, if 100% people comply with the lock down provisions, then it should happen instantly. But that may not be the case. So, let's say enough people comply with the provisions so that the effective reproduction number drops down to around 0.90 and the infection is deemed unsustainable. Then, we should expect a slow regression to zero cases eventually. But how slow or how fast? It cannot be as rapid as the cases upswung in the first place. Let's say we study the trend of new cases for days 30 through 37 (in our case this would be something like we study from March 28 through April 6) and try to predict when we can expect the new cases will drop down to zero, so the total case count will steady with no new cases anymore. 

Here is the figure:

![Postlockdown](https://i.imgur.com/9THX8LE.jpg)


As can be seen, that countdown to zero is much slower although there does seem that most days, we will not see many new cases above 15 as does seem to be the case. We assumed an R of 0.70; what happens if we settle for a more conservative R, such as 0.90 which may be a little more realistic?

See the next figure:

![Postlockdown conservative](https://i.imgur.com/ZbhMPg9.jpg)


```{r}
pred_post_lockdown_90 <- project(nz_inc[29:35], R = 0.90,
                              si = serial_interval,
                              model = "poisson",
                              n_days = 30,
                              n_sim = 5)

pred_post <- as.data.frame(pred_post_lockdown_90)

plot(pred_post_lockdown_90) +
  theme_classic() +
  ggsave("postlockdown90.jpg") 
  

```

This gives us a more realistic and perhaps a little more aligned to the numbers we are seeing. The numbers are 20 and then drop down to 10 and five etc on days. We are hopeful but if you study the trajectory, you will see that we are not going to see zero new cases anytime in the middle of April even with effective reproduction number of 0.90

Hence the message is that we need to be patient and continue to practice the conditions of hand and respiratory hygiene, and maintain safe distance, even as the "lock down" is lifted. A rush at this stage is not warranted. 

