# Fitting SIR model with 2019-nCoV data doesn't conververge

I am trying to calculate the basic reproduction number $$R_0$$ of the new 2019-nCoV virus by fitting a SIR model to the current data. My code is based on https://arxiv.org/pdf/1605.01931.pdf, p. 11ff:

library(deSolve)
library(RColorBrewer)

#https://en.wikipedia.org/wiki/Timeline_of_the_2019%E2%80%9320_Wuhan_coronavirus_outbreak#Cases_Chronology_in_Mainland_China
Infected <- c(45, 62, 121, 198, 291, 440, 571, 830, 1287, 1975, 2744, 4515)
day <- 0:(length(Infected)-1)
N <- 1400000000 #pop of china
init <- c(S = N-1, I = 1, R = 0)
plot(day, Infected)


SIR <- function(time, state, parameters) {
par <- as.list(c(state, parameters))
with(par, { dS <- -beta * S * I
dI <- beta * S * I - gamma * I
dR <- gamma * I
list(c(dS, dI, dR))
})
}

names(parameters) <- c("beta", "gamma")
out <- ode(y = init, times = day, func = SIR, parms = parameters)
fit <- out[ , 3]
}

lower = c(0, 0)
upper = c(0.1, 0.5)

set.seed(12)
Opt <- optim(c(0.001, 0.4), RSS.SIR, method = "L-BFGS-B", lower = lower, upper = upper)
Opt$message ## [1] "NEW_X" Opt_par <- Opt$par
names(Opt_par) <- c("beta", "gamma")
Opt_par
##      beta     gamma
## 0.0000000 0.4438188

t <- seq(0, 100, length = 100)
fit <- data.frame(ode(y = init, times = t, func = SIR, parms = Opt_par))
col <- brewer.pal(4, "GnBu")[-1]
matplot(fit$time, fit[ , 2:4], type = "l", xlab = "Day", ylab = "Number of subjects", lwd = 2, lty = 1, col = col) points(day, Infected) legend("right", c("Susceptibles", "Infecteds", "Recovereds"), lty = 1, lwd = 2, col = col, inset = 0.05)  R0 <- N * Opt_par[1] / Opt_par[2] names(R0) <- "R0" R0 ## R0 ## 0  I also tried fitting with GAs (as in the paper), also to no avail. My question Am I making any mistakes or is there just not enough data yet? Or is the SIR model too simple? I would appreciate suggestions on how to change the code so that I get some sensible numbers out of it. Addendum I wrote a blog post based on the final model and current data: Epidemiology: How contagious is Novel Coronavirus (2019-nCoV)? • You may want to try a SEIR model (Susceptible-Exposed-Infectious-Recovered), because the nCov virus has a significant incubation period. – dain Jan 28 '20 at 21:18 • @dain: Thank you... how has the SIR model to be enhanced to become a SEIR model? Jan 29 '20 at 1:04 • @vonjd in your blog you write that your fit results in$R_0 \approx 2$. Note that this is only because the optim function terminates too early. Below I am showing a graph with a solution$R_0 \approx 1.005$, which makes a better fit (I have also added the code to get to it by using the itterative solver, which is the most robust). Feb 12 '20 at 15:52 ## 5 Answers There are several points that you can improve in the code ### Wrong boundary conditions Your model is fixed to I=1 for time zero. You can either changes this point to the observed value or add a parameter in the model that shifts the time accordingly. init <- c(S = N-1, I = 1, R = 0) # should be init <- c(S = N-Infected[1], I = Infected[1], R = 0)  ### Unequal parameter scales As other people have noted the equation $$I' = \beta \cdot S \cdot I - \gamma \cdot I$$ has a very large value for $$S \cdot I$$ this makes that the value of the parameter $$\beta$$ very small and the algorithm which checks whether the step sizes in the iterations reach some point will not vary the steps in $$\beta$$ and $$\gamma$$ equally (the changes in $$\beta$$ will have a much larger effect than changes in $$\gamma$$). You can change scale in the call to the optim function to correct for these differences in size (and checking the hessian allows you to see whether it works a bit). This is done by using a control parameter. In addition you might want to solve the function in segregated steps making the optimization of the two parameters independent from each others (see more here: How to deal with unstable estimates during curve fitting? this is also done in the code below, and the result is much better convergence; although still you reach the limits of your lower and upper bounds) Opt <- optim(c(2*coefficients(mod)[2]/N, coefficients(mod)[2]), RSS.SIR, method = "L-BFGS-B", lower = lower, upper = upper, hessian = TRUE, control = list(parscale = c(1/N,1),factr = 1))  more intuitive might be to scale the parameter in the function (note the term beta/N in place of beta) SIR <- function(time, state, parameters) { par <- as.list(c(state, parameters)) with(par, { dS <- -beta/N * S * I dI <- beta/N * S * I - gamma * I dR <- gamma * I list(c(dS, dI, dR)) }) }  ### Starting condition Because the value of $$S$$ is in the beginning more or less constant (namely $$S \approx N$$) the expression for the infected in the beginning can be solved as a single equation: $$I' \approx (\beta \cdot N - \gamma) \cdot I$$ So you can find a starting condition using an initial exponential fit: # get a good starting condition mod <- nls(Infected ~ a*exp(b*day), start = list(a = Infected[1], b = log(Infected[2]/Infected[1])))  ### Unstable, correlation between $$\beta$$ and $$\gamma$$ There is a bit of ambiguity how to choose $$\beta$$ and $$\gamma$$ for the starting condition. This will also make the outcome of your analysis not so stable. The error in the individual parameters $$\beta$$ and $$\gamma$$ will be very large because many pairs of $$\beta$$ and $$\gamma$$ will give a more or less similarly low RSS. The plot below is for the solution $$\beta = 0.8310849; \gamma = 0.4137507$$ However the adjusted Opt_par value $$\beta = 0.8310849-0.2; \gamma = 0.4137507-0.2$$ works just as well: # Using a different parameterization The optim function allows you to read out the hessian > Opt <- optim(optimsstart, RSS.SIR, method = "L-BFGS-B", lower = lower, upper = upper, + hessian = TRUE) > Opt$hessian
b
b  7371274104 -7371294772
-7371294772  7371315619


