deSolve中不同参数的不同时间间隔的值

Question

我正在尝试使用deSolve解决SVEIR（敏感，接种疫苗，暴露，感染和移除）模型。爆发从第8天开始（通过在易感人群中导入指标病例）。用于捕获此我利用一个事件（通过添加值一（1）在时间t的状态变量（1）= 8。

# Model's parameters 

parms <- c(beta=1.29, 
     betaE=0.25, 
     betaI=1, 
     betaV=0.0, 
     sigma=0.5, 
     gama=0.2, 
     delta=1/365, 
     m=0.000046, 
     r=0.000052, 
     kapa=1.857/10000,    
     alpha=0.00643, 
     thita=1/365, 
     f=0.002)  
dt <- seq(0,50,0.25)  

inits <- c(S=14900, V=0, E=0, I=0, R=0)  
N <- sum(inits) 

eventdat <- data.frame(var = c("I"),time = c(8), 
        value = c(1), method = c("add")) 
eventdat 

#The SVEIR model 

SVEIR <- function(t, x, parms){ 

with(as.list(c(parms,x)),{ 
dS <- - beta*betaE*E*(S/N) - beta*betaI*I*(S/N) - f*S - m*S +delta*R + thita*V + r*N 
dV <- - beta*betaE*betaV*E*(V/N) - beta*betaI*betaV*I*(V/N) - m*V - thita*V + f*S 
dE <- + beta*betaE*E*(S/N) + beta*betaI*I*(S/N) + beta*betaE*betaV*E*(V/N) + beta*betaI*betaV*I*(V/N) - (m + kapa + sigma)*E 
dI <- + sigma*E - (m + alpha + gama)*I 
dR <- kapa*E + gama*I - m*R - delta*R  
der <- c(dS, dV, dE, dI, dR) 
list(der)  
}) 

} 

library(deSolve) 

out <- as.data.frame(lsoda(inits, dt, SVEIR, parms=parms, events = list(data = eventdat))) 

# Plotting the output 

attach(out) 

matplot(x = out[,1], y = out[,-1], type = "l", lwd = 2, 
    lty = "solid", col = c("red", "blue", "black", "green", "darkgreen"), 
    xlab = "time", ylab = "y", main = "SVEIR model") 

legend("bottomright", col = c("red", "blue", "black", "green", "darkgreen"), 
    legend = c("S", "V", "E", "I", "R"), lwd = 2) 

除此之外，我希望我的模型也捕获。在一些参数所以，我一直在努力（不成功迄今为止）的变化，以我的功能“而”内集成或“for”循环，其考虑以下几点：

1. 的时间段在0 - 9之间我需要参数 betaV的值为0
2. 对于10之间的时间段 - 50我需要的 值参数betaV为0.002

我曾尝试使用一个事件，而是[R给我一个错误（我想我可以利用事件仅用于变量而不用于参数）。

任何想法如何处理这？

非常感谢，

汤姆

PS：（。Samsuzzoha等，2012），该模型是基于工作。

Answer 1

你的基本问题似乎是如何指定两个不同的值betaV取决于时间。你不能在功能上做到这一点，如：

#The SVEIR model 
SVEIR <- function(t, x, parms){ 
    with(as.list(c(parms,x)),{ 
    betaV <- ifelse(t<10,betaV,0.002) # adjust betaV based on value of t 
    dS <- - beta*betaE*E*(S/N) - beta*betaI*I*(S/N) - f*S - m*S +delta*R + thita*V + r*N 
    dV <- - beta*betaE*betaV*E*(V/N) - beta*betaI*betaV*I*(V/N) - m*V - thita*V + f*S 
    dE <- + beta*betaE*E*(S/N) + beta*betaI*I*(S/N) + beta*betaE*betaV*E*(V/N) + beta*betaI*betaV*I*(V/N) - (m + kapa + sigma)*E 
    dI <- + sigma*E - (m + alpha + gama)*I 
    dR <- kapa*E + gama*I - m*R - delta*R  
    der <- c(dS, dV, dE, dI, dR) 
    list(der)  
    }) 

请注意，你的问题实际上并不指定betaV上9 < t < 10的值，所以我估计在10

截止当我运行这个与betaV = 0.002 (t>10)，没有明显的差异在输出。如果我将betaV设置为1或10的值为t > 10，则V(t)对于较大的t和S, E, I and R被抑制为较低的时间。这听起来正确吗？