0

I would like to do a standard SEIR model for the COVID pandemic on R. I have number of cases, deaths and daily cases (incidence). Adopting a Bayesian framework, I would like to calibrate and forecast the model using simple MCMC. Having consulted with the literature, I have tried to input some code below but I'm stuck with defining the MCMC (here, post.samp): nothing is working; R runs it without producing an output. Would anyone be able to help refine the code or offer some code? Thank you for your help in advance! 

 library(deSolve)
odefun <-function(t,state,parameters){
  with(as.list(c(state, parameters)),{
    beta <- parameters[1]
    sigma <- parameters[2]
    gamma <- parameters[3]

    x<- state

    dx <- rep(0,6)
    dx[1] <- -beta*x[1]*(x[4] + x[5]) #susceptible
    dx[2] <- beta*x[1]*(x[4] + x[5]) - sigma*x[2] #newly infected but latent
    dx[3] <- sigma*x[2] - sigma*x[3] #late stage latent
    dx[4] <- sigma*x[3] - gamma*x[4] #newly infectious
    dx[5] <- gamma*x[4] - gamma*x[5] #late infectious
    dx[6] <- gamma*x[5] #recovered
    return(list(dx))
  })}

library(date)

dates = as.Date(c("2020-3-2", "2020-3-3", "2020-3-4", "2020-3-5", "2020-3-6", "2020-3-7", "2020-3-8", "2020-3-9", "2020-3-10", "2020-3-11", "2020-3-12",
                  "2020-3-13", "2020-3-14", "2020-3-15", "2020-3-16", "2020-3-17", "2020-3-18", "2020-3-19", "2020-3-20", "2020-3-21", "2020-3-22",
                  "2020-3-23", "2020-3-24", "2020-3-25", "2020-3-26", "2020-3-27", "2020-3-28", "2020-3-29", "2020-3-30", "2020-3-31", "2020-4-1",
                  "2020-4-2", "2020-4-3", "2020-4-4", "2020-4-5", "2020-4-6", "2020-4-7", "2020-4-8", "2020-4-9", "2020-4-10", "2020-4-11", "2020-4-12",
                  "2020-4-13", "2020-4-14", "2020-4-15", "2020-4-16", "2020-4-17", "2020-4-18", "2020-4-19", "2020-4-20", "2020-4-21", "2020-4-22", 
                  "2020-4-23", "2020-4-24", "2020-4-25", "2020-4-26", "2020-4-27", "2020-4-28", "2020-4-29", "2020-4-30", "2020-5-1", "2020-5-2",
                  "2020-5-3", "2020-5-4", "2020-5-5", "2020-5-6", "2020-5-7", "2020-5-8", "2020-5-9", "2020-5-10", "2020-5-11", "2020-5-12", "2020-5-13",
                  "2020-5-14", "2020-5-15", "2020-5-16", "2020-5-17", "2020-5-18"))

epi <- data.frame(row.names = c("Morocco"), 
                  startdate =  as.Date(c("2020-3-2")),
                  interventiondate=as.Date(c("2020-3-20")),
                  N=c(3.7e7))

InfectedC <- c(1,1,1,2,2,2,2,2,3,6,6,8,18,28,37,44,54,63,79,96,115,143,170,225,275,345,
               402,479,556,617,654,708,791,919,1021,1120,1184,1275, 1374, 
               1448, 1545, 1661,1763,1888,2024,2283,2564,2685,2855,3046,3209,
               3446, 3568, 3758, 3897, 4065, 4120, 4252, 4321, 4423, 4569, 4729, 
               4903, 5053, 5219, 5408, 5548, 5711, 5910, 6063, 6281, 6418, 6516, 6607, 6652,
               6741, 6870, 6952) #C = cumulative

Infected <-c(1,0,0,1,0,0,0,0,1,3,0,2,10,10,9,7,10,9,16,17,19,28,27,55,50,70,57,77,77,
             61,37,54,89,122,102,99,64,91,99,74,97,116,102,125,136,259, 281, 121, 170, 
             191, 163, 237, 122, 190, 139, 168,55, 132, 69, 102, 146, 170, 174, 150, 166, 
             189, 140, 163, 199, 153, 218, 137, 98, 91, 45, 89, 139, 82)
# Infected = daily infections
deaths <- c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 6, 10, 
            23, 25, 26, 33, 36, 39, 44, 48, 59, 70, 80, 90, 93, 97, 107, 111, 118,126,126,
            127,130,135,137,141,143,145,149, 155, 158, 159, 161, 162, 165, 168, 170, 171, 
            173, 174, 179, 181, 183, 183, 186, 186, 188, 188, 188, 188, 190, 190, 192, 192, 192)