The hessian can be related to the variance of the parameters (In R, given an output from optim with a hessian matrix, how to calculate parameter confidence intervals using the hessian matrix?). But note that for this purpose you need the Hessian of the log likelihood which is not the same as the the RSS (it differs by a factor, see the code below).

Based on this you can see that the estimate of the sample variance of the parameters is very large (which means that your results/estimates are not very accurate). But also note that the error is a lot correlated. This means that you can change the parameters such that the outcome is not very correlated. Some example parameterization would be:

$$\begin{array}{} c &=& \beta - \gamma \\ R_0 &=& \frac{\beta}{\gamma} \end{array}$$

such that the old equations (note a scaling by 1/N is used):

$$\begin{array}{rccl} S^\prime &=& - \beta \frac{S}{N}& I\\ I^\prime &=& (\beta \frac{S}{N}-\gamma)& I\\ R^\prime &=& \gamma &I \end{array}$$

become

$$\begin{array}{rccl} S^\prime &=& -c\frac{R_0}{R_0-1} \frac{S}{N}& I&\\ I^\prime &=& c\frac{(S/N) R_0 - 1}{R_0-1} &I& \underbrace{\approx c I}_{\text{for t=0 when S/N \approx 1}}\\ R^\prime &=& c \frac{1}{R_0-1}& I& \end{array}$$

which is especially appealing since you get this approximate $$I^\prime = cI$$ for the beginning. This will make you see that you are basically estimating the first part which is approximately exponential growth. You will be able to very accurately determine the growth parameter, $$c = \beta - \gamma$$. However, $$\beta$$ and $$\gamma$$, or $$R_0$$, can not be easily determined.

In the code below a simulation is made with the same value $$c=\beta - \gamma$$ but with different values for $$R_0 = \beta / \gamma$$. You can see that the data is not capable to allow us differentiate which different scenario's (which different $$R_0$$) we are dealing with (and we would need more information, e.g. the locations of each infected individual and trying to see how the infection spread out).

