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

[胡艺]

Thanks for your detailed explanations, Thomas, and that is really helpful.
I am still struggling with fitting the model to data. I ran the following
code, but it failed. Any responses are appreciated.

Yi

require(deSolve)

n<-3
moverate = matrix(nr = n, nc = n, 0.008)
moverate[1,2] <- 0.005
moverate[1,3] <- 0.007
moverate[2,3] <- 0.006
diag(moverate) <- 0

b1<-approxfun(x=c(0,365),y=c(3,10),rule = 2)
b2<-approxfun(x=c(0,365),y=c(4,20),rule = 2)
b3<-approxfun(x=c(0,365),y=c(2,15),rule = 2)
b<-cbind(b1(0:365),b2(0:365),b3(0:365))

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

  a0 <- pars[1:n]
  a1 <- pars[(n+1):(2*n)]
  alpha <- pars[2*n+1]
  gamma <- pars[2*n+2]

  birth = b[time+1,] ###observed birth number for each subpopulation
  beta = a0*(1+a1*cos(2*pi*time/365)) ## transmission rate for each
subpopulation

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

}

pars  <- c(a0  = c(2,0.8,3.2),    # guess the contact rate
           a1  = c(0.02,0.01,0.1),
           alpha  = 1/4 ,   # incubation rate
           gamma = 1/8)    # recovery rate

yini  <- c(S = c(999,1000,1000),
           E = c(0,0,0),
           I = c(1,0,0),
           R = c(0,0,0))

times <- seq(0, 365, by = 1)

out <- ode(yini, times, trans, pars)
plot(out)

require(FME)

#to arbitarily simulate observed I1, I2, and I3, because only infected
individuals can be observed in epidemiological study###
sa1<-sample(1:365,100)
sa2<-sample(1:365,100)
sa3<-sample(1:365,100)

x <- seq(1, 365, 1)
y <- dnorm(x, mean=80, sd=18)*10000
plot(x, y, type="l", lwd=1)
sub1<-data.frame(time=x,infect=y)[sa1[order(sa1)],]

x <- seq(1, 365, 1)
y <- dnorm(x, mean=100, sd=20)*10000
plot(x, y, type="l", lwd=1)
sub2<-data.frame(time=x,infect=y)[sa2[order(sa2)],]

x <- seq(1, 365, 1)
y <- dnorm(x, mean=50, sd=15)*10000
plot(x, y, type="l", lwd=1)
sub3<-data.frame(time=x,infect=y)[sa3[order(sa3)],]

obs<-data.frame(name=c(rep("I1",50),rep("I2",50),rep("I3",50)),
                time=c(sub1$time,sub2$time,sub3$time),
                val = c(sub1$infect,sub2$infect,sub3$infect))

plot(obs$time,obs$val,type="l")

## simulate and plot the model and the data
plot(out, obs=obs, obspar=list(pch=15, col="red"))

cost <- function(p) {
  yy<-p[1:(4*n)]  # the initial values
  pp<-p[(4*n+1):(6*n+2)]
  out  <-  ode(yy, times, trans, pp)
  modCost(out, obs, weight = "none",y="val")
}

p<-c(S = c(1000,1020,1000),
     E = c(0,0,0),
     I = c(1,0,0),
     R = c(0,0,0),a0  = c(1,1.2,2.2),
     a1 = c(0.01,0.02,0.03),
     alpha  = 1/3.5 ,
     gamma = 1/6)

## Note: initial parameters taken from above, may be adjusted manually
fit <- modFit(f = cost, p = p)
summary(fit)

## Now plot original and fitted models and data
out1 <- ode(yini, times, trans, pars)
out2 <- ode(coef(fit)[1:(4*n)], times, trans, coef(fit)[(4*n+1):(6*n+2)])




On Thu, Sep 27, 2018 at 4:25 PM Thomas Petzoldt <
thomas.petzoldt using tu-dresden.de> wrote:

> Hi,
>
> this is a common task for the dynamic modeler, its solution depends on
> the model and the data. Package 'simecol' contains two convenience
> functions approxTime and approxTime1, but for small small models like
> yours its easiest to write just three approxfun'ctions. You can put them
> in a list, if you wish. approxfun is quite efficient, given that the
> interpolation tables are not too big.
>
> For larger models, one may consider to write an own interpolation
> function. This works well and can be faster than approx, appproxfun or
> approxTime, if data are given in equidistant resolution, because you can
> then directly calculate the position of grid points without searching.
>
> As an alternative, you can also supply a Fortran or C-model to deSolve.
> Then deSolve provides a built-in interpolator.
>
> Hope it helps,
>
> Thomas
>
> Am 27.09.2018 um 03:20 schrieb 胡艺:
> > Thomas, thanks for your studying files and they are very useful. As for
> the
> > above problem, there are three subpopulations in my model and I want each
> > to have different birth rate. If I just use  b1(time) in the model, there
> > is only one birth rate for one subpopulation at each integration time,
> > right? I'm wondering whether there is a way to let approxfun() to return
> > three birth rates (a vector) at each integration time. Thanks.
> >
> > Yi
> >
> > On Wed, Sep 26, 2018 at 4:44 PM Thomas Petzoldt <
> > thomas.petzoldt using tu-dresden.de> wrote:
> >
> >> Hi,
> >>
> >> Note that approxfun uses a special feature of R: It returns a function.
> >>
> >> If you create:
> >>
> >> b1<-approxfun(x=c(0,365),y=c(3,6),rule = 2)
> >>
> >> then just use
> >>
> >> b1(time)
> >>
> >> in your model.
> >>
> >>
> >> Please read the ?approxfun help page again. For examples how to use it
> >> with desolve, have a look at:
> >>
> >> http://limno-live.hydro.tu-dresden.de/user2017/
> >>
> >> and in particular:
> >>
> >>
> >>
> http://limno-live.hydro.tu-dresden.de/user2017/2.1-forcing+events/desolve-forcing.html#(5)
> >>
> >> Thomas
> >>
> >>
> >> --
> >> http://tu-dresden.de/Members/thomas.petzoldt
>
>

	[[alternative HTML version deleted]]