[R-sig-dyn-mod] metapopulation model with random effect in parameters
胡艺
@tevenhuyi @ending from gm@il@com
Wed Sep 26 10:29:52 CEST 2018
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 using 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 using 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)
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
[[alternative HTML version deleted]]
More information about the R-sig-dynamic-models
mailing list