Early analysis of Wuhan nCoV-2019 outbreak in Wuhan, China.

This R package implements an ODE-based model of the novel coronavirus outbreak in Wuhan, China. It presents a simulator and likelihood function assuming Poisson-distributed increments in the number of new cases in Wuhan, in the rest of China via the airline network, and to the rest of the world.

Data required:

• china_cases daily case reports in all Chinese cities (see data(package='wuhan'))
• world_cases daily case reports from other countries (see data(package='wuhan'))
• K daily numbers of passengers going between cities in China via airline network, available from OAG Traffic Analyzer
• W daily numbers of passengers going between Chinese cities and other countries via airline network, available from OAG Traffic Analyzer
• china_population the population size in each Chinese city (see data(package='wuhan'))

Parameters:

• beta the human-human basic transmission rate
• gamma the removal rate (inverse of infectious period)
• I0W the number of initial infectives in Wuhan
• phi the case ascertainment rate in Wuhan

To use the package, assume the following workflow in R:

# Load required packages
> install.packages('devtools')
> devtools::install_git('https://github.com/chrism0dwk/wuhan.git')
> library(wuhan)

# Instantiate ODE model, simulate up to day 22.
> simulator = NetworkODEModel(N=china_population, K=K, init_loc='Wuhan', alpha=1/4, max_t=22) 

# Instantiate LogLikelihood function
> llik = LogLikelihood(y=china_cases[,1:22], z=world_cases[,1:22], N=N, K=K, W=W, sim_fun=simulator)

# Find MLEs using optimisation
> par_init = c(0.4, 0.142857142857143, 1, 0.5)  # Starting point
> fit = optim(log(par_init), llik, control=list(fnscale=-1))
> p_hat = fit$par

Asymptotic assumptions for confidence intervals fail in our case, since the parameter space is highly non-orthogonal. Confidence intervals are therefore calculated using parametric bootstrap. p_hat is calculated on the log scale (logit scale for the phi parameter), so needs to be transformed first:

> p_hat[1:3] = exp(p_hat[1:3])
> p_hat[4] = invlogit(p_hat[4])

The samples can then be drawn by bootstrap, for which a computing cluster is highly recommended (thanks Lancaster University HEC facility!).

> samples = bootstrap(p_hat, K, W, alpha=1/4, max_t=22, n_samples=1000)

Since the airline connectivity matrices are not included in this package, samples from the parameters (for 4 different values of the latent period $1/\alpha$) are provided as in-build datasets. See data(package='wuhan').