It is interesting that several articles already pretend to have reasonable estimates of $$R_0$$. For instance this preprint Novel coronavirus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions (https://doi.org/10.1101/2020.01.23.20018549)

# Some code:

####
####
####

library(deSolve)
library(RColorBrewer)

#https://en.wikipedia.org/wiki/Timeline_of_the_2019%E2%80%9320_Wuhan_coronavirus_outbreak#Cases_Chronology_in_Mainland_China
Infected <- c(45, 62, 121, 198, 291, 440, 571, 830, 1287, 1975, 2744, 4515)
day <- 0:(length(Infected)-1)
N <- 1400000000 #pop of china

###edit 1: use different boundary condiotion
###init <- c(S = N-1, I = 1, R = 0)
init <- c(S = N-Infected[1], I = Infected[1], R = 0)
plot(day, Infected)

SIR <- function(time, state, parameters) {
par <- as.list(c(state, parameters))
####edit 2; use equally scaled variables
with(par, { dS <- -beta * (S/N) * I
dI <- beta * (S/N) * I - gamma * I
dR <- gamma * I
list(c(dS, dI, dR))
})
}

SIR2 <- function(time, state, parameters) {
par <- as.list(c(state, parameters))
####
#### use as change of variables variable
#### const = (beta-gamma)
#### delta = gamma/beta
#### R0 = beta/gamma > 1
####
#### beta-gamma = beta*(1-delta)
#### beta-gamma = beta*(1-1/R0)
#### gamma = beta/R0
with(par, {
beta  <- const/(1-1/R0)
gamma <- const/(R0-1)
dS <- -(beta * (S/N)      ) * I
dI <-  (beta * (S/N)-gamma) * I
dR <-  (             gamma) * I
list(c(dS, dI, dR))
})
}

names(parameters) <- c("const", "R0")
out <- ode(y = init, times = day, func = SIR2, parms = parameters)
fit <- out[ , 3]
}

### plotting different values R0

# use the ordinary exponential model to determine const = beta - gamma
const <- coef(mod)[2]

names(parameters) <- c("beta", "gamma")
out <- ode(y = init, times = day, func = SIR, parms = parameters)
fit <- out[ , 3]
}

lower = c(0, 0)
upper = c(1, 1)  ###adjust limit because different scale 1/N

### edit: get a good starting condition
mod <- nls(Infected ~ a*exp(b*day),
start = list(a = Infected[1],
b = log(Infected[2]/Infected[1])))
optimsstart <- c(2,1)*coef(mod)[2]

set.seed(12)
Opt <- optim(optimsstart, RSS.SIR, method = "L-BFGS-B", lower = lower, upper = upper,
hessian = TRUE)
Opt

### estimated covariance matrix of coefficients
### note the large error, but also strong correlation (nearly 1)
## note scaling with estimate of sigma because we need to use Hessian of loglikelihood
sigest <- sqrt(Opt$$value/(length(Infected)-1)) solve(1/(2*sigest^2)*Opt$$hessian)

####
####  using alternative parameters
####  for this we use the function SIR2
####

optimsstart <- c(coef(mod)[2],5)
lower = c(0, 1)
upper = c(1, 10^3)  ### adjust limit because we use R0 now which should be >1

set.seed(12)
Opt2 <- optim(optimsstart, RSS.SIR2, method = "L-BFGS-B",lower=lower, upper=upper,
hessian = TRUE, control = list(maxit = 1000,
parscale = c(10^-3,1)))
Opt2

# now the estimated variance of the 1st parameter is small
# the 2nd parameter is still with large variance
#
# thus we can predict beta - gamma very well
# this beta - gamma is the initial growth coefficient
# but the individual values of beta and gamma are not very well known
#
# also note that hessian is not at the MLE since we hit the lower boundary
#
sigest <- sqrt(Opt2$$value/(length(Infected)-1)) solve(1/(2*sigest^2)*Opt2$$hessian)

#### We can also estimated variance by
#### Monte Carlo estimation
##
## assuming data to be distributed as mean +/- q mean
## with q such that mean RSS = 52030
##
##
##

