# required libraries library(deSolve) library(shape) # for plotting arrows library(progress) # for drawing the progress bar #################################### ## ## The basic model (which does not result in the exact solution) ## Adapted from the previous blogpost but with small adaptations ## #################################### # the data in China #https://en.wikipedia.org/wiki/Timeline_of_the_2019%E2%80%9320_Wuhan_coronavirus_outbreak#Cases_Chronology_in_Mainland_China #https://en.wikipedia.org/wiki/Template:2019%E2%80%9320_coronavirus_pandemic_data/China_medical_cases IpRpD <- 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,59804,63851,66492,68500,70548,72436,74185,74576,75465,76288,76936,77150,77658,78064,78497,78824,79251,79824,80026,80151,80270,80409,80552,80651,80695,80735,80754,80778,80793,80813,80824,80844,80860,80881,80894,80928,80967,81008) R <- c(15,19, 24, 25, 25, 25, 25, 34, 38, 49, 51, 60, 103, 124, 171, 243, 328, 475, 632, 892, 1153, 1540, 2050, 2649, 3281, 3996, 4749, 5911, 6723, 8096, 9419,10844,12552,14376,16155,18264,20659,22888,24734,27323,29745,32495,36117,39002,41625,44462,47204,49856,52045,53726,55404,57065,58600,59897,61475,62793,64111,65541,66911,67749,68679,69601,70420,71150,71740) 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, 1259, 1380, 1523, 1665, 1770, 1868, 2004, 2118, 2236, 2345, 2442, 2592, 2663, 2715, 2744, 2788, 2835, 2870, 2912, 2943, 2981, 3012, 3042, 3070, 3097, 3119, 3136, 3158, 3169, 3176, 3189, 3199, 3213, 3226, 3237, 3245, 3248, 3255) # we take the first 19 to make the analysis relate to the previous blog post Infected <- (IpRpD)[1:19] Day <- 1:(length(Infected)) N <- 1400000000 # population of mainland china # ODE equation used for fitting # # I have removed the R(t) in comparison # to the function used in the odler blogpost # because we are not gonna use that value # also we have anyway: R(t) = N(0) - N(t) - I(t) SIR <- function(time, state, parameters) { par <- as.list(c(state, parameters)) with(par, { dS <- -beta/N * I * S dI <- beta/N * I * S - gamma * I list(c(dS, dI)) }) } # # cost function to be optimized in the fitting # RSS <- function(parameters) { names(parameters) <- c("beta", "gamma") out <- ode(y = init, times = Day, func = SIR, parms = parameters) fitInfected <- out[,3] # fitInfected <- N-out[,2] # this would be a better comparison since the data is not the number of Infectious people sum((Infected - fitInfected)^2) } # starting condition init <- c(S = N-Infected[1], I = Infected[1]) # init <- c(S = N-Infected[1], I = Infected[1]-R[1]-D[1]) use this starting condition when applying the different line in the RSS function # performing the fit Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, 1)) # optimize with some sensible conditions Opt$message ## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" Opt_par <- setNames(Opt$par, c("beta", "gamma")) Opt_par ## beta gamma ## 0.693018 0.306982 # plotting the result t <- 1:70 # time in days fit <- data.frame(ode(y = init, times = t, func = SIR, parms = Opt_par)) plot(Day,Infected, xlim = c(0,30), ylim = c(0,3*10^4)) lines(t,fit[,3]) ########################### ## ## Alternative model which provides a better fit ## ############################ # We transform the equations and instead of # parameters beta and gamma # we use parameters # # K = beta-gamma # R0 = beta/gamma # # or # # beta = K * R0/(R0-1) # gamma = K * 1/(R0-1) # # then the equations become # # dS = I * K * (-S/N * R0)/(R0-1) # dI = I * K * ( S/N * R0 - 1)/(R0-1) # note in the beginning, S/N = 1 # then in the start you get this approximate exponential growth # dI = I * K * (1) SIR2 <- function(time, state, parameters) { par <- as.list(c(state, parameters)) with(par, { dS <- I * K * (-S/N * R0/(R0-1)) dI <- I * K * ( S/N * R0/(R0-1) - 1/(R0-1)) list(c(dS, dI)) }) } RSS2 <- function(parameters) { names(parameters) <- c("K", "R0") out <- ode(y = init, times = Day, func = SIR2, parms = parameters) fitInfected <- out[,3] #fitInfected <- N-out[,2] sum((Infected - fitInfected)^2) } ### Two functions RSS to do the optimization in a nested way ### ### This nesting requires a lot more computational power ### However, it makes that we have to worry less about the different scale ### of the parameters Infected_MC <- Infected SIRMC2 <- function(R0,K) { parameters <- c(K=K, R0=R0) out <- ode(y = init, times = Day, func = SIR2, parms = parameters) fitInfected <- out[,3] #fitInfected <- N-out[,2] RSS <- sum((Infected_MC - fitInfected)^2) return(RSS) } SIRMC <- function(K) { optimize(SIRMC2, lower=1,upper=10^5,K=K, tol = .Machine$double.eps)$objective } # wrapper to optimize and return estimated values getOptim <- function() { opt1 <- optimize(SIRMC,lower=0,upper=1, tol = .Machine$double.eps) opt2 <- optimize(SIRMC2, lower=1,upper=10^5,K=opt1$minimum, tol = .Machine$double.eps) return(list(RSS=opt2$objective,K=opt1$minimum,R0=opt2$minimum)) } # starting condition #init <- c(S = N-Infected[1], I = Infected[1]-R[1]-D[1]) init <- c(S = N-Infected[1], I = Infected[1]) # performing the fit # starting K=0.3, R0 = 2 Opt2 <- optim(c(0.3, 3), RSS2, method = "L-BFGS-B", hessian = TRUE, control = list(parscale = c(10^0,10^0), factr = 1)) Opt2 Opt3 <- getOptim() Opt3 Opt_par2 <- setNames(Opt2$par, c("K", "R0")) Opt_par3 <- setNames(Opt3[2:3], c("K", "R0")) # plotting the result t <- seq(1,70,0.1) # time in days fit1 <- data.frame(ode(y = init, times = t, func = SIR , parms = Opt_par)) fit2 <- data.frame(ode(y = init, times = t, func = SIR2, parms = Opt_par2)) fit3 <- data.frame(ode(y = init, times = t, func = SIR2, parms = Opt_par3)) plot(Day,Infected, xlim = c(0,30), ylim = c(0,3*10^4), log = "", xaxt = "n", main = "Infected in China (including Recovered and Death)", xlab = "", ylab = "number infected") lines(t, fit3[,3], col = 1) lines(t, fit2[,3], col = 4, lty = 2) lines(t, fit1[,3], col = 2, lty = 3) axis(1, at = 1:30, labels = rep("",30), tck = -0.01) axis(1, at = c(1,8,15,22), labels = c("Jan 16", "Jan 23", "Jan 30", "Feb 6")) text(t[183]+2,fit1[183,3]+1800,"old optim fit",pos=4, col=2) text(t[183]+2,fit1[183,3],expression(R[0] == 2.07),pos=4, col=2) text(t[183]+2,fit1[183,3]-1400,expression(RSS == 74.3 %*% 10^6),pos=4, col=2) text(t[220]+3,fit2[220,3]+3200,"new optim fit",pos=3, col=4) text(t[220]+3,fit2[220,3]+1400,expression(R[0] == 1.0054626),pos=3, col=4) text(t[220]+3,fit2[220,3],expression(RSS == 6.5 %*% 10^6),pos=3, col=4) text(t[240]-3,fit3[240,3],"nested algorithm",pos=1, col=1) text(t[240]-3,fit3[240,3]+700-2500,expression(R[0] == 1.005332),pos=1, col=1) text(t[240]-3,fit3[240,3]-700-2500,expression(RSS == 5.9 %*% 10^6),pos=1, col=1) x1 <- t[240]-3; x2 <- t[225]; y1 <- fit3[240,3]; y2 <- fit3[225,3] Arrows(x1,y1,x1+(x2-x1)*0.65,y1+(y2-y1)*0.65, col = 1) x1 <- t[220]+2; x2 <- t[227]; y1 <- fit2[220,3]; y2 <- fit2[227,3] Arrows(x1,y1,x1+(x2-x1)*0.6,y1+(y2-y1)*0.6, col = 4) x1 <- t[183]+2; x2 <- t[183]; y1 <- fit1[183,3]; y2 <- fit1[183,3] Arrows(x1,y1,x1+(x2-x1)*0.6,y1+(y2-y1)*0.6, col = 2) ################# ## ## Plotting the cost function and how ## the optimize approaches the minimum ## ################ # use this value to select between the different figures of the blogpost figure = 2 # compute loss function on a grid n=101 K <- seq(0.25,0.45,length.out=n) R0 <- seq(1.001,1.03,length.out=n) #be sure not to use R0= 1 this gives a devision by 0 if (figure == 2) { K <- seq(0.34920,0.34924,length.out=n) R0 <- seq(1.8,2.01,length.out=n) } z <- matrix(rep(0,n^2),n) for (i in 1:n) { for(j in 1:n) { z[i,j] <- log(RSS2(c(K[i],R0[j]))) } } # levels for plotting levels <- log(10^seq(6,12,0.5)) key <- head(c(t(matrix(c(10^c(6:12),rep("",7)),7))),-1) if (figure == 2) { levels <- seq(18.12349,18.12354,0.000005) key <- head(c(t(matrix(c(seq(18.12349,18.12354,0.00001),rep("",5)),5))),-1) } key # colours for plotting colours <- function(n) {hsv(c(seq(0.15,0.7,length.out=n),0), c(seq(0.2,0.4,length.out=n),0), c(seq(1,1,length.out=n),0.9))} # empty plot ylims = c(1,1.03) if (figure == 2) { ylims = c(1.8,2.01) } plot(-1000,-1000, xlab="K",ylab=expression(R[0]), xlim=range(K), ylim = ylims #ylim=range(R0) ) title("RSS as function of parameters \n and path of optimization algorithm",cex.main=1) # add contours .filled.contour(K,(R0),z, col=colours(length(levels)), levels=levels) contour(K,(R0),z,add=1, levels=levels, labels = key) # compute path # start value new=c(0.3, 1.025) if (figure == 2) { new = c(0.34920,2) # with this starting condition convergence will be more difficult } # make lots of small steps for (i in 1:20) { ### safe old value old <- new ### compute 1 step OptTemp <- optim(new, RSS2, method = "L-BFGS-B", lower = c(0,1.00001), hessian = TRUE, control = list(parscale = c(10^0,10^0), factr = 1, maxit=1)) OptTemp new <- OptTemp$par ### plot path dist <- (old[1]-new[1])^2+(old[2]-new[2])^2 lines(c(old[1],new[1]),(c(old[2],new[2])), col=2) points(new[1],new[2], col=2, bg = 2, pch = 21, cex = 0.7) if (dist>0.01^2) { step <- new-old Arrows(old[1], old[2], old[1]+step[1]*0.85,old[2]+step[2]*0.85, col=2, arr.type="curved", arr.length=0.25, lty=1) } } OptTemp # make lots of small steps again (but different parscale) for (i in 1:20) { ### safe old value old <- new ### compute 1 step OptTemp <- optim(new, RSS2, method = "L-BFGS-B", lower = c(0,1.00001), hessian = TRUE, control = list(parscale = c(10^-4,10^-4), factr = 1, maxit=1)) # alternatively use: ndeps = c(10^-7,10^-7), parscale = c(1,1) OptTemp new <- OptTemp$par ### plot path dist <- (old[1]-new[1])^2+(old[2]-new[2])^2 lines(c(old[1],new[1]),(c(old[2],new[2])), col=1) points(new[1],new[2], col=1, bg = 1, pch = 21, cex = 0.7) if (dist>0.01^2) { step <- new-old Arrows(old[1], old[2], old[1]+step[1]*0.85,old[2]+step[2]*0.85, col=1, arr.type="curved", arr.length=0.25, lty=1) } } OptTemp if (figure == 1) { legend(0.255,1.003, c("path with parscale 1","path with parscale 0.0001"), col = c(2,1), pch = 21, lty=1, cex = 0.9, pt.bg = c(2,1) , yjust = 0) } if (figure == 2) { legend(0.349201,1.81, c("path with parscale 1","path with parscale 0.0001"), col = c(2,1), pch = 21, lty=1, cex = 0.9, pt.bg = c(2,1) , yjust = 0) } #################### ## ## Graph with various values of R0 ## ####################### # starting condition #init <- c(S = N-Infected[1], I = Infected[1]-R[1]-D[1]) init <- c(S = N-Infected[1], I = Infected[1]) Infected_MC <- Infected SIRMC3 <- function(R0,K) { parameters <- c(K=K, R0=R0) out <- ode(y = init, times = Day, func = SIR2, parms = parameters) fitInfected <- out[,3] #fitInfected <- N-out[,2] RSS <- sum((Infected_MC - fitInfected)^2) return(RSS) } plot(Day,Infected, xlim = c(0,80), ylim = c(1,10^9), log = "y", xaxt = "n", main = "scenario's for different R0", xlab = "", ylab = "number infected") axis(1, at = 1:30, labels = rep("",30), tck = -0.01) axis(1, at = c(1,8,15,22), labels = c("Jan 16", "Jan 23", "Jan 30", "Feb 6")) for (i in 1:10) { R0 <- c(1.005,1.01,1.05,1.1,1.2,1.5,2,2.5,4,20)[i] K <- optimize(SIRMC3, lower=0,upper=1,R0=R0, tol = .Machine$double.eps)$minimum parameters <- c(K,R0) xd <- seq(1,60,0.01) if (i == 1) { xd <- seq(1,40,0.01) } if (i == 2) { xd <- seq(1,50,0.01) } out <- ode(y = init, times = xd, func = SIR2, parms = parameters) lines(xd,out[,3]) text(tail(xd,1),tail(out[,3],1),bquote(R[0] == .(R0)), pos =4) } ##################### ## ## computing variance with Hessian ## ################### ### The output of optim will store valures for RSS and the hessian mod <- optim(c(0.3, 1.04), RSS2, method = "L-BFGS-B", hessian = TRUE, control = list(parscale = c(10^-4,10^-4), factr = 1)) # unbiased estimate of standard deviation # we divide by n-p # where n is the number of data points # and p is the number of estimated parameters sigma_estimate <- sqrt(mod$value/(length(Infected)-2)) # compute the inverse of the hessian # The hessian = the second order partial derivative of the objective function # in our case this is the RSS # we multiply by 1/(2 * sigma^2) covpar <- solve(1/(2*sigma_estimate^2)*mod$hessian) covpar # [,1] [,2] #[1,] 1.236666e-05 -2.349611e-07 #[2,] -2.349611e-07 9.175736e-09 ## the marginal variance of R0 is then approximately at least ## covpar[2,2]^0.5 ##################### ## ## computing variance with Monte Carlo method ## ################### set.seed(1) # modeled data that we use to repeatedly generate data with noise Infected_MC <- Infected Day <- c(1:length(Infected)) modMC <- getOptim() beginpar <- as.numeric(modMC[2:3]) #getOptim is the nested fitting model defined earlier names(beginpar) <- c("K", "R0") beginpar modInfected <- data.frame(ode(y = init, times = Day, func = SIR2, parms = beginpar))$I # use RSS to make a relative noise level errrate <- modMC$RSS/sum(Infected) ### generate random data with ### the mean of modInfected and ### Noise based on errrate ### Note: this assumes heterogeneous noise ### the Fisher Matrix assumed homogeneous noise (which is probabily less correct) par <- c(0,0) # this will store the parameter outcome and grow using rbind for (i in 1:100) { Infected_MC <- rnorm(length(modInfected),modInfected,(modInfected*errrate)^0.5) OptMC <- getOptim() par <- rbind(par,c(OptMC$K,OptMC$R0)) } par <- par[-1,] # strip the first c(0,0) entry which(par[,2]>2) # some fits went wrong par <- par[-c(36,55),] layout(matrix(1:2,1)) plot(par, xlab = "K",ylab=expression(R[0])) title("Monte Carlo simulation") cov(par) # [,1] [,2] #[1,] 4.461422e-04 -1.709287e-05 #[2,] -1.709287e-05 1.599932e-06 beta <- par[,1]*par[,2]/(par[,2]-1) gamma <- par[,1]/(par[,2]-1) plot(beta,gamma, xlab = expression(beta), ylab = expression(gamma)) ########################### ## ## The alternative model 2 ## ############################ ## ## 1 Now we make the starting number of infected I(0) and population N(0) ## fitting parameters ## (the population N(0) is fitted indirectly by scaling the outcome by a factor ## note that the equations allow scaling all variables by a fixed factor) ## dS = I * K * (-S/N * R0)/(R0-1) ## dI = I * K * ( S/N * R0 - 1)/(R0-1) ## d(f*S) = (f*I) * K * (-(f*S)/(f*N) * R0)/(R0-1) ## d(f*I) = (f*I) * K * ( (f*S)/(f*N) * R0 - 1)/(R0-1) ## 2 Also we consider the 'Infected' from the data cumulatively ## The data infected is not the infected from the model but ## the number of discovered people ## 3 We introduce a factor to scale the number of computed infected with the observed infected ## This scaling corrects for the idea that not all infectious/infected people are reported ## A 'problem' with this scaling is that it also relates to a different population size N(0) ## So this fitting will not allow to determine the effective population size and the over reporting independently Infected <- IpRpD Day <- 1:(length(Infected)) N <- 1400000000 # population of mainland china Ir <- Infected[1] Days <- 1:length(Infected) #### Some exploratory stuff to estimate starting values: #### This is based on the differential equation in #### Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic # model and of the SIR model with equal death and birth rates # from Harko, Lobo and Mak #x0 = N exp(beta/gamma N3) #ln(x0) = ln(N) + beta/gamma N3 #S/x0 = u #ln(u) = ln(S/N) - beta/gamma N3 #dS/dt = x0/phi = S(beta (N3-N) -gamma ln(S/N) + beta S) # function to estimate parameters based on linear least squares fit # and the *linear* relation between dS and S, S ln(S) and S^2 RSSdiff <- function(N) { corrInfected <- Infected * c(rep(1.22,27),rep(1,38)) dS <- -diff(corrInfected) S <- (N-corrInfected[-1]) moddif <- lm(dS ~ 0+I(S)+I(S*log(S/N))+I(S*S)) sum(moddif$residuals^2) # RSS #moddif$coefficients[3]/moddif$coefficients[2]*N # beta } RSSdiff <- Vectorize(RSSdiff) N <- 10:200*10^4 plot(N,RSSdiff(N),log="") corrInfected <- Infected * c(rep(1.22,27),rep(1,38)) N <- 320000 dS <- diff(corrInfected) S <- (N-corrInfected[-1]) Sp <- S Slog <- S*log(S/N) S2 <- S*S moddif <- lm(dS ~ 0+Sp+Slog+S2) moddif sum(moddif$residuals)^2 beta <- moddif$coefficients[3] gamma <- -moddif$coefficients[2] Rstart <- moddif$coefficients[1]/beta+N # this is an unrealistic Rstart = -130 beta <- beta*N #scaling due to different use of beta in Harko et al beta/gamma -(beta - gamma) Rstart plot(diff(corrInfected)) lines(predict(moddif)) # matrix used for computing the scaling # we use a different sacaling after the 27th day when reporting changed X <- cbind(c(rep(1,27),rep(0,max(Days)-27)), c(rep(0,27),rep(1,max(Days)-27))) SIR4 <- function(time, state, parameters) { par <- as.list(c(state, parameters)) with(par, { dS <- I * K * (-S/Nm * R0/(R0-1)) dI <- I * K * ( S/Nm * R0/(R0-1) - 1/(R0-1)) list(c(dS, dI)) }) } parameters <- c(0.4,4) Nm <- 10^5 #using a smaller population than entire China Istart = 30 # inner optimization loop # here we optimize K and R0 for a given Nm and Istart RSS4 <- function(parameters,Nm,Istart) { parametersx <- c(parameters, Nm, Istart) names(parametersx) <- c("K", "R0", "Nm", "Istart") start <- c(parametersx[3]-Infected[1],parametersx[4]) names(start) <- c("S","I") out <- ode(y = start, times = Day, func = SIR4, parms = parametersx[c(1,2,3)]) fitInfected <- parametersx[3]-out[,2] ## perform additional scaling (you could see this as yet another nested optimization) Xm <- X*fitInfected # the limit for a and b is because we should not scale such that the infected become more than obeserved # we introduced I(0) because of uncertainty in R(0) modscaled <- nls(Infected ~ a*Xm[,1]+ b* Xm[,2], start = list(a=0.5*Infected[1]/Istart,b=0.5*Infected[1]/Istart), upper = c(1,2)*Infected[1]/Istart,algorithm= "port") sum((Infected - predict(modscaled))^2) } # outer optimization loop # here we optimize Istart # for each Istart we optimize K and R0 # this is very slow but much more robust RSS4b <- function(parameters) { Istart <- parameters[1] RSS4start <- c(0.3,2.5) OptInner <- optim(RSS4start, RSS4, method = "L-BFGS-B", Istart = Istart, Nm = Nm, lower = c(0.1,1.001), upper = c(1,20), hessian = TRUE, control = list(parscale = c(10^-6,10^-6), factr = 1, maxit= 10^4)) OptInner$value } ### Doing the outer optimization loop by Optimize OptOuter <- optimize(RSS4b,interval = c(5,45) ) OptOuter ### Doing the outer optimization loop by manual grid search n=30 Nm <- 10^5 #using a smaller population than entire China rangeIstart <- 10.83284 + seq(-5,+30,length.out=n) zRSS <- c(NA,length(rangeIstart)) # empty container for values to be computed pb <- progress_bar$new(total = n) for (i in 1:n) { pb$tick() zRSS[i] <- RSS4b(rangeIstart[i]) } # plotting RSS as function of I(0) plot(rangeIstart,zRSS, log="y", yaxt = "n", xlab = "I(0)", ylab = "RSS") title("optimal RSS as function of parameters I(0)",cex.main=1) ### plot a fit with the data for some selected parameters # Doing model with Nm = 10^5 and Istart = 12 # But multiplied by a factor # In this way we end up with a scaling in the end of 1 # This is a bit exchangable, correcting for underreporting and for population size is similar Nm <- 10^5*4.155 Istart <- 10.83284*4.155 RSS4start <- c(0.24,2.5) Opt4 <- optim(RSS4start, RSS4, method = "L-BFGS-B", Nm = Nm, Istart = Istart, lower = c(0.1,1.001), upper = c(1,20), hessian = TRUE, control = list(parscale = c(10^-5,10^-4), ndeps = rep(10^-6,2), factr = 1000, maxit= 10^4, trace = 2 )) Opt4 # this may give an abnormal termination # but that could be due to the small factr parameter # Retrieve optimal parameters and compute curve parameters <- c(Opt4$par,Nm,Istart) names(parameters) <- c("K", "R0", "Nm","Istart") start <- c(parameters[3]-Infected[1],parameters[4]) names(start) <- c("S","I") result <- ode(y = start, times = Day, func = SIR4, parms = parameters[c(1,2,3)]) parameters # scaling fitInfected <- Nm-result[,2] Xm <- X*fitInfected modscaled <- lm(Infected ~ 0 + Xm) modscaled modscaled plot(Infected/Nm*100, log = "", xlim = c(0,70), ylim = c(0,100), xlab = "day", ylab = "Total people that got sick [%]", xaxt="n") axis(1, at = c(1,8,15,22,29,36,43), labels = c("Jan 16", "Jan 23", "Jan 30", "Feb 6", "Feb 13", "Feb20", "Feb 27"), cex.axis=0.8) axis(1, at = c(1:50), labels = rep("",50), cex.axis=0.8, tck=-0.01) lines((predict(modscaled))/Nm*100) title("modeling with Ntot = 4.155*10^5, I(0) = 45, K = 0.217721, R0 = 1.087", cex.main = 1)