N <- 3.603e7 # population of Morocco according to WB 2019 estimate

case <- "Morocco"

#choose a long or short run type
run_type <- "full" 
# full for long test
#run_type <- "test" for short

startdate <- epi[case,1]
interventiondate <- epi[case,2]
N <- epi[case,3]

#parameters for error calculation
report_err <- 0.2
model_err <- 0.05

#prior mean for parameters in order
par_pri <- c(4.5,2.5,-15,0.0075,3,1.)

#prior sd
par_sd <- c(.5,0.5,15,0.00125,1,.5)

print(case)

obs <- cbind(dates,deaths) 

  #code to run model and evaluate cost (log likelihood) vs observations. Including prior.

  modelcost <- function(params,obs){ 

    latent_period <- max(.5,min(params[1],10)) #bound with 0.1 and 10
    infectious_period <- max(.5,min(params[2],10)) #bound with 0.1 and 10
    i0 <- max(0.,min(exp(params[3]),.01)) #bound with 0. and 0.01 NB this one is logarithmic!
    death <- max(0.001,min(params[4],0.05)) #bound with 0.1 and 5%

    R0 <- max(1.,min(params[5],10)) #bound with 0.1 and 10 also not less than 1 in initial segment
    Rt <- max(0.01,min(params[6],10)) #bound with 0.1 and 10

    #prior mean for parameters in order
    par_pri <- c(4.5,2,-15,0.0075,3,1.)

    #prior sd
    par_sd <- c(.5,0.1,15,0.0003,1,.5)

    ###testing differnet priors    
    #prior mean for parameters in order
    par_pri <- c(4.5,2.5,-15,0.0075,3,1.)

    #prior sd
    par_sd <- c(.5,0.5,15,0.00125,1,.5)

  }


library(MCMCpack)


burn<-3000
runlength<-5000


if(run_type == "test"){
  burn<-000
  runlength<-500
}

set.seed(42) #reproducibility!
obs_save <- obs
num_obs <- length(obs_save[,2])


vals <- c("2020-03-22","2020-03-29","2020-04-05","2020-04-12","2020-04-19","2020-04-26","2020-05-03","2020-05-10",0) #
theta.start <- c(5.1,2.5,-15,0.0075,2.72,0.76)
post.samp <- MCMCmetrop1R(modelcost, theta.init=theta.start,#par_pri,
                          obs=obs,thin=1, mcmc=runlength, burnin=burn,
                          verbose=500, logfun=TRUE)

r0_mean <- mean(post.samp)
r0_sd <- sd(post.samp)
rt_mean <- mean(post.samp)
rt_sd <- sd(post.samp)

print(c(rt_mean,rt_mean-1.96*rt_sd,rt_mean+1.96*rt_sd))


nsamp <- length(post.samp)
rts <- sort(post.samp)

analyse_ensemble(post.samp)
run.obj <- run_ensemble(post.samp,500,120)
plot_result(run.obj,obs,post.samp,case,range=c(10,120))

Errors:

Error in optim(theta.init.0, maxfun, control = optim.control, lower = optim.lower,  : objective function in optim evaluates to length 6 not 1

Error in if (dates[tomorrow] > interventiondate) { : argument is of length zero

1 Answers1

0

dates is not defined in your code. When you try to cbind vectors of different lengths (here, deaths and dates), you have a problem.

Limey
  • 10,234
  • 2
  • 12
  • 32
  • I don't know. It's your data, not mine! :) Seriously though, it looks like you have a series of counts of deaths and infections. Presumably, there is a date associated with each. Those dates would be a logical candidate. But like I said, it's your data. I don't have the context to tell you. – Limey Jun 04 '20 at 11:59
  • Could have a look at the amended code and see what needs improving or what else is missing - thank you! –  Jun 05 '20 at 14:51
  • With all due respect, Zak, your code does not provide [a simple self-contained example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example/5963610#5963610). It's 266 lines long, including a function (`centile`) whch is nearly 200 lines long, contains several local function definitions, loads various libraries and yet is never called. Your code as supplied runs no analyses and produces no output. You have more fundamental problems to solve before you even consider trying to fit or debug your model. – Limey Jun 05 '20 at 16:40
  • Should I remove the centile function if it doesn't contribute to anything? I have cut some of the functions (see code above) which weren't necessary but one error appears. Could you help, pls? You state fundamental problems, which are graver than one's I've highlighted, but I'm not sure what they are and how to solve them. Could you elaborate further pls. Thank you! –  Jun 05 '20 at 18:32