[R-sig-dyn-mod] metapopulation model with random effect in parameters
胡艺
@tevenhuyi @ending from gm@il@com
Sat Sep 29 11:02:28 CEST 2018
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]]
More information about the R-sig-dynamic-models
mailing list