deSolve包可以参数包括矩阵吗？

I'm trying to code an SEIR model that is age-stratified; that is, in my differential equations I have a parameter for mass action that is the sum of beta*(proportion infected)*(number susceptible) over 20 age classes. The transmission coefficient (beta) is calculated from a contact matrix. The contact matrix has 20 columns and rows which represent the age classes (rows=person i, columns=person j), and contains the probability of contact between two people in any age classes. I designed it and read it into R. My problem is I don't know how (or if) I can use a matrix inside my parameters with deSolve. This code below I wrote doesn't work, I believe because of the matrix / I get this error:

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

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

Before I fool with it too much, I'd like to know if it is possible to use a matrix like this as a parameter for this model.

在我愚弄它之前，我想知道是否可以使用像这样的矩阵作为这个模型的参数。

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))

Also, if I print the parameters, this is what beta looks like:

此外，如果我打印参数，这就是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

....through $beta.V20. So I think it's creating 20 vectors, each with 20 arguments...I think each is a row of the original matrix 'mat' multiplied by the constant 0.04? However, when I multiply mat*0.04 outside "parameters" I get the expected matrix. I'm struggling a bit with how to realize these equations using deSolve and would appreciate any advice on whether it's possible. Thanks in advance.

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

1 个解决方案

#1

The error occurs in this line :

此行中出现错误：

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

non-numeric argument to binary operator means you try to multiply a function for example by a numeric. You can reproduce it by I*1

二元运算符的非数字参数意味着您尝试将函数乘以例如数字。您可以通过I * 1重现它

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

Here, R can' find beta , and beta is interpreted as the Special Functions of Mathematics, so the error. You need to define parameters as

在这里，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)

and

和

 ## 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)

I think you can even rewrite your function as:

我想你甚至可以改写你的功能：

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
  })
}

This will correct the problem above but the ODE will not work because the number of derivatives returned by SEIR0 must be equal to the length of the initial conditions xstart vector(4 here).

这将解决上述问题，但ODE将无法工作，因为SEIR0返回的导数必须等于初始条件xstart向量的长度（此处为4）。

I suggest for example :

我建议例如：

  res <- c(dS=mean(dS),dE=mean(dE),dI=dI,dR=dR)
  list(res)

#1