3

我正在尝试编写一个按年龄分层的 SEIR 模型；也就是说，在我的微分方程中，我有一个质量作用参数，它是 20 个年龄段的 beta*（感染比例）*（易感人数）之和。传输系数 (beta) 由接触矩阵计算得出。接触矩阵有 20 列和代表年龄类别的行（行 = 人 i，列 = 人 j），并包含任何年龄类别的两个人之间的接触概率。我设计了它并将其读入 R。我的问题是我不知道如何（或是否）可以在 deSolve 的参数中使用矩阵。我写的下面这段代码不起作用，我相信是因为矩阵/我得到了这个错误：

Error in beta * S : non-numeric argument to binary operator

在我胡思乱想之前，我想知道是否可以使用这样的矩阵作为该模型的参数。

mat <-as.matrix(read.csv("H:/IBS 796R/contactmatrix.csv", header=F))

times <- seq(0,20,by=1/52)
parameters <- c(mu=0,v=1/75,N=1,p=0,delta=2.4,beta=mat*0.04,sigma=1/8,gamma=1/15)
xstart <- c(S=0.06,E=0,I=0.001,R=0)

SEIR0 <- function(t,x,parameters){
    S=x[1]
    E=x[2]
    I=x[3]
    R=x[4]
    with(as.list(parameters), {
        dS=v*S -beta*S*I/N -delta*S
        dE=beta*S*1/N -E*(sigma+delta)
        dI=sigma*E -I*(gamma+delta)
        dR=gamma*I-delta*R
        res=c(dS,dE,dI,dR)
        list(res)
    })
}

out <- as.data.frame(lsoda(xstart,times,SEIR0,parameters))

另外，如果我打印参数，这就是 beta 的样子：

$beta.V1
 [1] 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04
[12] 4e-04 4e-04 8e-03 8e-03 8e-03 8e-03 8e-03 8e-03 8e-03

$beta.V2
 [1] 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04 4e-04
[12] 4e-04 4e-04 8e-03 8e-03 8e-03 8e-03 8e-03 8e-03 8e-03

....通过 $beta.V20。所以我认为它创建了 20 个向量，每个向量都有 20 个参数……我认为每个向量都是原始矩阵“mat”的一行乘以常数 0.04？但是，当我在“参数”之外乘以 mat*0.04 时，我得到了预期的矩阵。我正在为如何使用 deSolve 实现这些方程而苦苦挣扎，并希望得到任何关于它是否可能的建议。提前致谢。

4

1 回答 1

3

错误发生在这一行：

dS=v*S -beta*S*I/N -delta*S

non-numeric argument to binary operator意味着您尝试将一个函数例如乘以一个数字。您可以通过以下方式重现它I*1

Error in I * 1 : non-numeric argument to binary operator`

在这里，R 找不到 beta ，而 beta 被解释为数学的特殊函数，所以错误。您需要将参数定义为

# a list 
list(mu=0,v=1/75,N=1,p=0,delta=2.4,beta=mat*0.04,sigma=1/8,gamma=1/15)

和

 ## you get a vector mu,N,p,delta,beta1,bet2,...  
c(mu=0,v=1/75,N=1,p=0,delta=2.4,beta=mat*0.04,sigma=1/8,gamma=1/15)

我认为你甚至可以将你的函数重写为：

SEIR0 <- function(t,x,parameters){
  with(as.list(c(parameters, x)), {
    dS = v*S -beta*S*I/N -delta*S    ## matrix
    dE = beta*S*1/N -E*(sigma+delta) ## matrix
    dI = sigma*E -I*(gamma+delta)
    dR = gamma*I-delta*R
    res = c(dS,dE,dI,dR)
    list(res)                        ## different of the structure of xstart
  })
}

这将纠正上述问题，但 ODE 将不起作用，因为 SEIR0 返回的导数数量必须等于初始条件xstart向量的长度（此处为 4）。

我建议例如：

  res <- c(dS=mean(dS),dE=mean(dE),dI=dI,dR=dR)
  list(res)
于 2013-03-11T06:01:17.873 回答