deSolve软件包参数可以包含矩阵吗? [英] deSolve package can parameters include a matrix?

问题描述

我正在尝试编写年龄分层的SEIR模型；也就是说，在我的微分方程式中，我有一个针对大众行为的参数，它是20个年龄段的beta *(感染比例)*(易感人数)的总和.根据接触矩阵计算透射系数(β).接触矩阵有20个列和行，分别代表年龄段(行=人i，列=人j)，并且包含两个年龄段中两个人之间的接触概率.我设计了它并将其读入R.我的问题是我不知道如何(或者是否)可以在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:

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

此外，如果我打印参数，则beta如下所示:

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

$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来实现这些方程式方面有些挣扎，并且希望就是否有可能提供任何建议.预先感谢.

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

推荐答案

此行中发生错误:

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

non-numeric argument to binary operator表示您尝试将一个函数乘以一个数字.您可以通过I*1

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

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

在这里，R可以找到beta，并且beta被解释为数学的特殊函数，因此是错误的. 您需要将参数定义为

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

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

我认为您甚至可以将函数重写为:

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

这将解决上面的问题，但是ODE将不起作用，因为SEIR0返回的导数的数量必须等于初始条件xstart vector(此处为4)的长度.

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

我建议例如:

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

这篇关于deSolve软件包参数可以包含矩阵吗?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！