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))
})
}
RSS.SIR2 <- function(parameters) {
names(parameters) <- c("const", "R0")
out <- ode(y = init, times = day, func = SIR2, parms = parameters)
fit <- out[ , 3]
RSS <- sum((Infected - fit)^2)
return(RSS)
}
### plotting different values R0
# use the ordinary exponential model to determine const = beta - gamma
const <- coef(mod)[2]
RSS.SIR <- function(parameters) {
names(parameters) <- c("beta", "gamma")
out <- ode(y = init, times = day, func = SIR, parms = parameters)
fit <- out[ , 3]
RSS <- sum((Infected - fit)^2)
return(RSS)
}
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
RSS.SIRMC2 <- function(const,R0) {
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
}
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() {
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))
}
# doing the nested model to get RSS
Infected_MC <- Infected
modnested <- getOptim()
rss <- sapply(seq(0.3,0.5,0.01),
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+(modnested$R0==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).
optimfunction 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). $\endgroup$ – Sextus Empiricus Feb 12 '20 at 15:52