### Two functions RSS to do the optimization in a nested way
parameters <- c(const=const, R0=R0)
out <- ode(y = init, times = day, func = SIR2, parms = parameters)
fit <- out[ , 3]
}
optimize(RSS.SIRMC2, lower=1,upper=10^5,const=const)$objective } getOptim <- function() { opt1 <- optimize(RSS.SIRMC,lower=0,upper=1) opt2 <- optimize(RSS.SIRMC2, lower=1,upper=10^5,const=opt1$$minimum) return(list(RSS=opt2$$objective,const=opt1$$minimum,R0=opt2$$minimum)) } # modeled data that we use to repeatedly generate data with noise Opt_par <- Opt2$$par names(Opt_par) <- c("const", "R0") modInfected <- data.frame(ode(y = init, times = day, func = SIR2, parms = Opt_par))$$I # doing the nested model to get RSS set.seed(1) Infected_MC <- Infected modnested <- getOptim() errrate <- modnested$RSS/sum(Infected)

par <- c(0,0)
for (i in 1:100) {
Infected_MC <- rnorm(length(modInfected),modInfected,(modInfected*errrate)^0.5)
OptMC <- getOptim()
par <- rbind(par,c(OptMC$$const,OptMC$$R0))
}
par <- par[-1,]

plot(par, xlab = "const",ylab="R0",ylim=c(1,1))
title("Monte Carlo simulation")
cov(par)

###conclusion: the parameter R0 can not be reliably estimated

##### End of Monte Carlo estimation

### plotting different values R0

# use the ordinary exponential model to determine const = beta - gamma
const <- coef(mod)[2]
R0 <- 1.1

# graph
plot(-100,-100, xlim=c(0,80), ylim = c(1,N), log="y",
ylab = "infected", xlab = "days", yaxt = "n")
axis(2, las=2, at=10^c(0:9),
labels=c(expression(1),
expression(10^1),
expression(10^2),
expression(10^3),
expression(10^4),
expression(10^5),
expression(10^6),
expression(10^7),
expression(10^8),
expression(10^9)))
axis(2, at=rep(c(2:9),9)*rep(10^c(0:8),each=8), labels=rep("",8*9),tck=-0.02)
title(bquote(paste("scenario's for different ", R[0])), cex.main = 1)

# time
t <- seq(0,60,0.1)

# plot model with different R0
for (R0 in c(1.1,1.2,1.5,2,3,5,10)) {
fit <- data.frame(ode(y = init, times = t, func = SIR2, parms = c(const,R0)))$I lines(t,fit) text(t[601],fit[601], bquote(paste(R[0], " = ",.(R0))), cex=0.7,pos=4) } # plot observations points(day,Infected)  ### How is R0 estimated? The graph above (which is repeated below) showed that there is not much variation in the number of 'infected' as a function of $$R_0$$, and the data of the number of infected people is not providing much information about $$R_0$$ (except whether or not it is above or below zero). However, for the SIR model there is a large variation in the number of recovered or the ratio infected/recovered. This is shown in the image below where the model is plotted not only for the number of infected people but also for the number of recovered people. It is such information (as well additional data like detailed information where and when the people got infected and with whom they had contact) that allows the estimate of $$R_0$$. # Update In your blog article you write that the fit is leading to a value of $$R_0 \approx 2$$. However that is not the correct solution. You find this value only because the optim is terminating early when it has found a good enough solution and the improvements for given stepsize of the vector $$\beta, \gamma$$ are getting small. When you use the nested optimization then you will find a more precise solution with a $$R_0$$ very close to 1. We see this value $$R_0 \approx 1$$ because that is how the (wrong) model is able to get this change in the growth rate into the curve. ### #### #### library(deSolve) library(RColorBrewer) #https://en.wikipedia.org/wiki/Timeline_of_the_2019%E2%80%9320_Wuhan_coronavirus_outbreak#Cases_Chronology_in_Mainland_China Infected <- c(45,62,121,198,291,440,571,830,1287,1975, 2744,4515,5974,7711,9692,11791,14380,17205,20440) #Infected <- c(45,62,121,198,291,440,571,830,1287,1975, # 2744,4515,5974,7711,9692,11791,14380,17205,20440, # 24324,28018,31161,34546,37198,40171,42638,44653) day <- 0:(length(Infected)-1) N <- 1400000000 #pop of china init <- c(S = N-Infected[1], I = Infected[1], R = 0) # model function SIR2 <- function(time, state, parameters) { par <- as.list(c(state, parameters)) with(par, { beta <- const/(1-1/R0) gamma <- const/(R0-1) dS <- -(beta * (S/N) ) * I dI <- (beta * (S/N)-gamma) * I dR <- ( gamma) * I list(c(dS, dI, dR)) }) } ### Two functions RSS to do the optimization in a nested way RSS.SIRMC2 <- function(R0,const) { parameters <- c(const=const, R0=R0) out <- ode(y = init, times = day, func = SIR2, parms = parameters) fit <- out[ , 3] RSS <- sum((Infected_MC - fit)^2) return(RSS) } RSS.SIRMC <- function(const) { optimize(RSS.SIRMC2, lower=1,upper=10^5,const=const)$objective
}

