7
$\begingroup$

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)

enter image description here

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))
  })
}

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(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)

enter image description here

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)?

$\endgroup$
  • 2
    $\begingroup$ You may want to try a SEIR model (Susceptible-Exposed-Infectious-Recovered), because the nCov virus has a significant incubation period. $\endgroup$ – dain Jan 28 at 21:18
  • $\begingroup$ @dain: Thank you... how has the SIR model to be enhanced to become a SEIR model? $\endgroup$ – vonjd Jan 29 at 1:04
  • 1
    $\begingroup$ @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). $\endgroup$ – Sextus Empiricus Feb 12 at 15:52
14
$\begingroup$

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 $

example solution

However the adjusted Opt_par value $\beta = 0.8310849-0.2; \gamma = 0.4137507-0.2$ works just as well:

example variation


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)

fitting with different R_0


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$.

exmple

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.

fit close to 1

###
####
####

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).

$\endgroup$
  • 1
    $\begingroup$ @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. $\endgroup$ – Sextus Empiricus Jan 28 at 16:49
  • 1
    $\begingroup$ 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. $\endgroup$ – Sextus Empiricus Jan 28 at 17:20
  • 1
    $\begingroup$ 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. $\endgroup$ – Sextus Empiricus Jan 28 at 17:59
  • 1
    $\begingroup$ The plot on the right is log-scale. $\endgroup$ – Sextus Empiricus Jan 28 at 18:02
  • 1
    $\begingroup$ 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. $\endgroup$ – Sextus Empiricus Feb 12 at 17:19
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you... any suggestions to make this model work? $\endgroup$ – vonjd Jan 28 at 13:29
  • 1
    $\begingroup$ @vonjd See edits above $\endgroup$ – S. Catterall Reinstate Monica Jan 28 at 14:44
  • $\begingroup$ Thank you again... so you have to change that in the first two equations, right? $\endgroup$ – vonjd Jan 28 at 14:53
  • 1
    $\begingroup$ @vonjd Yes, see further edits $\endgroup$ – S. Catterall Reinstate Monica Jan 28 at 14:57
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you... could you please be more specific how to change the code? $\endgroup$ – vonjd Jan 28 at 16:30
  • 1
    $\begingroup$ @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. $\endgroup$ – Demetri Pananos Jan 28 at 16:33
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Warm welcome to you. Is it possible to add a brief summary of the paper - just a few lines, for ease of future reference. $\endgroup$ – user2696565 Feb 9 at 9:27
1
$\begingroup$

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.

New contributor
Filip Floegel is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
  • $\begingroup$ 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/… $\endgroup$ – Sextus Empiricus 2 days ago
  • $\begingroup$ Maybe this answer would have been better as a comment. $\endgroup$ – Sextus Empiricus 2 days ago

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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