Thank you for your help on the deterministic model once again, but it appears I still need help concerning my model.
I am trying to do particle filtering on the model below. I made some changes to include Euler multinomial distribution instead of the chained Binomial model after going through your work on Measles. The conditional likelihood goes to 0 very fast, and I am unable to get the standard error for the Monte-Carlo estimate for the loglikelihood. I get the value NA with a warning message "In max(x) : no non-missing arguments to max; returning -Inf" . I initially thought it has to do with my rSim function, but I still arrive at the same problem after changing it to the new model below.
rSim <- Csnippet("
double lambda;
double rate[6], trans[6];
//transmission rate
double beta1 = R0 * g1; //transmission rate for asymptomatic
double beta2 = R0 * g2; //transmission rate for mild
// expected force of infection (iota is the imported infection)
//lambda = (beta1 * pow(Ia + Ip + iota, alpha) + beta2 * Im) / N;
lambda = ((beta1*(Ia + Ip)) + beta2*Im)/ N;
//S
double dN_SE = rbinom(S, 1 - exp(-lambda*dt)); //movement form S to E
//E
double dN_EIp = rbinom(E, 1 - exp(-k*dt)); //movement from E to Ip
//Ip
rate[0] = tau1;
rate[1] = tau2;
reulermultinom(2, Ip, &rate[0], dt, &trans[0]);
double dN_IpIa = trans[0]; //rate of developing mild infection
double dN_IpIm = trans[1]; //rate of aymptomatic infection
//Ia
double dN_IaR = rbinom(Ia, 1 - exp(-g1*dt)); //rate of recovery form aymptomatic infection
//Im
rate[2] = g1;
rate[3] = delta;
reulermultinom(2, Im, &rate[2], dt, &trans[2]);
double dN_ImR = trans[2]; //rate of recovery form mild infection
double dN_ImIC = trans[3]; //rate of developing severe symptoms
//IC
rate[4] = g2;
rate[5] = phi;
reulermultinom(2, IC, &rate[4], dt, &trans[4]);
double dN_ICR = trans[4]; //rate of recovery at isolation centers
double dN_ICD = trans[5]; //rate of deaths at isolation centers
// Balance the equations
S -= dN_SE;
E += dN_SE - dN_EIp;
Ip += dN_EIp - dN_IpIa - dN_IpIm;
Ia += dN_IpIa - dN_IaR;
Im += dN_IpIm - dN_ImIC - dN_ImR;
IC += dN_ImIC - dN_ICR - dN_ICD;
R += dN_IaR + dN_ImR + dN_ICR;
D += dN_ICD;
H += dN_IaR + dN_ImR + dN_ICR + dN_ICD; //true incidence
I += dN_EIp + dN_IpIa + dN_IpIm + dN_ImIC; //total infected
//Rprintf(\"%lg %lg %lg %lg %lg \\n\",lambda, dN_SE, H, E, Ip);
N_SE += dN_SE;
N_EIp += dN_EIp; //number transition from E to Ip
N_IpIa += dN_IpIa; //number transition from Ip to Ia
N_IpIm += dN_IpIm; //number transition from Ip to Im
N_IaR += dN_IaR; //number transition from Ia to R
N_ImIC += dN_ImIC; // number of transition from Im to IC
N_ImR += dN_ImR; //number of transitions from I to R
N_ICR += dN_ICR; // number of transitions from IC to R
N_ICD += dN_ICD; //number of transitions from IC to D
")
#initial values simulator
rInit <- Csnippet("
double m = N /(S_0 + E_0 + Ip_0 + Ia_0 + Im_0 + IC_0 + R_0 + D_0);
S = nearbyint(m*S_0);
E = nearbyint(m*E_0);
Ip = nearbyint(m*Ip_0);
Ia = nearbyint(m*Ia_0);
Im = nearbyint(m*Im_0);
IC = nearbyint(m*IC_0);
R = nearbyint(m*R_0);
D = nearbyint(m*D_0);
H = 0;
I = 0;
N_SE = 0;
N_EIp = 0;
N_IpIa = 0;
N_IpIm = 0;
N_IaR = 0;
N_ImIC = 0;
N_ImR = 0;
N_ICR = 0;
N_ICD = 0;
")
#random generator of the measurement process
rmeas <- Csnippet("
if (tetha > 0.0) {
Cases = rnbinom_mu(1.0/tetha, H);
} else {
Cases = rpois(H);
}
"
)
rmeas2 <- Csnippet("
if (tetha > 0.0) {
Cases = rnbinom_mu(1.0/tetha, I);
} else {
Cases = rpois(I);
}
"
)
## density of the measurement process
dmeas <- Csnippet("
double f;
if(ISNA(Cases)) {
lik = (give_log) ? 0 : 1;
} else {
if (tetha > 0.0) {
f = dnbinom_mu(Cases, 1.0/tetha, H, give_log);
} else {
f = dpois(Cases, H, give_log);
}
}
lik = (give_log) ? f : exp(f);
Rprintf(\"%lg %lg %lg \\n\",f, H, Cases);
lik;
")
#for the deterministic skeleton
dmeas2 <- Csnippet("
double f;
if(ISNA(Cases)) {
lik = (give_log) ? 0 : 1;
} else {
if (tetha > 0.0) {
f = dnbinom_mu(nearbyint(Cases), 1.0/tetha, I, give_log);
} else {
f = dpois(nearbyint(Cases), I, give_log);
}
}
lik = (give_log) ? f : exp(f);
Rprintf(\"%lg %lg \\n\",f, I);
lik;
")
#Deterministic skeleton
skel <- Csnippet("
double lambda;
// transmission rate
double beta1 = R0 * g1; //transmission rate for asymptomatic
double beta2 = R0 * g2; //transmission rate for mild
// expected force of infection (iota is the imported infection)
//lambda = (beta1 * pow(Ia + Ip + iota, alpha) + beta2 * Im) / N;
lambda = ((beta1*(Ia + Ip)) + beta2*Im)/ 20000000;
//suceptible class
DS = -lambda * S;
//exposed class
DE = (lambda * S) - (k * E);
//presymptomatic infectious class
DIp = (k * E) - (tau1 * Ip + tau2 * Ip);
//asymptomatic infection class
DIa = (tau1 * Ip) - (g1 * Ia);
//mild infection class
DIm = (tau2 * Im) - (g1*Im + delta*Im);
//isolation center
DIC = (delta*Im) - (g2*IC + phi*IC);
//Recovery
DR = (g1*Ia) + (g1*Im) + (g2*IC);
//deaths
DD = phi*IC;
//true incidence
DH = g1*Ia + g1*Im + g2*IC + (phi*IC);
DI = (k * E) + (tau1*Ip) + (tau2 * Ip);
DN_SE = lambda * S;
DN_EIp = k * E;
DN_IpIa = tau1 * Ip;
DN_IpIm = tau2 * Ip;
DN_IaR = g1 * Ia;
DN_ImIC = delta*Im;
DN_ImR = g1*Im;
DN_ICR = g2*IC;
DN_ICD = phi*IC;
")
### pomp formulation
# Population size
N <- 20000000
covid_model <- function(state = unique(naija_covid_data2$State)[1:37],
covid_data = naija_covid_data2,
timestep = 1/24, pop = N, type = c("I", "H")) {
state <- match.arg(state)
type = match.arg(type)
dt <- timestep
globs <- paste0("static int N = ", pop, ";")
state_covid_data <- covid_data %>%
filter(State == state) %>%
mutate(Cases = as.numeric(Cases),
Deaths = as.numeric(Deaths),
Recovered = as.numeric(Recovered),
Date = as.Date(Date, format = "%d/%m/%Y")) %>%
#filter(!is.na(Cases)) %>%
mutate(dt = c(0, diff(Date)),
Day = cumsum(c(1, diff(Date)))
)
state_covid_data %>%
dplyr::select(Day, Cases) %>%
pomp(
times = "Day",
t0 = 0,
obsnames = c("Cases"),
globals = globs,
rprocess = euler(rSim, delta.t = dt),
rinit = rInit,
rmeasure = if(type == "H") rmeas else rmeas2,
dmeasure = if(type == "H") dmeas else dmeas2,
skeleton = vectorfield(skel),
accumvars = c("H", "I","N_SE", "N_EIp", "N_IpIa", "N_IpIm",
"N_IaR", "N_ImIC", "N_ImR", "N_ICR", "N_ICD"),
statenames = c("S", "E", "Ip", "Ia", "Im", "IC", "R", "D",
"H", "I","N_SE", "N_EIp", "N_IpIa", "N_IpIm",
"N_IaR", "N_ImIC", "N_ImR", "N_ICR", "N_ICD"),
paramnames = c("R0", "k", "tau1", "tau2", "delta",
"g1", "g2", "phi", "tetha", "S_0", "E_0",
"Ip_0", "Im_0", "Ia_0",
"IC_0", "R_0", "D_0"
),
partrans = parameter_trans(
log = c("R0", "k", "g1", "g2", "tau1", "tau2", "delta", "phi", "tetha"),
baycentric = c("S_0", "E_0", "Ip_0", "Im_0", "Ia_0"
, "IC_0", "R_0", "D_0"
)
)
) -> dis_model
}
# parameter estimation
#death rate function
drates <- function(X_death, Time_death, mu) {
(1 - mu*Time_death) * (X_death / Time_death)
}
# # recovery rate function
recrates <- function(X_death, Time_rec, mu) {
(1 - mu*Time_rec) * (1 - X_death) / Time_rec
}
cfr <- 0.028
T_death <- 7 #time to death in ic
T_ic_rec <- 10 #time to recover in ic
T_as_rec <- 3.5 #time to recover aymptotic
T_md_rec <- 3.5 #time to recover mild
T_mic <- 1.5 #time before develop developing severe infection
g1 <- recrates(cfr, T_as_rec, 0)
#g2 <- recrates(cfr, T_md_rec, 0)
g2 <- recrates(cfr, T_as_rec, 0)
g3 <- recrates(cfr, T_ic_rec, 0)
#sigma <- recrates(cfr, T_q, mu_Nig)
phi <- drates(cfr, T_death, 0)
params <- c(R0 = 1.75, #basic reproduction number (to estimate)
k = 1/3, #latency period (to estimate)
tau1 = 1/3.5, #rate of developing aymptomatic infection (to estimate)
tau2 = 1/3.5, #rate of developing aymptomatic infection (to estimate)
delta = 1/2, #rate of developing severe infection (estimate)
g1 = 1/7, #rate of recovery for asymptomatic infections (4.5 days) (to estimate)
g2 = 1/10, #recovery rate in isolation centers (fix)
phi = 1/7, #death rate in isolation centers (fix)
tetha = 0.0, #overdispersion parameter (to estimate)
E_0 = 8,
Ip_0 = 0,
Ia_0 = 1,
Im_0 = 0,
IC_0 = 0,
R_0 = 0,
D_0 = 0,
S_0 = N - 11
)
Lagos_covid_model <- covid_model(state = "Lagos",
covid_data = naija_covid_data2,
timestep = 1/24,
type = "H"
)
#####test particle filter
Lagos_covid_model %>% pfilter(Np = 100, params = params) -> pf_lagos
pf_lagos %>% logLik() %>% logmeanexp(se = TRUE)
plot(pf_lagos)
kingaa
commented
4 years ago
Thank you for your help on the deterministic model once again, but it appears I still need help concerning my model. I am trying to do particle filtering on the model below. I made some changes to include Euler multinomial distribution instead of the chained Binomial model after going through your work on Measles. The conditional likelihood goes to 0 very fast, and I am unable to get the standard error for the Monte-Carlo estimate for the loglikelihood. I get the value NA with a warning message "In max(x) : no non-missing arguments to max; returning -Inf" . I initially thought it has to do with my rSim function, but I still arrive at the same problem after changing it to the new model below.