# wrapper to optimize and return estimated values
getOptim <- function() {
opt2 <- optimize(RSS.SIRMC2, lower=1,upper=10^5,const=opt1$$minimum) return(list(RSS=opt2$$objective,const=opt1$$minimum,R0=opt2$$minimum))
}

# doing the nested model to get RSS
Infected_MC <- Infected
modnested <- getOptim()

FUN = function(x) optimize(RSS.SIRMC2, lower=1,upper=10^5,const=x)$objective) plot(seq(0.3,0.5,0.01),rss) optimize(RSS.SIRMC2, lower=1,upper=10^5,const=0.35) # view modnested ### plotting different values R0 const <- modnested$$const R0 <- modnested$$R0 # graph plot(-100,-100, xlim=c(0,80), ylim = c(1,6*10^4), log="", ylab = "infected", xlab = "days") title(bquote(paste("scenario's for different ", R[0])), cex.main = 1) ### this is what your beta and gamma from the blog beta = 0.6746089 gamma = 0.3253912 fit <- data.frame(ode(y = init, times = t, func = SIR, parms = c(beta,gamma)))$I
lines(t,fit,col=3)

# plot model with different R0
t <- seq(0,50,0.1)
for (R0 in c(modnested$$R0,1.07,1.08,1.09,1.1,1.11)) { fit <- data.frame(ode(y = init, times = t, func = SIR2, parms = c(const,R0)))I lines(t,fit,col=1+(modnestedR0==R0)) text(t[501],fit[501], bquote(paste(R[0], " = ",.(R0))), cex=0.7,pos=4,col=1+(modnested$$R0==R0))
}

# plot observations
points(day,Infected, cex = 0.7)


If we use the relation between recovered and infected people $$R^\prime = c (R_0-1)^{-1} I$$ then we also see the opposite, namely a large $$R_0$$ of around 18:

I <- c(45,62,121,198,291,440,571,830,1287,1975,2744,4515,5974,7711,9692,11791,14380,17205,20440, 24324,28018,31161,34546,37198,40171,42638,44653)
D <- c(2,2,2,3,6,9,17,25,41,56,80,106,132,170,213,259,304,361,425,490,563,637,722,811,908,1016,1113)
R <- c(12,15,19,25,25,25,25,34,38,49,51,60,103,124,171,243,328,475,632,892,1153,1540,2050,2649,3281,3996,4749)
A <- I-D-R

plot(A[-27],diff(R+D))
mod <- lm(diff(R+D) ~ A[-27])


giving:

> const
[1] 0.3577354
> const/mod$coefficients[2]+1 A[-27] 17.87653  This is a restriction of the SIR model which models $$R_0 = \frac{\beta}{\gamma}$$ where $$\frac{1}{\gamma}$$ is the period how long somebody is sick (time from Infected to Recovered) but that may not need to be the time that somebody is infectious. In addition, the compartment models is limited since the age of patients (how long one has been sick) is not taken into account and each age should be considered as a separate compartment. But in any case. If the numbers from wikipedia are meaningfull (they may be doubted) then only 2% of the active/infected recover daily, and thus the $$\gamma$$ parameter seems to be small (no matter what model you use). • @vondj the coefficients(mod)[2] is the exponential factor that relates to $$I' \approx \underbrace{(\beta \cdot N - \gamma)}_{\text{factor in exponential fit}} \cdot I$$ So you have$\text{factor} = \beta \cdot N - \gamma$And I chose here$\beta = 2\text{factor}/N$and$\gamma = \text{factor}$. The parameter coefficients(mod)[1] is the starting point (the height) of the exponential curve. Jan 28 '20 at 16:49 • extra note: I like to scale the parameters in the called function instead of via the parscale option. The parscale option does not scale the output of the hessian, which is a bit annoying. Jan 28 '20 at 17:20 • But keep in mind that those two parameters individually have a large error. It will be something like a linear combination of the two (something close to$\beta-\gamma$or$\beta/N -\gamma $, depending on the scaling) which will be approximated with high accuracy. Jan 28 '20 at 17:59 • The plot on the right is log-scale. Jan 28 '20 at 18:02 • Related to my last addition (including the recovery cases in the calculation) a very quick way to assess$R_0$is to use $$-\frac {S^\prime}{R^\prime} = R_0 \frac {S}{N} \approx R_0$$ and note that the rate at which people become infected is much larger than the rate at which people recover. So based on those data$R_0\$ should be considered high. This requires no fitting. Feb 12 '20 at 17:19

