[R-sig-dyn-mod] metapopulation model with random effect in parameters

[胡艺]

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