我想为R上的COVID大流行建立标准的SEIR模型。我有很多病例，死亡病例和日常病例（发病率）。我采用贝叶斯框架，希望使用简单的MCMC来校准和预测模型。在查阅了文献之后，我试图在下面输入一些代码，但是我不理解错误。任何人都可以帮助完善代码或提供一些代码吗？提前谢谢你的帮助！
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)
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]
print(case)
obs <- cbind(dates,deaths)
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.27,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[,5])
r0_sd <- sd(post.samp[,5])
rt_mean <- mean(post.samp[,6])
rt_sd <- sd(post.samp[,6])
print(c(rt_mean,rt_mean-1.96*rt_sd,rt_mean+1.96*rt_sd))
nsamp <- length(post.samp[,6])
rts <- sort(post.samp[,6])
analyse_ensemble(post.samp)
run.obj <- run_ensemble(post.samp,500,120)
plot_result(run.obj,obs,post.samp,case,range=c(10,120))
错误：
Warning message:
In cbind(dates, deaths) :
number of rows of result is not a multiple of vector length (arg 2)
Hessian from call to optim() not negative definite.
Sampling (as specified) cannot proceed.
Error: Check data and fun() and call MCMCmetrop1R() again.
dates未在您的代码中定义。当您尝试使用cbind个不同长度的矢量（此处为deaths和dates）时，您会遇到问题。