kingaa
commented
5 years ago The problem is that your deterministic skeleton is a vectorfield, not a map. That is, your deterministic model evolves in continuous time. See the following, which updates your code to pomp2.
set.seed(594709947L)
library(tidyverse)
theme_set(theme_bw())
library(pomp2)
library(readxl)
library(tgp)
demogra <- read_xls("demogra.xls")
data_1 <- read_excel("japan_data.xlsx",sheet=1)
skel<- "
double rate[6];
double term[6];
double Beta;
Beta = betam*(1 + betas * cos(M_2PI*(t+pha)));
rate[0] = birthrate1;
rate[1] = Beta * (I + iota)/pop;
rate[2] = deathrate1;
rate[3] = gamma;
rate[4] = deathrate1;
rate[5] = deathrate1;
term[0]=rate[0];
term[1]=rate[1]*S;
term[2]=rate[2]*S;
term[3]=rate[3]*I;
term[4]=rate[4]*I;
term[5]=rate[5]*R;
// for (int k = 0; k < 6; k++) Rprintf(\"%lg \",rate[k]);
// Rprintf(\"\\n\");
DS= term[0] - term[1] - term[2];
DI= term[1] - term[3] - term[4];
DR= term[3] - term[5];
DH= term[1];"
data_1 %>%
mutate(
date=as.Date(sprintf("%04d-%02d-01",year,month)),
time=julian(date,origin=as.Date("2000-01-01"))/365.25+2000
) %>%
select(time,case="Coxsackievirus A4") %>%
filter(time<2015) %>%
pomp(times = "time", t0=2000-1/12,
skeleton=vectorfield((Csnippet(skel))),
covar = covariate_table(
time=with(demogra,seq(from=min(year)-1/12,to=max(year)+1/12,by=1/12)),
pop=with(demogra,predict(smooth.spline(x=year,y=pop),x=time)$y),
birthrate1=with(demogra,predict(smooth.spline(x=year,y=birth,spar = 0.2),x=time)$y),
deathrate=with(demogra,predict(smooth.spline(x=year,y=death, spar = 0.2),x=time)$y),
deathrate1=deathrate/pop,
times="time"
),
# covar = covariate_table(covar1,times="time"),
accumvars = c("H"),
statenames = c("S", "I", "R", "H"),
paramnames = c("gamma", "rho", "sigma",
"S.0", "I.0", "R.0",
"betas" , "betam", "iota","pha"),
rinit = function (S.0, I.0, R.0, pop, ...) {
s <- pop/(S.0+I.0+R.0)
c(S=round(s*S.0),I=round(s*I.0),R=round(s*R.0),H=0)
},
params = c(iota = 92, betas = 0.25, betam = 1706,
gamma = 52, pha=0.7,
rho = 1.94e-04, sigma =0.0000001,
S.0 = 0.0254, I.0 = 1.48e-04, R.0 = 0.9745)
) -> model1
estpar_name <- c("iota","betas","betam","rho","pha","S.0","I.0","R.0")
rbind(
iota = c(0,1000),
betas = c(0,1),
betam = c(1,10000),
rho = c(0,0.3),
pha=c(0,1),
S.0=c(0,0.1),
I.0=c(0,0.1)
) -> para
lhs(40000,para) %>%
as.data.frame() %>%
setNames(rownames(para)) %>%
mutate(
R.0 = 1-S.0-I.0,
gamma = 52,
sigma = 1e-5
) -> estpar
model1 %>%
trajectory(params=t(estpar[1:20,]),format="d",maxstep=0.001) %>%
filter(time>2000) %>%
ggplot(aes(x=time,y=H,group=.id))+
geom_line()This shows that the trajectory calculations are free of NaN and NA. From here, to do trajectory matching, you construct an objective function using traj_objfun (adding parameter transformations if needed). This is passed to the optimizer of your choice to do the optimization.
minjuner
When estimating the SIR model parameters, I first set a parameter range, then use the trajectory matching to select some of the desirable parameters, but there is an error in running traj.match: Error in is.nloptr(ret) : objective in x0 returns NA.
Here is my code and data file.zip