3

年齢層別のSEIRモデルをコーディングしようとしています。つまり、私の微分方程式には、20歳以上のベータ*（感染率）*（影響を受けやすい数）の合計である質量作用のパラメーターがあります。透過係数（ベータ）は、接触マトリックスから計算されます。接触マトリックスには、年齢クラス（rows = person i、columns = person j）を表す20の列と行があり、任意の年齢クラスの2人の接触の確率が含まれています。私はそれを設計して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.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はベータを見つけることができ、ベータは数学の特殊関数として解釈されるため、エラーが発生します。パラメータを次のように定義する必要があります

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

xstartこれにより上記の問題は修正されますが、SEIR0によって返される導関数の数は初期条件ベクトルの長さ（ここでは4）と等しくなければならないため、ODEは機能しません。

私は例えば提案します：

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