You might be experiencing numerical issues due to the very large population size $$N$$, which will force the estimate of $$\beta$$ to be very close to zero. You could re-parameterise the model as \begin{align} {\mathrm d S \over \mathrm d t} &= -\beta {S I / N}\\[1.5ex] {\mathrm d I \over \mathrm d t} &= \beta {S I / N} - \gamma I \\[1.5ex] {\mathrm d R \over \mathrm d t} &= \gamma I \\ \end{align}

This will make the estimate of $$\beta$$ larger so hopefully you'll get something more sensible out of the optimisation.

In this context the SIR model is useful but it only gives a very crude fit to these data (it assumes that the whole population of China mixes homogenously). It's perhaps not too bad as a first attempt at analysis. Ideally you would want some kind of spatial or network model that would better reflect the true contact structure in the population. For example, a metapopulation model as described in Program 7.2 and the accompanying book (Modeling Infectious Diseases in Humans and Animals, Keeling & Rohani). However this approach would require much more work and also some data on the population structure. An approximate alternative could be to replace the $$I$$ in $$\beta SI/N$$ (in both of the first two equations) with $$I^\delta$$ where $$\delta$$, which is probably $$<1$$, is a third parameter to be estimated. Such a model tries to capture the fact that the force of infection on a susceptible increases less than linearly with the number of infecteds $$I$$, while avoiding specification of an explicit population structure. For more details on this approach, see e.g. Hochberg, Non-linear transmission rates and the dynamics of infectious disease, Journal of Theoretical Biology 153:301-321.

• Thank you... any suggestions to make this model work? Jan 28 '20 at 13:29
• @vonjd See edits above Jan 28 '20 at 14:44
• Thank you again... so you have to change that in the first two equations, right? Jan 28 '20 at 14:53
• @vonjd Yes, see further edits Jan 28 '20 at 14:57

Because the population of china is so huge, the parameters will be very small.

Since we are in the early days of the infection, and because N is so big, then $$S(t)I(t)/N \ll 1$$. It could me more reasonable to assume that at this stage of the infection, the number of infected people is approximately exponential, and fit a much simpler model.

• Thank you... could you please be more specific how to change the code? Jan 28 '20 at 16:30
• @vonjd I'm thinking about how to do it myself. I don't think the code is a problem, I think you will have to rethink the model. Jan 28 '20 at 16:33

This is only marginally related to the detailed coding discussion, but seems highly relevant to the original question concerning modeling of the current 2019-nCoV epidemic. Please see arxiv:2002.00418v1 (paper at https://arxiv.org/pdf/2002.00418v1.pdf ) for a delayed diff equation system ~5 component model, with parameter estimation and predictions using dde23 in MatLab. These are compared to daily published reports of confirmed cases, number cured, etc. To me, it is quite worthy of discussion, refinement, and updating. It concludes that there is a bifurcation in the solution space dependent upon the efficacy of isolation, thus explaining the strong public health measures recently taken, which have a fair chance of success so far.

• Warm welcome to you. Is it possible to add a brief summary of the paper - just a few lines, for ease of future reference. Feb 9 '20 at 9:27

what do you think about putting the initial number of infectious as an addition parameter in the optimization problem otherwise the fitting need to start with the initial condition.

• That works. And not just the initial number of infectious people, but also the initial number of susceptible people should be an additional fitting parameter. See here in a follow up on the question blog.ephorie.de/… Mar 31 '20 at 18:43
• Maybe this answer would have been better as a comment. Mar 31 '20 at 18:50