Thomas and all, my further question is that how can I use the forcing
function as a vector? In my above coding, I would like to interpolate the
birth rate for each subpopulation but encounter a problem. I've come up
with a solution, that is, the same way just like we deal with signals in
beta (i.e., to define a external birth), but I'm wondering if I can use the
approxfun directly within the trans function. I ran the following code, but
it ends with "Error in lsoda(y, times, func, parms, ...) : REAL() can
only be applied to a 'numeric', not a 'list'".
thanks
Yi
my updated coding is as follows:
require(deSolve)
n<-3
moverate = matrix(nr = n, nc = n, 0.08)
moverate[1,2] <- 0.05
moverate[1,3] <- 0.07
moverate[2,3] <- 0.06
diag(moverate) <- 0
b1<-approxfun(x=c(0,365),y=c(3,6),rule = 2)
b2<-approxfun(x=c(0,365),y=c(4,7),rule = 2)
b3<-approxfun(x=c(0,365),y=c(2,5),rule = 2)
e<-function(x) c(b1(x),b2(x),b3(x))
trans <- function(time, state, pars) {
S <- state[1:n]
E <- state[(n+1):(2*n)]
I <- state[(2*n+1):(3*n)]
R <- state[(3*n+1):(4*n)]
N <- S + E + I + R
birth <- pars[1:n]
a0 <- pars[(n+1):(2*n)]
alpha <- pars[2*n+1]
gamma <- pars[2*n+2]
*beta = a0+sigl[time+1,] *
* birth = e*
dS <- birth(time)-beta*S*I/N
dE <- beta*S*I/N-E*alpha
dI <- E*alpha-I*gamma
dR <- I*gamma
FluxS <- (S*moverate) # described as a matrix[i,j], flux of susceptibles
from subpop i to subpop j
FluxE <- (E*moverate)
FluxI <- (I*moverate) # described as a matrix[i,j], flux of infectives
from subpop i to subpop j
FluxR <- (R*moverate)
dS <- dS - (rowSums(FluxS) - colSums(FluxS))
dE <- dE - (rowSums(FluxE) - colSums(FluxE))
dI <- dI - (rowSums(FluxI) - colSums(FluxI))
dR <- dR - (rowSums(FluxR) - colSums(FluxR))
return(list (c(dS, dE, dI, dR),beta,birth=birth))
}
pars <- c(birth = c(5,4,6), # birth rate
a0 = runif(3,1,3), # guess the contact rate
alpha = 1/4 , # incubation rate
gamma = 1/8) # recovery rate
yini <- c(S = c(998,999,1000),
E = c(1,0,0),
I = c(1,1,0),
R = c(0,0,0))
times <- seq(0, 365, by = 1)
sigl<-matrix(rnorm(n*366,sd=0.5),ncol = n)
out <- ode(yini, times, trans, pars)
plot(out)
On Tue, Sep 25, 2018 at 10:18 AM ?? <stevenhuyi at gmail.com> wrote:
Thanks very much for your detailed explanations and help in the coding, Thomas. I've gained a lot. Yi On Mon, Sep 24, 2018 at 11:45 PM Thomas Petzoldt < thomas.petzoldt at tu-dresden.de> wrote:
Hello Yi,
a derivative with a built-in random component is non-reproducible. It is
not surprising, that an automatic step size integration algorithm fails
under such conditions. More precisely:
(1) a differential equation model is continuous by definition. It has (in
theory) no time steps. In practice, it is distinguished between "internal
time steps" used by the solver to approximate the continuous process with
sufficient precision and "external time steps", where the user wants to get
return values from the simulation.
(2) if you include a random term in the model call, it is unclear at
which time points it is called.
I see two possible solutions:
A) approximate the integration with the Euler method that has fixed steps.
B) define an external "signal" of random values that is passed to the
model and then accessed by the model in a fixed interval. The example below
uses a "1 time unit" interval, but it is of course also possible to use
something else. Please take in mind to restrict the maximum integration
step to the external interval by setting appropriate external time steps or
the ode-argument hmax.
See modified example below,
Thomas
On 23.09.2018 at 15:09 ?? wrote:
Dear all,
I would like to incorporate random effect in transmission rate (i.e.,
beta) for each subpopulation when setting a metapopulation SEIR model. I
used the code in my modelling
*beta<- a0+rnorm(1, sd=0.1)*
but it turns out that the random part is exactly the same for each beta,
I have tried to use
*beta<- a0+rnorm(3, sd=0.1)*
but the program fails. How can I get a set of betas, each of which has
different random part?
Any comments will be appreciated.
Yi
[original example omitted]
Here a suggestion that works technically, but without warranty that the
model itself is correct.
require(deSolve)
n<-3
moverate = matrix(nr = n, nc = n, 0.08)
moverate[1,2] <- 0.05
moverate[1,3] <- 0.07
moverate[2,3] <- 0.06
##for (i in 1:n){moverate[i,i]=0}
## easier:
diag(moverate) <- 0
trans <- function(time, state, pars, sigbeta) {
S <- state[1:n]
E <- state[(n+1):(2*n)]
I <- state[(2*n+1):(3*n)]
R <- state[(3*n+1):(4*n)]
N <- S + E + I + R
## the following can be made simpler if pars would be a list
birth <- pars[1:n]
a0 <- pars[(n+1):(2*n)]
alpha <- pars[2*n+1]
gamma <- pars[2*n+2]
beta<- a0 + sigbeta[floor(time)+1, ] # <----- access signal
dS <- birth-beta*S*I/N
dE <- beta*S*I/N-E*alpha
dI <- E*alpha-I*gamma
dR <- I*gamma
FluxS <- (S*moverate) # described as a matrix[i,j], flux of
susceptibles from subpop i to subpop j
FluxE <- (E*moverate)
FluxI <- (I*moverate) # described as a matrix[i,j], flux of infectives
from subpop i to subpop j
FluxR <- (R*moverate)
dS <- dS - (rowSums(FluxS) - colSums(FluxS))
dE <- dE - (rowSums(FluxE) - colSums(FluxE))
dI <- dI - (rowSums(FluxI) - colSums(FluxI))
dR <- dR - (rowSums(FluxR) - colSums(FluxR))
return(list (c(dS, dE, dI, dR), beta))
}
pars <- c(birth = c(5,4,6), # birth rate
a0 = runif(3,1,3), # guess the contact rate
alpha = 1/4 , # incubation rate
gamma = 1/8) # recovery rate
yini <- c(S = c(998,999,1000),
E = c(1,0,0),
I = c(1,1,0),
R = c(0,0,0))
times <- seq(0, 365, by = 1)
sigbeta <- matrix(rnorm(n * 366, sd=0.1), ncol=n) # <----- external
signal
out <- ode(yini, times, trans, pars, sigbeta=sigbeta, hmax=1)
plot(out)